On equilibria stability in an epidemiological SIR model with recovery-dependent infection rate

ABSTRACT

We consider an epidemiological SIR model with an infection rate depending on the recovered population. We establish sufficient conditions for existence, uniqueness, and stability (local and global) of endemic equilibria and consider also the stability of the disease-free equilibrium. We show that, in contrast with classical SIR models, a system with a recovery-dependent infection rate can have multiple endemic stable equilibria (multistability) and multiple unstable saddle points of equilibria. We establish conditions for the occurrence of these phenomena and illustrate the results with some examples.

Keywords:
SIR epidemiological model; recovery-dependent infection rate; endemic equilibria; multistability

RESUMO

Neste artigo é considerado um modelo epidemiológico SIR cuja taxa de infeção depende da população removida. São estabelecidas condições suficientes para a existência, unicidade e estabilidade (local e global) de equilíbrios endêmicos e é considerada também a estabilidade das soluções de equilíbrio livres da infecção. É mostrado que, a diferença do modelo SIR clássico, um sistema com taxa de infeção que depende dos removidos, pode ter múltiplos equilíbrios endêmicos estáveis (multiestabilidade) assim como múltiplos equilíbrios instáveis tipo sela. São apresentadas condições para a ocorrência destes fenômenos e os resultados são ilustrados com alguns exemplos.

Palavras-chave:
modelo epidemiológico SIR; taxa de infeção variável; equilíbrios endêmicos; multiestabilidade

1 INTRODUCTION

Compartmental models, and particularly SIR models, have been extensively used for mathematical modeling of infectious diseases within a population ^{4), (14}.

The main idea behind SIR models is to consider that a population N is divided into three disjoint categories or compartments: susceptible individuals, infected individuals, and recovered or deceased individuals, denoted by S, I, and R, respectively, so that $N=S+I+R$. Depending on the modeling approach, the variables S, I, and R are considered to be the absolute numbers of individuals in each group or the proportion of individuals relative to the total population. In this work, we consider this latter approach.

Within these considerations, an epidemiological SIR model with vital dynamics and constant population can be stated as

$$\begin{gathered} \frac{dS}{dt}=\mu -\beta SI-\mu S \\ \frac{dI}{dt}=\beta SI-\mu I-\gamma I \\ \frac{dR}{dt}=\gamma I-\mu R \end{gathered}$$ (1.1)

with $S\left(0\right)+I\left(0\right)+R\left(0\right)=1$. The positive real numbers µ, β, and γ can be interpreted as birthmortality rate, infection rate, and recovery rate, respectively. For more details about SIR models see for example ¹⁴. Note that from (1.1), we can obtain

$\frac{dN}{dt}=\frac{dS}{dt}+\frac{dI}{dt}+\frac{dR}{dt}=\mu (1-N),$

and since $N\left(t\right)\equiv 1$ is the only solution of this equation, satisfying $N\left(0\right)=1$, we can consider $N\left(t\right)=S\left(t\right)+I\left(t\right)+R\left(t\right)=1$ for all t.

Letting

$\tau =t\mu ;\tilde{\beta }=\frac{\beta }{\mu };\tilde{\gamma }=\frac{\gamma }{\mu };k=1+\tilde{\gamma };\text{and}{R}_0=\frac{\beta }{\mu +\gamma }=\frac{\tilde{\beta }}{1+\tilde{\gamma }}=\frac{\tilde{\beta }}{k},$

we obtain a redimensionalized version of (1.1):

$$\begin{gathered} \frac{dS}{d\tau }=1-k{R}_{0}SI-S \\ \frac{dI}{d\tau }=k{R}_{0}SI-kI \\ \frac{dR}{d\tau }=(k-1)I-R, \end{gathered}$$

with $S\left(0\right)+I\left(0\right)+R\left(0\right)=1$.

Note that the parameters $\tilde{\beta },\tilde{\gamma },k$ and R ₀ are all positive real numbers and in particular, $k>1$. The parameter R ₀ is called the basic reproduction number and its fundamental role in the description of the equilibria stability in the classical SIR model is well known ¹⁴. R ₀ can be interpreted as the number of cases one case generates, on average, in an uninfected population. It represents a measure of the effectiveness of the infection.

Several generalizations and modifications of the SIR model have been proposed by other authors, particularly considering non-constant epidemiological rates (see, ^{1),(7), (8), (12), (15), (19), (21), (22}). These kinds of considerations have been recognized as necessary features to model more realistic epidemic situations, like the interaction between human behavior and disease dynamics ^{13), (21}.

Consider, for example, the population behavior with respect to some possible anti-infection measures (like vaccination, quarantine or sexual precautions). The propagation of the disease can be affected by changes in the population behavior and, in the same way, the risk perception and state of the disease can influence the behavior of the population related to anti-infection measures ⁶. Recent measles outbreaks, for example, are considered to be a direct consequence of the increasing number of unvaccinated children, due to parental behavior and beliefs ²³. In the case of antiretroviral therapy (ART) for HIV, patients under successful ART have lower morbidity and mortality rates and, in many cases, their viral load becomes so low that the patient can be considered almost recovered. Massive scale-up of this successful ART has been considered one of the possible causes of an increase in the practice of sexual risky behaviors and as consequence, an increase in the number of HIV cases and other sexually transmitted diseases ^{3), (11), (24}. Most recently, in the context of the COVID-19 pandemic, human behavior has played a fundamental role in the diseases dynamics ^{2), (5), (10} and, at the same time, it has become evident that the increase in the number of infected and death cases changed the way population and policymakers embrace anti-infection strategies ^{16), (17), (20}.

In all the situations described above, infection rates changed during the evolution of the disease and, in the last two cases, these changes can be considered to be related to changes in the epidemiological variables S, I, or R. In the present paper, we are interested in the stability of equilibria in situations where the infection rate changes depending on the recovered/removed population.

Hence, we propose the following generalized SIR model with a recovery-dependent infection rate:

$$\begin{gathered} \frac{dS}{d\tau }=1-f\left(R\right)SI-S \\ \frac{dI}{d\tau }=f\left(R\right)SI-kI \\ \frac{dR}{d\tau }=(k-1)I-R, \end{gathered}$$ (1.2)

where f is a positive function of R, generalizing the infection rate β, and $S\left(0\right)+I\left(0\right)+R\left(0\right)=1$. The function f can be interpreted as a quantification of the effect on the infection rate, produced by control strategies that depend on the size of R, or, since $S+I+R=1$, that depend on the susceptible and infected population simultaneously.

It is worth noting that ¹⁵ is considered a deterministic model similar to (1.2), and, although the main focus was on the stochastic version, a recurrent solution to the model was obtained when the recovery-dependent infection rate is considered piecewise constant.

Our work focuses on the stability and multistability features of the equilibrium solutions of model (1.2).

The article is organized as follows: In Section 2 we develop a two-dimensional simplified model equivalent to (1.2) and, additional conditions on f, we prove several interesting results, including the non-existence of non-constant positive periodic solutions. In Section 3, the disease-free equilibrium is considered and two results about its local and global stability are established. The results of this section generalize well-known results for the classical SIR model.

Section 4 considers endemic equilibrium points. First, we define an auxiliary function g and establish sufficient conditions for the existence of endemic equilibrium points in terms of f and g. Later, we consider the local stability of endemic equilibrium points and we illustrate conditions for the occurrence of multiple locally stable endemic equilibrium points (multistability). Finally, we consider conditions for the uniqueness and global stability of an endemic equilibrium point. Final comments and concluding remarks are presented in Section 5.

2 SIMPLIFIED MODEL

In this section, we develop a simplified two-dimensional model equivalent to (1.2). In the following lemma, we show that model (1.2) is well defined in the sense that, for all solutions, the conditions $S,I,R\in [0,1]$ and $S+I+R=1$ are preserved under the dynamics described in model (1.2).

Lemma 2.1.The set$\Omega =\{S\ge 0,I\ge 0,R\ge 0\mathit{\text{and}}S+I+R=1\}$is positively invariant under (1.2).

Proof. First, we consider the behavior of the solutions with some initial condition equal to 0.

If $S\left(0\right)=0$, then $\frac{dS}{d\tau }\left(0\right)=1>0$. If $I\left(0\right)=0$, then $\frac{dI}{d\tau }\left(0\right)=0$. If $R\left(0\right)=0$, then $\frac{dR}{d\tau }\left(0\right)=(k-1)I\left(0\right)\ge 0$ since $k>1$ and we consider $I\left(0\right)\ge 0$. This proves the positive invariance of the positive octant.

Consider now $N\left(t\right)=S\left(t\right)+I\left(t\right)+R\left(t\right)$. From (1.2) we have that

$\frac{dN}{d\tau }=\frac{dS}{d\tau }+\frac{dI}{d\tau }+\frac{dR}{d\tau }=1-S-I-R=1-N.$

Since $N\left(0\right)=1$, then the solution of the above ordinary differential equation is $N\left(t\right)=1$. That is, $S\left(t\right)+I\left(t\right)+R\left(t\right)=1$ for all $t\ge 0$. □

Lemma 2.1 implies also that the solutions are bounded and, as a further consequence, we can consider $S=1-I-R$ to obtain the following simplified model:

$$\begin{gathered} \frac{dI}{d\tau }=I\left[f\right(R)(1-I-R)-k] \\ \frac{dR}{d\tau }=(k-1)I-R. \end{gathered}$$ (2.1)

The study of the equilibrium points of (1.2) will be done through the study of the simplified model (2.1). Hence, it will be relevant to consider the associated Jacobian matrix given by:

$J(I,R)=\left[\begin{array}{cc}f\left(R\right)(1-I-R)-k-If\left(R\right)& I\left[\frac{df}{dR}·(1-I-R)-f\left(R\right)\right]\\ (k-1)& -1\end{array}\right],$ (2.2)

valid when $\frac{df}{dR}$ is well defined. In fact, assuming some additional conditions on f, we can obtain the following useful result.

Lemma 2.2.Let f be a positive function, continuously differentiable on ℝ. The model given by (2.1) does not have non-constant periodic solutions with$0<I\left(t\right)$for all t.

Proof. Consider ϕ given by $\varphi (I,R)=\frac{1}{I}$ we have that

$\frac{\partial }{\partial I}\left(\varphi (I,R)I\left[f\right(R)(1-I-R)-k]\right)+\frac{\partial }{\partial R}\left(\varphi (I,R)\left[\right(k-1)I-R]\right)=-f\left(R\right)-\frac{1}{I}<0,$ (2.3)

if $I>0$. From the Bendixson-Dulac criterion, it follows that the system does not have a nonconstant periodic solution lying entirely in any simply connected region of the upper-plane $I>0$, so (2.1) does not have non-constant periodic solutions with $0<I\left(t\right)$ for all t. □

3 DISEASE-FREE EQUILIBRIUM

Note that $(I^{*},R^{*})=(0,0)$ is an equilibrium point of (2.1) corresponding to a disease-free state. The following results generalize well-known results related to the stability of the disease-free equilibrium in the classical SIR model ¹⁴, considering $\frac{f\left(R\right)}{k}$ as a variable reproduction number R ₀.

Lemma 3.3.Let f be a positive function, continuously differentiable on ℝ. If$\frac{f\left(0\right)}{k}<1$, then (0, 0) is a locally stable equilibrium point of (2.1).If$\frac{f\left(0\right)}{k}>1$, then (0, 0) is a local saddle equilibrium point.

Proof. The results follow from (2.2), since $J(0,0)=\left[\begin{array}{cc}f\left(0\right)-k& 0\\ k-1& -1\end{array}\right]$ has eigenvalues equal to −1 and $f\left(0\right)-k$. □

Lemma 3.4.Let f be a positive function, continuously differentiable on ℝ. If$\frac{f\left(0\right)}{k}<1$and (0, 0) is the only equilibrium point of the model given by (2.1), then (0, 0) is globally stable.

Proof. Consider $Z=\{0\le I\le 1;0\le R\le 1;I+R\le 1\}$, and X any open set on the plane such that $Z\subset X$. Because of Lemma 2.1, any solution of (2.1) with initial conditions $u^0=\left(I\right(0),R(0\left)\right)$ on Z, remains bounded and the ω-limit of u ₀, ω(u ⁰), satisfies $\omega \left(u^0\right)\subset Z\subset X$. Because we are considering that (0, 0) is the only equilibrium point of (2.1), from the Poincaré-Bendixson Theorem it follows that, for any initial condition $u^0=\left(I\right(0),R(0\left)\right)\in Z$ we have that

1. ω(u⁰) is a periodic orbit, or,

2. $(0,0)\in \omega \left(u^0\right)$ .

If ω(u ⁰) is a periodic orbit, Lemma 2.2 implies that (2.1) does not have periodic non-constant orbits with $I\left(t\right)>0$ for all t. Therefore, if ω(u ⁰) is a periodic orbit, then the orbit must intercept the axis $I=0$. Equations (2.1) imply that, in this case, I remains equal to zero and $R\to 0$ therefore $\left(I\right(t),R(t\left)\right)\to (0,0)$, then because (0, 0) is locally stable by Lemma 3.3, every solution that gets close enough, converges to (0, 0), so, in fact, in this case also $\left(I\right(t),R(t\left)\right)\to (0,0)$. □

4 ENDEMIC EQUILIBRIUM

4.1 Characterization and Existence

Now we consider the possibility of an endemic equilibrium point (I*, R*), so $I^{*}>0$. Note that if $I^{*}>0$, then any endemic equilibrium point of (2.1) must satisfy the following equations:

$f\left(R^{*}\right)(1-I^{*}-R^{*})-k=0\text{and}(k-1)I^{*}-R^{*}=0;$ (4.1)

which can be rewritten in terms of R* as

$f\left(R^{*}\right)=\frac{k-1}{\frac{k-1}{k}-R^{*}}\text{and}{I}^{*}=\frac{1}{k-1}{R}^{*}.$ (4.2)

If we define the auxiliary function g by

$g\left(R\right)=\frac{k-1}{\frac{k-1}{k}-R},$ (4.3)

it is clear from (4.2) that for the existence of endemic equilibrium, it is necessary that f and g intercept. In fact, it is possible to completely characterize the endemic equilibrium points of (2.1) in terms of functions f and g.

Theorem 4.1.Let f be a positive function, differentiable on [0,1] and g defined as in (4.3). A point (I ^* , R ^* ) is an endemic equilibrium of (2.1) if and only if$R^{*}\in \left(0,\frac{k-1}{k}\right),I^{*}\in \left(0,\frac{1}{k}\right),I^{*}=\frac{1}{k-1}{R}^{*}$and,$f\left(R^{*}\right)=g\left(R^{*}\right)$.

Proof. The results follow from the equivalence between Eqs. (4.1) and (4.2), the fact that g(R) is positive only if $R<\frac{k-1}{k}$ and that $R^{*}\in \left(0,\frac{k-1}{k}\right)$ if and only if $I^{*}=\frac{1}{k-1}{R}^{*}\in \left(0,\frac{1}{k}\right)$ because $k>1$. □

Theorem 4.1 establishes that endemic equilibrium points occur if and only if the functions f and g intercept each other on $\left(0,\frac{k-1}{k}\right)$. The next corollary establishes a simple condition to ensure that this interception will occur.

Theorem 4.2.Let f be a positive function on ℝ, differentiable [0;1]. If$f\left(R\right)>g\left(R\right)$for some$R\in [0,\frac{k-1}{k})$then (2.1) has at least one endemic equilibrium point (I ^* ;R ^* ), with$R^{*}\in \left(R,\frac{k-1}{k}\right)$and$I^{*}\in \left(\frac{R}{k-1},\frac{1}{k}\right)$In particular, if$f\left(0\right)>k$, there exists at least one endemic equilibrium.

Proof. Consider the function $h=f-g$. Note that, because f and g are continuous on $[0,\frac{k-1}{k})$, h is also continuous on $[0,\frac{k-1}{k})$. According to Theorem 4.1, to obtain the desired result, we must prove that h has at least one root on $(R,\frac{k-1}{k})$. If f$f\left(R\right)>g\left(R\right)$ for some $R\in [0,\frac{k-1}{k})$ then $h\left(R\right)>0$, and by the hypothesis on f and the definition of g we have

$\underset{R\to {\left(\frac{k-1}{k}\right)}^{-}}{\lim }h\left(R\right)=f\left(\frac{k-1}{k}\right)-\underset{R\to {\left(\frac{k-1}{k}\right)}^{-}}{\lim }g\left(R\right)=-\infty .$

From the Mean Value Theorem and the continuity of h, the previous statements imply that h has at least one root $R^{*}\in \left(R,\frac{k-1}{k}\right)$. By making $I^{*}=\frac{R^{*}}{k-1}$ we obtain the desired equilibrium point as (I ^∗ , R ^∗ ). The final statement follows from the fact that $g\left(0\right)=k$. □

4.2 Local Stability of Endemic Equilibrium

Theorems 4.1 and 4.2 establish conditions to verify the existence of endemic equilibrium points in terms of functions f and g. The next theorem shows that the relationship between the derivatives of f and g can be used to classify the local stability of the endemic equilibrium obtained.

Note that for all $R\ne \frac{k-1}{k}$, g satisfies $\frac{dg}{dR}=\frac{1}{k-1}{g}^2\left(R\right)$.

Theorem 4.3.Let f be a positive function, differentiable on [0, 1]; g defined as in (4.3); and (I ^* , R ^* ) an endemic equilibrium point of (2.1). If,

$\frac{df}{dR}\left(R^{*}\right)<\frac{dg}{dR}\left(R^{*}\right)\left(\mathit{\text{or}}<\frac{1}{k-1}{g}^2\left(R^{*}\right)\mathit{\text{or}}<\frac{1}{k-1}{f}^2\left(R^{*}\right)\right),$ (4.4)

then (I ^∗ , R ^∗ ) is a locally stable equilibrium point. If

$\frac{df}{dR}\left(R^{*}\right)>\frac{dg}{dR}\left(R^{*}\right)\left(\text{or}>\frac{1}{k-1}{g}^2\left(R^{*}\right)\text{or}>\frac{1}{k-1}{f}^2\left(R^{*}\right)\right),$ (4.5)

then (I ^∗ , R ^∗ ) is a locally saddle point.

Proof. If (I ^* , R ^* ) is an endemic equilibrium point of (2.1), then, from Theorem 4.1, we have that $f\left(R^{*}\right)=g\left(R^{*}\right)$. Because $\frac{dg}{dR}=\frac{1}{k-1}{g}^2\left(R\right)$, we have that $\frac{dg}{dR}\left(R^{*}\right)=\frac{1}{k-1}{g}^2\left(R^{*}\right)=\frac{1}{k-1}{f}^2\left(R^{*}\right)$.

Therefore, to obtain the desired result, we can use any of these three equivalent expressions. We will use $\frac{1}{k-1}{f}^2\left(R^{*}\right)$.

Since $k>0$, the first equation in (4.1) implies that $1-I^{*}-R^{*}\ne 0$ and $f\left(R^{*}\right)\ne 0$. Using (4.1), we have that the Jacobian matrix (2.2) evaluated on (I ^* , R ^* ) is given by:

$J(I^{*},R^{*})=\left[\begin{array}{cc}-I^{*}f\left(R^{*}\right)& I^{*}\left[\frac{df}{dR}\left(R^{*}\right)\frac{k}{f\left(R^{*}\right)}-f\left(R^{*}\right)\right]\\ (k-1)& -1\end{array}\right].$ (4.6)

By the Routh-Hurwitz criterion, in order to prove the local stability of (I ^* , R ^* ), it would be sufficient to show that the characteristic polynomial of the Jacobian matrix (4.6) has positive coefficients. Therefore, it would be sufficient to prove that:

$\text{(Trace):}{I}^{*}f\left(R^{*}\right)+1>0.$ (4.7)

$\text{(Determinant):}{I}^{*}f\left(R^{*}\right)-(k-1)I^{*}\left[\frac{df}{dR}\left(R^{*}\right)\frac{k}{f\left(R^{*}\right)}-f\left(R^{*}\right)\right]>0.$ (4.8)

Inequality (4.7) is satisfied because we are considering f as a positive function. As $I^{*}\ne 0$, the inequality given by (4.8) is equivalent to the following inequalities:

$$\begin{gathered} 0<f\left(R^{*}\right)-(k-1)\left[\frac{df}{dR}\left(R^{*}\right)\frac{k}{f\left(R^{*}\right)}-f\left(R^{*}\right)\right], \\ 0<kf\left(R^{*}\right)-(k-1)\frac{df}{dR}\left(R^{*}\right)\frac{k}{f\left(R^{*}\right)}, \\ (k-1)\frac{df}{dR}\left(R^{*}\right)\frac{k}{f\left(R^{*}\right)}<kf\left(R^{*}\right), \\ \frac{df}{dR}\left(R^{*}\right)<\frac{1}{k-1}{f}^2\left(R^{*}\right). \end{gathered}$$

Hence, if $\frac{df}{dR}\left(R^{*}\right)<\frac{1}{k-1}{f}^2\left(R^{*}\right)$, the characteristic polynomial of matrix (4.6) has only positive coefficients, which implies that both eigenvalues have a negative real part. So, (I ^* , R ^* ) is a locally stable equilibrium point.

Similarly, if $\frac{df}{dR}\left(R^{*}\right)>\frac{1}{k-1}{f}^2\left(R^{*}\right)$ then matrix (4.6) has a negative determinant, and its characteristic polynomial has the form ${\lambda }^2+b\lambda -c$ with $b,c>0$. This implies that matrix (4.6) has one positive and one negative real eigenvalue and, therefore, (I ^* , R ^* ) is a locally saddle point. □

Corollary 4.3.Let f be a positive function, differentiable on [0, 1] and let (I ^* , R ^* ) be an endemic equilibrium point of (2.1) such that$\frac{df}{dR}\left(R^{*}\right)\le 0$. Then (I ^* , R ^* ) is locally stable.

4.3 Multiple Endemic Equilibrium and local Multistability

Theorem 4.1 implies that multiple endemic equilibrium points occur if f and g have multiple interception points on $\left(0,\frac{k}{k-1}\right)$. The following result shows, that under some conditions, the existence of one endemic equilibrium implies the existence of another one.

Proposition 4.1.Let f be a positive function, differentiable on [0, 1]; g defined as in (4.3); and (I ^* , R ^* ) and endemic equilibrium of (2.1) such that$\frac{df}{dR}\left(R^{*}\right)\ne \frac{dg}{dR}\left(R^{*}\right)$. Then (I ^* , R ^* ) is locally stable or there exists another endemic equilibrium point$({\overline{I}}^{*},{\overline{R}}^{*})$with${\overline{R}}^{*}>R^{*}$.

Proof. Because (I ^* , R ^* ) is an endemic equilibrium of (2.1), we have $f\left(R^{*}\right)=g\left(R^{*}\right)$ and, because $\frac{df}{dR}\left(R^{*}\right)\ne \frac{dg}{dR}\left(R^{*}\right)$, then by Theorem 4.3 either is locally stable if $\frac{df}{dR}\left(R^{*}\right)<\frac{dg}{dR}\left(R^{*}\right)$ or, locally unstable if $\frac{df}{dR}\left(R^{*}\right)>\frac{dg}{dR}\left(R^{*}\right)$. In the last case, for values of R slightly bigger than R ^* , the function f will be greater than g so Theorem 4.2 implies the existence of at least one endemic equilibrium point $({\overline{I}}^{*},{\overline{R}}^{*})$ with ${\overline{R}}^{*}>R^{*}$. □

The following example illustrates a situation with multiple unstable and stable endemic equilibrium points.

Example 1.Let n be a fixed positive integer, $R_i^{*}=i\left(\frac{k-1}{k}\right)\frac{1}{2n}$for$i=0,1,\dots ,2n-1$and$g\left(R\right)=\frac{k-1}{\frac{k-1}{k}-R}$Let f be a positive and differentiable function on [0, 1] such that, $f\left(0\right)<k$and$f\left(R\right)=g\left(R\right)-\sin \left(2n\pi \frac{k}{k-1}R\right)$for$R\in [R_1^{*},R_{2n-1}^{*}]$(Figure 1). Note that

$f\left(R_i^{*}\right)=g\left(R_i^{*}\right)-\sin \left(i\pi \right)=g\left(R_i^{*}\right).$

Figure 1:
Multistability. Consider $n=5$ and $k=5$ in Example 1. Diamond markers correspond to saddle equilibrium points and circle markers correspond to locally stable equilibrium points. Functions f (dashed line) and g (solid line) are pictured on the Figure 1a. On Figure 1b, the partial phase plane $R\times I$ for (2.1) is pictured.

Therefore, Theorem 4.1 implies that, in this case, model (2.1) has at least 2n equilibrium points given by$\left(R_i^{*},\frac{1}{k-1}{R}_i^{*}\right)$. Furthermore,

$\frac{df}{dR}\left(R_i^{*}\right)=\frac{dg}{dR}\left(R_i^{*}\right)-\frac{2n\pi k}{k-1}\cos \left(i\pi \right).$

Hence, for$i=2,4,\dots ,2n-2$, we have$\frac{df}{dR}\left(R_i^{*}\right)=\frac{dg}{dR}\left(R_i^{*}\right)-\frac{2n\pi k}{k-1}<\frac{dg}{dR}\left(R_i^{*}\right)$and, from Theorem 4.3, these$n-1$equilibrium points are locally stable. On the other hand, for$i=1,3,\dots ,2n-1$, we have$\frac{df}{dR}\left(R_i^{*}\right)=\frac{dg}{dR}\left(R_i^{*}\right)+\frac{2n\pi k}{k-1}>\frac{dg}{dR}\left(R_i^{*}\right)$; which, by Theorem 4.3, implies that these n equilibrium points are locally unstable saddle points. Since$f\left(0\right)<k$, Lemma 3.3 implies the local stability of equilibrium point (0, 0). This alternation between stable and unstable saddle points is illustrated inFigure 1.

Note that in Example 1, the generalized variable reproduction number $\frac{f\left(R\right)}{k}$ may attain some values less than 1 and, in a similar fashion as in the backward bifurcation phenomenon ^{9), (25}, the disease-free equilibrium co-exists with several endemic locally stable equilibrium points. The multistability phenomenon, i.e., the coexistence of different stable equilibrium points for a given set of parameters, has been the focus of a lot of research in the applied dynamical systems community. In multistable systems, the asymptotic behavior depends crucially on the initial conditions. For an overview of instances of multistability across different areas, and an extensive list of references see ¹⁸.

4.4 Uniqueness of Endemic Equilibrium and Global Stability

In this final subsection, we consider the possibility of global stability for endemic equilibrium. Clearly, this is only possible if there exists only one locally stable endemic equilibrium. The next result establishes a sufficient condition to guarantee such a situation.

Proposition 4.2.Let f be a positive differentiable function on [0, 1] such that f is constant or f is non-increasing . If$\frac{f\left(0\right)}{k}<1$then (2.1) has no endemic equilibrium points. If$\frac{f\left(0\right)}{k}>1$then there exists a unique endemic equilibrium point for (2.1). Furthermore, the endemic equilibrium point is locally stable.

Proof.

Consider as in Theorem 4.2, that $h=f-g$ with $g\left(R\right)=\frac{k-1}{\frac{k-1}{k}-R}$. Note that because g is strictly increasing and f is constant or non-increasing, the function h is strictly decreasing. Note also that $g\left(0\right)=k$, so if $\frac{f\left(0\right)}{k}<1$ then $f\left(0\right)<g\left(0\right)$ and therefore $h\left(0\right)<0$. Because h is strictly decreasing, h has no roots. By Theorem 4.1, this means that there are no endemic equilibrium points for (2.1).

If $\frac{f\left(0\right)}{k}>1$, then by Theorem 4.2, the system (2.1) has at least one endemic equilibrium (I*, R*) and therefore h has at least one root $R^{*}\in \left(0,\frac{k-1}{k}\right)$. Since h is strictly decreasing, the root must be unique, so the endemic equilibrium is also unique. The local stability follows from Corollary 4.3. □

The conditions for f in Proposition 4.2, are not necessary for the uniqueness of a locally stable endemic equilibrium point, because it is possible to find a positive and strictly increasing function f that intercepts g only once, as illustrated in the following example.

Example 2.Let$f:[0,1]\to \mathrm{ℝ}$be given by$f\left(R\right)=k{R}^2+2k$. Note that f is a positive, differentiable, and strictly increasing function on [0, 1]. Furthermore,$f\left(R\right)=g\left(R\right)$has a unique solution$R^{*}\in \left[0,\frac{k-1}{k}\right)$, which corresponds to a locally stable endemic equilibrium point (Figure 2).

Figure 2:
Unique and stable endemic equilibrium point. Consider $k=5$ in Example 2. Functions f (dashed line) and g (solid line) are pictured on the Figure 2a. On Figure 2b, the partial phase plane $R\times I$ for (2.1) is pictured. The gray dot corresponds to the unique stable equilibrium point.

In the previous example, the partial phase-plane presented suggests that the endemic equilibrium is not only locally stable, but also globally stable for initial conditions with $I\left(0\right)>0$. This, in fact, is true as a consequence of the following theorem.

Theorem 4.4.Let f be a positive function and continuously differentiable on ℝ with$f\left(0\right)>k$. If model (2.1) has a unique endemic equilibrium point (I ^* , R ^* ) and$\frac{df\left(R^{*}\right)}{dR}\ne \frac{dg\left(R^{*}\right)}{dR}$, then (I ^* , R ^* ) is globally stable for$I\left(0\right)>0$.

Proof. Note that because we are considering a unique endemic equilibrium point (I ^* , R ^* ), f and g intercept only once at R ^* . Because $f\left(0\right)>k=g\left(0\right)$, the continuity of f and g imply that $f\left(R\right)>g\left(R\right)$ if $R<R^{*}$ and, if $R>R^{*}$, then $f\left(R\right)<g\left(R\right)$, otherwise Theorem 4.2 implies the existence of another interception point for some $R>R^{*}$. Hence,

$$\begin{gathered} \frac{df}{dR}\left(R^{*}\right)=\underset{h\to 0}{\lim }\frac{f(R^{*}+h)-f\left(R^{*}\right)}{h}=\underset{h\to 0}{\lim }\frac{f(R^{*}+h)-g\left(R^{*}\right)}{h} \\ \le \underset{h\to 0}{\lim }\frac{g(R^{*}+h)-g\left(R^{*}\right)}{h}=\frac{dg}{dR}\left(R^{*}\right). \end{gathered}$$

Since we are assuming that $\frac{df\left(R^{*}\right)}{dR}\ne \frac{dg\left(R^{*}\right)}{dR}$, the above inequality implies that $\frac{df\left(R^{*}\right)}{dR}<\frac{dg\left(R^{*}\right)}{dR}$. From Theorem 4.3, it follows that the unique endemic equilibrium point (I ^* , R ^* ) is locally stable.

We analyze now the global stability of (I ^* , R ^* ). Note first that, because $f\left(0\right)>k$, the Lemma 3.3 implies that the disease-free equilibrium (0, 0) is a saddle unstable point. We claim that the stable manifold associated with this disease-free saddle equilibrium point corresponds to the axis $I=0$. Note first that if $I\left(0\right)=0$, then Eqs. (2.1) implies that $I\left(t\right)=0$ for all $t>0$ and $R\left(t\right)=R\left(0\right)e^{-t}\to 0$ as $t\to \infty$. Note also that, from the continuity of f and the fact that $f\left(0\right)>k$, if I and R take small enough positive values, then from the first equation of (2.1) we have $\frac{dI}{d\tau }>0$. Therefore, if $\left(I\right(t),R(t\left)\right)\to (0,0)$, then I can not take values always strictly positive, which means that $I\left(t_0\right)=0$ for some t ₀. However, $(0,R(t_0\left)e^{-(t-t_0)}\right)$ is a solution passing through $\left(I\right(t_0),R(t_0\left)\right)=(0,R(t_0\left)\right)$.

Because the uniqueness of the solution of (2.1), it follows that $I\left(t\right)=0$ for all t. Hence, the stable manifold associated with (0, 0) is the axis $I=0$.

Consider now $Z=\{0\le I\le 1;0\le R\le 1;I+R\le 1\}$, and X any open set on the plane such that $Z\subset X$. Because Lemma 2.1, any solution of (2.1) with initial conditions $u^0=\left(I\right(0),R(0\left)\right)$ on Z remains bounded, and the ω−limit of u ₀, ω(u ⁰) satisfy $\omega \left(u^0\right)\subset Z\subset X$. From the Poincaré-Bendixson Theorem, one of the followings holds:

1. ω(u⁰) is a periodic orbit, or,

2. ω(u⁰) a connected set composed of a finite number of fixed points together with homoclinic and heteroclinic orbits connecting these, or,

3. ω(u⁰) consists of an equilibrium.

We claim that if $I\left(0\right)>0$ then $\omega \left(u^0\right)=(I^{*},R^{*})$. Assume first that $I\left(0\right)>0$ and ω(u ⁰) is a periodic orbit. By Lemma 2.2, the orbit ω(u ⁰) must intercept the axis $I=0$, but because this axis is the stable manifold of (0, 0), this implies that I remains equal to zero and $R\to 0$, so it is impossible for I to periodically return to $I\left(0\right)>0$. Hence, ω(u ⁰) can not be a periodic orbit.

We argue now that if $I\left(0\right)>0$, then the endemic equilibrium (I ^* , R ^* ) belongs to ω(u ⁰). If $(I^{*},R^{*})\notin \omega \left(u^0\right)$, then there are not heteroclinic orbits and, because the stable manifold of (0, 0) is the axis $I=0$, there are not homoclinic orbits either. Therefore, the only possibility is that $\omega \left(u^0\right)=\left\{\right(0,0\left)\right\}$. From the Bolzano-Weierstrass theorem and the compacity of Z, we would have that ${\lim }_{t\to \infty }u\left(t\right)={\lim }_{t\to \infty }\left(I\right(t),R(t\left)\right)=(0,0)$, but this implies that $I\left(0\right)=0$ which contradicts our hypothesis. Hence, if $I\left(0\right)>0$, then $(I^{*},R^{*})\in \omega \left(u^0\right)$.

Since the endemic equilibrium $(I^{*},R^{*})\in \omega \left(u^0\right)$ and is locally stable, every solution that gets close enough to it, converges to it. This implies that the disease-free equilibrium (0, 0) can not belong to ω(u ⁰), so we conclude that if $I\left(0\right)>0$ then ω(u ⁰) consists only of the endemic equilibrium (I ^* , R ^* ) and therefore all solutions with $I\left(0\right)>0$ converge to (I ^* , R ^* ). □

5 CONCLUSION

We considered a general SIR epidemiological model (1.2) with an infection rate f depending on the recovered population R. The main contribution of this paper is the determination of sufficient conditions for the existence, uniqueness (Theorems 4.1, 4.2 and Proposition 4.2), and stability of endemic equilibrium points (Theorems 4.3 and 4.4). These results are obtained in terms of f and the auxiliary function $g\left(R\right)=\frac{k-1}{\frac{k-1}{k}-R}$.

The auxiliary function g(R) can be considered a recovery-dependent threshold for the infection rate f(R). If this threshold is surpassed in some state (I, R), some of the consequences may be undesirable from an epidemiological point of view. Theorem 4.1 implies that when $f\left(R\right)>g\left(R\right)$, then there must exist an endemic equilibrium point and Proposition 4.1 implies that, in many situations, this equilibrium will be locally stable or it will lead to another endemic equilibrium point. On the other hand, the results in this paper imply that if one is able to control f to remain less than g for all R, then there are not endemic equilibrium points and the disease-free equilibrium will be globally stable.

The relationship between f and g generalizes the relationship between the constants R ₀ and 1 in the classical SIR model and shows that when considering variable parameters in epidemiological models, one could expect the appearance of variable thresholds relevant to the effective control of the diseases.

REFERENCES

• ¹
  M.E. Alexander & S.M. Moghadas. Bifurcation analysis of an SIRS epidemic model with generalized incidence. Journal on Applied Mathematics, 65(5) (2005), 1794-1816.
• ²
  C. Betsch. How behavioural science data helps mitigate the COVID-19 crisis. Nature Human Behaviour, (2020). doi:10.1038/s41562-020-0866-1.
  » https://doi.org/10.1038/s41562-020-0866-1
• ³
  M. Boily, G. Godin, M. Hogben, L. Sherr & F. Bastos. The impact of the transmission dynamics of the HIV/AIDS epidemic on sexual behaviour: A new hypothesis to explain recent increases in risk taking-behaviour among men who have sex with men. Medical Hypotheses, 65(2) (2005), 215-226.
• ⁴
  F. Brauer. Mathematical epidemiology: Past, present, and future. Infectious Disease Modelling, 2(2) (2017), 113-127.
• ⁵
  M. Chinazzi, J. Davis, M. Ajelli, C. Gioannini, M. Litvinova, S. Merler, A.P. Piontti, K. Mu, L. Rossi, K. Sun, C. Viboud, X. Xiong, H. Yu, M. Halloran, I. Longini & A. Vespignani. The effect of travel restrictions on the spread of the 2019 novel coronavirus (COVID-19) outbreak. Science, (2020). doi:10.1126/science.aba9757. URL https://science.sciencemag.org/content/early/2020/03/05/science.aba9757
  » https://doi.org/10.1126/science.aba9757» https://science.sciencemag.org/content/early/2020/03/05/science.aba9757
• ⁶
  R. Ferrer & W. Klein. Risk perceptions and health behavior. Curr Opin Psychol., 5 (2015), 85-89.
• ⁷
  D. Greenhalgh & R. Das. Modeling epidemics with variable contact rates. Theoretical population biology, 47(2) (1995), 129-179.
• ⁸
  S. Greenhalgh & T. Day. Time-varying and state-dependent recovery rates in epidemiological models. Infectious Disease Modelling , 2(4) (2017), 419-430.
• ⁹
  A. Gumel. Causes of backward bifurcations in some epidemiological models. Journal of Mathematical Analysis and Applications, 395(1) (2012), 355-365.
• ¹⁰
  M. Kraemer, C. Yang, B. Gutierrez, C. Wu, B. Klein, D. Pigott, L. Plessis, N. Faria, R. Li, W. Hanage, J. Brownstein, M. Layan, A. Vespignani , H. Tian, C. Dye, O.G. Pybus & S. Scarpino. The effect of human mobility and control measures on the COVID-19 epidemic in China. Science , (2020). doi:10.1126/science.abb4218. URL https://science.sciencemag.org/content/early/2020/03/25/science.abb4218
  » https://doi.org/10.1126/science.abb4218» https://science.sciencemag.org/content/early/2020/03/25/science.abb4218
• ¹¹
  I. Levy, Z. Mor, E. Anis, S. Maayan, E. Leshem, S. Pollack, M. Chowers, O. Mor, K. Riesenberg, Z. Sthoeger, D. Ram & Z. Grossman. Men Who Have Sex With Men, Risk Behavior, and HIV Infection: Integrative Analysis of Clinical, Epidemiological, and Laboratory Databases. Clinical Infectious Diseases, 52(11) (2011), 1363-1370.
• ¹²
  X. Liu & P. Stechlinski. Infectious disease models with time-varying parameters and general nonlinear incidence rate. Applied Mathematical Modelling, 36(5) (2012), 1974-1994.
• ¹³
  P. Manfredi & A. d’Onofrio. “Modeling the Interplay Between Human Behavior and the Spread of Infectious Diseases”. Springer (2013).
• ¹⁴
  M. Martcheva. “An introduction to mathematical epidemiology”, volume 61. Springer (2015).
• ¹⁵
  P. O’Neill. An epidemic model with removal-dependent infection rate. The Annals of Applied Probability, 7(1) (1997), 90-109.
• ¹⁶
  T. Parker-Pope. What You Can do About Coronavirus Right Now (2020). URL URL https://www. nytimes.com/interactive/2020/world/coronavirus-preparation-preparedness.html Last accessed on March 28, 2020.
  » https://www. nytimes.com/interactive/2020/world/coronavirus-preparation-preparedness.html
• ¹⁷
  G.P. Pisano, R. Sadun & M. Zanini. Lessons from Italy’s Response to Coronavirus (2020). URL URL https://hbr.org/2020/03/lessons-from-italys-response-to-coronavirus Last accessed on March 28, 2020.
  » https://hbr.org/2020/03/lessons-from-italys-response-to-coronavirus
• ¹⁸
  A. Pisarchik & U. Feudel. Control of multistability. Physics Reports, 540(4) (2014), 167 - 218. doi: https://doi.org/10.1016/j.physrep.2014.02.007. URL http://www.sciencedirect.com/science/ article/pii/S0370157314000453
  » https://doi.org/https://doi.org/10.1016/j.physrep.2014.02.007» http://www.sciencedirect.com/science/ article/pii/S0370157314000453
• ¹⁹
  G. Pulido, I. Barradas & B. Luna. Backward bifurcation for some general recovery functions. Mathematical Methods in the Applied Sciences, 40(5) (2017), 1505-1515.
• ²⁰
  A. Remuzzi & G. Remuzzi. COVID-19 and Italy: what next? The Lancet, (2020). doi:10.1016/ S0140-6736(20)30627-9.
  » https://doi.org/10.1016/ S0140-6736(20)30627-9
• ²¹
  M. Roberts, V. Andreasen, A. Lloyd & L. Pellis. Nine challenges for deterministic epidemic models. Epidemics, 10 (2015), 49-53.
• ²²
  H. Thieme & J. Yang. An endemic model with variable re-infection rate and applications to influenza. Mathematical biosciences, 180(1-2) (2002), 207-235.
• ²³
  WHO. Measles cases spike globally due to gaps in vaccination coverage (2018).URL URL https://www.who.int/news-room/detail/ 29-11-2018-measles-cases-spike-globally-due-to-gaps-in-vaccination-coverage Last accessed on March 03, 2020.
  » https://www.who.int/news-room/detail/ 29-11-2018-measles-cases-spike-globally-due-to-gaps-in-vaccination-coverage
• ²⁴
  J. Zaidi, E. Grapsa, F. Tanser, M. Newell & T. Barnighausen. Dramatic increase in HIV prevalence after scale-up of antiretroviral treatment. AIDS, 27(14) (2013), 2301-2305.
• ²⁵
  F. Zhang, T. Zhao, H. Liu & Y. Chen. Backward bifurcation in a stage-structured epidemic model. Applied Mathematics Letters, 89 (2019), 85-90.

Publication Dates

• Publication in this collection
  30 Nov 2020
• Date of issue
  Sep-Dec 2020

History

• Received
  08 Apr 2020
• Accepted
  22 June 2020