Mathematical Modeling and Stability Analysis of the COVID-19 Spread by Considering Quarantine and Hospitalize

The spread of COVID-19 model in this study, represents the interaction between seven subpopulations, that is namely susceptible (S), exposed (E), infected (I), quarantined-1 (Q₁), quarantined-2 (Q₂), hospitalized (H) and recovered (R). Hence, the total population at time t is represented by

$N(t)=S(t)+E(t)+I(t)+Q_1(t)+Q_2(t)+H(t)+R(t)$

In this model the migration rate is ignored. The transmission of COVID-19 is assumed from susceptible individuals (S) who are in contact with exposed individuals (E) and infected individuals (I) with the incidence rate $\beta_1$ and $\beta_2$​ respectively. It is assumed that the exposed individuals (E) are transmitted to infected individuals (I) and the quarantined exposed individuals (Q₁) is transmitted to quarantined infected individuals (Q₂) with the same average rate of v, also the infected individuals (I) and the quarantined infected (Q₂) are transmitted to hospitalized individuals (H) with the same transition rate of η. Hospitalized individuals (H) are transmitted to recovered individuals with rate of $\gamma$. The proposed COVID-19 disease transmission model is given in the scheme showed in Figure 1.

1.png

Figure 1. Schematic model of the spread of COVID-19

The dynamics of each compartment, which were used to formulate the mathematical model representing the COVID-19 spread, are described as follows. Changes in the susceptible individuals (S) cover: (1) increment due to the recruitment rate of A, (2) decrement due to infections from direct interactions between vulnerable individuals and COVID-19 carriers, and move to individuals E with rates $\beta_1$​ and $\beta_2$​, (3) decrement due to natural death by μ₁, Further, it increased due to the rate of development from Q₁ to S of $\varepsilon$.

The number of individuals in the exposed subpopulation (E) increased due to the interaction between susceptible individuals and COVID-19 carriers, causing E and I with rate $\beta_1$​ and $\beta_2$​ respectively. Further, it decreased due to the rate of development from E to Q₁ and I of q₁ and v respectively. Finally, there is a natural death of μ₁.

Changes in the number of the exposed individual who were quarantined (Q₁) cover: (1) increment due to the rate of development of exposed E to Q₁ by the rate of q₁, (2) decrement due to recovery with the rate of $\varepsilon$, (3) decrement due to the rate of transmission from Q₁ to Q₂  and natural death of v and μ₁respectively.

The number of the infected individuals (I) increased due to the rate of development from the E class to positive individuals I by the rate of v. The infected individuals by COVID-19 decreased due to the rate of change of infected individuals to quarantine Q₂ and hospitalized H by q₂ and η respectively, decrement because of the death rate due to COVID-19 infected and natural deaths with the rate of μ₂ and μ₁ respectively.

Changes in the number of the infected individual who were quarantined (Q₂) cover: (1) increment due to the rate of development of infected individuals I to Q₂ and the rate of transmission from Q₁ to Q₂ by the rate of q₂ and v respectively, (2) decrement because of the death rate due to COVID-19 infected and natural deaths with the rate of μ₂ nd μ₁ respectively. (3) decrement due to the rate of transmission from Q₂ to H of $\eta$.

The number of the hospitalized individuals (H) increased by subpopulation transmissions from Q₂ and I with the rate of $\eta$ Further, H reduced due to COVID-19 infected and natural deaths with the rate of μ₂ and μ₁ respectively. Finally, decrement due to the rate of transmission from H to R of $\gamma$. 

The number of the recovered individual (R) increased by population transmissions from H with the rate of $\gamma$ and reduced due to natural deaths by μ₁.

There are seven variables involved in the model, the dynamic of each variable is modeled as a differential equation. Below, we present the derivation of the differential equation of S only, where the model derivation for the other six variables was derived analogously. The dynamic of S at a time t with the time increment $\Delta t$ is modeled as follows:

$\begin{matrix}  S(t+\Delta t)-S(t)=(A+\varepsilon {{Q}_{1}}(t)-{{\beta }_{1}}S(t)I(t) \\   -{{\beta }_{2}}S(t)E(t)-{{\mu }_{1}}S(t))\Delta t. \\ \end{matrix}$

By taking the limit value of S for $\Delta t \rightarrow 0$, we derive the following differential equation:

$\lim _{\Delta t \rightarrow 0} \frac{S(t+\Delta t)-S(t)}{\Delta t}=\frac{d S}{d t}=A+\varepsilon Q_1(t)-\beta_1 S(t) I(t)-\beta_2 S(t) E(t)-\mu_1 S(t)$.

Based on the dynamics described above and Figure 1, the non-linear differential equation system for all seven variables as the complete mathematical model of the transmission of COVID-19 disease is obtained as follows

$\left\{\begin{array}{l}\frac{d S}{d t}=A+\varepsilon Q_1-\beta_1 S I-\beta_2 S E-\mu_1 S, \\ \frac{d E}{d t}=\beta_1 S I+\beta_2 S E-q_1 E-\mu_1 E-v E, \\ \frac{d I}{d t}=v E-q_2 I-\eta I-\left(\mu_1+\mu_2\right) I, \\ \frac{d Q_1}{d t}=q_1 E-\varepsilon Q_1-v Q_1-\mu_1 Q_1, \\ \frac{d Q_2}{d t}=q_2 I+v Q_1-\eta Q_2-\left(\mu_1+\mu_2\right) Q_2, \\ \frac{d H}{d t}=\eta I+\eta Q_2-\gamma H-\left(\mu_1+\mu_2\right) H, \\ \frac{d R}{d t}=\gamma H-\mu_1 R.\end{array}\right.$               (1)

with non-negative initial values

$S(0)=S_0, \quad E(0)=E_0, I(0)=I_0, Q_1(0)=Q_{1_0}, Q_2(0)=Q_{2_0}, H(0)=H_0, R(0)=R_0$

and all parameters are positive.

2.1 Positiveness and boundedness of solutions

In order to analyze whether the COVID-19 model as in (1) has epidemic significance, the first step is to prove that the variables are positive. This means that the solution of the system of equations with non-negative initial conditions must be positive for every t>0 and also model (1) represent the interaction between of human subpopulations, so the solutions are bounded. The next theorem has guaranteed the positiveness and boundedness of solutions of the system.

Theorem 1. Let S(0)≥0, E(0)≥0, I(0)≥0, Q₁(0)≥0, Q₂(0)≥0,, H(0)≥0, and R(0)≥0, then the solutions S(t), E(t), I(t), Q₁(t), Q₂(t), H(t),and R(t) of the model (1) are positive for every t>0.

Proof. Please see Appendix A for proof

Theorem 2. The feasible region of the model (1) is defined by $\Omega=\left\{\left(S(t), E(t), I(t), Q_1(t), Q_2(t), H(t), R(t)\right) \in\right.\left.\mathbb{R}_{+}^7: 0 \leq N(t) \leq \frac{A}{\mu_1}\right\}$ is positively invariant for the system (1).

Proof. Please see Appendix B for proof

2.2 The basic reproductive number

The basic reproductive number ($\mathfrak{N}_0$) is the average number of new infection cases in a population. If $\mathfrak{N}_0<1$, the spread of the disease can be controlled and will not become an epidemic [27, 28]. However, if $\mathfrak{N}_0<1$, then each infected individual will spread the disease to other individuals, so that it can lead to an epidemic [29, 30]. The reproduction number is found from the Next Generation Matrix (NGM) method, which is to build a matrix that generates the number of infected individuals [31] in this dynamic system model the infected compartments are the exposed (E) and infectious (I).

For example, $X=\left[\begin{array}{ll}E & I\end{array}\right]^T$ so that it can be written as

$\frac{d X}{d t}=F(X)-V(X)=\left[\begin{array}{c}\beta_1 S I+\beta_2 S E-q_1 E-\mu_1 E-v E \\ v E-q_2 I-\eta I-\left(\mu_1+\mu_2\right) I\end{array}\right]$    

where, F(x) is a matrix containing the initial infection rate, while V(x) is a matrix containing the rate of movement of the infected population, with the values of F(x) and V(x) being

$F(X)=\left[\begin{array}{l}F_1 \\ F_2\end{array}\right]=\left[\begin{array}{c}\beta_1 S I+\beta_2 S E \\ 0\end{array}\right]$,

and

$V(X)=\left[\begin{array}{l}V_1 \\ V_2\end{array}\right]=\left[\begin{array}{c}q_1 E+\mu_1 E+v E \\ q_2 I+\eta I+\left(\mu_1+\mu_2\right) I-v E\end{array}\right]$

Suppose F and V are Jacobian matrices of F(x) and V(x) and substitute the value of the variable with the equilibrium point $K^0\left(S^0, E^0, I^0, Q_1^0, Q_2^0\right)=\left(\frac{A}{\mu_1}, 0,0,0,0\right)$ then we get:

$F=\left[\begin{array}{ll}\frac{\partial F_1}{\partial E} & \frac{\partial F_1}{\partial I} \\ \frac{\partial F_2}{\partial E} & \frac{\partial F_2}{\partial I}\end{array}\right]=\left[\begin{array}{cc}\beta_2 S & \beta_1 S \\ 0 & 0\end{array}\right]=\left[\begin{array}{cc}\frac{\beta_2 A}{\mu_1} & \frac{\beta_1 A}{\mu_1} \\ 0 & 0\end{array}\right]$,

and

$V=\left[\begin{array}{cc}\frac{\partial V_1}{\partial E} & \frac{\partial V_1}{\partial I} \\ \frac{\partial V_2}{\partial E} & \frac{\partial V_2}{\partial I}\end{array}\right]=\left[\begin{array}{cc}q_1+\mu_1+v & 0 \\ -v & q_2+\eta+\mu_1+\mu_2\end{array}\right]$

Furthermore, to determine the inverse of the V method is used Next Generation Matrix (NGM)

$V^{-1}=\frac{1}{\operatorname{det}(V)} \operatorname{Adj}(V)$

$V^{-1}=\frac{1}{q_1+\eta+\mu_1+\mu_2}\left[\begin{array}{cc}q_1+\eta+\mu_1+\mu_2 & 0 \\ v & q_1+\mu_1+v\end{array}\right]$.

So that the inverse matrix V, as follows

$N G M=F(V)^{-1} N G M\left[\frac{\beta_2 A}{\mu_1\left(q_1+\mu_1+v\right)}+\frac{\beta_1 A v}{\mu_1\left(q_1+\mu_1+v\right)\left(q_2+\eta+\mu_1+\mu_2\right)} \frac{\beta_1 A}{\mu_1\left(q_2+\eta+\mu_1+\mu_2\right)}\right]$

The reproduction number is ($\mathfrak{R}_0$) obtained from the radius spectral of the NGM determined using a non-endemic equilibrium point. The reproduction number obtained are as follow

$\Re_0=\frac{\beta_2 A}{\mu_1\left(q_1+\mu_1+v\right)}+\frac{\beta_1 A v}{\mu_1\left(q_1+\mu_1+v\right)\left(q_2+\eta+\mu_1+\mu_2\right)}$.

The form $\mathfrak{R}_0$ above, it shows that the parameters q₁ and q₂ relating to quarantine and also η (transition rate from quarantined and infected individuals to hospitalized individuals) are inversely proportional to the value of the basic reproduction number which will determine the type of stability. The greater the value of $q_1, q_2$ and $\eta$, will impact the smaller the value of $\mathfrak{R}_0$. This parameter makes an important contribution to the occurrence of epidemics when the system has an endemic equilibrium and determines the speed in achieving COVID-19 disease-free equilibrium when $\mathfrak{R}_0$ is less than one.

Analysis is used to determine the behaviour around the equilibrium point. Non-endemic equilibrium is the case in which there are no individuals infected by COVID-19 in a population, denoted by

$K^0\left(S^0, E^0, I^0, Q_1^0, Q_2^0\right)=\left(\frac{A}{\mu_1}, 0,0,0,0\right)$.

The local stability analysis of the non-endemic equilibrium point for the system of Eq. (1) is expressed in the following theorem.

Theorem 3. If $\Re_0<1 $then the non-endemic equilibrium point K⁰ is locally asymptotically stable and $\Re_0<1$, then K⁰ is unstable.

Proof. Please see Appendix C for proof.

Further, we explore the global stability analysis of the endemic equilibrium point K* by using Lyapunov method. The conditions for global stability of K* is discuss as follows.

Theorem 4. Let the endemic equilibrium point K* of system (1) exists.

1. If $\mathfrak{R}_0>1$ then the endemic equilibrium point K* is globally asymptotically stable.

2. If $\mathfrak{R}_0>1$, then it K^* is unstable.

Proof. To prove the stability of the endemic equilibrium point, we will use the Lyapunov method. Assume that $\mathfrak{R}_0>1$ and let suitable Lyapunov function [32] in the form.

$\mathcal{L}(x)=\sum_{i=1}^n a_1\left(x_i-x_i^*-x_i^* \ln \frac{x_i}{x_i^*}\right)$,           (4)

with x=(S,E,I,Q₁,Q₂ ).

Adjusting to (4), a Lyapunov function is formed by $\mathcal{L}: \Omega \subset \mathbb{R}^5 \rightarrow \mathbb{R}$ with:

$\mathcal{L}\left(S, E, I, Q_1, Q_2\right)=\left(S-S^*-S^* \ln \frac{S}{S^*}\right)$

$+a_1\left(E-E^*-E^* \ln \frac{E}{E^*}\right)+a_2\left(I-I^*-I^* \ln \frac{I}{I^*}\right)$

$+a_3\left(Q_1-Q_1^*-Q_1^* \ln \frac{Q_1}{Q_1^*}\right)+a_4\left(Q_2-Q_2^*-Q_2^* \ln \frac{Q_2}{Q_2^*}\right)$,

where, $\forall\left(S, E, I, Q_1, Q_2, R\right) \in \Omega$ and $a_1, a_2, a_3$, and $a_4$ are real numbers. The function $\mathcal{L}$ is a Lyapunov function because it satisfies the definition of a Lyapunov function which will be shown as follows. The function $\mathcal{L}$ is continuous on Ω because the function $\mathcal{L}$ contains logarithms and has a continuous first partial derivative on Ω.

For any $K=\left(S, E, I, Q_1, Q_2\right) \in \Omega$ with $K \neq K^*$, then $\mathcal{L}(t)>0$ and if $K=K^*$, then $\mathcal{L}(t)=0$.

Next, we derived $\mathcal{L}(t)>0$ when $K \neq K^*$. Let $\frac{K}{K^*}=$ $\operatorname{aand} g(a)=K-K^*-K^* \ln \frac{K}{K^*}$, then

$g(a)=K^*\left(\frac{K}{K^*}-1-\ln \frac{K}{K^*}\right)=K^*(a-1-\ln a)$.

Note that the point $a=1$ is the minimum point of $g(a)$ with $g(1)=0$, because $g^{\prime}(1)=0$ and $g^{\prime \prime}(a)=\frac{1}{a^2}>0$. Thus it is obtained $g(a)=K-K^*-K^* \ln \frac{K}{K^*}>0$, for $K \neq K^*$.

Further, to show that the equilibrium point K^* is a global minimum point is done by obtaining the Hessian matrix at K^* as follows.

$H\left(K^*\right)=\left[\begin{array}{ccccc}\frac{\partial^2 \mathrm{~L}}{\partial S^2} & \frac{\partial^2 \mathrm{~L}}{\partial S \partial E} & \frac{\partial^2 \mathrm{~L}}{\partial S \partial I} & \frac{\partial^2 \mathrm{~L}}{\partial S \partial Q_1} & \frac{\partial^2 \mathrm{~L}}{\partial S \partial Q_2} \\ \frac{\partial^2 \mathrm{~L}}{\partial S \partial E} & \frac{\partial^2 \mathrm{~L}}{\partial E^2} & \frac{\partial^2 \mathrm{~L}}{\partial E \partial I} & \frac{\partial^2 \mathrm{~L}}{\partial E \partial Q_1} & \frac{\partial^2 \mathrm{~L}}{\partial E \partial Q_2} \\ \frac{\partial^2 \mathrm{~L}}{\partial S \partial I} & \frac{\partial^2 \mathrm{~L}}{\partial E \partial I} & \frac{\partial^2 \mathrm{~L}}{\partial I^2} & \frac{\partial^2 \mathrm{~L}}{\partial I \partial Q_1} & \frac{\partial^2 \mathrm{~L}}{\partial I \partial Q_2} \\ \frac{\partial^2 \mathrm{~L}}{\partial S \partial Q_1} & \frac{\partial^2 \mathrm{~L}}{\partial E \partial Q_1} & \frac{\partial^2 \mathrm{~L}}{\partial I \partial Q_1} & \frac{\partial^2 \mathrm{~L}}{\partial Q_1^2} & \frac{\partial^2 \mathrm{~L}}{\partial Q_1 \partial Q_2} \\ \frac{\partial^2 \mathrm{~L}}{\partial S \partial Q_2} & \frac{\partial^2 \mathrm{~L}}{\partial E \partial Q_2} & \frac{\partial^2 \mathrm{~L}}{\partial I \partial Q_2} & \frac{\partial^2 \mathrm{~L}}{\partial Q_1 \partial Q_2} & \frac{\partial^2 \mathrm{~L}}{\partial Q_2^2}\end{array}\right]=\left[\begin{array}{ccccc}\frac{1}{S^*} & 0 & 0 & 0 & 0 \\ 0 & \frac{a_1}{E^*} & 0 & 0 & 0 \\ 0 & 0 & \frac{a_2}{I^*} & 0 & 0 \\ 0 & 0 & 0 & \frac{a_3}{Q_1} & 0 \\ 0 & 0 & 0 & 0 & \frac{a_4}{Q_2}\end{array}\right]$.

Matrix $H\left(K^*\right)$ is positive due to $\operatorname{det}\left(H\left(K^*\right)\right)=$ $\frac{a_1 a_2 a_3}{s^* E^* I^* Q_1^* Q_2^*}>0$. Furthermore, when $K=K^*$ is obtained Lyapunov function $\mathcal{L}(t)$ as follows:

$\mathcal{L}\left(S, E, I, Q_1, Q_2\right)=\left(S-S^*-S^* \ln \frac{S}{S^*}\right)+a_1\left(E-E^*-E^* \ln \frac{E}{E^*}\right)+a_2\left(I-I^*-I^* \ln \frac{I}{I^*}\right)$

$+a_3\left(Q_1-Q_1^*-Q_1^* \ln \frac{Q_1}{Q_1^*}\right)+a_4\left(Q_2-Q_2^*-Q_2^* \ln \frac{Q_2}{Q_2^*}\right), =\left(S^* \ln (1)\right)+a_1\left(E^* \ln (1)\right)+a_2\left(I^* \ln (1)\right)+a_3\left(Q_1^* \ln (1)\right)+a_4\left(Q_2^* \ln (1)\right)=0$.

It is proved that $\mathcal{L}(t)>0$, if K≠K* and $\mathcal{L}(t)=0$ if K=K*, and K* are a global minimal.

The derivative of the function $\mathcal{L}$ with respect to t is:

$\begin{aligned} \frac{d \mathcal{L}}{d t} & =\frac{\partial \mathcal{L}}{\partial S} \frac{d S}{d t}+\frac{\partial \mathcal{L}}{\partial E} \frac{d E}{d t}+\frac{\partial \mathcal{L}}{\partial I} \frac{d I}{d t}+\frac{\partial \mathcal{L}}{\partial Q_1} \frac{d Q_1}{d t}+\frac{\partial \mathcal{L}}{\partial Q_2} \frac{d Q_2}{d t} \\ \frac{d \mathcal{L}}{d t} & =\left(1-\frac{S^*}{S}\right) \frac{d S}{d t}+a_1\left(1-\frac{E^*}{E}\right) \frac{d E}{d t}+a_2\left(1-\frac{I^*}{I}\right) \frac{d I}{d t} \\ & +a_3\left(1-\frac{Q_1^*}{Q_1}\right) \frac{d Q_1}{d t}+a_4\left(1-\frac{Q_2^*}{Q_2}\right) \frac{d Q_2}{d t}\end{aligned}$

$\frac{d \mathcal{L}}{d t}=\left(1-\frac{S^*}{S}\right)\left(A+\varepsilon Q_1-\beta_1 S I-\beta_2 S E-\mu_1 S\right)$

$+a_1\left(1-\frac{E^*}{E}\right)\left(\beta_1 S I+\beta_2 S E-B_1 E\right)$

$+a_2\left(1-\frac{I^*}{I}\right)\left(v E-B_2 I\right)+a_3\left(1-\frac{Q_1^*}{Q_1}\right)\left(q_1 E-B_3 Q_1\right)$

$+a_4\left(1-\frac{Q_2^*}{Q_2}\right)\left(q_2 I+v Q_1-B_4 Q_2\right)$,         (5)

where

$B_1=\left(q_1+\mu_1+v\right)$

$B_2=\left(q_2+\eta+\mu_1+\mu_2\right)$

$B_3=\left(\varepsilon+v+\mu_1\right)$

$B_4=\left(\eta+\mu_1+\mu_2\right)$.                  (6)

The relationship between $S^*, E *, I^*, Q_1^*$, and $Q_2^*$ with Eqns. (6) and (3) is as follows

$\left\{\begin{aligned} A+\varepsilon Q_1^* & =\beta_1 S^* I^*+\beta_2 S^* E^*+\mu_1 S^* \\ \beta_1 S^* I^*+\beta_2 S^* E^* & =B_1 E^* \\ v E^* & =B_2 I^* \\ q 1 E^* & =B_3 Q_1^* \\ q_2 I^*+v Q_1^* & =B_4 Q_2^*.\end{aligned}\right.$                  (7)

Thus, Eq. (5) becomes

$\frac{d \mathcal{L}}{d t}=\left[A+\varepsilon Q_1-\beta_1 S I-\beta_2 S E-\mu_1 S\right]$

$+\left[-A \frac{S^*}{S}-\varepsilon Q_1 \frac{S^*}{S}+\beta_1 I S^*+\beta_2 E S^*+\mu_1 S^*\right]$

$+a_1\left[\beta_1 S I+\beta_2 S E-B_1 E\right]+a_1\left[-\beta_1 E^* \frac{S I}{E}-\beta_2 E^* S+B_1 E^*\right]$

$+a_2\left[v E-B_2 I\right]+a_2\left[-v I^* \frac{E}{I}+B_2 I^*\right]+a_3\left[q_1 E-B_3 Q_1\right]$

$+a_3\left[-q_1 Q_1^* \frac{E}{Q_1}+B_3 Q_1^*\right]+a_4\left[q_2 I+v Q_1-B_4 Q_2\right]$

$+a_4\left[-q_2 Q_2^* \frac{I}{Q_2}-v Q_2^* \frac{Q_1}{Q_2}+B_4 Q_2^*\right]$,

which is equivalent to the following expression,

$\frac{d \mathcal{L}}{d t}=\left[A+\mu_1 S^*+a_1 B_1 E^*+a_2 B_2 I^*+a_3 B_3 Q_1^*+a_4 B_4 Q_2^*\right]$

$-a_1 \beta_2 E^* S-\mu_1 S-\beta_1 SI+a_1 \beta_1 S I+\beta_2 E S^*-\beta_2 S E+$

$a_1 \beta_2 S E-a_1 B_1 E+a_2 v E+a_3 q_1 E+\beta_1 I S^*-a_2 B_2 I+$

$a_4 q_2 I+\varepsilon Q_1-a_3 B_3 Q_1+a_4 v Q_1-a_4 B_4 Q_2-A \frac{S^*}{S}-$

$\varepsilon Q_1 \frac{S^*}{S}-a_1 \beta_1 E^* \frac{S I}{E}+a_2 vI^* \frac{E}{I}-a_3 q_1 Q_4 \frac{E}{Q_1}-$

$a_4 q_2 Q_2 \frac{I}{Q_2}-a_4 v Q_2 \frac{Q_1}{Q_2}$.         (8)

To $\quad$ simplify $\quad(8), \quad$ define $\quad\left(\frac{S^*}{S}, \frac{E^*}{E}, \frac{I^*}{I}, \frac{Q_1^*}{Q_1}, \frac{Q_2^*}{Q_2}\right)=$ $\left(y_0, y_1, y_2, y_3, y_4\right)$, and let $D=A+\mu_1 S^*+a_1 B_1 E^*+$ $a_2 B_2 I^*+a_3 B_3 Q_1^*+a_4 B_4 Q_2^*$. Thus, we have

$\frac{d \mathcal{L}}{d t}=D-\left(\mu_1+a_1 \beta_2 E^*\right) \frac{1}{y_0} S^*$

$+\left(1-a_1\right) \beta_1 S^* I^* \frac{1}{y_0 y_2}-\left(1-a_1\right) \beta_2 S^* E^* \frac{1}{y_0 y_1}$

$+\left(\beta_2 S^*-a_1 B_1+a_2 v+a_3 q_1\right) E^* \frac{1}{y_1}$

$+\left(\beta_1 S^*-a_2 B_2+a_4 q_2\right) I^* \frac{1}{y_2}+\left(\varepsilon-a_3 B_3+a_4 v\right) Q_1^* \frac{1}{y_3}$

$-a_4 B_4 Q_2^* \frac{1}{y_4}-A y_0-\varepsilon Q_1^* \frac{y_0}{y_3}-a_1 \beta S^* I^* \frac{y_1}{y_0 y_2}$

$+a_2 v E^* \frac{y_2}{y_1}-a_3 q_1 E^* \frac{y_3}{y_1}-a_4 q_2 I^* \frac{y_4}{y_3}-a_4 Q_1^* \frac{y_4}{y_3}$.        (9)

Construct the function set $\gamma$

$\Upsilon=\left\{y_0, \frac{1}{y_0}, \frac{1}{y_1}, \frac{1}{y_2}, \frac{1}{y_3}, \frac{1}{y_4}, \frac{1}{y_0 y_2}, \frac{1}{y_0 y_1}, \frac{y_1}{y_0 y_2}, \frac{y_2}{y_1}, \frac{y_3}{y_1}, \frac{y_4}{y_2}, \frac{y_4}{y_3}, \frac{y_0}{y_3}\right\}$.

There are two groups corresponding to ϒ such that the product of all functions within each group is unity. The two groups are

$\left\{y_0, \frac{1}{y_0}\right\}$ and $\left\{y_0, \frac{y_1}{y_0 y_2}, \frac{y_2}{y_1}\right\}$.

Furthermore, associating with this groups, we can define the function

$\frac{d \mathcal{L}}{d t}=b_1\left(2-y_0-\frac{1}{y_0}\right)+b_2\left(3-y_0-\frac{y_1}{y_0 y_2}-\frac{y_2}{y_1}\right)$,        (10)

where, the coefficients b₁, b₂ are the unknown. We will determine coefficients a₁, a₂, a₃, and a₄and b₁ and b₂ such that (9) is equal to (10). Taking those equations together, we derive

$\left\{\begin{array}{l}1-a_1=0 \\ A=b_1+b_2 \\ \mu_1+a_1 \beta_2 E^*=b_1 \\ \beta_2 S^*-a_1 B_1+a_2 v+a_3 q_1=0 \\ \beta_1 S^*-a_2 B_2+a_4 q_2=0 \\ \varepsilon-a_3 B_3+a_4 v=0 \\ a_4 B_4 Q_2^*=0 \\ \varepsilon Q_1^*=0 \\ a_1 \beta_1 S^* I^*=b_2 \\ a_2 v E=b_2 \\ a_3 q_1 E^*=0 \\ a_4 q_2 I^*=0 \\ a_4 v Q_1^*=0\end{array}\right.$            (11)

From (11), by using some basic algebraic manipulations, we find

$a_1=1, a_3=0, a_4=0, a_2=\frac{\beta_1 S^*}{B_2}, b_1=\mu_1+a_1 \beta_2 E^*=\mu_1+\beta_2 E^*, b_2=a_1 \beta_1 S^* I^*=\beta_1 S^* I^*.$             (12)

Subtituting (12) into (11), we have

$\frac{d \mathcal{L}}{d t}=\left(\mu_1+\beta_2 E^*\right)\left(2-y_0-\frac{1}{y_0}\right)+\left(\beta_1 S^* I^*\right)(3-\left.y_0-\frac{y_1}{y_0 y_2}-\frac{y_2}{y_1}\right) \leq 0$                (13)

From the principle of arithmetical inequalities and geometrical means, we obtain

$y_0+\frac{y_1}{y_0 y_2}+\frac{y_2}{y_1} \geq 3 \sqrt{y_0 \frac{y_1}{y_0 y_2} \frac{y_2}{y_1}}$

$y_0+\frac{y_1}{y_0 y_2}+\frac{y_2}{y_1} \geq 3$

$3-y_0-\frac{y_1}{y_0 y_2}-\frac{y_2}{y_1} \leq 0$   

and

$y_0+\frac{1}{y_0} \geq 2 \sqrt{y_0 \frac{1}{y_0}}$

$y_0+\frac{1}{y_0} \geq 2$

$2-y_0-\frac{1}{y_0} \leq 0$.

Thus, it is proved that $\frac{d \mathcal{L}}{d t} \leq 0$ and $\frac{d \mathcal{L}}{d t}=0$ if $y_0=1$ and $y_1=y_2$. Since we can construct the Lyapunov function, $\mathcal{L}>$ 0 and $\frac{d \mathcal{L}}{d t} \leq 0$ so that the endemic equilibrium point $\left(K^*\right)$ reach globally asymptotically stable.

We proposed in Section 2 the mathematical model that represents the dynamics of population in dealing with COVID-19 spread including a case study using data in Central Java Province, Indonesia presented in Section 5. In the case study we demonstrated a numerical simulation based on data variables (S,E,I,Q₁,Q₂,H,R) Central Java Province, Indonesia from 02 April 2021 to 5 October 2021 (about six months) to estimate the parameters. One may wonder whether these six-month data are sufficient. In fact, six months of data collection are enough because COVID-19 spread is not related to the annual cycle. For one year data, it does not significantly affect the transmission pattern of COVID-19 compared with six-month data. Besides, the analysis regarding the mathematical proofs and stability were completely observed in this study. For comparison purposes, some other published references also used data with less than or equal to six months of the observation in their models and have produced good simulation results (see e.g. [2, 4, 5, 11, 13, 15, 22, 24-26, 28]).

Numerical simulation results showed that the model has an endemic equilibrium point. In early stages, the number of individuals exposed and infected with COVID-19 increased, this is shown by the basic reproductive number value of 2.15, which means that one infected person can on average transmit the COVID-19 to two susceptible people. To reduce the reproduction number value, it can be done by increasing the average rate of exposed and infected individuals by applying self-quarantine and increasing the transition rate from quarantined and infected individuals to hospitalized individuals. This means that the more individuals who are quarantined, the number of infected individuals can be significantly reduced. Likewise, the more infected individuals who are admitted to the hospital, the faster the recovery. From the simulation results, by increasing quarantine and hospitalization, we found R₀=0.2192<1. This condition results in an asymptotically stable non-endemic equilibrium point. This indicates that the transmission of COVID-19 is decreasing.

Based on the formulated model (1) this study still leaves open problems for future research prospects related to how to control the spread of COVID 19. Control variables could be added in the model that plays as strategies to reduce the number of exposed individual, quarantine individual, and infected individuals in the population and to increase the number of recovered individuals. Control variables could involve number of vaccinated individuals, self-precaution (wearing medical masks, washing hands/using hand sanitizers and physical distancing) and treatment. Furthermore, control methods that suit the dynamical model could be implemented in calculating the optimal strategy in controlling the spread of COVID-19. This will indeed produce more complicated differential equations and therefore further analysis regarding the existence and uniqueness of the optimal control are needed.

This work is supported by Diponegoro University, Indonesia, under RPI Research Grant with contract number 225-26/UN7.6.1/PP/2022. This research was supported by the Laboratory of Computer Modelling Mathematics Department, Faculty of Science and Mathematics, Diponegoro University, Semarang, Indonesia.