An SIRS model with nonmonotone incidence and saturated treatment in a changing environment.

Nonmonotone incidence and saturated treatment are incorporated into an SIRS model under constant and changing environments. The nonmonotone incidence rate describes the psychological or inhibitory effect: when the number of the infected individuals exceeds a certain level, the infection function decreases. The saturated treatment function describes the effect of infected individuals being delayed for treatment due to the limitation of medical resources. In a constant environment, the model undergoes a sequence of bifurcations including backward bifurcation, degenerate Bogdanov-Takens bifurcation of codimension 3, degenerate Hopf bifurcation as the parameters vary, and the model exhibits rich dynamics such as bistability, tristability, multiple periodic orbits, and homoclinic orbits. Moreover, we provide some sufficient conditions to guarantee the global asymptotical stability of the disease-free equilibrium or the unique positive equilibrium. Our results indicate that there exist three critical values [Formula: see text] and [Formula: see text] for the treatment rate r: (i) when [Formula: see text], the disease will disappear; (ii) when [Formula: see text], the disease will persist. In a changing environment, the infective population starts along the stable disease-free state (or an endemic state) and surprisingly continues tracking the unstable disease-free state (or a limit cycle) when the system crosses a bifurcation point, and eventually tends to the stable endemic state (or the stable disease-free state). This transient tracking of the unstable disease-free state when [Formula: see text] predicts regime shifts that cause the delayed disease outbreak in a changing environment. Furthermore, the disease can disappear in advance (or belatedly) if the rate of environmental change is negative and large (or small). The transient dynamics of an infectious disease heavily depend on the initial infection number and rate or the speed of environmental change.

Introduction

Once again the COVID-19 pandemic raised the alarm for humans to prevent and control infectious diseases. A suitable epidemic model can provide helpful insights to understand, predict and control the transmission of an emerging infectious disease. In the epidemic model formulation, the incidence rate and the treatment rate play key roles in disease dynamics. We will incorporate a more realistic version of incidence and treatment into the classical susceptible-infectious-recovered-susceptible (SIRS) model with treatment:where the variables and parameters are listed in Table 1. Here the SIRS model assumes that new recruits are susceptible and recovered individuals have temporary immunity.

Table 1

Definitions of variables and parameters in model (1.1)

Variables	Description
S(t)	The susceptible population at time t
I(t)	The infected population at time t
R(t)	The recovered population at time t
f(I)	The infectious force
T(I)	The treatment function
b	The recruitment rate of the population
d	The per capita natural death rate of the population
μ	The per capita natural recovery rate of the infective individuals
γ	The per capita rate of recovered individuals who lose immunity and return to the susceptible class

Many variations of the incidence rate were proposed in the literature. The simplest one is a bilinear incidence rate (Kermack and McKendrick 1927) (where is the infection rate). Since it ignores the crowding effect of infected individuals, protection measures and awareness or inhibition, the bilinear incidence rate is not realistic enough for many infectious diseases. To describe the behavioral change and crowding effect of infected individuals, Capasso et al. (1977) and Capasso and Serio (1978) used a saturated incidence force to describe a “crowding effect” in modeling the cholera epidemics in Bari in 1973 (see Liu and Yang 2012; Yang et al. 2012; Zhou et al. 2014 for additional references). This function is a nonlinear increasing function in I but eventually tends to a saturation level as I increases to infinity. Notice that both the bilinear and saturated incidence rates are monotone. However, the incidence rate exhibits nonmonotonicity in some epidemic diseases, induced by “psychological” effect (Capasso and Serio 1978): the contact rate and the infection probability usually increase when a new infectious disease emerges because people have little knowledge about the disease, however, when more and more individuals are infected by the disease, people usually tend to reduce the number of contacts. To incorporate the psychological effect, Xiao and Ruan (2007) proposed the nonmonotone incidence rate , which is increasing when but decreasing when .

Treatment is a pivotal methodology in the control of an infectious disease. Different types of treatment rate were proposed by many researchers for different epidemic scenarios. Traditional epidemic models assume the linear treatment rate , which is reasonable given sufficient medical resources including vaccines, medicines, hospital beds, etc. However, the supply of these medical facilities is always limited once the disease spans rapidly and widely. For this reason, Wang and Ruan (2004) introduced the constant treatment function . Later Wang (2006) modified the constant treatment function to a piecewise-smooth treatment function , where r is a positive constant and is the infective level at which the health care system reaches its capacity. Furthermore, Zhang and Liu (2008) proposed a continuously differentiable treatment function , where is the cure rate, measures the extent of the effect of the infected population being delayed for treatment, and for small I, , whereas for large I, . Moreover, when , the saturated treatment function returns to the linear one. There have been many studies (Zhang and Liu 2008; Zhang and Suo 2010; Ghosh et al. 2021; Jana et al. 2016; Zhou and Meng 2012) that used the saturated treatment function to characterize limited medical resources.

Ghosh et al. (2021) proposed an SIR model with nonmonotone incidence and saturated treatment as well as the disease-induced death rate and vaccination. They obtained a necessary and sufficient condition for backward bifurcation, and investigated saddle-node and Hopf bifurcations. Meanwhile, subcritical Hopf bifurcation was exhibited by numerical simulation. Pan et al. (2022) formulated an SIRS model by considering nonmonotone incidence and piecewise-smooth treatment, and found that the model can possess Bogdanov-Takens and subcritical Hopf bifurcations. In this paper, we replace the piecewise-smooth treatment function by a saturated treatment function as below:which is a specific form of model (1.1). In addition, model (1.2) extends the model in Wang (2006) that assumed a bilinear incidence rate and the model in Xiao and Ruan (2007) that assumed no treatment.

The response of populations to environmental changes has been and remains a cutting-edge area in ecological modeling (Arumugam et al. 2020, 2021; Xiang et al. 2022). Following the same logic in studying environmental changes in ecology, modeling a changing environment is equally important in epidemic research. Linking transient dynamics under a changing environment to stable and unstable states in a constant environment provides novel insights into understanding disease transmission mechanisms and outbreaks. The infection rate highly depends on nonpharmaceutical interventions, weather conditions, holidays, gatherings, and many other factors. To describe the impact of environmental changes on the infection rate, we assume that the infection rate is a linear function of time t as the simplest dependence. We formulate the following model in a changing environment:where u is the speed of environmental change. From the fourth equation of system (1.3), we have ( is the initial value), which represents the possible directional environmental change.

In a constant environment (i.e., ), system (1.3) is reduced to system (1.2). We will study the dynamics of model (1.2) via rigorous bifurcation analysis. In a changing environment (i.e., ), we will study how the speed of environmental change regulates the dynamics of system (1.3). Meanwhile, we will compare and contrast dynamics in the changing environment with those obtained by bifurcation analysis in the constant environment. In addition, the long-term dynamics and persistence of epidemic systems are predicted under continuous environmental changes.

The remaining paper is organized as follows. In Sect. 2, we simplify system (1.2) into a two-dimensional system based on the fact that the total population size is constant, then we discuss the stability and types of disease-free and endemic equilibria. Furthermore, we study the bifurcations and global dynamics of system (2.2). In Sect. 3, we explore the impact of environmental change on the transient dynamics of system (1.3). We summarize the results and suggest future research directions in the last section.

Constant environment

Under a constant environment , we will study the dynamical behaviors of model (1.2) on the plane where the limit set of model (1.2) locates.

Model simplification

Model (1.2) assumes there is no disease-caused mortality, the total population size is hence completely determined by the vital dynamics in model (1.2). Therefore, the total population size eventually tends to a constant. This allows us to simplify system (1.2) into a two-dimensional system.

Lemma 2.1

The plane is an invariant manifold of system (1.2), which is attracting in the first octant.

Proof

Summing up the three equations in (1.2) and denoting , we haveIt is clear that is a solution and for any , the general solution isThuswhich implies the conclusion.

Lemma 2.1 implies that the limit set of system (1.2) is on the plane . Thus we focus on the reduced systemWe rescale system (2.1) by introducing then system (2.1) becomes (still denote by t)whereandIt is easy to see that the positive invariant and bounded region of system (2.2) is

Equilibria and their types

To find the equilibria of system (2.2) in , we setwhich yieldClearly, system (2.2) always has a disease-free equilibrium . By the next generation matrix method in Driessche and Watmough (2002), we find the basic reproduction number

Remark 2.1

From (2.3) and (2.6), using the original parameters, we have , it is obvious that decreases as r increases, which implies that we will overestimate the risk of the disease if we take in system (2.1).

We have the following results for the types of , which are important to determine the global dynamics in .

Theorem 2.2

The disease-free equilibrium of system (2.2) is The phase portraits are given in Fig. 1.

Fig. 1

Types of disease-free equilibrium of system (2.2). a A hyperbolic stable node if . b A hyperbolic saddle if . c A saddle-node with a stable parabolic sector in the right half plane if and ; d A stable degenerate node if and

a hyperbolic stable node if ;

a hyperbolic saddle if ;

a degenerate equilibrium if , more precisely,

when , is a saddle-node with a stable (or unstable) parabolic sector in the right half plane if (or ;

when , is a stable degenerate node.

The Jacobian matrix of system (2.2) at the equilibrium isIt is easy to see that the two eigenvalues of are and . Obviously, when , then is a hyperbolic stable node.

When , we have , . We first make the following transformationstill denote by t, and the Taylor expansions of system (2.2) around origin are as follows:whereBy Theorem 7.1 of Zhang et al. (1992), we know that is a saddle-node if .

If , then , by the center manifold theorem, we assume and substitute it into the first equation of (2.7). By using the second equation of system (2.7), we haveand substitute into the second equation of system (2.7), then the reduced equation restricted to the center manifold is described asNotice that , by Theorem 7.1 of Zhang et al. (1992), and we have made a time transformation , then is a stable degenerate node.

To find the endemic equilibria of system (2.2), we setthen an endemic equilibrium of system (2.2) is given by , where x is a positive real root of the cubic equation . Then, we discuss the existence of the positive real roots of (i.e., the existence of the endemic equilibria of (2.2)) according to the properties of (2.10). First of all, we notice thatandIf (i.e., ) and , all coefficients of f(x) are positive. Thus, has no positive root due to Descartes’ Rule of signs (or the facts and for all ).

If (i.e., ) and , has two rootsConsequently, (i.e., f(x) is an increase function) for , but (i.e., f(x) is a decreasing function) for . Combining these and the property (2.11), we can see that has at most two positive real roots (see Fig. 2).

Fig. 2

The positive real roots of when and : a no positive root. b A double positive root (i.e., ). c Two single positive roots ,

If (i.e., ) and , it is easy to see that has no positive root. If (i.e., ) and , then has a unique positive root.

If (i.e., ), whether or , only has a single positive root. In fact, from , we can get . If , then for (i.e., f(x) is increasing in ), and by (2.11), only has a positive root. If , then has two roots . Moreover, f(x) is a decreasing function for , then . Since f(x) is increasing in , then by (2.11), only has a positive root.

We further discuss the stability of endemic equilibria. The Jacobian matrix of system (2.2) at any equilibrium E(x, y) isFrom , we havethenIt implies that E(x, y) is a hyperbolic saddle if , an elementary equilibrium if , and a degenerate equilibrium if . We have the following results.

Lemma 2.3

Let f(x) and be given by (2.10) and (2.13), respectively. System (2.2) has at most two positive equilibria. Moreover,

when , we have

if , or and , then system (2.2) has no positive equilibrium;

if and , then system (2.2) has a unique positive equilibrium , which is degenerate and , ;

if and , then system (2.2) has two positive equilibria and , which are all elementary and is a hyperbolic saddle. , , ;

when , we have

if , then system (2.2) has no positive equilibrium;

if , then system (2.2) has a positive equilibrium ;

when , then system (2.2) has a positive equilibrium .

From (2.15) and Fig. 2, it is easy to see that , and , then and are all elementary equilibria and only is a hyperbolic saddle, and is a degenerate equilibrium.

We first consider case of Lemma 2.3. In this case, system (2.2) has a unique positive equilibrium , which is degenerate. From ==0, v and n can be expressed by m, p, s, q and asMoreover, from =0 and (2.17), m can be expressed by p, s, q and asLetwe have the following results.

Theorem 2.4

If , and the conditions in (2.4) and (2.17) are satisfied, then system (2.2) has a unique positive equilibrium . Moreover,

if , then is a saddle-node, which includes a stable parabolic sector (an unstable parabolic sector) if ();

if , then is a cusp. Moreover,

if , then is a cusp of codimension 2;

if , then is a cusp of codimension 3.

The Proof of Theorem 2.4 is given in Appendix A.

Next we discuss the case (I)(iii), (II)(ii) and (III) of Lemma 2.3. Let

Theorem 2.5

If equilibria and exist, then is always a hyperbolic saddle, and

when , is a hyperbolic stable node or focus,

when ,

is a hyperbolic unstable node or focus if , and ;

is a hyperbolic stable node or focus if ;

is a weak focus or center if , and .

From (2.15) and Fig. 2, it is easy to see that and , then and are all elementary equilibria and is a hyperbolic saddle. From (2.16), we haveIf , then ; if , it is easy to see that , and if , (the conditions and guarantee ), and , respectively, leading to the conclusions.

Bifurcation analysis

In this subsection, we discuss different kinds of bifurcations for system (2.2) in depth.

Saddle-node bifurcation

From Theorem 2.4, we know that the surfaceis a saddle-node bifurcation surface. When the parameters vary to cross the surface from one side to the other, the number of positive equilibria of system (2.2) change from zero to two, the saddle-node bifurcation yields two positive equilibria.

Backward bifurcation

From Wang (2006) (see also Lu et al. (2021) and Zhang and Liu (2008)), backward bifurcation is an interesting and crucial topic in epidemic models. For some classical epidemic models, the basic reproduction number serves as a threshold in the sense that a disease is persistent if and goes extinct if , where the transition from a disease-free equilibrium to an endemic equilibrium is called a forward bifurcation. While this bifurcation is backward if an endemic equilibrium occurs when . From Lemma 2.3, we know that system (2.2) can have one or two positive equilibria when , and . Then we have the following result.

Theorem 2.6

System (2.2) admits a backward bifurcation as crosses one if and .

Degenerate Bogdanov-Takens bifurcation of codimension three

The case (II)(ii) of Theorem 2.4 indicates that system (2.2) may exhibit a degenerate Boganov-Takens bifurcation of codimension 3 around . In this subsection, we want to make sure if such a bifurcation can be fully unfolded inside the class of system (2.2).

Letwhere , , , , , and are given in (2.17), (2.18), (2.19) and (A10) of Appendix A, respectively.

Firstly, we choose n, m and q as bifurcation parameters and obtain the following unfolding system:where and . If we can transform, by a series of near-identity transformations, the unfolding system (2.23) into the following versal unfolding of a cusp of codimension 3:whereand , then we can claim that system (2.23) (i.e., system (2.2)) undergoes a degenerate Bogdanov-Takens bifurcation of codimension three (see Dumortier et al. (1987) and Chow et al. (1994)). In fact, we have the following theorem.

Theorem 2.7

When and , system (2.2) has a unique positive equilibrium , which is a nilpotent cusp of codimension 3. If we choose n, m and q as bifurcation parameters, then system (2.2) can undergo a degenerate Bogdanov-Takens bifurcation of codimension 3 in a small neighborhood of . Hence, system (2.2) can exhibit the coexistence of a stable homoclinic loop and an unstable limit cycle, coexistence of two limit cycles (the inner one unstable), and a semi-stable limit cycle for different sets of parameters.

Proof

Following the procedure in Li et al. (2015), we use several steps to transform system (2.23) into the versal unfolding of a Bogdanov-Takens singularity (cusp case) of codimension three. We make the following transformations successively:where the expressions of , , , , , , , , and are given in supplementary materials. Then system (2.23) becomes (still denote by t)where has the property of (2.25). To save space, we omit the coefficient , and here. With the help of Mathematica software, we obtain thatsince , and , which have been shown in the Proof of Theorem 2.4, and , and are given in (A19) of Appendix A. System (2.26) is exactly in the form of system (2.24). Therefore, system (2.23) undergoes a degenereate Bogdanov-Takens of codimension 3 in a small neighborhood of .

Following Figs. 2-5 in Dumortier et al. (1987), we next describe the bifurcation diagram and phase portraits in system (2.26), the bifurcation diagram has the conical structure in , starting from =(0, 0, 0). It can be best shown by drawing its intersection with the half sphereTo clearly see the traces of the intersections, we draw the projections of traces onto the -plane in Fig. 3.

Fig. 3

The projection of bifurcation diagram for system (2.26) on S

In Fig. 3, both curves H and C are tangent to (the boundary of S) at the points and , and the two curves cross each other at point d. Along , the circle , excluding the points and , is a saddle-node bifurcation curve, while and correspond to Bogdanov-Takens bifurcation points of codimension 2; H is the Hopf bifurcation curve, on which is a Hopf bifurcation point of codimension 2; C is the homoclinic bifurcation curves, on which is a homoclinic bifurcation point of codimension 2; L is the saddle-node bifurcation curve of limit cycles, which is tangent to curves H and C at and , respectively.

We next present some numerical simulations plotted by Matcont to illustrate the existence of almost all bifurcations appearing in Fig. 3. Firstly, we obtain the numerical bifurcation diagram of system (2.2) in plane as shown in Fig. 4 by fixing , , , . While we could not numerically plot homoclinic bifurcation curve C, we will give evidences for its existence in Fig. 4. Indeed, we plot some phase portraits in Figs. 5, 6 and 7, where parameter values and corresponding dynamical behaviors are given in Tables 2, 3 and 4, respectively. The parameter values are taken from three parallel lines , and in Fig. 4, respectively.

Fig. 4

a Bifurcation diagram for system (2.2) in plane when , , , . BT and GH denote Bogdanov-Takens bifurcation point and degenerate Hopf bifurcation point, respectively. Blue, red and green curves denote saddle-node bifurcation SN, Hopf bifurcation H, saddle-node bifurcation L of limit cycles, respectively. b The local enlarged view of (a) (colour figure online)

Fig. 5

Phase portraits of system (2.2) with , , , and . a ; b ; c ; d ; e ; f . The detailed dynamical behaviors are described in Table 2

Fig. 6

Phase portraits of system (2.2) with , , , and . a ; b ; c ; d ; e ; f . The detailed dynamical behaviors are described in Table 3

Fig. 7

Phase portraits of system (2.2) with , , , and . a ; b ; c ; d . The detailed dynamical behaviors are described in Table 4

Table 2

, , , and

m	Positive equilibria and types	Closed orbits and homoclinic orbits
1.59	No	No(Fig. 5a )
1.57	E_1(saddle), E_2(unstable focus)	No (Fig. 5b)
1.566815	E_1(saddle), E_2(stable focus)	Unstable limit cycle (Fig. 5c)
1.566758405	E_1(saddle), E_2(stable focus)	Unstable limit cycle, a homoclinic orbit (Fig. 5d)
1.5667506	E_1(saddle), E_2(stable focus)	Unstable limit cycle, a stable limit cycle (Fig. 5e)
1.56	E_1(saddle), E_2(stable focus)	No (Fig. 5f)

Table 3

, , , and

m	Positive equilibria and types	Closed orbits and homoclinic orbits
=1.58	No	No (Fig. 6a)
1.569	E_1(saddle), E_2(unstable focus)	No (Fig. 6b)
1.563985	E_1(saddle), E_2(unstable focus)	A homoclinic orbit (Fig. 6c)
1.563956	E_1(saddle), E_2(unstable focus)	A stable limit cycle (Fig. 6d)
1.56389	E_1(saddle), E_2(stable focus)	A stable limit cycle, unstable limit cycle (Fig. 6e)
1.562	E_1(saddle), E_2(stable focus)	No (Fig. 6f)

Table 4

, , , and

m	Positive equilibria and types	Closed orbits and homoclinic orbits
1.56	E_1(saddle), E_2(unstable focus)	No (Fig. 7a)
1.555476	E_1(saddle), E_2(unstable focus)	A homoclinic orbit (Figure  7b)
1.55495	E_1(saddle), E_2(unstable focus)	A stable limit cycle (Fig. 7c)
1.55	E_1(saddle), E_2(stable focus)	No (Fig. 7d)

From Fig. 5 and Table 2, we can see that when m decreases, system (2.2) undergoes successively saddle-node bifurcation, subcritical Hopf bifurcation, attracting homoclinic bifurcation and saddle-node bifurcation of limit cycles.

From Fig. 6 and Table 3, as m decreases, system (2.2) undergoes successively saddle-node bifurcation, attraction homoclinic bifurcation, subcritical Hopf bifurcation and saddle-node bifurcation of limit cycles;

From Fig. 7 and Table 4, as m decreases, system (2.2) undergoes successively attraction homoclinic bifurcation, supercritical Hopf bifurcation.

From the phase portraits in Figs. 5, 6 and 7, we can infer that there exists a curve C in Fig. 4. As m decreases, the curve C is tangent to the curve SN at BT point, intersects the curve H at d point and finally intersects the curve L at point.

Hopf bifurcation

Now we discuss Hopf bifurcation around in system (2.2). Let

Theorem 2.8

When , , , and the conditions in (2.4) are satisfied, then is a weak focus with multiplicity one for system (2.2). Where is given in (2.20).

When , we have from . We make the following transformations successively:where and , then system (2.2) becomes as (still denote , and by x, y and t, respectively)where and () are given in supplementary materials.

Using the formal series method in Zhou et al. (2014) and Mathematica software, we obtain the first Lyapunov coefficientSince , then the denominator of is negative, therefore the sign of is determined by . If and , then , which leads to the conclusion.

There exists a set of parameter values that satisfies the conditions of Theorem 2.8 such that . On the other hand, there exists a set of parameter values that also satisfies the conditions of Theorem 2.8 such that . Therefore, there exists an open set in the parameter space such that , i.e.,And there exists another open set in the parameter space such that , i.e.,Summarizing the above discussion, we have the following results.

Theorem 2.9

If and , then the equilibrium of system (2.2) is a stable weak focus with multiplicity one. There exists a stable limit cycle arising from supercritical Hopf bifurcation (see Fig. 8a);

Fig. 8

a A stable limit cycle arising from supercritical Hopf bifurcation. b An unstable limit cycle arising from subcritical Hopf bifurcation

If and , then the equilibrium of system (2.2) is an unstable weak focus with multiplicity one. There exists an unstable limit cycle arising from subcritical Hopf bifurcation (see Fig. 8b).

Using the formal series method in Zhou et al. (2014) and Mathematica software, when and , we can obtain the second Lyapunov coefficient , which is given in supplementary materials.

Next, we present an example to show that the positive equilibrium of system (2.2) is a stable weak focus of multiplicity two, system (2.2) can undergo a degenerate Hopf bifurcation around and two limit cycles occur.

Theorem 2.10

When , system (2.2) has a positive equilibrium which is a stable weak focus of multiplicity two. There exist two limit cycles arising from degenerate Hopf bifurcation around , the repelling cycle is surrounded by an attracting one (see Fig. 9).

Fig. 9

Two limit cycles (the inner one is unstable) around in system (2.2)

We first fix , , and , from the expression of , we get , furthermore, from , and , we get , and , respectively. Therefore, the positive equilibrium is a stable weak focus of multiplicity two when . Next, we can check by Mathematica thatwhich means that and with respect to m, q have full rank 2. We first perturb q such that q decreases to , then becomes an unstable weak focus with multiplicity one, a stable limit cycle occurs around which is the outer limit cycle in Fig. 9. Secondly, we perturb m such that m increases , then becomes a stable hyperbolic focus, another unstable limit cycle occurs around , which is the inner limit cycle in Fig. 9.

Global dynamics of system (2.2)

In this subsection, we discuss the global asymptotical stability of the unique boundary equilibrium of system (2.2) (i.e., the disease-free euilibrium or the unique positive equilibrium of system (1.2)).

Theorem 2.11

The unique boundary equilibrium of system (2.2) (i.e., the disease-free euilibrium of system (1.2)) is globally asymptotically stable if one of the following conditions holds:

and (see Fig. 1a and c);

, and (see Fig. 10a).

Fig. 10

a of system (2.2) is globally asymptotically stable when , and . b of system (2.2) is globally asymptotically stable when , and

Lemma 2.1 implies that the stability of disease-free equilibrium of system (1.2) in the interior is equivalent to that of the equilibrium of system (2.2) in , thus we only need to discuss the stability of equilibrium of system (2.2) in .

On one hand, is a positive invariant and bounded region and is an invariant line. On the other hand, from Theorem 2.2, we know that system (2.2) has no positive equilibrium and only has a boundary equilibrium when the conditions in (i) or (ii) are satisfied. Thus, Poincar-Bendixson Theorem implies that of system (2.2) is globally asymptotically stable under the condition (i) or (ii).

Remark 2.2

Using the original parameters, we have ; , whereFrom Theorem 2.11, we can see that the disease will disappear for all positive initial populations if one of the following cases holds:

;

Theorem 2.12

If the conditions in (2.4) and hold, then system (2.2) does not have nontrivial periodic orbits in the interior of .

Taking a Dulac function , we havefor and , we can immediately getThus, we obtain the conclusion by Dulac’s criteria.

Theorem 2.13

The unique positive equilibrium is globally asymptotically stable if one of the following conditions is satisfied:

, and (see Fig. 10b);

and (see Fig. 1b).

From the Theorem 2.2, Lemma 2.3 and Theorem 2.12, we can see that (I) if , and , then the unique positive equilibrium of system (2.2) is globally asymptotically stable since the unique boundary equilibrium is a saddle-node with a stable parabolic sector lying in the left half plane of and is a positive invariant and bounded region. (II) if and , then the unique positive equilibrium is globally asymptotically stable since the unique boundary equilibrium is a hyperbolic saddle and is a positive invariant and bounded region. We use simulations to illustrate the results. We take a set of parameter values: , , , , , such that , and for case (I), as shown in Fig. 10b: is globally asymptotically stable. We take a set of parameter values: , , , , , such that and for case (II), as shown in Fig. 1b: is globally asymptotically stable.

Remark 2.3

Using the original parameters, we have , whereTheorem 2.13 indicates that the disease will persist for all positive initial populations if one of the following cases holds:

.

Remark 2.4

If the cure rate (i.e., no treatment), then and . From Theorem 2.2, Lemma 2.3 and Theorem 2.12, if , we have the following conclusions, which are the results in Xiao and Ruan (2007):

When , system (2.2) has a unique disease-free equilibrium , which is a global attractor in the first octant.

When , system (2.2) has two equilibria: a disease-free equilibrium and an endemic equilibrium , where is a global attractor in the interior of the first octant.

Changing environment

In this section, we study how the speed of environmental change regulates the dynamics of system (1.3). In detail, we first use the bifurcation software Matcont to plot the one-parameter bifurcation diagram in plane for system (1.2), then plot some “representative” trajectories (time series) for system (1.3) and include them into the bifurcation diagram (see Figs. 11, 12 and 13). In Figs. 11, 12 and 13, the stable and unstable steady states for system (1.2) are described by blue solid and dashed curves, respectively. The maximum and minimum values of stable and unstable oscillations are described by blue filled and open circles, respectively. The “representative” trajectories (time series) for system (1.3) are plotted by red or green solid curves.

Fig. 11

a Bifurcation diagram in plane for system (1.2), where , and . b The curve of . c Time series in model (1.3) with and the initial point: . Other parameters: , , , , , , ,

Fig. 12

a) Bifurcation diagram in plane for system (1.2) (, and ). b The curve of (). c Time series in model (1.3) with , different initial points: (green curve) and (red curve). d Time series in model (1.3) with (green curve) or (red curve), and the same initial point: . Other parameters: , , , , , , (colour figure online)

Fig. 13

a) Bifurcation diagram in plane for system (1.2) (, and ). b The curve of (). c Time series in model (1.3) with (green curve) or (red curve) and the same initial point: . d Time series in model (1.3) with and the initial point: . Other parameters: , , , , , , , (colour figure online)

In system (1.3), the infection rate k(t) continuously increases over time when , and continuously decreases when . According to the bifurcations and dynamics of system (1.2) given in the previous sections, we classify the effect of k in system (1.3) as the two cases.

Case (I): and system (1.2) has no endemic equilibrium.

In this case, the basic reproduction number is a threshold in the sense that a disease is persistent if and goes extinct if in system (1.2).

In Fig. 11a, we choose , , , , , , in system (1.2), then we get , and . We can see that when , system (1.2) only has a stable disease-free equilibrium ; when (i.e., ), system (1.2) undergoes transcritical bifurcation, the disease-free equilibrium becomes unstable and a stable endemic equilibrium occurs when . As k further increases, becomes unstable, system (1.2) exhibits supercritical Hopf bifurcation around the unique endemic equilibrium at or , thus system (1.2) has a stable limit cycle when . In Fig. 11b, we can see that increases as k increases, and if . In Fig. 11c, we choose in system (1.3), the infective population I(t) (red curve) of system (1.3) starts along the stable disease-free state of system (1.2), tracks the unstable disease-free state when (i.e., ), then tends to the stable oscillation and finally to the stable endemic state .

This transient tracking on the unstable disease-free state when predicts regime shifts that cause the delayed disease outbreak under environmental changes.

Case (II): and system (1.2) has endemic equilibria.

In this case, it is not the basic reproduction number but a subthreshold that determines whether the disease persists in system (1.2), i.e., the disease will disappear when and persist when in system (1.2).

Case (II.1): a limit cycle when .

In Fig. 12a, we choose , , , , , , in system (1.2), then we get , and . System (1.2) only has a stable disease-free equilibrium when , and exhibits backward bifurcation (i.e., saddle-node bifurcation) at , where two endemic equilibria and occur and exist when . Furthermore, system (1.2) exhibits subcritical Hopf bifurcation around at , which implies the coexistence of two stable equilibria , and an unstable limit cycle. When , becomes unstable, and system (1.2) has a unique endemic equilibrium which is stable. Figure 12b implies if , and if . In Fig. 12c, it is shown that the outcome of a disease spread for system (1.3) can heavily depend on the initial infection number I(0) and the initial infection rate under the same rate u of environmental changes. We can see that I(t) (green curve) of system (1.3) tracks the stable disease-free state , then the unstable disease-free state when (i.e., ), finally the stable endemic state ; while I(t) (red curve) of system (1.3) tracks the unstable oscillation, and finally the stable endemic state . In Fig. 12d, we choose indicating that the infection rate decreases over time in the changing environment (e.g., nonpharmaceutical interventions), the disease will eventually disappear. For , I(t) (red curve) of system (1.3) reaches the stable disease-free state earlier compared to that of system (1.2), while for , I(t) (green curve) has a delay in disappearance.

Case (II.2): two limit cycles when .

In Fig. 13a, we choose , , , , , , in system (1.2), and get , and . System (1.2) undergoes degenerate Hopf bifurcation around at , and exhibits the coexistence of tristability (two stable equilibria , and a stable big limit cycle) and an unstable small limit cycle (see Fig. 13a). In Fig. 13b, it is shown that if , and if . In Fig. 13c, we choose in system (1.3) indicating that the infection rate increases over time in the changing environment (e.g., more gatherings). I(t) (green curve) of system (1.3) for smaller u starts along the stable disease-free state of system (1.2), and then continues following the unstable disease-free state when (i.e., ), and finally tends to the stable endemic state in an oscillatory way. However, from the same initial state, I(t) (red curve) of system (1.3) for larger u tracks the unstable small cycle first, and then tends to the unique stable endemic equilibrium . We can observe that the rate u affects the transient dynamics of system (1.3). In Fig. 13d, for , I(t) (red curve) of system (1.3) tracks the unstable cycle, and finally tends to the disease-free state . Moreover, I(t) reaches in advance the stable disease-free state , where the corresponding (i.e., ).

Concluding remarks

In epidemiology, the incidence rate and effective treatment are two most important covariates in determining the disease transmission and control. To explore the significance of the awareness factors, crowing effect, limited medical resources, and intervention policies on the dynamics of emerging infectious diseases, we proposed and studied an SIRS model with nonmonotone incidence and saturated treatment. Firstly, we simplified the model into a two-dimensional system (2.2). Then we analyzed the stability and types of disease-free and endemic equilibria as well as different bifurcations of system (2.2) in detail. We found that system (2.2) could have at most two positive equilibria when . The model exhibits rich dynamics, such as bistability (a disease-free equilibrium and an endemic equilibrium), tristability (a disease-free equilibrium, an endemic equilibrium and a big limit cycle), homoclinic orbits, the coexistence of a stable homoclinic loop and an unstable limit cycle, the coexistence of two limit cycles (the inner one unstable and the outer stable), and a semi-stable limit cycle for different sets of parameters. In particular, as parameters vary, the model undergoes saddle-node bifurcation, backward bifurcation, degenerate Bogdanov-Takens bifurcation of codimension three, Hopf bifurcation and degenerated Hopf bifurcation with codimension at least two. Numerical simulations are provided to demonstrate our results.

In Pan et al. (2022), Pan et al. consider an SIRS model with the same nonmonotone incidence rate and a piecewise-smooth treatment rate, where the piecewise-smooth treatment rate describes the situation where the community has limited medical resources, treatment rises linearly with I until the treatment capacity is reached, after which constant treatment (i.e., the maximum treatment) is taken. Their results indicate that there exists a critical value for the infective level at which the health care system reaches its capacity such that: (i) When , separates disease persistence from disease eradication. (ii) When , their model can exhibit multiple endemic equilibria, periodic orbits, homoclinic orbits, Bogdanov-Takens bifurcations, and subcritical Hopf bifurcation. In this paper, for system (2.2) with saturated treatment rate, we showed that there exist three critical values and for the treatment rate r: (i) when , the disease will disappear; (ii) when , the disease will persist. Compared the results in Pan et al. (2022) with ours, system (2.2) can exhibit more complex bifurcation phenomena: degenerate Bogdanov-Takens bifurcation of codimension three and degenerated Hopf bifurcation with codimension at least two. While the system in Pan et al. (2022) can have three endemic equilibria. On the other hand, compared with no treatment in Xiao and Ruan (2007), our results can cover their ones.

The researchers in infectious disease modeling suggested and studied many different incidence functions and treatment functions according to different scenarios. Some incorporated vaccination and even boosted vaccination, which are clearly significant in mitigating and control disease epidemics. The model formulated in this paper only considers treatment control, while vaccination and its boosted version will be mechanistically modeled in our next work. In addition, for mathematical simplicity we assumed no disease-caused mortality in our model. This assumption is reasonable for most diseases with a low death probability due to infection such as flu and COVID-19. It would also be interesting to explore complicated dynamics of model (1.1) by applying other incidence rate functions such as (Xiao and Zhou 2006; Lu et al. 2021) and (Lu et al. 2019).

Linking transient dynamics driven by environmental change to key states in a constant environment can be insightful in epidemiological studies. The transmissibility depends on nonpharmaceutical interventions, weather conditions, etc. Instead of assuming the infection rate as an explicit function of time, future research efforts should focus on the mechanistic modeling of transmissibility versus social behaviors, policies, and environmental factors.