Mathematical analysis of a COVID-19 model with different types of quarantine and isolation

1.   Introduction

The COVID-19 pandemic first initiated in December 2019 in the city of Wuhan, China ^[1], and since then, it has been spreading worldwide. The coronavirus disease 19 is caused by severe acute respiratory syndrome coronavirus 2 (SARS- CoV-2) ^[2,3]. On February 24, 2020, the Oman Ministry Of Health registered the first two COVID-19 cases for Omani women associated with travel. On January 30, 2020, the World Health Organization (WHO) declared the outbreak a public health emergency of international concern ^[4]. On March 11, 2020, it characterized the outbreak as a pandemic ^[5].

Governments around the world regularly respond to the announcements made by the WHO, and different medical protocols are followed depending on the epidemiological status. On the same day of declaring COVID-19 as a pandemic, the Omani government responded and His Majesty The Sultan issued his royal orders to set up a supreme committee for dealing with the developments resulting from the spread of COVID-19 ^[6]. Since then, the Supreme Committee has made many decisions and regulations to help monitor and protect the country's healthcare system. As is the case in many countries, the main methods adopted for controlling the spread of COVID-19 have been a lockdown, quarantine, isolation and creating awareness among individuals to adopt social distancing and maintain personal hygiene. By the end of December 2020, some types of vaccines against COVID-19 had been approved and distributed ^[7]. However, authorities still recommend that people practice non-pharmaceutical intervention strategies (such as social distancing and wearing masks) to avoid and fight against the spread of COVID-19. The timeline chart in Figure 1 summarizes the events mentioned above.

Many researchers proposed mathematical models to study the effectiveness of the various strategies and interventions in controlling the spread of COVID-19 before any effective vaccine had been developed. Aronna, Guglielmi and Moschen ^[8] proposed a compartmental model for the dynamics of COVID-19. They took into account the presence of asymptomatic infections, and also the main policies that had been adopted by several countries to fight the disease, these being: isolation, quarantine and testing. They modeled isolation by separating the population into two groups: one composed of essential workers that keep working during the pandemic, and the other group consisting of people that were recommended to stay at home. They showed that people that remain in isolation significantly reduce their probability of contagion, so risk groups should be recommended to maintain a low contact rate during the course of the epidemic.

Džiugys et al.^[9] developed a model under quarantine conditions and methods to estimate quarantine effectiveness. It is based on the daily growth rate of new infections. Using collected epidemiological data of the COVID-19 pandemic, they found that the daily growth rate of new infections has a tendency to decrease linearly when the quarantine is imposed until it reaches a constant value.

Varghese et al. ^[10] formulated a deterministic compartmental model including various stages of infection to study the spreading of COVID-19; they estimated the model parameters by fitting the model with the reported data in Oman. They proved the model's steady state, stability, and final pandemic size mathematically. In addition, they performed a sensitivity analysis to identify the key model parameters and found that contact with symptomatic moderately infected is more significant in spreading the disease.

For more compartmental models based on ordinary differential equations we refer the reader to ^{[11,12,13,14,15,16,17]}. For example, the authors of ^[11] developed a mathematical model to investigate the role of the diagnosis rate in the transmission dynamics of COVID-19 together with the combined effects of quarantine and isolation. They concluded that quarantine and diagnosis, especially for the presymptomatic humans, play the most crucial role in controlling the disease transmission.

The researchers in ^[18] found that if isolation can be implemented effectively (high efficacy and coverage), then the quarantine of people suspected of contracting COVID-19 may not be necessary. Many studies with mathematical models based on fractional differential equations also investigated COVID-19 transmission dynamics. Oud et al. ^[19] considered a fractional model to explore the transmission dynamics and possible control of the COVID-19 pandemic in Pakistan. Initially, they examined a model without optimal control variables. It provided a good fit to the reported cases and then estimated the model parameters using the non-linear least square curving fitting approach. They further reformulated the model by adding two time-dependent control variables. They observed that the most effective strategy to minimize the disease is the implementation of maintaining strict social distancing and contact tracing to quarantine the exposed people. Other studies can be found in References ^{[20,21,22,23]}. The researches in ^[24] found that fractional models better predict future dynamics and give closer estimates to the real data than integer-order models.

In this paper, we propose and analyze a mathematical model to study the roles of different types of quarantine and isolation in COVID-19 transmission dynamics. In particular, we will consider short- and long-term quarantines and home and institutional isolations, taking Oman as a case study and utilizing its COVID-19 data for model parameter estimations. Although the literature is rich in compartmental models, and the effects of quarantine and isolation on COVID-19 spread are well established, our model attempts to reflect the situation during the period of time when data were collected for fitting. The model aims to investigate the roles of different types of quarantine and isolation on the dynamics of COVID-19. Hence our main contribution is the inclusion of corresponding classes of quarantine and isolation to capture the effects of the related control strategies that were implemented by Oman's health authorities on disease dynamics.

The paper is organized as follows. The model formulation is described in Section 2. The model analysis includes the positivity, boundedness of the solution, control reproduction number and existence and stability analysis of disease-free and endemic equilibrium points, which are discussed in Sections 3 and 4. Furthermore, bifurcation analysis and sensitivity analysis are also presented. The fitting of the model parameters and numerical analysis of the model using estimated parameters are given in Section 5 to illustrate the effects of quarantine and isolation on the spread of the disease. Finally, the conclusion is given in Section 6.

2.   Model formulation

The formulation of this model is based on classifying the human population $N(t)$ according to their health status into 11 classes. These classes are susceptible $S(t)$, exposed $E(t)$, asymptomatic $A(t)$, presymptomatic $P(t)$, symptomatic $I(t)$, quarantined for a short time (at home or institutionally) $Q_S(t)$, quarantined for a long time (mostly at home) $Q_L(t)$, non-tested home-isolated $J_N(t)$, tested home isolated $J_{H}(t)$, tested isolated at health institutions $J_{I}(t)$ and recovered $R(t)$, thus,

2.1.   Model assumptions and description

Susceptible individuals, $S(t)$, are those at risk of becoming infected. They are recruited into the population at a constant rate $\Lambda$. They acquire COVID-19 by direct contact with asymptomatic, presymptomatic, symptomatic, non-tested and tested home-isolated at rates of $\beta_{{A}}$, $\beta_{{P}}$, $\beta_{{I}}$, $\beta_{{N}}$ and $\beta_{{H}}$, respectively. It is assumed that there is no transmission from the institutional isolation class. This class is assumed to be completely isolated from the susceptible group.

Following governmental regulations regarding travel restrictions, it is assumed that some susceptible individuals undergo a short quarantine, $Q_S(t)$, at a rate of $\rho_{{S}}$ and leave the class by either natural death, returning back to the susceptible class at a rate of $\lambda_S$ after 14 days or by moving to long quarantine, $Q_L(t)$, at a rate of $\alpha_{S}$. Long quarantine is considered a behavioural practice adhered to by some people due to their awareness or fear. Individuals in the long quarantine class leave by natural death or return back to the susceptible class at a rate of $\lambda_L$.

All newly infected individuals enter the exposed class $E(t)$, in which they remain non-infectious during their latent period. After this period, they become infectious at different rates $\lambda_A, \lambda_P$ or $\lambda_I$ as asymptomatic, presymptomatic or symptomatic, respectively.

Asymptomatic individuals, $A(t)$, are infected and contagious but do not show any symptoms until they recover. It is assumed that they leave their class by a natural death, isolation at home for non-tested individuals at a rate of $\rho_{{A}}$, getting diagnosed and isolated at home (health institutions) at a rate of $\varepsilon_A (\theta_A)$ or recovery at a rate of $\gamma_A$.

Presymptomatic infectees, $P(t)$, currently have no symptoms but will develop symptoms after some time. It is assumed that they leave their class due to natural death, tested positive and hence are isolated at homes/health institutions at a rate of $\varepsilon_P /\theta_P$ or move to the infected class with symptoms at a rate of $\sigma_P$ when symptoms are onset. As a result of contact-tracing or human behaviour, some presymptomatic individuals move to the non-tested home isolated class at a rate of $\rho_{{P}}$. It is assumed that there is no disease-related death in the presymptomatic class.

Individuals in the symptomatic class, $I(t)$, are those with mild COVID-19 symptoms but have not been tested. It is assumed that when diagnosed or symptoms become severe, medical care is needed; some are moved to tested home isolation at a rate of $\theta_{I}$, and some are moved to institutional isolation at a rate of $\varepsilon_{I}$, while the remaining stay in the symptomatic class until they recover at a rate of $\gamma_{I}$. No disease-related death is assumed in this class, but natural death may occur.

Non-tested home-isolated individuals, $J_N(t)$, are those individuals with suspected infection, with or without symptoms, who have not been tested. They come from the exposed class at a rate of $\rho_E$, from the asymptomatic class at a rate of $\rho_A$ and from the presymptomatic class at a rate of $\rho_P$. It is assumed that some people have COVID-19 symptoms but do not wish to obtain a diagnosis. Instead, they prefer to isolate themselves in their home. Also, due to identification by contact tracing, some are instructed to isolate themselves at home, and they are usually tracked by electronic track-bands to make sure that they adhere to isolation. Another group tends to self-isolate themselves at home as human behaviour. All of these categories are classified in the $J_N(t)$ compartment.

Depending on their symptoms or random testing, it is assumed that non-tested home isolated individuals leave their class by transferring to tested home (institutional) isolation at a rate of $\varepsilon_{{N}}$($\theta_{{N}}$). Some leave the class after recovery at a rate of $\gamma_{{N}}$ or due to natural death.

The tested home-isolated infectees ($J_H(t)$) leave their class by natural death, disease-related death at a rate of $\delta_H$, moving to institutional isolation at a rate of $\theta_H$ for health care after their symptoms become more severe or recovery at a rate of $\gamma_{H}$.

It is assumed that individuals isolated at health institutions, $J_I(t)$, leave this class due to natural death, recovery at a rate of $\gamma_{D}$, disease-related death at a rate of $\delta_D$ or moving to home isolation at a rate of $\varepsilon_D$ after they have received the required health care. The above description is summarized in the flow chart in Figure 2.

2.2.   Model equations

Based on the above description, the proposed mathematical model is given by the following set of differential equations:

where, $N'(t) = \Lambda - \delta_H J_H(t) - \delta_D J_I(t)- \mu N(t)$. Model state variables are defined in Table 1, and the parameters are defined in Table 2. For simplicity of the calculations, we denote the following:

3.   Positivity and boundedness

In this section, we study the positivity and boundedness of solutions of the system (2.1) subject to the positive initial conditions

where

Theorem 1. Solutions of the model (2.1) with the initial conditions (3.1) are positive for all $t \ge 0.$

Proof Consider the first equation of the model (2.1):

and rewrite the force of the infection term in the form

Then

Let $t_1 = {\rm{sup}}\left(t > 0 \mid {\textbf{U}}(t) > 0 \right)$. We have

Similarly, one can prove that the remaining components are positive at $t = t_1$. Hence, $t_1$ cannot be a supremum. Therefore, all solutions remain positive for all time $t\ge 0$.

Theorem 2. All solutions $S(t), \: E(t), \: A(t), \: P(t), \:I(t), \:Q_S(t), \:Q_L(t), \: J_N(t), \:J_H(t), \:J_I(t) \; {\rm{ and }}\; \:R(t)$ of the model (2.1) with the given initial conditions (3.1) are bounded.

Proof From Theorem 1, solutions of the model (2.1) with the positive initial conditions are positive for $t \ge 0$.

We have that $N = S \: + \: E \: + \: A \: + \:P\: + \: I \: + \: Q_S \: + \: Q_L \: + \: J_N \: + \: J_H \: + \: J_I \: + \:R.$

${\rm{Thus, }} \: \frac{dN}{dt} = \Lambda - \delta_H J_H(t) - \delta_D J_I(t) - \mu N(t) \:\le \: \Lambda - \mu N(t)$.

${\rm{Therefore, }}\; \lim\limits_{t \rightarrow \infty} \: {\rm{sup}} \: N(t) \le \frac{\Lambda}{\mu}, \: \:$ and hence the solutions are bounded. The feasible invariant region is given by

$\varOmega = \left\lbrace \left(S(t), E(t), A(t), P(t), I(t), Q_S(t), Q_L(t), J_N(t), J_H(t), J_I(t), R(t)\right) \in \mathbb{R}_+^{11}: N(t) \le \frac{\Lambda}{\mu} \right\rbrace$.

4.   Disease free equilibrium and control reproduction number

The disease free equilibrium (DFE) of the model is given by

where

Hence

The control reproduction number, denoted by $R_c$, is the number of secondary infections caused by a single infective in a totally susceptible population with control measures in place. It is calculated by using the next generation matrix (NGM) method ^[25] as follows. First, the matrix of new infection is found to be

where, $X = N-Q_{{S}}-Q_{{L}}- J_{I},$ and it has a Jacobian, evaluated at the DFE $E_0$ and $N^* = S^* + {Q_S}^* + {Q_L}^*$ that is given by

Then, the matrix of transition is given by

and its Jacobian is

where $a$, $b$, $c$, $d$, $e$ and $f$ are as defined earlier after the model equations.

Now, NGM is given by

where

The control reproduction number is the spectral radius of the NGM, i.e., $R_c = \rho(FV^{-1})$ and it is given by

where

Here, $R_A, \: R_P, R_I, R_N,$ and $R_H$ represent the contribution from asymptomatic, presymptomatic, symptomatic, non-tested home isolated and tested home isolated individuals, respectively.

Note that the reproductive number $(R_c)$ is independent of short and long quarantines parameters. In fact, according to the assumptions in formulating the model, short and long quarantines have no direct impact on the control reproductive number ($R_c$). However, they play an essential and considerable role in separating susceptible humans from infected individuals and reducing the contact rate between them. Hence both types of quarantines contribute indirectly to controlling the spread of COVID-19.

4.1.   Stability analysis of the disease free equilibrium

Studying the local stability of the DFE using linearization is complicated since the characteristic equation is of order 11. Therefore, global stability is checked instead by using the Lyapunov function method as described in ^[26].

Set

where $\mathcal{F}, \: F, \: \mathcal{V}, \; {\rm{ and }}\; V$ as previously defined in (4.1), (4.2), (4.3) and (4.4), respectively, $x^T = [E \quad A \quad P \quad I \quad J_N \quad J_H \quad J_I] \in \mathbb{R}_+^{7}$ is the disease compartment, and $y^T = [S \quad Q_S \quad Q_L \quad R] \in \mathbb{R}_+^{4}$ is the non-disease compartment.

Then, the disease compartment can be written as $x' = (F-V)x - f(x, y)$ with $f(0, y) = 0$.

Let $w^T \ge 0$ be the left eigenvector of the non-negative matrix $V^{-1}F$ corresponding to the eigenvalue $R_c$. The vector is found to be $w^T = [0 \quad w_A \quad w_P \quad w_I \quad w_N \quad w_H \quad 0]$, where all of the entries are non-negative.

The following theorem provides a general method to construct a Lyapunov function for the system (2.1).

Theorem 4.1. ^[26] Let $F$, $V$ and $f(x, y)$ be as defined above. If $f(x, y) \geq 0$ in $\varOmega \subset \mathbb{R}_+^{11}$, $F \geq 0,$ $V^{-1} \geq 0$ and $R_c\leq 1$, then the function $Q = w^T V^{-1} x$ is a Lyapunov function for the model (2.1) on $\varOmega$.

Proof We need to prove that the system has a Lyapunov function $Q$ and that $Q' \leq 0$ on $\varOmega \subset \mathbb{R}_+^{11}$.

We have

Now, one can show that $1- \frac{S}{N - Q_S - Q_L - J_I} \ge 0$ as follows:

Differentiating the Lyapunov function gives

since $w^T \geq 0$, $V^{-1} \geq 0$ and $f(x, y) \geq 0.$ Hence $Q$ is a Lyapunov function for the model (2.1) on $\varOmega \subset\mathbb{R}_+^{11}$; thus, the DFE is asymptotically stable.

Furthermore, $Q' = 0$ only at the DFE ($E_0$). Hence the largest invariant set contained in the set $\varOmega$ is reduced to the DFE. Since we are in a compact positively invariant set, by the LaSalle's invariance principle ^[27], the DFE is globally asymptotically stable in $\varOmega$.

4.2.   Existence of endemic equilibrium

In this section, we investigate the existence of an endemic equilibrium by putting the right-hand side of the system (2.1) equal to zero and solving the equations simultaneously. We get

where

Note that $N'(t) = \Lambda - \delta_H J_H(t) - \delta_D J_I(t) - \mu N(t)$. Therefore, the equilibrium state $N'(t) = 0$ implies that $\hat{N}(t) = \frac{1}{\mu} \left(\Lambda - \delta_H \hat{J_H}(t) - \delta_D \hat{J_I}(t) \right)$.

To find a relation between $\hat{S}$ and $\hat{E}$, the first two equations of the model can be written as follows:

where

Solving Eq (4.6) for $\hat{E}$ gives

Solving Eq (4.7) also for $\hat{E}$ gives

Solving Eqs (4.8) and (4.9) together gives

Substituting the expression of $\hat{S}$ from (4.10) into (4.9) gives

Note that $a_1 = \mu \left(a_2 +1\right)$; therefore, $\hat{S}$ and $\hat{E}$ can be written in the form

Note that the expressions $a, \: m, \: a_2$ and: $T$ are positive. One can also show that $\left(a \: m - T\right)$ is positive. Now, we check the signs of both $\hat{E}$ and $\hat{S}$ for the following cases:

(Ⅰ) $R_c > 1$: the denominators of $\hat{E}$ and $\hat{S}$ can be shown to be positive as follows:

Clearly, both of the expressions of $\hat{E}$ and $\hat{S}$ are positive, resulting in an endemic equilibrium $E_e = \left(\hat{S}, \hat{E}, \hat{A}, \hat{P}, \hat{I}, \hat{Q_S}, \hat{Q_L}, \hat{J_N}, \hat{J_H}, \hat{J_I}, \hat{R}\right)$.

(Ⅱ) $R_c = 1$: one can easily verify that the expressions $(4.12)$ and $(4.13)$ are reduced to the DFE ($E_0$) since $a \: m \left(R_c + a_2\right) - \left(a_2 + 1\right)T > 0.$

(Ⅲ) $R_c < 1$: $\hat{S}$ and $\hat{E}$ are always of different signs as it is clear from the expressions (4.12) and (4.13). Therefore, there is no endemic equilibrium.

4.3.   Bifurcation analysis

This section uses the center manifold theorem, as described in ^[28], to check if any type of bifurcation exists.

To apply the theorem, we rewrite the system (2.1) in the following formulation:

Take $\beta_{{I}} = \phi$ to be the bifurcation parameter and $\beta_{{I}}^* = \phi^*$ to be the corresponding bifurcation value at $R_c = 1$. Then, $R_c$ can be written as

where $K = R_A +R_P + R_N + R_H$.

At $R_c = 1$, the bifurcation value is

System (4.14) has the following Jacobian at the DFE $E_0$ and $\beta_I = \phi^*$:

One can prove that the constant term of the characteristic equation of the above matrix is zero. Therefore, at least one of the eigenvalues is zero. The right and left eigenvectors ($w$ and $v$, respectively) corresponding to the simple eigenvalue $0$ are: $w = [w_1 \qquad w_2 \qquad... \qquad w_{11}]^T$, and $v = [v_1 \qquad v_2 \qquad... \qquad v_{11}]^T$, where,

and

We compute the following second-order derivatives:

where $f_i, i = 1, 2, \cdot \cdot \cdot \: 9$ denotes the right-hand side of the equation number $i$ of the system (4.15). Then, we have the following second-order derivatives:

${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_2 \partial x_9 }} = { \frac {\beta_{{H}}}{S^*} }$, ${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_2 \partial x_8 }} = { \frac {\beta_{{N}}}{S^*} }$, ${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_2 \partial x_5 }} = { \frac {\phi}{S^*} }$, ${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_2 \partial x_4 }} = { \frac {\beta_{{P}}}{S^*} }$, ${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_2 \partial x_3 }} = { \frac {\beta_A}{S^*} }$,

${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_3 \partial x_{11} }} = { \frac {\beta_{{A}}}{S^*} }$, ${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_3 \partial x_9 }} = { \frac {\beta_{{A}} +\beta_{{H}}}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_3 \partial x_8 }} = { \frac {\beta_{{A}} + \beta_{{N}}}{S^*}}$,

${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_3 \partial x_5 }} = { \frac {\beta_{{A}} +\phi}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_3 \partial x_4 }} = { \frac {\beta_{{A}} + \beta_{{P}}}{S^*}}$, ${ \frac{ {d }^{2}{\it f_1}}{d {x_3}^2 }} = { \frac {2\beta_{{A}}}{S^*}}$,

${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_4 \partial x_{11} }} = { \frac {\beta_{{P}}}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_4 \partial x_9 }} = { \frac {\beta_{{P}} + \beta_{{H}}}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_4 \partial x_8 }} = { \frac {\beta_{{P}} + \beta_{{N}}}{S^*}}$,

${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_4 \partial x_5 }} = { \frac {\beta_{{P}}+\phi}{S^*}}$, ${ \frac{ {d }^{2}{\it f_1}}{d {x_4}^2 } } = { \frac {2\beta_{{P}}}{S^*}}$,

${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_5 \partial x_{11} }} = { \frac {\phi}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_5 \partial x_9 }} = { \frac {\phi+ \beta_{{H}}}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_5 \partial x_8 }} = { \frac {\phi+\beta_N}{S^*}}$,

${ \frac{ {d }^{2}{\it f_1}}{d {x_5}^2 }} = { \frac {2\phi}{S^*}}$,

${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_8 \partial x_{11} }} = { \frac {\beta_{{N}}}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_8 \partial x_9 }} = { \frac {\beta_{{N}}+\beta_{{H}}}{S^*}}$, ${ \frac{ {d }^{2}{\it f_1}}{d {x_8}^2 }} = { \frac {2\beta_{{N}}}{S^*}}$,

${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_9 \partial x_{11} }} = { \frac {\beta_{{H}}}{S^*}}$, ${ \frac{ {d }^{2}{\it f_1}}{d {x_9}^2 }} = { \frac {2\beta_{{H}}}{S^*}}$,

${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_2 \partial x_9 }} = { \frac {-\beta_{{H}}}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_2 \partial x_8 }} = { \frac {-\beta_{{N}}}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_2 \partial x_5 }} = { \frac {-\phi}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_2 \partial x_4 }} = { \frac {-\beta_{{P}}}{S^*}}$,

${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_2 \partial x_3 }} = { \frac {-\beta_A}{S^*}}$,

${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_3 \partial x_{11} }} = { \frac {-\beta_A}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_3 \partial x_9 }} = { \frac {-\beta_A - \beta_{{H}}}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_3 \partial x_8 }} = { \frac {-\beta_A - \beta_{{N}}}{S^*}}$,

${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_3 \partial x_5 }} = { \frac {-\beta_A - \phi}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_3 \partial x_4 }} = { \frac {-\beta_A - \beta_{{P}}}{S^*}}$, ${ \frac{ {d }^{2}{\it f_2}}{d {x_3}^2 }} = { \frac {-2\beta_A}{S^*}}$,

${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_4 \partial x_{11} }} = { \frac {-\beta_{{P}}}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_4 \partial x_{9} }} = { \frac {-\beta_{{P}} - \beta_{H}}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_4 \partial x_{8} }} = { \frac {-\beta_{{P}} - \beta_{N}}{S^*}}$,

${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_4 \partial x_{5} }} = { \frac {-\beta_{{P}} - \phi}{S^*}}$, ${ \frac{ {d }^{2}{\it f_2}}{d {x_4}^2 }} = { \frac {-2\beta_{{P}}}{S^*}}$,

${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_5 \partial x_{11} }} = { \frac {- \phi}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_5 \partial x_{9} }} = { \frac {-\phi - \beta_{H}}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_5 \partial x_{8} }} = { \frac {-\phi - \beta_{N}}{S^*}}$,

${ \frac{ {d }^{2}{\it f_2}}{d {x_5}^2 }} = { \frac {-2\phi}{S^*}}$,

${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_8 \partial x_{11} }} = { \frac {- \beta_{N}}{S^*}}$, ${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_8 \partial x_{9} }} = { \frac {-\beta_{{N}} - \beta_{H}}{S^*}}$, ${ \frac{ {d }^{2}{\it f_2}}{d {x_8}^2 }} = { \frac {-2\beta_{{N}}}{S^*}}$,

${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_9 \partial x_{11} }} = { \frac {- \beta_{H}}{S^*}}$, ${ \frac{ {d }^{2}{\it f_2}}{d {x_9}^2 }} = { \frac {-2\beta_{{H}}}{S^*}}$,

${ \frac{ {\partial }^{2}{\it f_1}}{\partial x_5 \partial \phi }} = -1$, ${ \frac{ {\partial }^{2}{\it f_2}}{\partial x_5 \partial \phi }} = 1$.

The rest of the second-order derivatives are all zero.

Calculating the coefficients $\widetilde{a}$ and $\widetilde{b}$ as defined in Theorem 4.1 ^[28] by Castillo Chavez and Song gives

Obviously, $\widetilde{a} < 0$ and $\widetilde{b} > 0$. Therefore, the system has a forward transcritical bifurcation.

Remark: Since the system has a forward transcritical bifurcation, the endemic equilibrium is locally stable when $R_c > 1$ and unstable when $R_c < 1$.

Numerical simulation of the model (2.1) using the fitted parameters also illustrated that the system undergoes transcritical bifurcation. The MatCont package ^[29] was used to study the model's behaviour using the transmission rate from symptomatic individuals $\beta_{{I}}$ as the bifurcation parameter. Figure 3 shows how the stability of the equilibrium state changes as $\beta_{{I}}$ varies. It is clear that before the bifurcation point (where $\beta_{I}^* \approx 0.287$ and $R_c = 1$), the DFE is stable, and then the stability state bifurcates to the endemic equilibrium branch after the bifurcation parameter passes the bifurcation point $\beta_{I}^*$. This indicates that for a high transmission rate $\beta_{{I}} > \beta_{I}^*$, the disease will persists in the population.

It is worth noticing that the expression derived for the bifurcation point in (4.15) matches exactly with its numerical value ($\beta_{I}^* \approx 0.28707$) when calculated using the fixed and fitted parameters.

5.   Fitting of model parameters and numerical results

5.1.   Fitting of model parameters

The proposed model (2.1) incorporates 36 parameters; some of them are available in literature and the remaining will be estimated using the daily reported cases for Oman. Data were collected for the period from August 4, 2020, to October 31, 2020 ^[30]. This period of time has been chosen for fitting because no vaccine was introduced at that time, so the effects of quarantine and isolation, in Oman, are better investigated. We employed the non-linear least square curve fitting approach ^[31] to fit the model to cumulative cases of the 89 days by using the Matlab optimization function fminsearchbnd. To improve the fitting and get the best of it, we set the initial guess as found in the literature for those parameters that are already available in references, as shown in Table 3. We also set lower and upper bounds for each parameter depending on the estimation of some biological and precautionary measures criteria. The fixed parameters shown in Table 3 are based on Oman population statistics and a previous study ^[10] about COVID-19 in Oman. The remaining 24 parameters were fitted and are listed in Table 3. The best fit curve of the model to the reported data is shown in Figure 4.

The graphs in Figure 5 represent the predicted situation of the pandemic for months following the three months used to estimate the parameters. Figure 5(a) shows a good prediction until the end of December 2020. However, the forecast is not fine when using more following months, as shown in Figure 5(b). This is because a new wave of the disease began at the end of December 2020.

Using the fitted parameters, the approximate value of the reproduction number is $R_c\approx1.202$ and the contribution from each transmission route is presented in the table below

The significant contribution to disease transmission comes from symptomatic and presymptomatic classes. These results agree with the findings in the review studies ^[35,36] and the references therein.

5.2.   Sensitivity analysis

This section is devoted to the use of sensitivity analysis to study the influence of different models' parameters on the spread of COVID-19. The sensitivity index ($\Upsilon_{\varphi}^{R_{c}}$) was measured when varying the reproduction number ($R_c$) with respect to the models' parameters by using the formula: $\Upsilon_{\varphi}^{R_{c}} = \frac{\partial R_c}{\partial \varphi } \frac{\varphi}{R_c}$.

The results are presented in Figure 6 and Table 5. Note that these results are based on using the set of estimated parameters in Table 3 as the baseline.

It is clear that the transmission rate from symptomatic individuals ($\beta_I$) has the most impact on disease spreading, followed by the transmission rate from presymptomatic ($\beta_P$), then the rate at which exposed become presymptomatic ($\lambda_P$). Higher rates of $\lambda_P$ indicate that it takes a shorter time for exposed ones to be infectious as preysmyptomatic, leading to higher transmission. On the other hand, the rate at which symptomatics are tested and moved to home isolation ($\varepsilon_{ I }$) has the best influence on reducing the spread of disease, followed by the rate at which the exposed are home-isolated without testing ($\rho_E$) due to identification by contact tracing or behavioral self-isolation. In addition, the rate at which the presymptomatic class develops symptoms ($\sigma_P$) has a good influence on reducing the spread. The shorter it takes them to become symptomatic, the better it is to recognize them and control the spread. These findings can help policymakers in taking necessary action to prevent the disease from wide-spreading, keeping in mind the limitation of the study.

5.3.   Numerical simulation

This section presents some numerical solutions for the model (2.1) based on the fitted parameters and ODE MATLAB solver ode45. The focus will be on the effects of different parameters on the model dynamics.

In all graphs below, the black curve represents the fitted curve, while the others are used to illustrate the dynamics when a specific parameter is varied or controlled. The peaks of symptomatic and presymptomatic cases of the fitted curves are 1200 and 545 cases, respectively, in all graphs.

As shown in Figure 7, the transmission rates from presymptomatic ($\beta_P$) and symptomatic ($\beta_I$) classes have the most considerable impact on spreading the disease, while the transmission rates from asymptomatic ($\beta_A$) and non-tested home isolated ($\beta_N$) have the least impact. The transmission rate from the tested home isolated group ($\beta_{{H}}$) has a moderate effect.

To compare the strength of transmission from different groups, we take, for example, Figure 7(a), (b), which show the effects of the transmission rates from asymptomatic and presymptomatic, respectively. Increasing or decreasing the transmission rate from asymptomatic ($\beta_A$) by a factor of 5 has a negligible effect on the number of infected cases. On the other hand, increasing or decreasing the transmission rate from presymptomatic people ($\beta_{P}$) by a significantly smaller factor of 1.1 dramatically impacts the number of infected cases. This effect is clear from the curves in Figure 7(b), where disease prevalence increases significantly with a slight increase in $\beta_P$. The behaviour mentioned above for $\beta_P$ applies to the transmission rate from the symptomatic group ($\beta_I$) when it is varied by a factor of 1.05, as shown in Figure 7(c). Figure 7(d), (e) illustrates the dynamics when varying the transmission rates from non-tested and tested home isolated groups, respectively, by a factor of 2.

It is interesting to note that controlling the transmission rates from classes of strong influence can flatten the curve of infection, but the time it takes to reach the peak remains almost the same.

These simulations can hint at concentrating efforts to limit the disease spread by controlling the parameters of high impacts. In this current case, it is helpful to identify presymptomatic individuals by contact tracing and increasing testing intensity.

Figure 8 shows the effects of short and long quarantines on disease dynamics. Although both types of quarantine have no direct influence on the reproductive number, they still considerably impact the disease transmission dynamics. As the more susceptible are quarantined for a short or long time (at rates $\rho_{S}$ and $\alpha_L$), the better it is for controlling the spread, and the peak of the infection curve will decrease. This behavior is clear from Figure 8(a) and (b). On the other hand, when more people leave their short or long quarantines (at rates $\lambda_S$ and $\lambda_L$) return to the susceptible class, more infections occur, and the peak becomes higher, as shown in Figure 8(c) and (d). The current simulation shows that moving from short to long quarantine helps with reducing the infection when the rate of moving ($\alpha_S$) increases, as one can conclude from Figure 8(e).

The curves in all of the graphs are plotted when each parameter is changed by a factor of 1.5. Long quarantine has more influence on the spread of the disease, as it can be seen from Figure 8(b), (d).

Note that quarantine helps to lessen the number of infections and flatten the curve, though it speeds up the occurrence of the peak of infection.

The acts of instituting full or partial lockdown policies from time to time and putting a significant portion of communities under quarantine can mitigate disease spreading.

Figure 9 shows the effects of isolating symptomatic infectees institutionally or at home on disease spreading. It is clear from the graphs that increasing the isolation rates helps with reducing the intensity of infection and flattening the curve. From Figure 9(b), varying the rate of tested home isolation ($\varepsilon_I$) by a factor of 1.2 has a great influence on the disease the dynamics compared to dynamics resulting from varying the rate of institutional isolation ($\theta_I$) by the same factor as shown in Figure 9(a).

Policymakers can utilize such observations, keeping in mind the limitation of the study, by enhancing home isolation over institutional isolation, where the latter is more costly and less effective. One may wonder how institutional isolation is less effective than home isolation, given that the first is more systematic and organized in a way that contact between isolated individuals is not likely to happen. The answer is that the number of isolations (whether at hospitals or other isolation centers) is limited and very small compared to the number at homes. Therefore, putting more effort into controlling and regulating home isolation is more economical and will help reduce infected cases.

Figure 10 illustrates the effects of isolating presymptomatic individuals on disease transmission dynamics. It can be concluded from the graphs that increasing the rates of isolation of presymptomatic individuals results in decreasing the overall infected cases. However, the presymptomatic group's challenge lies in the difficulty of identifying them, since they do not show any symptoms. Thus, increasing the intensity of random testing and contact tracing and promoting self-isolation are essential to control the spread from this group.

When comparing the three types of isolation, note that changing the rate of non-tested home isolation ($\rho_P$) by a relatively small factor of 1.1, as shown in Figure 10(a), results in a considerable variations in the number of infections. On the other side, to get a comparable effect for the rates of institutional isolation ($\theta_P$) and tested home isolation ($\varepsilon_P$), it is required to vary them by a higher factor of 2 as shown in Figure 10(b), (c). Hence, it is very crucial to deal with the presymptomatic group who are not identified through random testing. This mainly depends on human behaviour, so raising people's awareness on taking precautions and isolating themselves once there is a suspection of infection will definitely help containing the disease.

Figure 11 demonstrates the effects of isolating exposed individuals at home without testing on the number of symptomatic (Figure 11(a)) and presymptomatic (Figure 11(b)) cases. The isolation of exposed individuals could be due to contact tracing, travel regulations or a human behavior in response to awareness programs. The impact of varying the isolation rate $\rho_E$ by a factor of 1.2 can be clearly seen in Figure 11, where the isolation rate is altered only by a factor of 1.2, but the variation in the maximum numbers of infections for both symptomatic and presymptomatic individuals is significant. Therefore, employing programs targeted to enhance human behavior toward dealing with the pandemic and utilizing technology for contact tracing will help control the fast spread of the disease.

6.   Conclusions

A mathematical model for COVID-19 with different types of quarantine and isolation has been proposed, namely, short- and long-term quarantine, home isolation with or without testing and institutional isolation. The model has been fully analyzed to study the effects of the different types of quarantine and isolation on the disease transmission dynamics. The positivity and boundedness of solutions to the proposed model have been discussed. Model parameters were estimated by using the Oman dataset, and the best fit curve of the model to the reported data was illustrated graphically. The control reproduction number, $R_c$, was calculated using the NGM method and expressed as a sum of terms representing the contribution from all transmission routes. Using the obtained parameter values, the control reproduction number was estimated to be $R_c = 1.202$, with major contributions coming from symptomatic and presymptomatic transmissions. Sensitivity analysis for $R_c$ (to model parameters) was carried out, and it was found that the symptomatic transmission rate, $\beta_I$, has the strongest positive impact on $R_c$, followed by the presymptomatic transmission rate, $\beta_P$. On the other hand, the home isolation rate for tested symptomatic individuals, $\varepsilon_I$, was found to have the strongest negative impact on $R_c$, followed by the home isolation rate for exposed individuals, $\rho_E$. The global asymptotic stability of the DFE was shown by using an appropriate Lyapunov function together with LaSalle's invariance principle whenever $R_c \le 1$. The endemic equilibrium was found to exist if $R_c > 1$, and its local stability was deduced from the bifurcation analysis, since the system was found to have a forward transcritical bifurcation as a result of using the center manifold theorem and taking $\beta_I$ to be the bifurcation parameter. The bifurcation results have been demonstrated graphically as well. Finally, the effects of some model parameters related to the different types of quarantine and isolation on the disease transmission dynamics were investigated numerically, and the results have been demonstrated graphically. The numerical results show that quarantine reduces the number of infected individuals and flattens the curve of infections, i.e., slowing the spread of the disease, with long quarantine having more impact on the disease transmission dynamics than short quarantine. Out of all of the different types of isolation, home isolation for tested symptomatic, non-tested exposed and non-tested presymptomatic individuals has the most significant impact on the disease transmission dynamics. These results agree very well with the obtained sensitivity results. Hence, adopting regulations for home quarantine and isolation and enhancing the culture of self-isolation through awareness programs are essential in controlling the spread of COVID-19.

This paper provides a deep insight into COVID-19 dynamics when many types of quarantine and isolations are employed. Vaccination can play an important role in controlling transmission and, accordingly, impacts the regulations of quarantine and isolation. However, it is not considered in this study. We will discuss this in future work.

Conflict of interest

The authors declare that there is no conflict of interest.