Global Stability Analysis of Two-Stage Quarantine-Isolation Model with Holling Type II Incidence Function

Abstract

A new two-stage model for assessing the effect of basic control measures, quarantine and isolation, on a general disease transmission dynamic in a population is designed and rigorously analyzed. The model uses the Holling II incidence function for the infection rate. First, the basic reproduction number (${\mathcal{R}}_0$) is determined. The model has both locally and globally asymptotically stable disease-free equilibrium whenever ${\mathcal{R}}_0<1$. If ${\mathcal{R}}_0>1,$ then the disease is shown to be uniformly persistent. The model has a unique endemic equilibrium when ${\mathcal{R}}_0>1.$ A nonlinear Lyapunov function is used in conjunction with LaSalle Invariance Principle to show that the endemic equilibrium is globally asymptotically stable for a special case.

1. Introduction

Over the decades, quarantine (of individuals suspected of being exposed to a communicable disease) and isolation (of individuals with disease symptoms) have been widely used to control the spread of numerous communicable diseases, such as pandemic influenza, cholera, Ebola, Severe Acute Respiratory Syndrome (SARS), and most recently swine influenza pandemic [1,2,3,4,5,6,7,8,9]. Numerous mathematical models have been studying the effect of quarantine and isolation in combatting the spread of the diseases (see, for instance, refs. [1,2,4,5,6,7,8,10,11,12,13] and the references therein). In the aforementioned studies, mass action or standard incidence functions were used in the modeling of the transmission dynamics of the diseases. In this study, another nonlinear incidence function (called the Holling type II incidence function) will be used in the modeling of the transmission dynamics of a general disease. The Holling type II incidence function is given by ${\displaystyle g\left(I\right)=\frac{\beta I}{1+\alpha I},}$ with $\alpha >0$, where I is the number of infectious individuals and $\beta$ is the effective contact rate (the average number of contacts sufficient for transmitting infection). The incidence function $g\left(I\right)$ was first used in the study of the cholera epidemic in Bari, Italy by Capasso and Serio [14]. The reason for using the Holling type II incidence functional comes from the information that the number of effective contacts between susceptible individuals and infective individuals may saturate at very high levels due to behavioral changes or due to crowding of infective people taken by the people in reaction to the severity of the disease [15,16]. It is well known that some infectious diseases, such as influenza [17] and HIV [18], have multiple disease (infection) stages in their transmission dynamics.

The main purpose of this study is to offer a deep qualitative analysis of a new two-stage model for the transmission dynamics of a disease that can be controlled by using quarantine and isolation, where the Holling type II incidence function is used.

The paper is organized as follows. The formulation of the model is given in Section 2. The local and global asymptotic stability of the disease-free equilibrium (DFE) is analyzed in Section 3. The existence of the endemic equilibrium is provided in Section 4. Global stability proof for the endemic equilibrium for the special case is also analyzed using a nonlinear Lyapunov function.

2. Model Formulation

The total population at time t, denoted by $N\left(t\right)$ is sub-divided into ten compartments of susceptible $\left(S\right(t\left)\right),$ exposed (with two stages $(E_1\left(t\right)E_2\left(t\right)$), infectious individuals (with two stages $(I_1\left(t\right)I_2\left(t\right)$), Isolated individuals (with two stages $H_1\left(t\right)H_2\left(t\right)$), and recovered $\left(R\right(t\left)\right)$ individuals, so that

$$N\left(t\right)=S\left(t\right)+E_1\left(t\right)+E_2\left(t\right)+Q_1\left(t\right)+Q_2\left(t\right)+I_1\left(t\right)+I_2\left(t\right)+H_1\left(t\right)+H_2\left(t\right)+R\left(t\right).$$

The model is given by the following system of nonlinear differential equations

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{dS}{dt}& =\mathsf{\Pi }-\left(\lambda \right(t)+\mu )S\left(t\right),\hfill \\ \hfill \frac{d{E}_1}{dt}& =\lambda \left(t\right)S\left(t\right)-(a_1+b_1+\mu )E_1\left(t\right),\hfill \\ \hfill \frac{d{E}_2}{dt}& =a_1{E}_1-(a_2+b_2+\mu )E_2\left(t\right),\hfill \\ \hfill \frac{d{Q}_1}{dt}& =b_1{E}_1-(c_1+\mu )Q_1,\hfill \\ \hfill \frac{d{Q}_2}{dt}& =c_1{Q}_1+b_2{E}_2-(c_2+\mu )Q_2,\hfill \\ \hfill \frac{d{I}_1}{dt}& =a_2{E}_2\left(t\right)-(d_1+e_1+{\delta }_1+\mu )I_1\left(t\right),\hfill \\ \hfill \frac{d{I}_2}{dt}& =d_1{I}_1\left(t\right)-(e_2+{\gamma }_1+{\delta }_2+\mu )I_2\left(t\right),\hfill \\ \hfill \frac{d{H}_1}{dt}& =c_2{Q}_2+e_1{I}_1-(f_1+{\delta }_3+\mu )H_1,\hfill \\ \hfill \frac{d{H}_2}{dt}& =f_1{H}_1+e_2{I}_2-({\gamma }_2+{\delta }_4+\mu )H_2,\hfill \\ \hfill \frac{dR}{dlt}& ={\gamma }_1{I}_2\left(t\right)+{\gamma }_2{H}_2\left(t\right)-\mu R\left(t\right),\hfill \end{array}\end{array}$$

(1)

where $\lambda \left(t\right)$ is the infection rate given by

$$\lambda \left(t\right)=\beta \left(\frac{I_1}{1+{\alpha }_1{I}_1}+\frac{\eta I_2}{1+{\alpha }_2{I}_2}\right).$$

(2)

In (2), $\beta$ represents the effective contact rate, where $0<\eta <1$ is a parameter that accounts for the reduction in disease transmission given by infectious individuals ($I_1$) in comparison to infectious individuals in the $I_2$ stage.

Susceptible people $\left(S\right)$ is increased by the recruitment of individuals into the population, at a rate $\mathsf{\Pi }$. This class is decreased by infection (with the rate of $\lambda$). Furthermore, this population is decreased by natural death (at a rate $\mu$; populations in all classes are assumed to have the same natural death rate).

Exposed individuals in stage 1 ($E_1$) are generated with the rate of $\lambda$ and reduced by progression to the next exposed stage ($E_2$; at a rate $a_1$) and quarantine (at a rate $b_1$). Exposed individuals in stage 2 are generated at the rate $a_1$. This population is decreased by the development of clinical symptoms of the disease (at a rate $a_2$) and quarantine (at a rate $b_2$).

The class of quarantined individuals in stage 1 is increased by quarantine of exposed people in stage $E_1$ (at the rate $b_1$) and it is reduced by progression to the second quarantined stage (at a rate $c_1$). Similarly, quarantined people in stage 2 are increased by the quarantine of exposed people in the second stage (at the rate $b_2$) and the progression of quarantined people from the first stage into the second stage (at the rate $c_1$). It is decreased by hospitalization (at a rate $c_2$).

The infectious people in stage 1 are increased when exposed people in the second stage develop symptoms (at the rate $a_2$). It is reduced by progression to the second infectious stage (at a rate $d_1$), hospitalization (isolation) (at a rate $e_1$) and disease-induced death (at a rate ${\delta }_1$). The population of infectious class in the second stage is generated by progression of individuals in the first stage (at a rate $d_1$). It is reduced by isolation (at a rate $e_2$), recovery (at a rate ${\gamma }_1$) and disease-induced death (at a rate ${\delta }_2$).

The population of Isolated individuals in the first stage is increased by the hospitalization of infectious people in stage 1 (at the rate $e_1$) and quarantined individuals in the second stage (at the rate $c_2$). It is decreased by progression to the second Isolated stage (at a rate $f_1$), and disease-induced death (at a rate ${\delta }_3$). The population of Isolated individuals in the second stage is generated by the progression of Isolated individuals from the first stage into the second one (at the rate $f_1$). It is decreased by recovery (a rate ${\gamma }_2$) and disease-induced death (at a rate ${\delta }_4$).

Finally, the recovered individuals is increased by the recovery of infectious individuals and hospitalization individuals (at the rates ${\gamma }_1$ and ${\gamma }_2$, respectively). It is reduced by natural death (at the rate $\mu$). (A flow diagram of the model is depicted in Figure 1. The associated variables and parameters are described in

Table 1. Description of variables and parameters of the model (1).

Variable	Description
$S\left(t\right)$	Population of susceptible individuals
$E_1\left(t\right)$	Population of exposed individuals on the first exposed stage
$E_2\left(t\right)$	Population of exposed individuals on the second exposed stage
$Q_1\left(t\right)$	Population of quarantined individuals on the first quarantined stage
$Q_2\left(t\right)$	Population of quarantined individuals on the second quarantined stage
$I_1\left(t\right)$	Population of infected individuals on the first infectious stage
$I_2\left(t\right)$	Population of infected individuals on the second infectious stage
$H_1\left(t\right)$	Population of Isolated individuals on the first Isolated stage
$H_2\left(t\right)$	Population of Isolated individuals on the second Isolated stage
$R\left(t\right)$	Population of recovered individuals
Parameter	Description
$\mathsf{\Pi }$	Recruitment rate
$\beta$	Effective contact rate
$a_1$	Progression rate from the first exposed stage to the second one
$a_2$	Progression rate to first infectious class from exposed individuals
	in the second stage
$b_1$	Quarantine rate of exposed individuals on the first exposed stage
$b_2$	Quarantine rate of exposed individuals on the second exposed stage
$c_1$	Progression rate from the first quarantined stage to the second one
$c_2$	Progression rate to first Isolated class from quarantined individuals
	in the second stage
$d_1$	Progression rate from the first infectious stage to the second one
$e_1$	Hospitalization rate of infectious individuals on the first infectious
$e_2$	Hospitalization rate of infectious individuals on the second infectious
$f_1$	Progression rate from the first Isolated stage to the second one
${\gamma }_1$	Recovery rate of infectious individuals in the second stage
${\gamma }_2$	Recovery rate of Isolated individuals in the second stage
${\delta }_1$	Disease-induced death rate of the first infectious stage
${\delta }_2$	Disease-induced death rate of the second infectious stage
${\delta }_3$	Disease-induced death rate of the first Isolated stage
${\delta }_4$	Disease-induced death rate of the second Isolated stage
$\mu$	Natural death rate

):

It should be noted the model (1) is different by the basic model considered in [19] by

(a)

Using a Holling type incidence function to model the infection rate (the standard incidence function was used in [19])

(b)

Considering two stages for the infectious compartments (Exposed, infected, quarantined, and isolated compartments)

2.1. Preliminaries and Basic Properties

Since the model (1) for human populations, all its parameters are non-negative. Furthermore, the following non-negativity result holds.

Theorem 1.

All variables of the model (1) are non-negative for all $t>0$. This mean, the solutions of system (1) with positive initial conditions will remain positive for all time $t>0$.

Proof. 

Let

$$t_1=sup\{t>0:S>0,E_1>0,E_2>0,Q_1>0,Q_2>0,I_1>0,I_2>0,H_1>0,H_2>0,R>0\in [0,t]\}.$$

Hence, $t_1>0.$ From the first equation of the system (1) it follows that

$$\frac{d}{dt}\left[S\left(t\right)exp\left\{\mu t+{\int }_0^t\lambda \left(\tau \right)d\tau \right\}\right]=\mathsf{\Pi }exp\left\{\mu t+{\int }_0^t\lambda \left(\tau \right)d\tau \right\}.$$

which gives,

$$S\left(t_1\right)exp\left\{\mu t_1+{\int }_0^{t_1}\lambda \left(\tau \right)d\tau \right\}-S\left(0\right)={\int }_0^{t_1}\mathsf{\Pi }exp\left\{\mu y+{\int }_0^y\lambda \left(\tau \right)d\tau \right\}dy,$$

hence,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill S\left(t\right)& =S\left(0\right)exp\left\{-\mu t_1-{\int }_0^{t_1}\lambda \left(\tau \right)d\tau \right\}\hfill \\ & \left[exp\left\{-\mu t_1-{\int }_0^{t_1}\lambda \left(\tau \right)d\tau \right\}\right]{\int }_0^{t_1}\mathsf{\Pi }exp\left\{\mu y+{\int }_0^y\lambda \left(\tau \right)d\tau \right\}dy>0.\hfill \end{array}\end{array}$$

In the same way, it can be shown that $E_1>0,E_2>0,Q_1>0,Q_2>0,I_1>0,I_2>0,H_1>0,$$H_2>0$ and $R>0$ for all time $t>0$. □

Lemma 1.

The closed set

$$\mathcal{D}=\left\{(S,E_1,E_2,Q_1,Q_2,I_1,I_2,H_1,H_2,R)\in {\mathbb{R}}_{+}^{10}:S+E_1+E_2+Q_1+Q_2+I_1+I_2+H_1+H_2+R\le \frac{\mathsf{\Pi }}{\mu }\right\}$$

is positively invariant.

Proof. 

Adding all the equations of the model (1) gives,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{dN}{dt}& =\mathsf{\Pi }-\mu N-({\delta }_1{I}_1+{\delta }_2{I}_2+{\delta }_3{H}_1+{\delta }_4{H}_2).\hfill \end{array}\end{array}$$

(3)

It follows that ${\displaystyle \frac{dN}{dt}\le \mathsf{\Pi }-\mu N}$, thus ${\displaystyle \frac{dN}{dt}\le 0}$ provided that ${\displaystyle N\ge \frac{\mathsf{\Pi }}{\mu }}$. By using standard comparison theorem [20] it can be shown that $N\le {\displaystyle N\left(0\right)e^{-\mu t}+\frac{\mathsf{\Pi }}{\mu }(1-e^{-\mu t}})$. In particular, ${\displaystyle N\left(t\right)\le \frac{\mathsf{\Pi }}{\mu }}$ if ${\displaystyle N\left(0\right)\le \frac{\mathsf{\Pi }}{\mu }.}$ Thus, the region $\mathcal{D}$ is positively invariant. Furthermore, if ${\displaystyle N\left(0\right)>\frac{\mathsf{\Pi }}{\mu },}$ then either the solution enters $\mathcal{D}$ in finite time, or $N\left(t\right)$ approaches ${\displaystyle \frac{\mathsf{\Pi }}{\mu }}$ asymptotically. Hence, the region $\mathcal{D}$ attracts all solutions in ${\mathbb{R}}_{+}^{10}$. □

Since the region $\mathcal{D}$ is positively invariant, it is sufficient to consider the dynamics of the flow generated by the model (1) in $\mathcal{D}$, where the usual existence, uniqueness, continuation results hold for the system [21].

Next-Generation Method

Suppose that the population is divided into n compartments, with $m<n$ infected compartments. At time $t,$ let $x_i\left(t\right)$ be the number of infected individuals in the $i^{th}$ infected class such that

$$\frac{d{x}_i}{dt}=F_i\left(x\right)-V_i\left(x\right),\mathrm{with}{V}_i=V_i^{-}\left(x\right)-V_i^{+}\left(x\right)\mathrm{for}i=1,2,\dots ,m,$$

(4)

where $F_i\left(x\right)$ represents the rate of appearance of new infections in class i, $V_i+\left(x\right)$ represents the rate of transfer of individuals into class i by all other means, and $V_i^{-}\left(x\right)$ represents the rate of transfer of individuals out of class i. System can be rewritten as follows

$$\dot{X}=F\left(X\right)-V\left(X\right),$$

(5)

with, $F\left(X\right)={(F_1,F_2,\dots ,F_m)}^T$ and $V\left(X\right)={(V_1,V_2,\dots ,V_m)}^T.$

Lemma 2.

(van den Driessche and Watmough [22]). If $\overline{x}$ is a DFE of (5), then the derivatives $D\mathrm{F}\left(\overline{x}\right)$ and $D\mathrm{V}\left(\overline{x}\right)$ are partitioned as

$$D\mathrm{F}\left(\overline{x}\right)=\left(\begin{array}{cc}F& 0\\ 0& 0\end{array}\right),D\mathrm{V}\left(\overline{x}\right)=\left(\begin{array}{cc}V& 0\\ J_3& J_4\end{array}\right),$$

where F and V are the $m\times m$ matrices defined by,

$$F=\left[\begin{array}{c}{\displaystyle \frac{\partial {\mathrm{F}}_i}{\partial x_j}\left(\overline{x}\right)}\end{array}\right]andV=\left[\begin{array}{c}{\displaystyle \frac{\partial {\mathrm{V}}_i}{\partial x_j}\left(\overline{x}\right)}\end{array}\right]with1\le i,j\le m.$$

Furthermore, F is non-negative, V is a non-singular $M-$matrix and $J_3,J_4$ are matrices associated with the transition terms of the model, and all eigenvalues of $J_4$ have positive real parts.

Now, the next-generation matrix is given by ${\displaystyle F{V}^{-1}}$ and the spectral radius (the largest eigenvalue) of ${\displaystyle F{V}^{-1}}$ is the basic reproduction number of the model (5) [22].

3. Stability of DFE

3.1. Local Stability

The DFE of the model (1) is given by

$${\mathcal{E}}_0=(S^{*},E_1^{*},E_2^{*},Q_1^{*},Q_2^{*},I_1^{*},I_2^{*},H_1^{*},H_2^{*},R^{*})=(\mathsf{\Pi }/\mu ,0,0,0,0,0,0,0,0,0).$$

(6)

The next-generation operator method [22,23] will be used to analyze the stability of ${\mathcal{E}}_0$. Using the same notation in the previous section, the non-negative matrix, F and the M-matrix, V are given by

$$F=\left(\begin{array}{cccccccc}0& 0& 0& 0& {\displaystyle \frac{\beta \mathsf{\Pi }}{\mu }}& {\displaystyle \frac{\eta \beta \mathsf{\Pi }}{\mu }}& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right),$$

and,

$$V=\left(\begin{array}{cccccccc}{k}_1& 0& 0& 0& 0& 0& 0& 0\\ -a_1& k_2& 0& 0& 0& 0& 0& 0\\ -b_1& 0& k_3& 0& 0& 0& 0& 0\\ 0& -b_2& -c_1& k_4& 0& 0& 0& 0\\ 0& -a_2& 0& 0& k_5& 0& 0& 0\\ 0& 0& 0& 0& -d_1& k_6& 0& 0\\ 0& 0& 0& -c_2& -e_1& 0& k_7& 0\\ 0& 0& 0& 0& 0& -e_2& -f_1& k_8\end{array}\right).$$

The control reproduction number [24,25], denoted by ${\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)$ is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{R}}_0& {\displaystyle =\frac{\beta \mathsf{\Pi }{a}_1{a}_2(k_6+\eta d_1)}{\mu k_1{k}_2{k}_5{k}_6},}\hfill \end{array}\end{array}$$

where,

$$\begin{array}{c}\hfill \begin{array}{c}{k}_1=\mu +a_1+b_1,k_2=\mu +a_2+b_2,k_3=\mu +c_1,k_4=\mu +c_2,k_5=\mu +{\delta }_1+d_1+e_1,\hfill \\ k_6=\mu +{\delta }_2+{\gamma }_1+e_2,k_7=\mu +{\delta }_3+f_1,k_8=\mu +{\delta }_4+{\gamma }_2.\hfill \end{array}\end{array}$$

The following result is established by using Theorem 2 in [22].

Lemma 3.

The model (1) has a locally asymptotically stable (LAS) DFE whenever ${\mathcal{R}}_0<1$. Moreover, the DFE of this model is unstable if ${\mathcal{R}}_0>1$.

The average number of new infections generated by a single infectious individual in a population is measured by the quantity ${\mathcal{R}}_0$. The epidemiological implication of Lemma 3 is that the disease dies out from the population (when ${\mathcal{R}}_0<1$) if the initial sizes of the sub-populations of the model are in the basin of attraction of the DFE ($E_0$). To make sure that disease dies out form the population regardless of the initial sizes of sub-populations, it is necessary to show that the DFE is globally asymptotically stable (GAS) if ${\mathcal{R}}_0<1$. This is established below.

3.2. Global Stability of DFE

Theorem 2.

The model (1) has GAS DFE, given by (6), in $\mathcal{D}$ whenever ${\mathcal{R}}_0\le 1$.

Proof. 

Define the following Lyapunov function:

$$\mathcal{F}=\left(\frac{a_1{a}_2(k_6+\eta d_1)}{\eta k_1{k}_2{k}_5}\right)E_1+\left(\frac{a_1(k_6+\eta d_1)}{\eta k_2{k}_5}\right)E_2+\left(\frac{k_6+\eta d_1}{\eta k_5}\right)I_1+I_2,$$

differentiate $\mathcal{F}$ with respect to t gives

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \dot{\mathcal{F}}& =\left(\frac{a_1{a}_2(k_6+\eta d_1)}{\eta k_1{k}_2{k}_5}\right){\dot{E}}_1+\left(\frac{a_2(k_6+\eta d_1)}{\eta k_2{k}_5}\right){\dot{E}}_2+\left(\frac{k_6+\eta d_1}{\eta k_5}\right){\dot{I}}_1+{\dot{I}}_2\hfill \\ & =\left(\frac{a_1{a}_2(k_6+\eta d_1)}{\eta k_1{k}_2{k}_5}\right)\left[\beta S\left(\frac{I_1}{1+{\alpha }_1{I}_1}+\frac{\eta I_2}{1+{\alpha }_2{I}_2}\right)-k_1{E}_1\right]+\left(\frac{a_2(k_6+\eta d_1)}{\eta k_2{k}_5}\right)\left[a_1{E}_1-k_2{E}_2\right]\hfill \\ & +\left(\frac{k_6+\eta d_1}{\eta k_5}\right)\left[a_2{E}_2-k_5{I}_1\right]+d_1{I}_1-k_6{I}_2\hfill \\ & \le \left(\frac{a_1{a}_2(k_6+\eta d_1)}{\eta k_1{k}_2{k}_5}\right)\left[\beta \frac{\mathsf{\Pi }}{\mu }(I_1+\eta I_2)-k_1{E}_1\right]+\left(\frac{a_2(k_6+\eta d_1)}{\eta k_2{k}_5}\right)\left[a_1{E}_1-k_2{E}_2\right]\hfill \\ & +\left(\frac{k_6+\eta d_1}{\eta k_5}\right)\left[a_2{E}_2-k_5{I}_1\right]+d_1{I}_1-k_6{I}_2\hfill \\ & =\left(\frac{\beta \mathsf{\Pi }{a}_1{a}_2(k_6+\eta d_1)}{\mu \eta k_1{k}_2{k}_5}\right)(I_1+\eta I_2)+\left(d_1-\frac{k_6+\eta d_1}{\eta }\right)I_1-k_6{I}_2\hfill \\ & =\frac{k_6}{\eta }({\mathcal{R}}_0-1)(I_1+\eta I_2)\hfill \end{array}\end{array}$$

Since all the variables and the parameters of the model (1) are non-negative, it follows that $\dot{\mathcal{F}}\le 0$ for ${\mathcal{R}}_0\le 1$ with $\dot{\mathcal{F}}=0$ if and only if $E=I_u=I_e=0$. Thus, $\mathcal{F}$ defined a Lyapunov function on $\mathcal{D}$. Hence,

$$(E_1,E_2,I_1,I_2)\to (0,0,0,0)\mathrm{as}t\to \infty .$$

(7)

It can be easily shown that $(Q_1,Q_2,H_1,H_2,R)\to (0,0,0,0,0)$ and ${\displaystyle S\to \frac{\mathsf{\Pi }}{\mu }}$ as $t\to \infty .$ Furthermore, the region $\mathcal{D}$ is an invariant and attracting set of $R_{+}^{10}$, and the largest compact invariant set in $(S,E_1,E_2,Q_1,Q_2,I_1,I_2,H_1,H_2,R)\in \mathcal{D}:\mathcal{F}=0$ is the singleton $\left\{{\mathcal{E}}_0\right\}$. Thus, by Invariance Principle [26], every solution of the system (1), and initial conditions in $R_{+}^{10}$, approaches the DFE (${\mathcal{E}}_0$) as $t\to \infty$ whenever ${\mathcal{R}}_0<1$. □

The above result shows that the disease dies out from the population if the reproduction number of the model is less than one. The epidemiological implication of the above theorem is that the use of isolation and quarantine can lead to elimination of the disease if both controls can keep the threshold quantity, ${\mathcal{R}}_0$, to a value less than unity (i.e., The condition ${\mathcal{R}}_0<1$ is sufficient and necessary for the elimination of the disease ). Figure 2 illustrate numerical results obtained by simulating the model (1) using various initial conditions for the case ${\mathcal{R}}_0<1.$ Its clear that the solutions are converged to the DFE.

4. Existence and Stability for Endemic Equilibrium Point

4.1. Persistence of the Disease

The persistence of the disease in the population will be investigated below. The model system (1) is said to be uniformly persistent if there exists a constant c such that any solution $(S\left(t\right),E_1\left(t\right),E_2\left(t\right),Q_1\left(t\right),Q_2\left(t\right),I_1\left(t\right),I_2\left(t\right),H_1\left(t\right),H_2\left(t\right),R\left(t\right))$ satisfies ([27,28]):

$$\begin{array}{c}\hfill \begin{array}{c}\underset{t\to \infty }{lim\; inf}S\left(t\right)\ge c,\underset{t\to \infty }{lim\; inf}{E}_1\left(t\right)\ge c,\underset{t\to \infty }{lim\; inf}{E}_2\left(t\right)\ge c,\underset{t\to \infty }{lim\; inf}{Q}_1\left(t\right)\ge c,\hfill \\ \underset{t\to \infty }{lim\; inf}{Q}_2\left(t\right)\ge c,\underset{t\to \infty }{lim\; inf}{I}_1\left(t\right)\ge c,\underset{t\to \infty }{lim\; inf}{I}_2\left(t\right)\ge c,\underset{t\to \infty }{lim\; inf}{H}_1\left(t\right)\ge c,\hfill \\ \underset{t\to \infty }{lim\; inf}{H}_2\left(t\right)\ge c,\underset{t\to \infty }{lim\; inf}R\left(t\right)\ge c,\hfill \end{array}\end{array}$$

provided that $(S\left(0\right),E_1\left(0\right),E_2\left(0\right),Q_1\left(0\right),Q_2\left(0\right),I_1\left(0\right),I_2\left(0\right),H_1\left(0\right),H_2\left(0\right),R\left(0\right))\in \mathcal{D}$.

Theorem 3.

The model (1) is uniformly persistent in $\mathcal{D}$ if and only if ${\mathcal{R}}_0>1$.

Proof. 

The proof of the above theorem follows from using the same approach given in [29] to prove Proposition 3.3 of [29], which is applying a uniform persistence theorem in [27] and noting that the DFE of the model (1) is unstable whenever ${\mathcal{R}}_0>1$ (Lemma 3). □

Whenever ${\mathcal{R}}_0>1$ its clear (from Theorem (3)) that the model (1) is uniformly persistent. Moreover using Theorem 2.8.6 in [30] and Theorem D.3 in [20] gives the model (1) has at least one endemic equilibrium in $\mathcal{D}$. Hence, the following Lemma is concluded.

Lemma 4.

System (1) has at least one endemic equilibrium provided that ${\mathcal{R}}_0>1.$

The uniqueness of this equilibrium will be analyzed in the coming subsection.

4.2. Uniqueness of Endemic Equilibrium Point (EEP)

Let,

$${\mathcal{E}}_1=(S^{**},E_1^{**},E_2^{**},Q_1^{**},Q_2^{**},I_1^{**},I_2^{**},H_1^{**},H_2^{**},R^{**})$$

represents any arbitrary EEP of the model (1). Furthermore, define

$${\lambda }^{**}=\frac{\beta I_1^{**}}{1+{\alpha }_1{I}_1^{**}}+\frac{\beta \eta I_2^{**}}{1+{\alpha }_2{I}_2^{**}}$$

(8)

(the force of infection of the model (1) at steady-state). It follows, by solving the equations in (1) at steady-state that

$$\begin{array}{c}\hfill \begin{array}{c}{S}^{**}=\frac{\mathsf{\Pi }}{{\lambda }^{**}+\mu },E_1^{**}=\frac{\mathsf{\Pi }{\lambda }^{**}}{({\lambda }^{**}+\mu )k_1},E_2^{**}=\frac{a_1\mathsf{\Pi }{\lambda }^{**}}{({\lambda }^{**}+\mu )k_1{k}_2},\hfill \\ Q_1^{**}=\frac{b_1\mathsf{\Pi }{\lambda }^{**}}{({\lambda }^{**}+\mu )k_1{k}_3},Q_2^{**}=\frac{(a_1{b}_2{k}_3+b_1{c}_1{k}_2)\mathsf{\Pi }{\lambda }^{**}}{({\lambda }^{**}+\mu )k_1{k}_2{k}_3{k}_4},I_1^{**}=\frac{a_1{a}_2\mathsf{\Pi }{\lambda }^{**}}{({\lambda }^{**}+\mu )k_1{k}_2{k}_5},\hfill \\ I_2^{**}=\frac{a_1{a}_2{d}_1\mathsf{\Pi }{\lambda }^{**}}{({\lambda }^{**}+\mu )k_1{k}_2{k}_5{k}_6}.\hfill \end{array}\end{array}$$

(9)

Substituting $I_1^{**}$ and $I_2^{**}$ in (9) into (8) gives the following quadratic equation (in terms of ${\lambda }^{**}$):

$$M_0{\left({\lambda }^{**}\right)}^2+M_1{\lambda }^{**}+M_2=0,$$

(10)

with,

$$\begin{array}{c}\hfill \begin{array}{c}{M}_0=(\mathsf{\Pi }{a}_1{a}_2{\alpha }_1+k_1{k}_2{k}_5)(\mathsf{\Pi }{a}_1{a}_2{d}_1{\alpha }_2+k_1{k}_2{k}_5{k}_6),\hfill \\ M_1=-{\mathsf{\Pi }}^2{a}_1^2{a}_2^2{\alpha }_1{d}_1\eta \beta -{\mathsf{\Pi }}^2{a}_1^2{a}_2^2{\alpha }_2{d}_1\beta -\mathsf{\Pi }{a}_1{a}_2{k}_1{k}_2{k}_5{k}_6\beta -\mathsf{\Pi }{a}_1{a}_2{d}_1{k}_1{k}_2{k}_5\eta \beta \hfill \\ +\mathsf{\Pi }{a}_1{a}_2{k}_1{k}_2{k}_5{k}_6{\alpha }_1\mu +\mathsf{\Pi }{a}_1{a}_2{d}_1{k}_1{k}_2{k}_5{\alpha }_2\mu +2\mu k_1^2{k}_2^2{k}_5^2{k}_6,\hfill \\ M_2={\mu }^2{k}_1^2{k}_2^2{k}_5^2{k}_6(1-{\mathcal{R}}_0).\hfill \end{array}\end{array}$$

By solving for ${\lambda }^{**}$ in (10) and substituting the positive values of ${\lambda }^{**}$ into the expressions in (9) the endemic equilibria of the model (1) can then be obtained. It should be noted that $M_0>0$ and $M_2<0$ whenever ${\mathcal{R}}_0>1.$ Thus, by using the Descartes Rule of Signs on the quadratic Equation (10), the following result is established.

Lemma 5.

The model (1) has a unique endemic (positive) equilibrium, given by ${\mathcal{E}}_1$ whenever ${\mathcal{R}}_0>1.$

4.3. Global Stability for Endemic Equilibrium

In this section, the global stability of the endemic equilibrium of the model (1) is given for the special case where the associated disease-induced mortality in all classes is negligible (so that ${\delta }_1={\delta }_2={\delta }_3={\delta }_4=0$).

Define

$${\mathcal{R}}_s={\mathcal{R}}_0{|}_{{\delta }_1={\delta }_2={\delta }_3={\delta }_4=0}$$

$${\mathcal{D}}_0=\left\{(S,E_1,E_2,Q_1,Q_2,I_1,I_2,H_1,H_2,R)\in \mathcal{D}:E_1=E_2=I_1=I_2=H_1=H_2=R=0\right\}.$$

We claim the following result,

Theorem 4.

The endemic equilibrium of the model (1) with ${\delta }_1={\delta }_2={\delta }_3={\delta }_4=0$ is GAS in $\mathcal{D}\backslash {\mathcal{D}}_0$ if ${\mathcal{R}}_s>1$.

Proof. 

Let ${\mathcal{R}}_s>1$, such that the endemic equilibrium exists. Furthermore, define the following Lyapunov function:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{F}=& \frac{1}{2}\left[(S-S^{**})+(E_1-E_1^{**})+(E_1-E_1^{**})+(Q_1-Q_1^{**})+(Q_2-Q_2^{**})\right.\hfill \\ & +{\left.(I_1-I_1^{**})+(I_2-I_2^{**})+(H_1-H_1^{**})+(H_2-H_2^{**})+(R-R^{**})\right]}^2\hfill \end{array}\end{array}$$

with Lyapunov derivative

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \dot{\mathcal{F}}=& \left[(S-S^{**})+(E_1-E_1^{**})+(E_2-E_2^{**})+(Q_1-Q_1^{**})+(Q_2-Q_2^{**})\right.\hfill \\ & +\left.(I_1-I_1^{**})+(I_2-I_2^{**})+(H_1-H_1^{**})+(H_2-H_2^{**})+(R-R^{**})\right]\dot{N}\left(t\right).\hfill \end{array}\end{array}$$

Since ${\displaystyle S^{**}+E_1^{**}+E_2^{**}+Q_1^{**}+Q_2^{**}+I_1^{**}+I_2^{**}+H_1^{**}+H_2^{**}+R^{**}=\frac{\mathsf{\Pi }}{\mu }}$ and $\dot{N}\left(t\right)=\mathsf{\Pi }-\mu N\left(t\right)$ it follows that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \dot{\mathcal{F}}& =\left[N\left(t\right)-\frac{\mathsf{\Pi }}{\mu }\right]\left[\mathsf{\Pi }-\mu N\left(t\right)\right]\hfill \\ & =\frac{1}{\mu }\left[\mu N\left(t\right)-\mathsf{\Pi }\right]\left[\mathsf{\Pi }-\mu N\left(t\right)\right]\hfill \\ & =\frac{-1}{\mu }{\left[\mathsf{\Pi }-\mu N\left(t\right)\right]}^2.\hfill \end{array}\end{array}$$

□

It follows that $\dot{\mathcal{F}}\le 0$ for ${\mathcal{R}}_s>1$ with $\dot{\mathcal{F}}=0$ if and only if $S=S^{**},E_1=E_1^{**},E_2=E_2^{**},$$Q_1=Q_1^{**},Q_2=Q_2^{**},I_1=I_1^{**},I_2=I_2^{**},H_1=H_1^{**},H_2=H_2^{**},$ and $R=R^{**}$. Hence, $\mathcal{F}$ is a Lyapunov function on $\mathcal{D}\backslash {\mathcal{D}}_0.$ Thus, $S\left(t\right)\to S^{**},$$E_1\left(t\right)\to E_1^{**},$$E_2\left(t\right)\to E_2^{**},$$Q_1\left(t\right)\to Q_1^{**},$$Q_2\left(t\right)\to Q_2^{**},$$I_1\left(t\right)\to I_1^{**},$$I_2\left(t\right)\to I_2^{**},$$H_1\left(t\right)\to H_1^{**},$$H_2\left(t\right)\to H_2^{**},$ and $R\left(t\right)\to R^{**}$ as $t\to \infty$. The proof is concluded as in the proof of Theorem 2. Thus, the unique endemic equilibrium of the model (1) with ${\delta }_1={\delta }_2={\delta }_3={\delta }_4=0$ is GAS in $\mathcal{D}\backslash {\mathcal{D}}_0$ whenever $R_s>1$. The epidemiological implication of the above result is that the disease will persist in the community (with the use of isolation and quarantine) if threshold quantity (${\mathcal{R}}_s$) exceeds unity. Numerical simulation results, done in Figure 3 show convergence to the EEP for the case when ${\mathcal{R}}_0>1.$

5. Conclusions

In this paper, a new two stages quarantine/isolation model with a nonlinear incidence rate is designed and rigorously analyzed. The model, which consists of ten mutually exclusive epidemiological compartments, uses the Holling type II incidence function for the infection rate. Some of the theoretical findings of the study are the following:

(i)

The model (1) has a locally asymptotically stable DFE if the associated reproduction number (${\mathcal{R}}_0$) is less than one.

(ii)

The model (1) has a GAS whenever ${\mathcal{R}}_0<1$.

(iii)

System (1) is uniformly persistent in $\mathcal{D}$ if and only if the reproduction number exceeds unity.

(iv)

The model has a unique endemic equilibrium whenever ${\mathcal{R}}_0>1$.

(v)

The unique endemic equilibrium of the model is shown to be GAS for a special case.