SIR Model with Dependent Infectivity and Death Rates

1. Introduction

We develop, analyze and simulate a new extended version of the SIR model for the spread of a communicable disease, in which both infectivity and death rates are system variables. The motivation for including system dependent infectivity and death rates comes from observations on the evolution of the COVID-19 epidemic. It seems that by allowing for such expansion of the basic SIR model, a more realistic and more accurate description is provided for certain long lasting diseases. The model here is generic and needs to be updated with infectivity and death rate equations that are tailored to each specific disease.

This work expands on the work “A modified SIR model with dependent infectivity and hospitalization rates” [1]. That paper constructs, analyzes and simulates the SIR-IH model, a modified SIR model in which the infectivity rate and hospitalization ratio are system variables. Using this idea, our model has system-dependent infectivity and death rates. This concept was already considered in “A mathematical model of the COVID-19 pandemic dynamics with dependent variable infection rate: Application to the Republic of Korea,” [2].

A recent work, [3], constructed an epidemiological model for cancer in which, to account for the adaptive nature of the tumor’s PD-L1 expression, the authors introduced an additional ordinary differential equation governing the dynamics of the state variable $\epsilon$, which describes the rate of change in the tumor PD-L1 expression level over time.

A basic and comprehensive review of the SIR-type compartmental models can be found in the survey by Heathcote, “The Mathematics of Infectious Diseases,” [4]. We note that the literature on various aspects of SIR models is very large, and in addition to [4], many references can be found in [10,11] and in Wikipedia [5].

Here, we extend the SIR model into the SIR-ID model by adding the infectivity and death rates as system variables. It is known, see, e.g., [6,7,8] and the references therein, that the death rate for COVID-19 did vary over time, even for a single variant. Factors such as vaccination, new treatments (antiviral and steroids), and post-infection immunity evolved alongside the virus, altering the infectivity rate and death rate of infected individuals. These rates reflect the complex interaction between a variant’s virulence and the body’s and health responses. For example, vaccination significantly reduced severe outcomes and deaths, even for highly transmissible variants such as the Omicron [8]. Development and deployment of antivirals (such as Paxlovid) lowered death rates [7]. Repeated infections built up population-level immunity, making subsequent waves less deadly [9]. These are just examples of many factors that affected the death and infectivity rates for COVID-19. These observations apply to most other contagious diseases, and similar factors should affect the infectivity rate, which motivate our general model.

Following this introduction, we derive the model in detail in Section 2. There, the compartmental flow chart and the model equations can be found. The model consists of five rate equations, the basic equation of the SIR model for the Susceptibles, Infected and Recovered, and the addition of two equations for the death rate due to the disease and the infectivity rate. As noted above, the novelty lies in the introduction of the differential equations for the disease related death rate $\tilde{\alpha }$ and the infectivity rate $\tilde{\beta }$, which are treated as unknowns. However, these variables are constrained to avoid unreasonable values that the differential equations may produce. Although each specific disease has its own dynamics, here for simplicity, we assume simple structures for the rate equations for these rates. Related explanations can be found in Remarks 2.1 and 2.3.

The positivity and boundedness of the solutions, which makes them biologically viable, and the global existence of the solution, is described in Section 3. The steady states of the system and their stability are studied in Section 4.

Section 5 describes a simple explicit finite difference algorithm for the computer simulations of the model and shows two typical solutions. The first depicts the Disease Free Equilibrium (DFE), with system constants that make it stable and attracting. The second simulation is in the case when the DFE is unstable, and then the Endemic Equilibrium (EE) is stable and attracting. Furthermore, the simulations show how the system approaches to these steady states. To obtain confidence in the simulations, Section 6 presents a numerical argument for the convergence of the numerical method. The conclusions of the work and some ideas for future study can be found in Section 7.

2. The Model

We present a general SIR model of disease transmission, modified so it includes dependent rates of infectivity and death. We use the terminology of a generic disease in a large city in which the disease is spreading. Modifications to the model for a specific disease and location need to be tailored to the disease. We follow, in part, the presentation and ideas in [1].

We deal with whole populations, do not consider their spatial spread, and assume that the populations are large enough to justify the use of ordinary differential equations (ODEs). As noted in the introduction, we describe the dynamics of three populations consisting of individuals: susceptible (S), infected (I), and recovered (R). These are functions of time, and the rate of change of each population is described in the model. For simplicity, we consider the time unit to be a day. The novelty in this work is the assumption that the disease-induced death rate function $\left(\alpha =\alpha \left(t\right)\right)$ and the infectivity rate function $\left(\beta =\beta \left(t\right)\right)$ are dependent variables. The infectivity rate function $\beta$ is the probability that one susceptible person who comes into contact with one infected becomes sick. This measures how rapidly the disease may spread.

We assume that the death rate function $\alpha$ and the infectivity rate function $\beta$ are only mildly related. We use a differential equation for each one that depends on the infected and recovered populations. Moreover, we assume that the fraction $\rho$ of those who have recovered can lose their immunity and be reinfected.

The compartmental structure of the model, see e.g., [5], is depicted in Figure 1.

The susceptible population $S=S\left(t\right)$ consists of all individuals who are healthy and can become sick. The infected population $I=I\left(t\right)$ consists of individuals infected with the disease. Individuals in $I\left(t\right)$ may or may not present clinical symptoms of the disease. The population $R=R\left(t\right)$ denotes those who recovered from the disease.

We denote by $P\left(t\right)$ the net number of susceptible individuals that are added to the city per day, either by birth or outside arrivals. We let $\mu (1/\mathrm{day})$ denote the death rate in the absence of the disease, and then the death rate for the infected is $\mu +\alpha (1/\mathrm{day})$.

The disease recovery constant is $\sigma$ and the number of recovered (per day) is $\sigma I$. We let $\rho \in (0,1]$ denote the fraction of recovered that become susceptible again. For the sake of simplicity, we assume that $\sigma$ and $\rho$ are constants. It is straightforward to make them dependent variables, however, at the expense of making the model more complex.

Remark 1.

The equations for the changes in α and β are likely to be of the general form

$${\displaystyle \frac{d\alpha }{dt}}={\mathrm{\Psi }}_{\alpha }(I,R,\alpha ,\beta ),{\displaystyle \frac{d\beta }{dt}}={\mathrm{\Psi }}_{\beta }(I,R,\alpha ,\beta ),$$

(1)

for given functions ${\mathrm{\Psi }}_{\alpha }$ and ${\mathrm{\Psi }}_{\beta }$, which need to be established for each disease, based on field data. These are likely to be set inclusions, which we will comment on in Remark 2.2. They are also likely to depend on the disease history.

In this work, for simplicity and lack of specific information, we assume that

$${\mathrm{\Psi }}_{\alpha }(I,R,\alpha ,\beta )={\delta }^{+}I-{\delta }^{-}\alpha R,{\mathrm{\Psi }}_{\beta }(I,R,\alpha ,\beta )={\gamma }^{+}I-{\gamma }^{-}\beta R.$$

Here, ${\delta }^{+},{\delta }^{-},{\gamma }^{+},{\gamma }^{-}$ are positive constants, to be determined from field data. The dependence on R is due to population immunity from prior infections. The dependence on I seems to be natural.

Since $\alpha$ and $\beta$ are dependent variables, we need to make sure that their values remain within reasonable intervals, no matter what the populations are. To that end, let $0<{\alpha }_{-}<{\alpha }_{+}\le 1$ and $0<{\beta }_{-}<{\beta }_{+}\le 1$ be the limits for the allowed values of $\alpha$ and $\beta$, respectively, and we want to guarantee that for $0\le t\le T$,

$$\alpha \left(t\right)\in [{\alpha }_{-},{\alpha }_{+}],\beta \left(t\right)\in [{\beta }_{-},{\beta }_{+}].$$

(2)

To that end, we introduce the restriction or cutoff operators ${\mathrm{\Phi }}_{\alpha },{\mathrm{\Phi }}_{\beta }$, acting on $s\in \mathbb{R}$, as follows,

$${\mathrm{\Phi }}_{\alpha }\left(s\right)=\left\{\begin{array}{cc}{\alpha }_{-},\hfill & s\le {\alpha }_{-},\hfill \\ s,\hfill & {\alpha }_{-}<s<{\alpha }_{+},\hfill \\ {\alpha }_{+}\hfill & s\ge {\alpha }_{+}.\hfill \end{array}\right.{\mathrm{\Phi }}_{\beta }\left(s\right)=\left\{\begin{array}{cc}{\beta }_{-},\hfill & s\le {\beta }_{-},\hfill \\ s,\hfill & {\beta }_{-}<s<{\beta }_{+},\hfill \\ {\beta }_{+}\hfill & s\ge {\beta }_{+},\hfill \end{array}\right.$$

(3)

Both operators are bounded, piecewise linear and Lipschitz continuous, and are used in the model below.

The total living population, at time t, is given by

$$N\left(t\right)=S\left(t\right)+I\left(t\right)+R\left(t\right),$$

(4)

and is not assumed to be constant.

These assumptions lead to the following mathematical model for SIR-ID, a modified SIR disease dynamical system with system-dependent infectivity and death rates.

Model 1.

SIR-ID Find five functions $(S,I,R,\tilde{\alpha },\tilde{\beta })$, defined on $[0,T]$, such that the following equations hold, for $0<t\le T$:

$$\begin{array}{cc}\hfill {\displaystyle \frac{dS}{dt}}& =P-\beta I{\displaystyle \frac{S}{N}}-\mu S+\rho R,\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill {\displaystyle \frac{dI}{dt}}& =\beta I{\displaystyle \frac{S}{N}}-(\mu +\alpha +\sigma )I,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\displaystyle \frac{dR}{dt}}& =\sigma I-\mu R-\rho R,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d\tilde{\alpha }}{dt}}& ={\delta }^{+}I-{\delta }^{-}\tilde{\alpha }R\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d\tilde{\beta }}{dt}}& ={\gamma }^{+}I-{\gamma }^{-}\tilde{\beta }R,\hfill \end{array}$$

(9)

where

$$\alpha \left(t\right)={\mathrm{\Phi }}_{\alpha }\left(\widetilde{\alpha \left(t\right)}\right),\beta \left(t\right)={\mathrm{\Phi }}_{\beta }\left(\widetilde{\beta \left(t\right)}\right)$$

(10)

and N is given in (4). The initial conditions are

$$S\left(0\right)=S_0,I\left(0\right)=I_0,R\left(0\right)=R_{init},$$

$$\alpha \left(0\right)={\alpha }_0,\beta \left(0\right)={\beta }_0.$$

(11)

Here, $\tilde{\alpha }$ and $\tilde{\beta }$ are the unrestricted auxiliary dependent variables. Everywhere below we refer to the restricted variables $\alpha$ and $\beta$ as the disease-induced death and infectivity rates, respectively. We denote by $R_{init}$ the initial value of R, since in mathematical epidemiology $R_0$ usually denotes the basic stability number, or the basic reproduction number, see [4].

Furthermore, $S_0>0$ is the initial population before the outbreak of the disease, and $I_0\ge 0$, $R_{init}\ge 0$ are the initial populations at $t=0$. Additionally, $\alpha \left(0\right)={\alpha }_0\in [{\alpha }_{-},{\alpha }_{+}]$ is the initial disease-induced death rate and $\beta \left(0\right)={\beta }_0\in [{\beta }_{-},{\beta }_{+}]$ is the initial infectivity rate.

In practice, one typically assumes that initially $S_0=N\left(0\right)>0$, so that at the outbreak the population is only of susceptibles, and all the other populations vanish. However, for the sake of generality, we allow for them to be nonnegative but satisfying $N\left(0\right)=S_0+I_0+R_{init}$, so that (4) holds at $t=0$.

Next, we discuss equations (8) and (9). We assume that the infectivity of the disease increases with the number of infected, I, as the disease vector has more options to grow and possibly mutate, and decreases with the number of recovered, as more individuals are immune. As noted above, we assume that the relationship is linear, which is very likely a severe restriction. It may be of interest to study the question of the form of ${\mathrm{\Phi }}_{\alpha }$ and ${\mathrm{\Phi }}_{\beta }$ in the future, as it deserves a thorough investigation.

We note that (8) and (9) can be integrated, thus,

$$\widetilde{\alpha \left(t\right)}={\alpha }_0+{\int }_0^t\left({\delta }^{+}I\left(s\right)-{\delta }^{-}\widetilde{\alpha \left(s\right)}R\left(s\right)\right)ds,$$

$$\widetilde{\beta \left(t\right)}={\beta }_0+{\int }_0^t\left({\gamma }^{+}I\left(s\right)-{\gamma }^{-}\widetilde{\beta \left(s\right)}R\left(s\right)\right)ds,$$

which makes the system into a history dependent system of integral equations, however, we do not pursue it here.

The following remark can be found in [1].

Remark 2.

Instead of using the restriction operators ${\mathrm{\Phi }}_{\alpha },{\mathrm{\Phi }}_{\beta }$, we may use the notion of the subdifferential. To that end, we let $I_{\alpha },I_{\beta }$ be the indicator functions of $[{\alpha }_{-},{\alpha }_{+}]$ and $[{\beta }_{-},{\beta }_{+}]$, that is, $I_{\alpha }\left(\alpha \right)=0$ if $\alpha \in [{\alpha }_{-},{\alpha }_{+}]$ and $I_{\alpha }=+\infty$ otherwise, and similarly, $I_{\beta }\left(\beta \right)=0$ if $\beta \in [{\beta }_{-},{\beta }_{+}]$ and $I_{\beta }=+\infty$ otherwise. Let $\partial I_{\alpha },\partial I_{\beta }$ be the subdifferentials,

$$\partial I_{\alpha }\left(s\right)=\left\{\begin{array}{cc}(-\infty ,0],\hfill & s={\alpha }_{-},\hfill \\ 0\hfill & {\alpha }_{-}<s<{\alpha }_{+},\hfill \\ \infty )\hfill & s={\alpha }_{+},\hfill \end{array}\right.$$

and similarly for $\partial I_{\beta }$, with β replacing α.

Then, instead of using the restriction operators, we may use the set-valued inclusions for α and β, which are

$${\displaystyle \frac{d\alpha }{dt}}\in {\mathrm{\Psi }}_{\alpha }(I,R,\alpha ,\beta )-\partial I_{\alpha }\left(\alpha \right),{\displaystyle \frac{d\beta }{dt}}\in {\mathrm{\Psi }}_{\beta }(I,R,\alpha ,\beta )-\partial I_{\beta }\left(\beta \right).$$

The subdifferentials guarantee that the variables cannot leave the designated intervals.

• A summary of the definitions of the model parameters is given in Table 1.

3. Positivity, Boundedness and Existence

This section establishes the non-negativity, boundedness, and existence of global solutions to the model. These are important factors in ensuring that the model makes sense both biologically and mathematically. We note that the local existence is guaranteed by the general theory of dynamical systems, since the system is Lipschitz continuous.

Outside of this section, we take the recruitment rate of susceptible individuals (denoted P) to be constant with respect to time. However, for theoretical purposes, we allow P to be a bounded and nonnegative function of time with upper bound $P_{\max }$.

Theorem 2.

(Positivity and boundedness) Let $S_0>0,I_0\ge 0,$ and $R_{init}\ge 0$. If $z\left(t\right)=\left(S\right(t),I(t),R(t\left)\right)$ is any solution to the system, then all of the components of z are nonnegative and uniformly bounded on each interval [0,T] on which they exist, for $0<T\le \infty$. Then, for $0\le t\le T$,

$$0\le S\left(t\right),I\left(t\right),R\left(t\right)\le N_0+{\displaystyle \frac{P_{max}}{\mu }}$$

where $0\le P\left(t\right)\le P_{\max }$.

Proof. 

We first assume, in addition, that $I_0>0$ and $R_{init}>0$, and because of the continuity of the solutions, there exists a time interval $[0,T^{*}]$, with $T^{*}>0$, such that $S\left(t\right),I\left(t\right),R\left(t\right)>0$ for $0\le t<T^{*}$. Let $0<t<T<T^{*}$. Then, since $0<I,N$ and $I/N<1$, it follows from (5) and (2) that

$${\displaystyle \frac{dS}{dt}}=P-{\displaystyle \frac{\beta IS}{N}}-\mu S+\rho R\ge -\left({\displaystyle \frac{\beta I}{N}}+\mu \right)S\ge -({\beta }_{+}+\mu )S.$$

Using the separation of variables and integration on $[0,t]$ together with the initial conditions, we obtain

$$S\left(t\right)\ge S_0{e}^{-({\beta }_{+}+\mu )t}>0,$$

since $S_0>0$. Thus, the solution is positive for all $0<t\le T$. Similarly, from (5) and (6) we have

$$I\left(t\right)\ge I_0{e}^{-(\mu +{\alpha }_{+}+\sigma )t}\ge 0,$$

$$R\left(t\right)\ge R_{init}{e}^{-(\mu +\rho )t}\ge 0.$$

Consequently, the solution is nonnegative on $[0,T]$, for all T for which the solution exists. We note that the last two estimates are positive, unless $I_0=R_{init}=0$. These inequalities allows us to remove the restriction $T<T^{*}$ and, also using continuity, we may assume that $I_0,R_{init}\ge 0$. This means that the system will be nonnegative on every interval on which it exists. Moreover, we have $S\left(t\right)+I\left(t\right)+R\left(t\right)=N\left(t\right)>0$.

Next, we show the boundedness of the solutions. By adding equations (5), (6), and (7) and using (4), we find

$${\displaystyle \frac{dN}{dt}}=P-\mu N-\alpha I.$$

Hence,

$${\displaystyle \frac{dN}{dt}}\le P_{\max }-\mu N.$$

Next, we rearrange and and use the usual manipulations to get:

$$\begin{array}{c}\hfill 0<N\left(t\right)\le N_0+{\displaystyle \frac{P_{max}}{\mu }},\end{array}$$

showing $N\left(t\right)$ is bounded from above. Since N is bounded and each of its components (S, I, R) are nonnegative, we have that all of these components are bounded as well. This completes the proof of boundedness of the populations.

Furthermore, because of the cutoff functions ${\mathrm{\Phi }}_{\alpha }$ and ${\mathrm{\Phi }}_{\beta }$, we know that $\alpha ,\beta$ are nonnegative and bounded. However, for the sake of completeness, we show that $\widetilde{\alpha \left(t\right)}$ and $\widetilde{\beta \left(t\right)}$ are also bounded and positive. First, we have the initial conditions of ${\tilde{\alpha }}_0,{\tilde{\beta }}_0>0$. Using equation (8) and the upper bound for R, we obtain

$$\begin{array}{cc}\hfill {\displaystyle \frac{d\tilde{\alpha }}{dt}}& \ge -{\delta }^{-}\tilde{\alpha }R\ge -{\delta }^{-}\tilde{\alpha }(N_0+P_{max}/\mu ).\hfill \end{array}$$

Again, using separation of variables, integration on $[0,t]$, and initial conditions, we have

$$\tilde{\alpha }\left(t\right)\ge {\tilde{\alpha }}_0{e}^{-{\delta }^{-}(N_0+P_{max}/\mu )t}>0.$$

Thus, $\tilde{\alpha }\left(t\right)>0$ on the interval $[0,t]$, since ${\tilde{\alpha }}_0>0$, and we extend this to $[0,T]$. Now, for an upper bound, we take the same equation (8) and find

$${\displaystyle \frac{d\tilde{\alpha }}{dt}}\le {\delta }^{+}I\le {\delta }^{+}(N_0+P_{max}/\mu ).$$

Again using separation of variables and integrating both sides on $[0,t]$ gives, for $0\le t\le T$,

$$\tilde{\alpha }\left(t\right)\le {\tilde{\alpha }}_0+{\delta }^{+}(N_0+P_{max}/\mu )T.$$

So, we have $\tilde{\alpha }\left(t\right)$ bounded from above. We follow similar steps for $\tilde{\beta }\left(t\right)$ and get the following two inequalities:

$$\begin{array}{cc}\hfill 0\le \tilde{\alpha }\left(t\right)& \le {\tilde{\alpha }}_0+{\delta }^{+}(N_0+P_{max}/\mu )T,\hfill \\ \hfill 0\le \tilde{\beta }\left(t\right)& \le {\tilde{\beta }}_0+{\gamma }^{+}(N_0+P_{max}/\mu )T.\hfill \end{array}$$

Thus, these variables are also nonnegative and bounded on every finite time interval $[0,T]$.    □

In addition to proving that the solutions are positive and bounded, we need to establish their global existence. Using the estimates above, having bounded parameters, and knowing that the right-hand sides of the system are Lipschitz, we infer the existence of a unique solution from standard results for dynamical system (see [12]).

Let $x=x\left(t\right)={(S,I,R,\tilde{\alpha },\tilde{\beta })}^T$, then the model can be written in the condensed form

$${\displaystyle \frac{d\mathit{x}}{dt}}=\mathit{F}(\mathit{x},t),$$

(12)

where

$$\mathit{F}=\left[\begin{array}{c}P-\beta I{\displaystyle \frac{S}{N}}-\mu S+\rho R\\ \beta I{\displaystyle \frac{S}{N}}-(\mu +\alpha +\sigma )I\\ \sigma I-\mu R-\rho R\\ {\delta }^{+}I-{\delta }^{-}\tilde{\alpha }R\\ {\gamma }^{+}I-{\gamma }^{-}\tilde{\beta }R\end{array}\right].$$

(13)

Here, $\alpha \left(t\right)={\mathrm{\Phi }}_{\alpha }\left(\widetilde{\alpha \left(t\right)}\right),\beta \left(t\right)={\mathrm{\Phi }}_{\beta }\left(\widetilde{\beta \left(t\right)}\right)$. It follows from the estimates in Theorem 1 that, for $0<T<\infty$,

$$\left|\right|\mathit{x}\left(t\right)\left|\right|\le L(N_0,P,T),$$

where L depends only on $N_0,P,T$, and the system parameters.

Remark 3.

We note that the bound $L=L(N_0,P,T)$ depends on T via the dependence of $\tilde{\alpha }$ and $\tilde{\beta }$ on T. This restricts the existence of the solution $x\left(t\right)$ to every finite time interval $[0,T]$ for $T<\infty$. However, even when $\tilde{\alpha }$ or $\tilde{\beta }$ tend to infinity, α and β remain bounded and, therefore, the solution $\left(S\right(t),I(t),R(t\left)\right)$ exists on the interval $[0,\infty )$.

• We summarize our results as follows.

Theorem 3.

(Existence and uniqueness) Let $S_0>0,I_0\ge 0$, $R_{init}\ge 0$,  ${\alpha }_0\in [{\alpha }_{-},{\alpha }_{+}]$, ${\beta }_0\in [{\beta }_{-},{\beta }_{+}]$, and let all the coefficients be positive constants. Then, there exists a unique solution to Model 1 on every time interval $[0,T]$ for $T<\infty$. Furthermore, the solution $\left(S\right(t),I(t),R(t\left)\right)$ exists on the interval $[0,\infty )$.

Proof. 

It is straightforward to show that given bounded coefficients and functions (see Theorem 1), the system is locally Lipschitz continuous. Then, it follows from the estimates above that it is Lipschitz on every finite time interval. Therefore, the result follows from standard results for dynamical systems (see e.g., [12]).    □

4. Steady States

This section studies the equilibrium points of the system and their stability. We have the Disease-Free Equilibrium (DFE) and the Endemic Equilibrium (EE). We assume that $P>0$ is constant, and denote the steady states as: $\overline{S},\overline{I},\overline{R},\overline{\alpha },\overline{\beta }$. To that end, we compute the Jacobian and then the eigenvalues of the system, thus establishing the stability of the system at a given state.

To compute the Jacobian, we find in turns the partial derivatives with respect to S, I, and R in equations (5), (6), and (7). We note that since $N=S+I+R$:

$${\displaystyle \frac{\partial (S/N)}{\partial S}}={\displaystyle \frac{I+R}{N^2}},{\displaystyle \frac{\partial (I/N)}{\partial I}}={\displaystyle \frac{S+R}{N^2}},{\displaystyle \frac{1/N}{\partial R}}=-{\displaystyle \frac{1}{N^2}}.$$

The Jacobian of the populations is the following.

$$J(S,I,R)=\left(\begin{array}{ccc}-\beta {\displaystyle \frac{I(I+R)}{N^2}}-\mu & -\beta {\displaystyle \frac{S(S+R)}{N^2}}& \beta {\displaystyle \frac{IS}{N^2}}+\rho \\ \beta {\displaystyle \frac{I(I+R)}{N^2}}& \beta {\displaystyle \frac{S(S+R)}{N^2}}-(\mu +\alpha +\sigma )& -\beta {\displaystyle \frac{IS}{N^2}}\\ 0& \sigma & -\mu -\rho \end{array}\right)$$

Next, we consider the equilibrium points and their stability.

4.1. Disease Free Equilibrium

We begin with the DFE, which is $(\overline{S},0,0)$, where $\overline{S}=\overline{N}=P/\mu$. We evaluate the Jacobian at the DFE,

$$J(\overline{S},0,0)=\left(\begin{array}{ccc}-\mu & -\overline{\beta }& \rho \\ 0& \overline{\beta }-(\mu +\overline{\alpha }+\sigma )& 0\\ 0& \sigma & -\mu -\rho \end{array}\right).$$

The eigenvalues are

$${\lambda }_1=-\mu ,{\lambda }_2=-\mu -\rho ,{\lambda }_3=\overline{\beta }-(\mu +\overline{\alpha }+\sigma ).$$

In order for the disease-free state to be stable and attracting, the eigenvalues must be negative. We conclude that:

Proposition 1.

The DFE, $(\overline{S}=P/\mu ,0,0)$, is stable and attracting (asymptotically stable) when

$$\overline{\beta }<(\mu +\overline{\alpha }+\sigma ).$$

(14)

It is unstable when

$$\overline{\beta }>(\mu +\overline{\alpha }+\sigma ).$$

(15)

• It follows, that somewhat artificially, we may write the basic stability number as

  $$R_0={\displaystyle \frac{\overline{\beta }}{\mu +\overline{\alpha }+\sigma }}$$

  Then, when $R_0<1$ the DFE is stable and attracting, and when $R_0>1$ it is unstable. As can be seen below, when the DFE is unstable the EE appears and is stable and attracting.

First, we look at a case such that a stable and attracting DFE is attained. Starting with initial conditions $S_0=40,000,I_0=2,000,$ and $R_{init}=200$, numerical simulations produced the following steady state values:

$$\overline{S}=41,926,\overline{I}=0,\overline{R}=0,\overline{\alpha }=0.0204,\overline{\beta }=0.1618,$$

and the eigenvalues are

$${\lambda }_1=-3.4\times {10}^{-5},{\lambda }_2=-0.7500,{\lambda }_3=-0.359.$$

By Proposition 4.1, we have that ${\lambda }_3<0,$ therefore the DFE is stable and attracting.

4.2. Endemic Equilibrium

We find the conditions for the existence of an endemic equilibrium state, and denote its populations by $\overline{z}=(\overline{S},\overline{I},\overline{R})$, where $\overline{I}>0$. Then, $\overline{z}$ is the solution of the system,

$$\begin{array}{cc}\hfill 0& =P-\left(\overline{\beta }{\displaystyle \frac{\overline{I}}{\overline{N}}}+\mu \right)\overline{S}+\rho \overline{R},\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill 0& =\overline{\beta }\overline{I}{\displaystyle \frac{\overline{S}}{\overline{N}}}-(\mu +\overline{\alpha }+\sigma )\overline{I},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill 0& =\sigma \overline{I}-\left(\mu +\rho \right)\overline{R}.\hfill \end{array}$$

(18)

Here, $\overline{N}=\overline{S}+\overline{I}+\overline{R}$.

Next, it follows from (8) and (9) that

$$0={\delta }^{+}\overline{I}-{\delta }^{-}\tilde{\alpha }\overline{R},0={\gamma }^{+}\overline{I}-{\gamma }^{-}\tilde{\beta }\overline{R},$$

hence,

$$\tilde{\alpha }={\displaystyle \frac{{\delta }^{+}\overline{I}}{{\delta }^{-}\overline{R}}},\tilde{\beta }={\displaystyle \frac{{\gamma }^{+}\overline{I}}{{\gamma }^{-}\overline{R}}}.$$

Then,

$$\overline{\alpha }={\mathrm{\Phi }}_{\alpha }\left(\tilde{\alpha }\right),\overline{\beta }={\mathrm{\Phi }}_{\beta }\left(\tilde{\beta }\right).$$

(19)

It follows from (18) that

$$\overline{R}={\displaystyle \frac{\sigma }{(\mu +\rho )}}\overline{I}.$$

Next, since $\overline{I}>0$, we get from (17),

$$\overline{\beta }{\displaystyle \frac{\overline{S}}{\overline{N}}}={\gamma }_0,$$

(20)

where, ${\gamma }_0=\mu +\overline{\alpha }+\sigma$. Next, we let

$${\gamma }_1={\displaystyle \frac{1}{\mu }}\left({\gamma }_0-{\displaystyle \frac{\rho \sigma }{\mu +\rho }}\right),$$

then, it follows from (16), that

$$\overline{S}={\displaystyle \frac{P}{\mu }}-{\gamma }_1\overline{I},$$

Furthermore,

$$\overline{N}=\overline{S}+\overline{I}+\overline{R}={\displaystyle \frac{P}{\mu }}-{\gamma }_2\overline{I},$$

(21)

where

$${\gamma }_2={\gamma }_1-1-{\displaystyle \frac{\sigma }{\mu +\rho }}.$$

Using these expressions in (20), we obtain

$$\left({\displaystyle \frac{P}{\mu }}-{\gamma }_1\overline{I}\right)={\displaystyle \frac{{\gamma }_0}{\overline{\beta }}}\left({\displaystyle \frac{P}{\mu }}-{\gamma }_2\overline{I}\right).$$

Therefore,

$$\overline{I}={\displaystyle \frac{P}{\mu {\gamma }_3}}\left(1-{\displaystyle \frac{{\gamma }_0}{\overline{\beta }}}\right),$$

(22)

where

$${\gamma }_3={\gamma }_1-{\displaystyle \frac{{\gamma }_0{\gamma }_2}{\overline{\beta }}}.$$

The Jacobian evaluated at the EE is,

$$J(\overline{S},\overline{I},\overline{R})=\left(\begin{array}{ccc}-\mu -\overline{\beta }\lambda {\gamma }_4{\left(\overline{I}\right)}^2& -\overline{\beta }\lambda {\gamma }_5{\gamma }_6& \rho +\overline{\beta }\lambda \overline{I}{\gamma }_5\\ \overline{\beta }\lambda {\gamma }_4{\left(\overline{I}\right)}^2& -{\gamma }_0+\overline{\beta }\lambda {\gamma }_5{\gamma }_6& -\overline{\beta }\lambda \overline{I}{\gamma }_5\\ 0& \sigma & -\mu -\rho \end{array}\right)$$

(23)

Here, using (21), we set

$$\begin{array}{cc}\hfill & \lambda ={\displaystyle \frac{1}{{\overline{N}}^2}}={\displaystyle \frac{1}{{\left({\displaystyle \frac{P}{\mu }}-{\gamma }_2\overline{I}\right)}^2}},{\gamma }_4=\left(1+{\displaystyle \frac{\sigma }{\mu +\rho }}\right)\hfill \\ \hfill & {\gamma }_5=\left({\displaystyle \frac{P}{\mu }}-{\gamma }_1\overline{I}\right),{\gamma }_6=\left({\displaystyle \frac{P}{\mu }}-({\gamma }_2+1)\overline{I}\right)\hfill \end{array}$$

Collecting these results yields the following result.

Result 4.

The endemic equilibrium $(\overline{S},\overline{I},\overline{R})$ exists only when (15) holds, and is given by

$$\overline{S}={\displaystyle \frac{P}{\mu }}-{\gamma }_1\overline{I},\overline{I}={\displaystyle \frac{P}{\mu {\gamma }_3}}\left(1-{\displaystyle \frac{{\gamma }_0}{\overline{\beta }}}\right),\overline{R}={\displaystyle \frac{\sigma }{(\mu +\rho )}}\overline{I},$$

(24)

where $\overline{\alpha }$ and $\overline{\beta }$ are given in (19).

We recall that ${\gamma }_0=\mu +\overline{\alpha }+\sigma$, therefore condition (15) shows that the DFE is unstable, and then (24) guaranties that $\overline{I}>0$ and then that $\overline{R}>0$.

Next, in our simulations we found a case in which the EE exists. Starting with the initial conditions $S_0=9,500,I_0=400,$ and $R_{init}=100$, numerical simulations produced the following EE steady state values:

$$\overline{S}=2,183,\overline{I}=310,\overline{R}=43,\overline{\alpha }=0.0107,\overline{\beta }=0.1345$$

Indeed, (15) holds, i.e., $R_0>1$. Then, the eigenvalues of $J(\overline{S},\overline{I},\overline{R})$, (23), were found numerically to be

$${\lambda }_1=-0.7477,{\lambda }_2=-0.0174,{\lambda }_3=-0.0014,$$

confirming that the EE is stable and attracting.

5. Simulations

To obtain insight into the model predictions, an explicit time-stepping algorithm was constructed, and the related program was written in MATLAB. This provided quantitative approximations of the model solutions. The code performed very well, and the runs typically took a few seconds on a laptop. We first present the numerical scheme and then two simulation results. The first deals with the DFE and the second with the EE. In the next section, we show the numerical convergence of the algorithm.

5.1. Numerical Scheme

The time interval $[0,T]$ is discretized into N subintervals of length $\Delta t=T/N$. Typically, in simulations $\Delta t={10}^{-4}$ is used. Denoting the value of a function f at time $t_n=n\Delta t$ by $f^n$, the discretized system, for $n=0,1,\dots N-1$, is as follows.

$$\begin{array}{ccc}\hfill S^{n+1}& =& S^n+\left(P^n-{\beta }^n{I}^n{\displaystyle \frac{S^n}{N^n}}-\mu S^n+\rho R^n\right)\Delta t,\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill I^{n+1}& =& I^n+\left({\beta }^n{I}^n{\displaystyle \frac{S^{n+1}}{N^n}}-(\mu +{\alpha }^n)I^n-\sigma I^n\right)\Delta t,\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill R^{n+1}& =& R^n+\left(\sigma I^{n+1}-\mu R^n-\rho R^n\right)\Delta t,\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill {\tilde{\alpha }}^{n+1}& =& {\tilde{\alpha }}^n+\left({\delta }^{+}{I}^{n+1}-{\delta }^{-}{\tilde{\alpha }}^n{R}^{n+1}\right)\Delta t.\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill {\tilde{\beta }}^{n+1}& =& {\tilde{\beta }}^n+\left({\gamma }^{+}{I}^{n+1}-{\gamma }^{-}{\tilde{\beta }}^n{R}^{n+1}\right)\Delta t,\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill {\alpha }^{n+1}& =& {\mathrm{\Phi }}_{\alpha }\left({\tilde{\alpha }}^{n+1}\right),{\beta }^{n+1}={\mathrm{\Phi }}_{\beta }\left({\tilde{\beta }}^{n+1}\right).\hfill \end{array}$$

(30)

Together with the initial conditions,

$$S^0=S_0,I^0=I_0,R^0=R_{init},{\tilde{\alpha }}^0={\alpha }_0,{\tilde{\beta }}^0={\beta }_0.$$

(31)

We note that in this algorithm, we use the updated values of the functions when they are available. Using the above discretization of the model equations, we have the following algorithm.

Algorithm 1 Explicit time-stepping

Set $S^0,I^0,R^0,{\alpha }^0$ and ${\beta }^0$ using (31)       for n = 0 to N-1 do Calculate, in this order, $S^{n+1},I^{n+1},R^{n+1},{\tilde{\alpha }}^{n+1}$ and ${\tilde{\beta }}^{n+1}$ according to (25)–(29) Calculate ${\alpha }^{n+1}$ and ${\beta }^{n+1}$ according to (30)      end for

5.2. DFE and EE Simulations

We show two examples of numerical simulations: one in which the Disease-Free Equilibrium is attained and the other one with an Endemic Equilibrium. These examples were produced using Algorithm 5.1 and the parameters in Table 2.

As noted in section 4.1, the DFE example (Figure 2, see p. 17) begins with 40,000 susceptibles, 2,000 infected, and 200 already recovered from the disease. Since it satisfies equation (14), it tends to the disease-free equilibrium. The results are shown in the graphs below. The infected and recovered populations decrease to zero after about 20 days, at which point we have 41,926 susceptibles making up the total population. We see that the infected population gradually decreases over time, indicating that infected individuals are dying or recovering faster than the rate at which susceptibles are becoming infected. The recovered population increases to almost 1,000 before decreasing to 0, showing that many of the previously infected individuals recovered from the disease. The susceptible population takes a small dip, showing that some individuals moved from the susceptible to the infected group before the disease dies off and the populations reach their DFE. In this simulation, the death rate from the disease initially jumps to the upper bound, but after 3 days it decreases until leveling off between the upper and lower bounds. The infectivity rate changes similarly over the same time period.

The simulations of the Endemic Equilibrium are shown in Figure 3. The initial populations are 9,500 susceptibles, 400 infected, and 200 recovered. Since $\overline{\beta }>(\mu +\overline{\alpha }+\sigma )$, the disease remains in the population as t tends to infinity. When $0\le t\le 300$, we see that both $\alpha$ and $\beta$ are increasing. This results in a significant decline in the susceptible population and a steep increase of the infected population. When $100\le t\le 500$, the infected and recovered populations increase until they reach their maximum values at around $t=500$. After this time, they decrease to their EE values of $\overline{I}\approx 310,$ and $\overline{R}\approx 43$. Both $\alpha$ and $\beta$ reach equilibrium at $t\approx 600$, while the populations $S,I$ and R stabilize when $t\approx 3500$.

6. Numerical Convergence

We show the convergence of the algorithm numerically. Since we do not have the exact solution in closed form, we assume that the numerical solution obtained on the finest refinement is the reference `exact’ solution. Thus, the exact or reference solution ($\tilde{S},\tilde{I},\tilde{R}$) is the one obtained using Algorithm 5.1 with $\Delta t={10}^{-6}$. The parameters are as in Table 2. We calculate the numerical solution $u_k$ with $\Delta t={10}^{-k}$, for $k=2,\dots ,5$, and calculate the ${\ell }^{\infty }$ norm of the error up to the final time $T=20$. The rate of convergence is calculated as

$$R_{conv}={\displaystyle \frac{1}{log\left(0.1\right)}}log\left({\displaystyle \frac{\parallel \tilde{u}-u_k{\parallel }_{{\ell }^{\infty }}}{\parallel \tilde{u}-u_{k-1}{\parallel }_{{\ell }^{\infty }}}}\right).$$

The results are presented in Table 3, where we observe that the subsequent numerical solutions converge to the reference solution, and the error decays at a rate of at least 1. This result provides substantial confidence in the numerical method and therefore, in the numerical solutions.

7. Conclusions

This work describes a new augmented version of the SIR model that mathematically studies the spread of a generic disease. Our SIR-ID model modifies the usual SIR model by incorporating system dependent variables for the rate of death, $\alpha$, and the rate of infectivity, $\beta$. The motivation behind our particular model comes from observing the COVID-19 pandemic. As noted in the Introduction, various factors cause changes in the death rate and similarly the infectivity rate, such as improved health care, vaccinations, and partial immunity to subsequent infections. Moreover, the limits of the healthcare system to treat infected populations, thus straining hospitals, is another factor.

To make the model applicable to a specific disease, the death and infectivity rate functions should be determined from its field data. For the model to make epidemiological sense, the differential equations for the death rate ($\tilde{\alpha }$) and the infectivity rate ($\tilde{\beta }$) are incorporated and modified so that they do not fall below or exceed certain realistic threshold values.

Following the detailed derivation of the model, we establish the positivity and boundedness of the solutions, which leads to the proof of their existence. This shows that the model makes sense both biologically and mathematically. Then, we study the system’s steady states and show that the disease-free equilibrium (DFE) is stable and attracting (asymptotically stable) when condition (4.2) holds. Following the DFE, we study the conditions for the emergence of the endemic equilibrium (EE). Similarly to the DFE, we find the eigenvalues of the Jacobian matrix at the EE. Then, we present a simple time-stepping algorithm for the numerical approximations of the solutions. The algorithm has been incorporated in a computer code and two sets of simulations performed. In the first, the system parameters are such that it converges, as $t\to \infty$, to the DFE, showing that in this case the DFE is stable and attracting (asymptotically stable). In the second simulation the system parameters are such that it converges to the EE state, since the DFE is unstable. Indeed, the simulations of the model support the DFE and EE, and show how the system approaches them. As noted above, condition (15) guarantees that when the DFE is unstable (24) guarantees that $\overline{I}>0$, and therefore $\overline{R}>0$, too, so that the EE is found to be stable and attracting.

We emphasize that in both simulations the values of $\alpha$ and $\beta$ eventually do not reach the constraints, allowing them to be determined by the system. This shows that in these cases imposing their values would lead to different results, which is one of the main motivations of this research.

Now that the basic structure and behavior of the model have been investigated, there are related issues that warrant further study. Indeed, there is some interest in the study of related extensions of our model, some of which we mention below.

One of the main open questions comes from the introduction of system dependent death and infectivity rates. We could expand these rates to model certain disease-informed functional forms of ${\mathrm{\Psi }}_{\alpha }$ and ${\mathrm{\Psi }}_{\beta }$. This would make it necessary to find the relevant information and see if we can fit the functions to the data. This would make our model applicable to disease-specific results. Moreover, it may be of practical interest to make the recovery rate $\sigma$ a state variable, too.

Since our analysis suggests the possibility of multiple steady state solutions, the question is of the existence and reachability of such steady states. Reachability refers to the question of given a steady state, how to find the initial conditions that guarantee that the time dependent solution converges to this steady state ias $t\to \infty$. To better understand the model, reachability would be an area of interest.

Another area of interest may be the extension of the model that allows for the sudden appearance of new variants of the disease agents, as was the case with COVID-19. This could be done by introducing randomness to some of the system coefficients.

Furthermore, our model could be used in optimization analysis, leading to applications in vaccination/treatment, general public health, or economic studies. Additional simulations are always of interest, together with a comparison to the given information one may have on current diseases.

Thus, this work is a contribution to the rapidly growing discipline of mathematical epidemiology.