Stability and Sensitivity Analysis of the COVID-19 Spread with Comorbid Diseases

Abstract

This research investigates a model of the spread of COVID-19 in Indonesia by paying attention to comorbid disease, self-quarantine, government-provided quarantine, and vaccination factors. The symmetrical aspects of the model are studied. The evaluation of the model reveals non-endemic and endemic equilibrium points and the basic reproduction number (BRN). We provide the local and global stability analysis of the equilibriums. According to the sensitivity analysis of the BRN, the key parameters impacting the spread of COVID-19 are the susceptible recruitment rate, contact rate, infection death rate, and probability of infected individuals having no comorbidities. In addition, we provide a sensitivity analysis to examine the effect of parameter changes in each subpopulation. We discovered that the natural death rate is the most sensitive parameter based on the sensitivity index after reaching equilibrium. Symmetry aspects appear in some of the visualizations of the model’s solution and the sensitivity of the BRN and parameters.

1. Introduction

COVID-19 is an illness caused by the severe acute respiratory syndrome coronavirus 2 infection (SARS-CoV-2). The illness can induce a range of respiratory system problems, from mild symptoms like the flu to severe lung infections like pneumonia [1]. Antiviral medicines and vaccines for COVID 19 prevention were identified in 2021. In order to combat COVID-19, health professionals and mathematicians can work together to analyze, evaluate, anticipate, and optimize the prevention strategies. For example, mathematical models can be utilized to analyze COVID-19 prevention with less effective vaccines [2], to study the vaccine and treatment for the disease [3,4], and to examine the optimal control to combat the spread of the disease [5,6,7,8,9].

There are many approaches for modeling the COVID-19 spread by considering various factors. Ahmed et al. [10] analyzed the COVID-19 model with a numerical analysis and logistic model. Arino and Portet [11] proposed a COVID-19 model by taking into account susceptible, latently infected, symptomatic and asymptomatic infectious, and recovered individuals. Arfan et al. [12] investigated a fractional dynamical system of COVID-19 based on a modified susceptible–exposed–infectious–recovered type with a nonsingular kernel derivative. Arfan et al. [13] examined the non-integer-order population dynamical model of the COVID-19 pandemic in relation to various values of the heavily impacted system parameter of immigration for both susceptible and infected populations. Sugiyanto and Abrori [14] studied the model of COVID-19 spread by considering the case of with and without lockdown. Khan et al. [15] considered the model of COVID-19 dynamics where the infected subpopulations have isolation and quarantine.

Several researchers have studied models of the COVID-19 spread in different parts of the world. Annas et al. [16] analyzed and provided numerical simulations for the model of COVID-19 spread of the susceptible–exposed–infectious–recovered type in Indonesia. Fitriani et al. [17] studied the dynamic of COVID-19 transmission in Central Java, Indonesia. Anggriani et al. [18] investigated the model of COVID-19 propagation in West Java, Indonesia. Okuonghae and Omame [19] reviewed the model of the COVID-19 contagion in Lagos, Nigeria. Kifle and Obsu [20] considered the model of COVID-19 infectious in Ethiopia between March 2020 and July 2021. Asempapa et al. [21] studied the model of COVID-19 spread in Brazil and South Africa among populations at risk of death through the provision of natural medicine preventatives.

Various attempts have been made to estimate the data of COVID-19 cases in order to predict its spread. Postnikov [22] studied the criteria for predicting the COVID-19 spread using the susceptible–infectious–recovered-type model. Using data provided in China, Macau, Hong Kong, and Taiwan, Ivorra et al. [23] determined the characteristics of the COVID-19 model which is used to predict the disease transmission in other countries and also in Wuhan, China, as used by Yang and Wang [24]. Barmparis and Tsironis [25] provided an estimation technique of COVID-19 distribution parameters using data from nine nations. Sreeramula and Rahardjo [26] presented an estimation of COVID-19 propagation parameters in Indonesia.

Modeling COVID-19 transmission in relation to other diseases can give a rich examination of the infectious dynamics characteristics. Mekonen et al. [27] evaluated the model for the spread of tuberculosis coinfected with COVID-19. Ssebuliba et al. [28] studied the COVID-19 contagion by considering the infected individuals having partially comorbid diseases. In [29,30], the COVID-19 spread was examined where the individuals with comorbid diseases had interventions of natural remedies and vaccination. Okyere and Ackora-Prah [31] analyzed the COVID-19 transmission dynamics by taking into account diabetes.

In the dynamic model of COVID-19 spread studied by Ssebuliba et al. [28], individuals infected with COVID-19 with or without comorbid diseases are shifted to the hospital (local government quarantine) without reference to the self-quarantine in their respective homes. On the basis of the Ssebuliba et al. [28] model, in this paper, we extend the model by taking into account the transfer rate from individuals infected with COVID-19 without comorbid disease to the self-quarantine individuals. This is important to note because, according to real data from Indonesia, individuals infected with COVID-19 without comorbid diseases are generally isolated in their respective homes because they can recover by consuming nutritious foods, taking vitamins, getting enough rest, exercising, and adhering to health protocols [1]. The symmetrical aspects for the dynamical system are analyzed, such as the non-endemic and endemic equilibrium points and their stability and also the basic reproduction number. We perform some numerical simulations to visualize the solution of the system, the sensitivity analysis of the basic reproduction number to study the stability of the system, and the sensitivity of the model’s parameters to study their effects on the solution of the system. The symmetry appears in some of these visualizations. Based on the sensitivity index after reaching equilibrium, we determine which parameter is the most sensitive.

2. The Spread of COVID-19 with Comorbid Diseases

Consider a population which is divided into six subpopulations. S denotes the susceptible subpopulation which consists of individuals who are still healthy but susceptible to infection of COVID-19. $I_h$ denotes the subpopulation which is infected with COVID-19 but has no comorbid diseases. $I_c$ denotes the subpopulation which is infected with COVID-19 and has a history of comorbid diseases, such as hypertension, diabetes mellitus, heart disease, respiratory disorders, chronic renal disease, neurological disorders, endocrine disorders, and liver. $Q_h$ represents the self-quarantine subpopulation which consists of persons infected with COVID-19 without comorbid diseases. $Q_H$ represents the subpopulation of COVID-19-infected persons with comorbid diseases who are isolated in a local government-supplied quarantine. R represents the recovered subpopulation which consists of individuals who have recovered from or are immune to COVID-19.

The following are the model’s underlying assumptions: (1) New recruits are introduced into the susceptible subpopulation. (2) Individuals from the susceptible subpopulation are immunized against COVID-19 and enter the recovered subpopulation. (3) Susceptible subpopulation who come in contact with COVID-19-infected individuals become infected. (4) COVID-19-infected individuals who do not have comorbid diseases are self-quarantined in their residences. (5) Individuals infected with COVID-19 and comorbid diseases are quarantined in a local government-provided facility. (6) The subpopulation died of natural causes. (7) The majority of individuals who received treatment or supplements healed. (8) Those infected with COVID-19 may perish. A schematic representation of the spread of COVID-19 by paying attention to comorbid diseases is depicted in Figure 1.

On the basis of the assumptions and the schematic diagram in Figure 1, the following system of differential equations is derived,

$$\left\{\begin{array}{cc}\hfill \frac{dS}{dt}& =\mathsf{\Lambda }-\frac{\beta S(I_h+I_c)}{N}-(p_1\nu +\mu )S\hfill \\ \hfill \frac{d{I}_h}{dt}& =\frac{p\beta S(I_h+I_c)}{N}-(\alpha +\mu )I_h\hfill \\ \hfill \frac{d{I}_c}{dt}& =\frac{(1-p)\beta S(I_h+I_c)}{N}-(d_1+\delta +\mu )I_c\hfill \\ \hfill \frac{d{Q}_h}{dt}& =\alpha I_h-(\gamma +\mu )Q_h\hfill \\ \hfill \frac{d{Q}_H}{dt}& =\delta I_c-(d_2+\eta +\mu )Q_H\hfill \\ \hfill \frac{dR}{dt}& =p_1\nu S+\gamma Q_h+\eta Q_H-\mu R\hfill \end{array}.\right.$$

(1)

The description of the parameters of system (1) are presented in

Table 1. Description of the parameters and their value.

Parameters	Description	Value	Source
$\mathsf{\Lambda }$	Recruitment rate of individuals entering the susceptible subpopulation	4589	[1]
p	Probability of infected individuals entering the $I_h$ subpopulation	0.65	Estimated
$\beta$	Individual contact rate of the susceptible subpopulation with the infected subpopulations	0.036	[14]
$\nu$	Individual rate of the susceptible subpopulation being vaccinated	0.1	[16]
$p_1$	Chances of vaccine effectiveness	0.55	Estimated
$\mu$	Natural death rate of each subpopulation	$\frac{1}{65\times 365}$	Estimated
$\alpha$	Individual transfer rates from the $I_h$ to $Q_h$ subpopulations	0.2143	Assumed
$d_1$	Death rate of individuals infected with COVID-19	$7.344\times {10}^{-7}$	[16]
$\delta$	Individual transfer rates from subpopulation $I_c$ to $Q_H$	0.0714	Assumed
$\gamma$	Recovery rate of individuals who are infected with COVID-19 but have no comorbidities	0.003065	Estimated
$d_2$	Death rate of individuals infected with COVID-19 and having comorbid diseases	0.0177	[1]
$\eta$	Recovery rate of individuals infected with COVID-19 and having comorbid diseases	0.00124	[1]

. Their value is derived in part from Indonesian journals, calculated based on data from the Ministry of Health of the Republic of Indonesia, and assumed in a minor part.

3. Model Analysis

The total population is bounded at any time $t>0$. The following lemmas show that the solution of system (1) is bounded.

Lemma 1. 

The feasible region Ω defined by $\mathsf{\Omega }=\{(S\left(t\right),I_h\left(t\right),I_c\left(t\right),Q_h\left(t\right),Q_H\left(t\right),R\left(t\right))\in {\mathbb{R}}_{+}^6|0<N\left(t\right)\le max\{N\left(0\right),\mathsf{\Lambda }/\mu \}$, with initial conditions $S\left(t\right)\ge 0,I_h\left(t\right)\ge 0,I_c\left(t\right)\ge 0,Q_h\left(t\right)\ge 0,Q_H\left(t\right)\ge 0,R\left(t\right)\ge 0$, is non-negative invariant and associated with Equation (1) for $t>0$.

Proof. 

Equation $dS/dt$ in (1), for all $t>0$, can show the change in the number of susceptible per time,

$$\frac{dS}{dt}\ge \frac{p\beta S(I_h+I_c)}{N}+\frac{(1-p)\beta S(I_h+I_c)}{N}+(\nu +\mu )S=-\left(\frac{\beta (I_h+I_c)}{N}+\nu +\mu \right)S.$$

Solving the differential equation yields

$$S\left(t\right)\ge S_{0}exp\left(-\left[\mu t+{\int }_0^t\frac{\beta (I_h\left(k\right)+I_c\left(k\right))}{N}dk\right]\right),$$

and therefore

$$\underset{t\to \infty }{lim}infS\left(t\right)\ge 0.$$

Similarly, it can be shown that $I_h\left(t\right)\ge 0,I_c\left(t\right)\ge 0,Q_h\left(t\right)\ge 0,Q_H\left(t\right)\ge 0,R\left(t\right)\ge 0$. When all initial conditions are non-negative, it can be shown that the solutions of system (1) are non-negative.    □

Lemma 2. 

Solutions of system (1) are bounded for all $t\in [0,t_1]$.

Proof. 

The change in total population per unit time is

$$\frac{dN}{dt}=\mathsf{\Lambda }-d_1{I}_c-d_2{Q}_H-\mu N.$$

Therefore, the total population is bounded, that is

$$0\le \underset{t\to \infty }{lim}supN\left(t\right)\le \frac{\mathsf{\Lambda }}{\mu }.$$

This indicates that all solutions of system (1) are bounded for all $t\in [0,t_1]$.    □

3.1. Non-Endemic Equilibrium Point

The necessary and sufficient conditions for obtaining a non-endemic equilibrium point are the number of changes in subpopulations per unit time is constant, and the number of individuals confirmed positive for COVID-19 and the number of individuals quarantined are both zero. From (1), by setting $\frac{dS}{dt}=\frac{d{I}_h}{dt}=\frac{d{I}_c}{dt}=\frac{d{Q}_h}{dt}=\frac{d{Q}_H}{dt}=\frac{dR}{dt}=0$ and $I_h=I_c=Q_h=Q_H=0$, we obtain the non-endemic equilibrium point as follows

$$E_1=\left(S_1^{*},I_{h1}^{*},I_{c1}^{*},Q_{h1}^{*},Q_{H1}^{*},R_1^{*}\right)=\left(\frac{\mathsf{\Lambda }}{\mu },0,0,0,0,\frac{\nu \mathsf{\Lambda }}{\mu (\nu +\mu )}\right).$$

(2)

The local stability of the non-endemic equilibrium is given below.

Theorem 1. 

The non-endemic equilibrium point $E_1$ is locally asymptotically is stable.

Proof. 

In the manner used by Diekmann and Heesterbeek [32], by substituting the non-endemic equilibrium point in the Jacobian matrix of system (1), we obtain

$$J\left(E_1\right)=\left[\begin{array}{cccccc}-p_1\nu -\mu & -k& -k& 0& 0& 0\\ 0& pk-\alpha -\mu & pk& 0& 0& 0\\ 0& (1-p)k& (1-p)k-d_1-\delta -\mu & 0& 0& 0\\ 0& \alpha & 0& -\gamma -\mu & 0& 0\\ 0& 0& \delta & 0& -d_2-\eta -\mu & 0\\ p_1\nu & 0& 0& \gamma & \eta & -\mu \end{array}\right],$$

with $k=\frac{\beta \mathsf{\Lambda }}{\mu N}$. Eigenvalues of the Jacobian matrix $J\left(E_1\right)$ are ${\lambda }_1=-\mu$, ${\lambda }_2=--p_1\nu -\mu$, ${\lambda }_3=-d_2-\eta -\mu$, ${\lambda }_4=-\gamma -\mu$, ${\lambda }_5=-\frac{1}{2}(\alpha +d_1+\delta -2\mu )+\frac{1}{2}k-\frac{1}{2}\sqrt{A}$, ${\lambda }_6=-\frac{1}{2}(\alpha +d_1+\delta -2\mu )+\frac{1}{2}k+\frac{1}{2}\sqrt{A}$, where

$$A=4kp(d_1+\delta -\alpha )+{(\alpha -d_1)}^2+2\alpha (k-\delta )+d_1(d_1+2\delta -2k)+{(\delta -k)}^2.$$

Because it can be observed that ${\lambda }_i<0$, $i=1,2,3,4,5$, and ${\lambda }_6$ will be negative if

$$\frac{1}{2}(kp+\alpha +d_1)+\mu >d+\frac{1}{2}\sqrt{A},$$

(3)

then the equilibrium point $E_1$ is asymptotically stable if (3) is fulfilled.    □

3.2. The Basic Reproduction Number

The basic reproduction number is an essential indicator for determining whether the rate of disease transmission is rising or decreasing. The basic reproduction number is usually denoted by $R_0$. The interpretation is that if $R_0>1$, then the disease is increasing, and if $R_0<1$, then the disease is decreasing or disappearing. Based on the system (1), the basic reproduction number can be determined using the next-generation matrix [33]. The Jacobian matrix, taking into account the susceptible subpopulation in contact with the infected, denoted by $J_G$, is given as follows

$$J_G=\frac{1}{N}\left[\begin{array}{ccc}\beta (I_h+I_c)& -\beta S& -\beta S\\ p\beta (I_h+I_c)& p\beta S& p\beta S\\ (1-p)(I_h+I_c)& (1-p)\beta S& (1-p)\beta S\end{array}\right].$$

Evaluating $J_G$ at the non-endemic equilibrium results

$$G_F:=J_G\left(E_1\right)=\left[\begin{array}{ccc}0& -\frac{\beta \mathsf{\Lambda }}{\mu N}& -\frac{\beta \mathsf{\Lambda }}{\mu N}\\ 0& \frac{p\beta \mathsf{\Lambda }}{\mu N}& \frac{p\beta \mathsf{\Lambda }}{\mu N}\\ 0& \frac{(1-p)\beta \mathsf{\Lambda }}{\mu N}& \frac{(1-p)\beta \mathsf{\Lambda }}{\mu N}\end{array}\right].$$

The Jacobian matrix of susceptible and infected subpopulations that enter and are contactless is

$$G_V=\left[\begin{array}{ccc}{p}_1\nu +\mu & 0& 0\\ 0& \alpha +\mu & 0\\ 0& 0& d_1+\delta +\mu \end{array}\right].$$

Eigenvalues of $G_F{G}_V^{-1}$ are ${\lambda }_1=0$, ${\lambda }_2=0$, ${\lambda }_3=\frac{\beta \mathsf{\Lambda }[p(d_1+\delta -\alpha )+\alpha +\mu ]}{\mu N(\alpha +\mu )(d_1+\delta +\mu )}$. According to [33], the basic reproduction number of the system (1) is

$$R_0=\frac{\beta \mathsf{\Lambda }[p(d_1+\delta -\alpha )+\alpha +\mu ]}{\mu N(\alpha +\mu )(d_1+\delta +\mu )}.$$

(4)

3.3. Endemic Equilibrium Point

The necessary conditions for obtaining the endemic equilibrium point of system (1) are $\frac{dS}{dt}=\frac{d{I}_h}{dt}=\frac{d{I}_c}{dt}=\frac{d{Q}_h}{dt}=\frac{d{Q}_H}{dt}=\frac{dR}{dt}=0$. The endemic equilibrium point is

$$E_2=\left(S_2^{*},I_{h2}^{*},I_{c2}^{*},Q_{h2}^{*},Q_{H2}^{*},R_2^{*}\right),$$

with

$$\begin{array}{cc}\hfill S_2^{*}& =\frac{k_2{k}_{3}N}{\beta [p(k_3-k_2)+k_2]},\hfill \\ \hfill I_{h2}^{*}& =\frac{p\beta \mathsf{\Lambda }(k_2-k_3)+k_2(k_1{k}_{3}N-\beta \mathsf{\Lambda })}{\beta k_2[p(k_2-k_3)-k_2]},\hfill \\ \hfill I_{c2}^{*}& =\frac{(1-p)\left[p\beta \mathsf{\Lambda }(k_2-k_3)+k_2(k_1{k}_{3}N-\beta \mathsf{\Lambda })\right]}{p\beta k_3[p(k_2-k_3)-k_2]},\hfill \\ \hfill Q_{h2}^{*}& =\frac{p\alpha \beta \mathsf{\Lambda }(k_2-k_3)+k_2(k_1{k}_{3}N-\beta \mathsf{\Lambda })}{\beta k_2{k}_4[p(k_2-k_3)-k_2]},\hfill \\ \hfill Q_{H2}^{*}& =\frac{(1-p)\delta \left[p\beta \mathsf{\Lambda }(k_2-k_3)+\delta k_2(k_1{k}_{3}N-\beta \mathsf{\Lambda })\right]}{p\beta k_3{k}_5[p(k_2-k_3)-k_2]},\hfill \\ \hfill R_2^{*}& =\frac{p_{1}r{S}_2^{*}+\gamma Q_{h2}^{*}+\eta Q_{H2}^{*}}{\mu },\hfill \end{array}$$

where $k_1=p_1\nu +\mu$, $k_2=\alpha +\mu$, $k_3=d_1+\delta +\mu$, $k_4=\gamma +\mu$, $k_5=d_2+\eta +\mu$.

The local stability of the endemic equilibrium is provided below.

Theorem 2. 

The endemic equilibrium point $E_2$ is locally asymptotically stable.

Proof. 

By substituting the endemic equilibrium point $E_2$ into the Jacobian matrix of system (1), we obtain

$$J\left(E_2\right)=\left[\begin{array}{cccccc}-{\beta }^{*}-p_1\nu -\mu & -\frac{\beta S_2^{*}}{N}& -\frac{\beta S_2^{*}}{N}& 0& 0& 0\\ p{\beta }^{*}& \frac{p\beta S_2^{*}}{N}-\alpha -\mu & \frac{p\beta S_2^{*}}{N}& 0& 0& 0\\ (1-p){\beta }^{*}& \frac{(1-p)\beta S_2^{*}}{N}& \frac{(1-p)\beta S_2^{*}}{N}-d_1-\delta -\mu & 0& 0& 0\\ 0& \alpha & 0& -\gamma -\mu & 0& 0\\ 0& 0& \delta & 0& -d_2-\eta -\mu & 0\\ p_1\nu & 0& 0& \gamma & \eta & -\mu \end{array}\right],$$

where ${\beta }^{*}=\frac{\beta (I_{h2}^{*}+I_{c2}^{*})}{N}$.

The characteristic polynomial of $J\left(E_2\right)$ is

$$\mathcal{P}\left(x\right)=\frac{1}{N^2}({\lambda }_1+\mu )(\lambda +\gamma +\mu )(\lambda +d_2+\eta +\mu )\mathcal{G}\left(\lambda \right),$$

with ${\lambda }_1=-\mu$, ${\lambda }_2=-\gamma -\mu$, ${\lambda }_3=-d_2-\eta -\mu$, and

$$\mathcal{G}\left(\lambda \right)={\lambda }^3+c_1{\lambda }^2+c_2\lambda +c_3,$$

(5)

where

$$\begin{array}{cc}\hfill c_1& =\frac{1}{N}\left([(p+1)k_6-p_1\nu +\alpha +k_3]N+\beta S_2^{*}(p-1)\right),\hfill \\ \hfill c_2& =(\left[-p{k}_6^2+((p_1\nu -d_1-\delta )p+k_2+k_3)k_6-\nu (\alpha +k_2+k_3)p_1+\alpha (\delta +d_1)-{\mu }^2\right]N^2\hfill \\ \hfill & S_2^{*}\left[(2-p)p{k}_6+(1-p)(p_1\nu -\alpha )\right]\beta N+p{\beta }^2{S_2^{*}}^2(p-1))\frac{\lambda }{N^2},\hfill \\ \hfill c_3& =-\frac{1}{N^2}[k_3{N}^2(k_1-k_6)(-p{k}_6+k_2)-S_2^{*}\left[(p-1){\mu }^2+-p^2{k}_6+k_1(p-1)\right]\mu \hfill \\ \hfill & -p^2\nu k_6{p}_1+\left[(p_1\nu -d_1-\delta )k_6+\alpha p_1\nu \right]p-\alpha p_1\nu \beta N-p(p-1){\beta }^2{S_2^{*}}^2{k}_1],\hfill \end{array}$$

with $k_6=\frac{\beta (I_{h2}^{*}+I_{c2}^{*})}{N}$.

To show that the equilibrium point $E_2$ is locally asymptotically stable, we must prove that all roots of the polynomial (5) are negative. Based on the Routh–Hurwitz criteria, all roots of (5) are negative if they are fulfilled $c_1>0,c_2>0,c_3>0$ and $c_1{c}_2-c_3>0$. Using the parameter values in Table 1 and the initial subpopulation values $S_2^{*},I_h{2}^{*},I_c{2}^{*}$. We obtain $c_1=0.218156$, $c_2=0.01161589$, $c_3=0.00067804$, and we can verify that $c_1{c}_2-c_3=0.00185604$. This shows that $c_1>0,c_2>0,c_3>0$, and $c_1{c}_2-c_3>0$. Therefore, the endemic equilibrium point $E_2$ is locally asymptotically stable.    □

3.4. Global Stability Analysis

Following [34], the global stability of equilibrium point $E_2$ is analyzed. The expression of $E_2$ can be written as

$$\begin{array}{cc}\hfill & \mathsf{\Lambda }=\frac{q{S}_2^{*}(I_{h2}^{*}+I_{c2}^{*})}{N}+(p_1\nu +\mu )S_2^{*},I_{h2}^{*}=\frac{S_2^{*}{I}_{c2}^{*}}{(\alpha +\mu )N-pq{S}_2^{*}},\hfill \\ \hfill & I_{c2}^{*}=\frac{(1-p)q{S}_2^{*}{I}_{h2}^{*}}{(d_1+\delta +\mu )N-(1-p)q{S}_2^{*}},(\gamma +\mu )Q_{h2}^{*}=\alpha I_{h2}^{*},\hfill \\ \hfill & (d_2+\eta +\mu )Q_{H2}^{*}=\delta I_{c2}^{*},\mu R_2^{*}=p_1\nu S_2^{*}+\gamma Q_{h2}^{*}+\eta Q_{H2}^{*},\hfill \end{array}$$

with q is multiplication of a constant by the parameter $\beta$.

The global stability of the endemic equilibrium is given below.

Theorem 3. 

If $R_0>1$, then the endemic equilibrium $E_2$ is globally asymptotically stable.

Proof. 

Define a Lyapunov function

$$\begin{array}{cc}\hfill L& ={\int }_{S_2^{*}}^S\left(1-\frac{S_2^{*}}{y}\right)dy+{\int }_{I_{h2}^{*}}^{I_{h2}}\left(1-\frac{I_{h2}^{*}}{y}\right)dy+{\int }_{I_{c2}^{*}}^{I_{c2}}\left(1-\frac{I_{c2}^{*}}{y}\right)dy\hfill \\ \hfill & +{\int }_{Q_{h2}^{*}}^{Q_{h2}}\left(1-\frac{Q_{h2}^{*}}{y}\right)dy+{\int }_{Q_{H2}^{*}}^{Q_{H2}}\left(1-\frac{Q_{H2}^{*}}{y}\right)dy+{\int }_{R_2^{*}}^R\left(1-\frac{R_2^{*}}{y}\right)dy.\hfill \end{array}$$

The first derivative of the Lyapunov function is given by

$$\begin{array}{cc}\hfill L^{\prime }& =\left(1-\frac{S_2^{*}}{S}\right)S_2^{\prime }+\left(1-\frac{I_{h2}^{*}}{I_{h2}}\right)I_{h2}^{\prime }+\left(1-\frac{I_{c2}^{*}}{I_{c2}}\right)I_{c2}^{\prime }+\left(1-\frac{Q_{h2}^{*}}{Q_{h2}}\right)Q_{h2}^{\prime }\hfill \\ \hfill & +\left(1-\frac{Q_{H2}^{*}}{Q_{H2}}\right)Q_{H2}^{\prime }+\left(1-\frac{R_2^{*}}{R}\right)R_2^{\prime },\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill \left(1-\frac{S_2^{*}}{S}\right)S_2^{\prime }& =\left(1-\frac{S_2^{*}}{S}\right)\left(\frac{q{S}_2^{*}(I_{h2}^{*}+I_{c2}^{*})}{N}+k_1{S}_2^{*}-\frac{\beta S(I_h+I_c)}{N}-k_{1}S\right)\hfill \\ \hfill & =(p_1\nu +\mu )S_2^{*}\left(2-\frac{S}{S_2^{*}}-\frac{S_2^{*}}{S}\right)+\left(1-\frac{S_2^{*}}{S}\right)\left[\frac{q{S}_2^{*}(I_{h2}^{*}+I_{c2}^{*})}{N}-\frac{\beta S(I_h+I_c)}{N}\right],\hfill \\ \hfill \left(1-\frac{I_{h2}^{*}}{I_h}\right)I_{h2}^{\prime }& =\left(1-\frac{I_{h2}^{*}}{I_h}\right)\left(\frac{p\beta S(I_h+I_c)}{N}-k_2{I}_h\right)\hfill \\ \hfill & =\left(1-\frac{I_{h2}^{*}}{I_h}\right)\left(\frac{p\beta S(S_2^{*}{I}_{c2}^{*}+[k_{2}N-pq{S}_2^{*}])I_c}{k_2{N}^2-pqN{S}_2^{*}}-\frac{k_2{S}_2^{*}{I}_{c2}^{*}}{k_{2}N-pq{S}_2^{*}}\right),\hfill \\ \hfill \left(1-\frac{I_{c2}^{*}}{I_c}\right)I_{c2}^{\prime }& =\left(1-\frac{I_{c2}^{*}}{I_c}\right)\left(\frac{(1-p)\beta S(I_h+I_c)}{N}-k_3{I}_c\right)\hfill \\ \hfill & =\left(1-\frac{I_{c2}^{*}}{I_c}\right)\left(\frac{(1-p)\beta S(I_h+I_{h2}^{*})}{N}-\frac{(1-p)q{k}_3{S}_2^{*}{I}_{h2}^{*}}{k_{3}N-(1-p)q{S}_2^{*}}\right),\hfill \\ \hfill \left(1-\frac{Q_{h2}^{*}}{Q_h}\right)Q_{h2}^{\prime }& =\left(1-\frac{Q_{h2}^{*}}{Q_h}\right)\left[\alpha I_h-(\gamma +\mu )Q_h\right]\hfill \\ \hfill & =\left(1-\frac{Q_{h2}^{*}}{Q_h}\right)\alpha \left(\frac{S_2^{*}{I}_{c2}^{*}}{k_{2}N-pq{S}_2^{*}}-I_{h2}^{*}\right),\hfill \\ \hfill \left(1-\frac{Q_{H2}^{*}}{Q_H}\right)Q_{H2}^{\prime }& =\left(1-\frac{Q_{H2}^{*}}{Q_H}\right)(\delta I_c-k_3{Q}_H)\hfill \\ \hfill & =\left(1-\frac{Q_{H2}^{*}}{Q_H}\right)\delta \left(\frac{(1-p)q{S}_2^{*}{I}_{h2}^{*}}{k_{3}N-(1-p)q{S}_2^{*}}-I_{c2}^{*}\right),\hfill \\ \hfill \left(1-\frac{R_2^{*}}{R}\right)R_2^{\prime }& =\left(1-\frac{R_2^{*}}{R}\right)\left[p_1\nu (S-S_2^{*})+\gamma (Q_h-Q_{h2}^{*})+\eta (Q_H-Q_{H2}^{*})\right].\hfill \end{array}$$

Then, we have

$$\begin{array}{cc}\hfill L^{\prime }& =\mu S_2^{*}\left(2-\frac{S}{S_2^{*}}-\frac{S_2^{*}}{S}\right)+\frac{\beta }{N}\left(3-\frac{S_2^{*}}{S}-\frac{I_h}{I_{h2}^{*}}-\frac{I_{h2}^{*}{Q}_h}{Q_{h2}^{*}{I}_h}-\frac{I_c}{I_{c2}^{*}}\left[\frac{S{Q}_{h2}^{*}}{S_2^{*}{Q}_h}-1\right]\right)\hfill \\ \hfill & +\frac{\beta }{N}\left(3-\frac{S_2^{*}}{S}-\frac{I_c}{I_{c2}^{*}}-\frac{I_{c2}^{*}{Q}_H}{Q_{H2}^{*}{I}_c}-\frac{I_h}{I_{h2}^{*}}\left[\frac{S{Q}_{H2}^{*}}{S_2^{*}{Q}_H}-1\right]\right).\hfill \end{array}$$

(6)

From (6), we have

$$\begin{array}{cc}\hfill & \mu S_2^{*}\left(2-\frac{S}{S_2^{*}}-\frac{S_2^{*}}{S}\right)\le 0,\hfill \\ \hfill & \frac{\beta }{N}\left(3-\frac{S_2^{*}}{S}-\frac{I_h}{I_{h2}^{*}}-\frac{I_{h2}^{*}{Q}_h}{Q_{h2}^{*}{I}_h}-\frac{I_c}{I_{c2}^{*}}\left[\frac{S{Q}_{h2}^{*}}{S_2^{*}{Q}_h}-1\right]\right)\le 0,\hfill \\ \hfill & \frac{\beta }{N}\left(3-\frac{S_2^{*}}{S}-\frac{I_c}{I_{c2}^{*}}-\frac{I_{c2}^{*}{Q}_H}{Q_{H2}^{*}{I}_c}-\frac{I_h}{I_{h2}^{*}}\left[\frac{S{Q}_{H2}^{*}}{S_2^{*}{Q}_H}-1\right]\right)\le 0.\hfill \end{array}$$

Thus, the largest invariant subset is obtained $L^{\prime }=0$. By Lasalle’s invariance principle [35], if $R_0>1$, then $E_2$ is globally asymptotically stable.    □

4. Result and Discussion

4.1. Numerical Solutions

The numerical simulations in this study use the parameter values in Table 1, and the number of initial subpopulations in 2020 of the spread of COVID-19 in Indonesia with the initial number of each subpopulation is $S\left(0\right)=267823870$, $I_h\left(0\right)=23601$, $I_c\left(0\right)=3842$, $Q_h\left(0\right)=8630$, $Q_H\left(0\right)=2034$, $R\left(0\right)=286$, and $N=268000000$. We solve the system (1) numerically by using the Runge–Kutta order 45 method; in Matlab, we use package function ode45. The graphs of the numerical solution of the system are given in Figure 2. We group the subpopulations S with R, $I_h$ with $I_c$, and $Q_h$ and $Q_H$, based on the value level and comparison meaning. Here, we can observe that a symmetrical aspect appears in the susceptible and recovered subpopulations.

4.2. Basic Reproduction Number Sensitivity Analysis

The basic reproduction number $R_0$ in (4) can be viewed as a function of parameters. So, we can visualize $R_0$ as a function of two parameters by making other parameters the constant and then plot the graph. We have plotted the contour plots of the $R_0$ function for every combination of two parameters that appeared in (4), and most of them give results with the contour of $R_0$ having a value below 1. Only from the combination of parameters $\mathsf{\Lambda }$, $\beta$, p, and $d_1$, there can be seen the spectrum of the $R_0$ value varying from below 1 up to above 1, as shown in Figure 3. We can conclude that the combination of these parameters has a big role in determining the stability of the system. In other words, these parameters have the most effect on the spread of COVID-19.

Geometrically, Figure 3 shows the contour plot of the basic reproduction number $R_0$ as a function of two parameters. The color bar placed on the right shows the value of $R_0$; in other words, it shows the height of the $R_0$’s surface. Because 1 is the critical point for $R_0$ that concludes the system’s stability, then the interpretation of Figure 3 can be performed by observing the effect of the changes in the parameter value on the $R_0$ value, whether it makes the $R_0$ approach or leave 1. Figure 3a–c successively show that the combination of parameters $\mathsf{\Lambda }$ and $\beta$, p and $\beta$, and $d_1$ and $\beta$ have a dominant influence on increasing or decreasing the spread of COVID-19, whereas Figure 3d–f show that the combination of parameters p and $\mathsf{\Lambda }$, $d_1$ and $\mathsf{\Lambda }$, and $d_1$ and p does not have a dominant effect on increasing or decreasing the spread of COVID-19.

4.3. Sensitivity Analysis of Parameters

A sensitivity analysis is performed to examine the effect of changes in the model’s parameters to the dynamics of the model’s compartments. Some papers have used this sensitivity technique in several research areas [36,37,38]. Suppose we have a vector of variables $\mathbf{X}={(S,I_h,I_c,Q_h,Q_H,R)}^T$, vector parameters $\mathbf{P}={(\mathsf{\Lambda },p,\beta ,\nu ,p_1,\mu ,\alpha ,d_1,\delta ,\gamma ,d_2,\eta )}^T$, and a vector of the right side of system (1) denoted by $\mathbf{F}={(\mathcal{F},\mathcal{G},\mathcal{H},\mathcal{I},\mathcal{J},\mathcal{K})}^T$, where $dS/dt=\mathcal{F}$, $d{I}_h/dt=\mathcal{G}$, $d{I}_c/dt=\mathcal{H}$, $d{Q}_h/dt=\mathcal{I}$, $d{Q}_H/dt=\mathcal{J}$, $dR/dt=\mathcal{K}$.

Define a sensitivity function $V=\frac{\partial \mathbf{X}}{\partial \mathbf{P}}$. Because $\mathbf{X}$ contains $\mathbf{P}$ and t, and $\mathbf{F}$ contains $\mathbf{X}$ and $\mathbf{P}$, then by performing total derivative for function V, we can derive a system of differential equations as follows

$$\frac{dV}{dt}=J\left(\mathbf{X}\right)V+\frac{\partial \mathbf{F}}{\partial \mathbf{P}},$$

(7)

where $J\left(\mathbf{X}\right)$ is the Jacobian matrix $6\times 6$. There are six variables and twelve parameters. The size of matrix V is $6\times 12$ and of matrix $\frac{\partial \mathbf{F}}{\partial \mathbf{P}}$ is $6\times 12$. Therefore, Equation (7) is a $6\times 12$ matrix equation. If we pull out the elements of this matrix equation, we have a system with seventy-two (as many as the number of elements of V) plus six (number of elements of $\mathbf{X}$) differential equations.

For further writings, we denote $\frac{\partial S}{\partial \mathsf{\Lambda }}=v_{\mathsf{\Lambda }}^S$, $\frac{\partial S}{\partial p}=v_p^S$, $\frac{\partial S}{\partial \beta }=v_{\beta }^S$, …, $\frac{\partial S}{\partial \eta }=v_{\eta }^S$, and so on for the other variables. By solving the system (7) numerically, we can obtain solutions for the sensitivity of parameters over time. This time-dependent sensitivity is presented in Figure 4. At any given time, we can observe the effect of each parameter’s changes to each variable. The positive value of the sensitivity means that if the parameter is raised, then the variable goes higher; otherwise, the negative value of the sensitivity means that if the parameter is raised, then the variable goes lower, and the zero value of the sensitivity means that the changes in the parameter’s value do not affect the variable at all.

Based on the sensitivity analysis in Figure 4a, if the value of the parameter $\mathsf{\Lambda }$ increases, then the number of subpopulations R and S increases. In Figure 4b,h,k, the increase in the values of the parameters p, $d_1$, or $d_2$ results in the increase in the number of subpopulation $Q_h$ and the decrease in the number of subpopulations $I_c$ or $Q_H$. In Figure 4c, if the value of parameter $\beta$ increases, then the number of subpopulations $Q_h$, $Q_H$, $I_c$, and $I_h$ increases, but the number of subpopulations R and S decreases. In Figure 4d,e, if the value of the parameters $\nu$ or $p_1$ increases, then the number of subpopulation R increases, but the number of subpopulation S decreases. In Figure 4f, if the value of parameter $\mu$ increases, then the number of subpopulations R and S decreases. In Figure 4g, the increase in the value of parameter $\alpha$ results in the increase in the number of subpopulations R, $I_c$, and S but the decrease in the number of subpopulation $I_h$. In Figure 4i, if the value of parameter $\delta$ increases, then the number of subpopulations R, $Q_H$, and S increases, but the number of subpopulations $I_c$ and $Q_H$ decreases. In Figure 4j,l, if the value of the parameters $\gamma$ or $\eta$ increases, then the number of subpopulation R increases, but the number of subpopulation $Q_h$ decreases.

From Figure 4, we can observe that some sensitivity has a positive value at some point in the first period, but at the next period, the sensitivity value becomes negative or tends to zero. From these sensitivity graphs, especially at the peak point, we can observe the time when the parameter changes have the most effect on the variable. Therefore, another interesting result can be derived by observing the sensitivity index after arriving at the equilibrium. This can be performed by taking the sensitivity values when the time is very long. By this method, we can determine which one of the parameters is the most sensitive to the model when converging to the equilibrium. In Figure 5, we present the sensitivity index of all the parameters for each variable. We give four different scales for depicting the sensitivity index, because some variables have values quite larger than other variables; thus, it makes the sensitivity index have a different level for each variable. From each panel of Figure 5, we can observe that parameter $\mu$, which stands for the natural death rate for each subpopulation, is the most sensitive parameter.

5. Discussion and Conclusions

This article explores the dynamic model of COVID-19 where infected individuals with comorbid disease are quarantined in hospitals provided by the local government and infected individuals without comorbid disease are self-isolated. Based on the results of the model analysis, the following are obtained: non-endemic and endemic equilibrium points; the basic reproduction number to determine whether the spread of COVID-19 is increasing or decreasing; the local asymptotic stability theorem of the non-endemic equilibrium point; and the local and global asymptotic stability theorems of the endemic equilibrium point.

According to the results of the sensitivity analysis of the basic reproduction number using a contour plot, the parameters for which their combination significantly impacts the system stability are the recruitment rate of individuals entering the susceptible subpopulation, probability of infected individuals that have no comorbid disease, contact rate between the susceptible subpopulation and the infected subpopulation, and death rate of individuals infected with COVID-19. In addition, the combination of the contact rate with the other three parameters produces a wider range of COVID-19 spread. Therefore, minimizing the contact rate is very good advice to be implemented in the government policy in order to reduce the disease contagion.

In addition, we provided a sensitivity analysis to examine the effect of the parameter changes in each subpopulation. From the visualization of the sensitivity function over time, we can observe the dynamical effect of each of the parameter’s changes. We found, based on the sensitivity index after achieving equilibrium, that the natural death rate is the most sensitive parameter. The implication of this result is that the estimation of the natural death rate parameter must be performed carefully in order to avoid a transgression in depicting the true dynamics of the COVID-19 spread.

We also found that the symmetry aspects appear in the visualization of the susceptible and recovered subpopulations over time, the image of some of the basic reproduction number sensitivity, and the depiction of some parameters’ sensitivity over time.