Analyzing the Asymptotic Behavior of an Extended SEIR Model with Vaccination for COVID-19

Abstract

Several research papers have attempted to describe the dynamics of COVID-19 based on systems of differential equations. These systems have taken into account quarantined or isolated cases, vaccinations, control measures, and demographic parameters, presenting propositions regarding theoretical results that often investigate the asymptotic behavior of the system. In this paper, we discuss issues that concern the theoretical results proposed in the paper “An Extended SEIR Model with Vaccination for Forecasting the COVID-19 Pandemic in Saudi Arabia Using an Ensemble Kalman Filter”. We propose detailed explanations regarding the resolution of these issues. Additionally, this paper focuses on extending the local stability analysis of the disease-free equilibrium, as presented in the aforementioned paper, while emphasizing the derivation of theorems that validate the global stability of both epidemic equilibria. Emphasis is placed on the basic reproduction number $R_0$, which determines the asymptotic behavior of the system. This index represents the expected number of secondary infections that are generated from an already infected case in a population where almost all individuals are susceptible. The derived propositions can inform health authorities about the long-term behavior of the phenomenon, potentially leading to more precise and efficient public measures. Finally, it is worth noting that the examined paper still presents an interesting epidemiological scheme, and the utilization of the Kalman filtering approach remains one of the state-of-the-art methods for modeling epidemic phenomena.

1. Introduction

The authors in [1] propose a compartmental model that contains a system of seven differential equations with the aim of describing the changing dynamics of the spread of COVID-19. The model divides the population into smaller parts, considering susceptible (S), exposed (E), infected (I), quarantined (Q), recovered (R), deceased (D), and vaccinated (V) cases. The first part of the aforementioned paper is dedicated to the proposal of theoretical results regarding the non-negativity of model’s states, the boundedness of the total population, and the existence and local stability of the disease-free equilibrium (DFE), based on the basic reproduction number, $R_0$.

However, there are certain issues regarding the presented proofs and formulas. In the present analysis, we aim to introduce the existing errors, providing detailed comments that rectify them. Moreover, emphasis is placed on proposing theorems concerning the global stability analysis of epidemic equilibria, in accordance with the above-mentioned scheme. In this way, we extend the theoretical results displayed in [1], while providing valuable information regarding the asymptotic behavior of the epidemiological system. This information can be employed for the establishment of more accurate measures that can facilitate the limiting of the virus’s transmission.

Typically, the spread of infectious diseases is explained through compartmental models, among which the SIR model—representing susceptible-infected-recovered—is the most recognized [2]. Consequently, numerous studies delve into the dynamics of COVID-19 relying on the SIR model or its adaptations like SIRS [3], SEIR, or the SEIRD models [4]. Furthermore, Malkov [5] proposed a deterministic SEIRS model that incorporates time-varying transmission rates for the description of the transmission of COVID-19. Other compartmental extensions can be found in [6,7,8,9,10,11]. All the abovementioned endeavors are based on systems of differential equations that can be numerically solved.

Several techniques have been proposed in the literature to establish numerically stable methodologies for solving systems of differential equations. Many papers employ Runge–Kutta methodologies, with the 4th-order Runge–Kutta being the most widely known [12,13]. Several extensions have been proposed in the articles of Kalogiratou and Monovasilis, which refer to two-derivative Runge–Kutta methods with optimal phase properties [14], optimized dispersion and dissipation error [15], and constant and frequency dependent coefficients [16]. Moreover, additional advanced techniques for solving systems of partial differential equations have been proposed in [17,18,19].

In summary, the present paper provides valuable corrections concerning the theoretical results displayed in [1] that pertain to the non-negativity and boundedness of a system of seven differential equations, which describe the transition of COVID-19 after the onset of the vaccination period. These modifications are crucial in validating the suitability of the epidemiological model for accurately describing the spread of COVID-19. More importantly, we provide novel properties regarding the global asymptotic stability of both the disease-free and endemic equilibria based on the values of the basic reproduction number ($R_0)$. These theoretical aspects are more crucial than the local stability analysis, offering insights into long-term behavior when the system approaches the aforementioned equilibria. Finally, a novel addition to the literature is the computation of the convergence rate to the endemic equilibria, offering a more comprehensive understanding of the system’s asymptotic behavior. Using real values for the basic reproduction number derived from experimental data, we can evaluate the severity of the phenomenon and validate previous predictions about the future course of the pandemic in the literature.

The rest of the article is structured as follows: In Section 2, we present a series of issues regarding the non-negativity and boundedness theorems that are proposed in [1], while Section 3 and Section 4 are dedicated to the rectification of issues concerning the local stability of the disease-free equilibrium and the existence and uniqueness of the endemic equilibrium, respectively. Finally, in Section 5, we present novel results regarding the global stability of the epidemic equilibria, while in Section 6, we conclude with the advantages of epidemiological modeling, emphasizing the main contribution of the present paper.

2. Non-Negativity of Model’s States and Boundedness of the Total Population

To begin, the authors in [1] have proposed an ODE system of seven equations to describe the transmission of COVID-19 after the opening of the vaccination period. As a result, the examined population has been split into seven compartments (classes) based on the state of the population’s members; Equation (1) displays the transitions between these classes, namely

$$\begin{gathered} \frac{dS\left(t\right)}{dt}=\Lambda -\beta S\left(t\right)I\left(t\right)-aS\left(t\right)-\mu S\left(t\right), \\ \frac{dE\left(t\right)}{dt}=\beta S\left(t\right)I\left(t\right)-\gamma E\left(t\right)+\sigma \beta V\left(t\right)I\left(t\right)-\mu E\left(t\right), \\ \frac{dI\left(t\right)}{dt}=\gamma E\left(t\right)-\delta I\left(t\right)-\mu I\left(t\right), \\ \frac{dQ\left(t\right)}{dt}=\delta I\left(t\right)-\left(1-\kappa \right)\lambda Q\left(t\right)-\kappa \rho Q\left(t\right)-\mu Q\left(t\right), \\ \frac{dR\left(t\right)}{dt}=\left(1-\kappa \right)\lambda Q\left(t\right)-\mu R\left(t\right), \\ \frac{dD\left(t\right)}{dt}=\kappa \rho Q\left(t\right), \\ \frac{dV\left(t\right)}{dt}=\alpha S\left(t\right)-\sigma \beta V\left(t\right)I\left(t\right)-\mu V\left(t\right), \end{gathered}$$

(1)

with non-negative initial conditions. In

Table 1. Parameter and state definition of the proposed SEIHCRDV model.

Symbol	Definition of Parameter/State
S	Susceptible
E	Exposed
I	Infectious
Q	Quarantined
R	Recovered
D	Deceased
V	Vaccinated
Λ	New births and new residents
$a$	Vaccination rate
$\beta$	Transmission rate
$\gamma$	Incubation rate
$\delta$	Infection rate
$\lambda$	Recovery rate
$\kappa$	Case fatality rate
$\mu$	Natural death rate
$\rho$	Death rate
$\sigma$	Vaccine inefficacy

, we present the definition of the system’s states and parameters.

In the first theorem of [1], the authors aim to prove the non-negativity of the system’s states based on the proposed system of differential equations. More specifically, an attempt to prove the non-negativity of the number of susceptible cases $S\left(t\right), \forall t\ge 0$, when $S_0>0$ is displayed. This attempt leads to

$$S\left(t\right)\ge S_0{e}^{-\mu t}\ge 0 .$$

(2)

However, this inequality does not seem to hold true when considering the first differential equation of the system. Equation

$$\frac{dS\left(t\right)}{dt}=\Lambda -\beta S\left(t\right)I\left(t\right)-aS\left(t\right)-\mu S\left(t\right)>-\mu S\left(t\right) ,$$

(3)

holds true only when $\Lambda >\beta S\left(t\right)I\left(t\right)+aS\left(t\right)$. It is evident that there are several instances of parameter selections where the aforementioned expression is not satisfied. We note that all system’s parameters are assumed to be positive constants, as they represent ingoing or outgoing transition rates of the system’s states. Therefore, a modification of (2) is required to lead to the desired outcome. Specifically, we take

$$\begin{gathered} \frac{dS\left(t\right)}{dt}=\Lambda -\beta S\left(t\right)I\left(t\right)-aS\left(t\right)-\mu S\left(t\right)>-\beta S\left(t\right)I\left(t\right)-aS\left(t\right)-\mu S\left(t\right) \\ =-\left(\beta {\rm I}\left(t\right)+a+\mu \right)S\left(t\right)\ge -\left(\beta \underset{t\in \left[0,\infty \right)}{\max }I\left(t\right)+a+\mu \right)S\left(t\right) . \end{gathered}$$

Using the infinity norm ${\Vert I\left(t\right)\Vert }_{\infty }=\underset{t\in \left[0,\infty \right)}{\max }I\left(t\right)$, we obtain

$$\frac{d\ln \left(S\left(t\right)\right)}{dt}\ge -\left(\beta {\Vert I\left(t\right)\Vert }_{\infty }+a+\mu \right).$$

Consequently, by integrating the above expression with respect to $t$, and substituting $t=0$, we result in

$$S\left(t\right)\ge S_0{e}^{-\left(\beta {\Vert I\left(t\right)\Vert }_{\infty }+a+\mu \right)t}\ge 0 , \forall t\ge 0 .$$

(4)

Notice that ${\Vert I\left(t\right)\Vert }_{\infty }<\infty$, as we refer to a finite population function, $N\left(t\right)$. We believe that this approach now rectifies the proof of Theorem 1. The utilization of the infinity norm can be employed for proving the non-negativity of vaccinated cases, too, while the proof for the remaining states is omitted for the sake of brevity.

In the second theorem of [1], the authors aim to prove the boundedness of the total population function $N\left(t\right), \forall t\ge 0$. They claim that since $N\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+Q\left(t\right)+R\left(t\right)+D\left(t\right)+V\left(t\right)$, for the derivative with respect to time, we have

$$\frac{dN\left(t\right)}{dt}=\frac{dS\left(t\right)}{dt}+\frac{dE\left(t\right)}{dt}+\frac{dI\left(t\right)}{dt}+\frac{dQ\left(t\right)}{dt}+\frac{dR\left(t\right)}{dt}+\frac{dD\left(t\right)}{dt}+\frac{dV\left(t\right)}{dt} ,$$

(5)

which leads to

$$\frac{dN\left(t\right)}{dt}=\Lambda -\mu {\rm N}\left(t\right), \forall t\ge 0 .$$

(6)

However, Equation (6) does not hold based on the proposed epidemiological system, which is presented in Equation (1). This derives from the inclusion of the deceased cases, $D\left(t\right)$, in the total population. The above expression should be rectified as

$$\frac{dN\left(t\right)}{dt}=\Lambda -\mu \left(S\left(t\right)+E\left(t\right)+I\left(t\right)+Q\left(t\right)+R\left(t\right)+V\left(t\right)\right)=\Lambda -\mu {\rm N}\left(t\right)+\mu D\left(t\right),$$

(7)

after the summation of all equations of the ODE system, as $\frac{dN\left(t\right)}{dt}=\frac{dS\left(t\right)}{dt}+\frac{dE\left(t\right)}{dt}+\frac{dI\left(t\right)}{dt}+\frac{dQ\left(t\right)}{dt}+\frac{dR\left(t\right)}{dt}+\frac{dD\left(t\right)}{dt}+\frac{dV\left(t\right)}{dt}.$ Thus, after moving $\mu {\rm N}\left(t\right)$ to the left side, we lead to

$$\frac{dN\left(t\right)}{dt}+\mu {\rm N}\left(t\right)=\Lambda +\mu D\left(t\right),$$

or

$$\frac{d{e}^{\mu t}N\left(t\right)}{dt}=\left(\Lambda +\mu D\left(t\right)\right)e^{\mu t}$$

and integrating with respect to $t$, we obtain

$$e^{\mu t}N\left(t\right)-N_0=\frac{\Lambda }{\mu }\left(e^{\mu t}-1\right)+\mu {{\displaystyle \int }}_0^{t}D\left(s\right)e^{\mu s}ds,$$

(8)

or

$$N\left(t\right)=\frac{\Lambda }{\mu }+\left(N_0-\frac{\Lambda }{\mu }\right)e^{-\mu t}+\mu e^{-\mu t}{{\displaystyle \int }}_0^{t}D\left(s\right)e^{\mu s}ds.$$

(9)

As a result, $N\left(t\right)$ is bounded if and only if ${{\displaystyle \int }}_0^{t}D\left(s\right)e^{\mu s}ds$ is bounded for all $t>0$.

Moreover, there is another major issue in the proof of Theorem 2 in [1]. The authors claim that $N\left(t\right)\le \frac{\Lambda }{\mu }$ for all $t>0$, regardless of the system’s parameters. According to Expression (6), which as we mentioned earlier is not true, the authors lead to

$$N\left(t\right)=\frac{\Lambda }{\mu }+\left({{\rm N}}_0-\frac{\Lambda }{\mu }\right)e^{-\mu t}, \forall t>0 .$$

(10)

Apparently, $N\left(t\right)\le \frac{\Lambda }{\mu }$ is not true for every parametric set. Based on Equation (10), this is valid for all $t>0$, only when $N_0<\frac{\Lambda }{\mu }$.

Finally based on (9) for $t\to \infty$, we lead to

$$\underset{t\to \infty }{\lim }N\left(t\right)=\frac{\Lambda }{\mu } ,$$

(11)

in case $\underset{t\to \infty }{\lim } e^{-\mu t}{{\displaystyle \int }}_0^{t}D\left(s\right)e^{\mu s}ds=0$.

As a result, it becomes evident that the authors’ proof for Theorem 2 does not validate the theorem’s statement.

3. Local Stability of the Disease-Free Equilibrium (DFE)

Moving on to Theorem 3, the authors claim that the disease-free equilibrium (DFE) is locally asymptotically stable if $R_0<1$ and unstable when $R_0>1$. First, we notice that during the computation of $R_0$, the vector $\mathit{W}\left(X\right)$ of Section 3.3 in [1] should be rectified to $\mathit{W}\left(X\right)={\left(\left(\mu +\delta \right){\rm E},-\gamma {\rm E}+\left(\mu +\delta \right){\rm I}\right)}^T$, since the number of infected cases is missing from the second component of the vector. This modification does not alter the final formula for $R_0,$ where

$$R_0=\frac{\beta \gamma \Lambda \left(\mu +\alpha \sigma \right)}{\mu \left(\mu +\gamma \right)\left(\mu +\delta \right)\left(\mu +\alpha \right)} .$$

(12)

To prove the local asymptotic stability of the DFE ${\mathit{X}}^0$, the respective Jacobian matrix of the epidemiological system is employed. Equilibrium ${\mathit{X}}^0$ is locally asymptotically stable when all six eigenvalues of the Jacobian $J\left({\mathit{X}}^0\right)$ are negative. So, it is claimed that there are two eigenvalues ${\lambda }_5, {\lambda }_6$, where

$${\lambda }_5=-\frac{1}{2}\left({\epsilon }_2+{\epsilon }_3+\sqrt{{\left({\epsilon }_2-{\epsilon }_3\right)}^2+4{\epsilon }_2{\epsilon }_3{R}_0}\right),$$

(13)

and

$${\lambda }_6=-\frac{1}{2}\left({\epsilon }_2+{\epsilon }_3-\sqrt{{\left({\epsilon }_2-{\epsilon }_3\right)}^2+4{\epsilon }_2{\epsilon }_3{R}_0}\right).$$

(14)

with ${\epsilon }_2=-\left(\alpha +\mu \right)$ and ${\epsilon }_3=-\left(\delta +\mu \right)$. While ${\lambda }_6$ is indeed smaller than 0 when $R_0<1$, the same does not hold true for Expression (13). More specifically, we have

$$\sqrt{{\left({\epsilon }_2-{\epsilon }_3\right)}^2+4{\epsilon }_2{\epsilon }_3{R}_0}<\sqrt{{\left({\epsilon }_2-{\epsilon }_3\right)}^2+4{\epsilon }_2{\epsilon }_3}=\left|{\epsilon }_2+{\epsilon }_3\right|=-\left({\epsilon }_2+{\epsilon }_3\right),$$

or

$${\epsilon }_2+{\epsilon }_3+\sqrt{{\left({\epsilon }_2-{\epsilon }_3\right)}^2+4{\epsilon }_2{\epsilon }_3{R}_0}<0 ,$$

or

$$-\frac{1}{2}\left({\epsilon }_2+{\epsilon }_3+\sqrt{{\left({\epsilon }_2-{\epsilon }_3\right)}^2+4{\epsilon }_2{\epsilon }_3{R}_0}\right)>0 ,$$

leading to ${\lambda }_5>0$. According to the above, the DFE becomes asymptotically unstable when $R_0<1$, which contradicts with the statement of Theorem 3. Moreover, this outcome opposes several analyses in literature [6,8,9,12,20,21,22,23,24,25,26]. On the other hand, when $R_0$ is greater than 1, the DFE becomes asymptotically unstable as ${\lambda }_5<0$ and ${\lambda }_6>0$.

The aforementioned issues possibly derive from the form of the Jacobin matrix $J\left({\mathit{X}}^0\right),$ which is presented in [1] (Equation (20), Section 3.3), as there are several mistakes concerning the signs of the elements that take place on the matrix diagonal. The authors’ proof for Theorem 3 does not validate the theorem’s statement.

4. Existence and Uniqueness of the Endemic Equilibrium

Following the theorem that concerns the local stability of the DFE, the authors emphasize the existence and uniqueness of an endemic equilibrium, denoted by ${\mathit{X}}^{*}.$ To begin with, the expression of the endemic equilibrium should be rectified to ${\mathit{X}}^{*}=\left(S^{*},E^{*},I^{*},Q^{*},R^{*},V^{*}\right)$, as the number of diseased cases—and the respective differential equation—are excluded from the determination of the equilibrium.

At the first part of the proof, the authors in [1] describe the components of ${\mathit{X}}^{*}$ with respect to $I^{*}$ after adding the second and sixth equation of the system evaluated on the endemic equilibrium. In Section 3.4 of [1], the authors use the notations ${\epsilon }_1=\mu +\alpha$,${\epsilon }_3=\mu +\delta$ and ${\epsilon }_4=\mu +\lambda \left(1-\kappa \right)+\kappa \rho$, and lead to expression

$$V^{*}=\frac{\Lambda \beta \gamma I+\Lambda \alpha \gamma -{\epsilon }_2{\epsilon }_3\left(\beta I+{\epsilon }_1\right)I}{\mu \gamma \left(\beta {\rm I}+{\epsilon }_1\right)} ,$$

(15)

which represents the number of vaccinated cases when the system has entered the endemic equilibrium. The endemic equilibrium is obtained after setting all derivatives of system (1) to zero. After solving with respect to $V^{*},$ we reach Expression (15).

First, the $I$ symbols should be replaced with $I^{*},$ as the system has to be evaluated at the endemic equilibrium ${\mathit{X}}^{*}$ to describe $V^{*}$. Also, considering the notation in Section 3.3 of [1], where ${\mathsf{\epsilon }}_2=-\left(\gamma +\mu \right)$, the minus sign of the numerator must be replaced with a plus sign.

Afterwards, the statement that $a_2$ is always positive and $a_0$ is negative when $R_0>1$ contradicts Formula (26) in that paper. It should be emphasized that the opposite behavior holds for these two quantities, namely $a_2<0$ and $a_0>0$.

Finally, we notice that an alternative, simpler formula can be derived for the number of vaccinated cases at the equilibrium. Using the sixth equation of the system evaluated at ${\mathit{X}}^{*},$ we culminate in $V^{*}=\frac{a{S}^{*}}{\mu +\sigma \beta {{\rm I}}^{*}}=\frac{\alpha \Lambda }{\left(\mu +\sigma \beta {{\rm I}}^{*}\right)\left(\beta {{\rm I}}^{*}+\mu +\alpha \right)} .$

5. Global Stability Analysis of Epidemic Equilibria

At this point we emphasize the extension of the results concerning the stability analysis of epidemic equilibria. Global stability analysis provides information about the behavior of a system across its entire state space. Therefore, it determines the stability of the system for all initial conditions. Thus, it offers a comprehensive view of the system’s behavior in contrast to the local stability analysis, which can provide insights only around the equilibria.

Theorem 1. 

The DFE${\mathit{X}}^0,$ is globally asymptotically stable if and only if $R_0<1$.

Proof of Theorem 1. 

Based on the LaSalle’s invariance principle, we choose a Lyapunov function $L\left(t\right)$ that is positive semidefinite in the feasible region $\Omega =\{\left(S,E,I,Q,R,V\right) | S, E,I,Q,R,V\ge 0\}={ℝ}_{+}^6$, while its derivative is negative definite in the same region. Let

$$L\left(t\right)=\frac{1}{2}\left[{\left(S-S^0\right)}^2+E^2+I^2+Q^2+R^2+{\left(V-V^0\right)}^2\right]=\frac{1}{2}{\mathit{X}}^T\mathit{X}\ge 0,$$

(16)

where $\mathit{X}={\left(S-S^0, E, I, Q, R, V-V^0\right)}^T$, as $E^0=I^0=Q^0=R^0=0$. Function $L\left(t\right)$ becomes 0 only on ${\mathit{X}}^0$. For the derivate, we obtain

$$\begin{array}{ccc}\hfill \frac{dL}{dt}& =& \left(S-S^0\right)\frac{dS}{dt}+E\frac{dE}{dt}+I\frac{dI}{dt}+Q\frac{dQ}{dt}+R\frac{dR}{dt}+\left(V-V^0\right)\frac{dV}{dt}\hfill \\ & =& \left(S-S^0\right)\left[\Lambda -\beta \left(S-S^0\right)I-\beta S^{0}I-\left(\alpha +\mu \right)\left(S-S^0\right)-\left(\alpha +\mu \right)S^0\right]\hfill \\ & & +E\left[-\left(\gamma +\mu \right)E+\beta \left(S-S^0\right)I+\beta S^{0}I+\sigma \beta (V-V^0){\rm I}+\sigma \beta V^{0}I\right]\hfill \\ & & +I\left[\gamma E-\left(\delta +\mu \right)I\right]+Q\left[\delta {\rm I}-\left(\left(1-\kappa \right)\lambda +\kappa \rho +\mu \right)Q\right]+R\left[\left(1-\kappa \right)\lambda Q-\mu R\right]\hfill \\ & & +\left(V-V^0\right)\left[a\left(S-S^0\right)+a{S}^0-\mu \left(V-V^0\right)-\mu V^0-\sigma \beta (V-V^0\right)I+\sigma \beta V^{0}I]\\ & =& -\beta {\left(S-S^0\right)}^2{\rm I}-\beta S^0\left(S-S^0\right)I-\left(\alpha +\mu \right){\left(S-S^0\right)}^2\hfill \\ & & -\left(\gamma +\mu \right){{\rm E}}^2+\beta \left(S-S^0\right){\rm E}I+\sigma \beta \left(V-V^0\right){\rm E}I+\beta S^0{\rm E}I+\sigma \beta V^0{\rm E}I\hfill \\ & & +\gamma {\rm I}E-\left(\delta +\mu \right){{\rm I}}^2+\delta QI-\left(\left(1-\kappa \right)\lambda +\kappa \rho +\mu \right)Q^2+\left(1-\kappa \right)\lambda RQ-\mu R^2\hfill \\ & & +\alpha \left(V-V^0\right)\left(S-S^0\right)-\mu {\left(V-V^0\right)}^2-\sigma \beta {\rm I}{\left(V-V^0\right)}^2+\sigma \beta V^0\left(V-V^0\right)I\hfill \\ & =& {\mathit{X}}^T\mathit{A}\mathit{X}={\mathit{X}}^T{\mathit{A}}_1\mathit{X}+{\mathit{X}}^T{\mathit{A}}_2\mathit{X}, \hfill \end{array}$$

(17)

where

$${\mathit{A}}_1=\left(\begin{array}{cccccc}-\left(\alpha +\mu \right)& 0& -\beta S^0& 0& 0& 0\\ 0& -\left(\gamma +\mu \right)& \beta \left(S^0+\sigma V^0\right)& 0& 0& 0\\ 0& \gamma & -\left(\delta +\mu \right)& 0& 0& 0\\ 0& 0& \delta & -\left(\left(1-\kappa \right)\lambda +\kappa \rho +\mu \right)& 0& 0\\ 0& 0& 0& \left(1-\kappa \right)\lambda & -\mu & 0\\ \alpha & 0& -\sigma \beta V^0& 0& 0& -\mu \end{array}\right),$$

and

$${\mathit{A}}_2=\left(\begin{array}{ccccccc}-\beta {\rm I}& 0& 0& 0& 0& 0& 0\\ \beta {\rm I}& 0& 0& 0& 0& 0& \sigma \beta {\rm I}\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& -\sigma \beta {\rm I}\end{array}\right),$$

where ${\mathit{A}}_2$ is a polynomial matrix. We have dropped the $\left(t\right)$ notation for the sake of simplicity. Furthermore, we notice that $\Lambda -\left(\alpha +\mu \right)S^0=0$ and $a{S}^0-\mu V^0=0$. Quantities $\beta \left(S-S^0\right){\rm I}E, \sigma \beta \left(V-V^0\right){\rm I}{\rm E}$, $-\sigma \beta {\rm I}{\left(V-V^0\right)}^2$ and $-\beta {\left(S-S^0\right)}^2{\rm I}$ are matched with ${\mathit{A}}_2.$

The characteristic polynomial of the 6 × 6 matrix ${\mathit{A}}_1$ is

$$\begin{gathered} p\left(x\right)=\det \left(x{\mathit{I}}_6-{\mathit{A}}_1\right)=\left(x+\alpha +\mu \right){\left(x+\mu \right)}^2\left(x+\left(\left(1-\kappa \right)\lambda +\kappa \rho +\mu \right)\right) \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\left(x^2+\left(\gamma +\delta +2\mu \right)x+\left(\mu +\delta \right)\left(\mu +\gamma \right)-\beta \gamma \left(S^0+\sigma V^0\right)\right)=0, \end{gathered}$$

(18)

leading to 4 negative eigenvalues. Now, for the second-order polynomial in (18), we implement the 2nd-order Routh–Hurwitz criterion, where the roots of the polynomial, lay on the left-hand side of the complex plane when coefficients $\left(\gamma +\delta +2\mu \right)$ and $\left(\mu +\delta \right)\left(\mu +\gamma \right)-\beta \gamma \left(S^0+\sigma V^0\right)$ are both positive. Apparently, $\left(\gamma +\delta +2\mu \right)>0$. Then the Routh–Hurwitz criterion is satisfied when $\left(\mu +\delta \right)\left(\mu +\gamma \right)-\beta \gamma \left(S^0+\sigma V^0\right)>0,$ and based on the $R_0$ formula displayed in [1], this inequality is true if and only if $R_0<1$.

Ultimately, the eigenvalues of ${\mathit{A}}_2$ are all non-positive due to $I\left(t\right)$ being non-negative, (non-negativity of system’s states). In parallel, the eigenvalues of ${\mathit{A}}_1$ are all negative if and only if $R_0<1$. According to the above observations we get that ${\mathit{X}}^T{\mathit{A}}_1\mathit{X}<0,$ and ${\mathit{X}}^T{\mathit{A}}_2\mathit{X}\le 0$. To summarize, $\frac{dL}{dt}={\mathit{X}}^T\mathit{A}\mathit{X}<0$, if and only if $R_0<1$, which proves the global asymptotic stability of the DFE. □

Theorem 2. 

The endemic equilibrium ${\mathit{X}}^{*},$ is globally asymptotically stable when $R_0>1$.

Proof of Theorem 2. 

Following a similar approach to that of the previous theorem, we note that according to [1] the endemic equilibrium exists only when $R_0$ is greater than 1. We choose the quadratic Lyapunov function,

$$L\left(t\right)=\frac{1}{2}\left[{\left(S-S^{*}\right)}^2+{\left(E-E^{*}\right)}^2+{\left(I-I^{*}\right)}^2+{\left(Q-Q^{*}\right)}^2+{\left(R-R^{*}\right)}^2+{\left(V-V^{*}\right)}^2\right]=\frac{1}{2}{\mathit{Y}}^T\mathit{Y}\ge 0,$$

(19)

where $\mathit{Y}={\left(S-S^{*}, E-E^{*}, I-I^{*}, Q-Q^{*}, R-R^{*}, V-V^{*}\right)}^T$, and $L\left({\mathit{X}}^{*}\right)=0$. Employing the second equation of the proposed system of differential equations evaluated on the endemic equilibrium, we get

$$\beta S^{*}{I}^{*}+\sigma \beta V^{*}{I}^{*}-\gamma E^{*}-\mu E^{*}=0 ,$$

or

$$\beta \left(S^{*}+\sigma V^{*}\right)=\frac{\left(\gamma +\mu \right){{\rm E}}^{*}}{{{\rm I}}^{*}}=\frac{\left(\gamma +\mu \right)\left(\delta +\mu \right)}{\gamma }$$

(20)

.

For the derivative of the selected Lyapunov function, we have

$$\frac{dL}{dt}=\left(S-S^{*}\right)\frac{dS}{dt}+\left(E-E^{*}\right)\frac{dE}{dt}+\left(I-I^{*}\right)\frac{dI}{dt}+\left(Q-Q^{*}\right)\frac{dQ}{dt}+\left(R-R^{*}\right)\frac{dR}{dt}+\left(V-V^{*}\right)\frac{dV}{dt},$$

(21)

and after some algebraic manipulations we obtain

$$\begin{gathered} \hspace{1em}\hspace{1em}\hspace{1em}\begin{array}{cc}\hfill \frac{dL}{dt}=& \left(S-S^{*}\right)\left[-\left(a+\mu \right)\left(S-S^{*}\right)-\beta \left(S-S^{*}\right)I-\beta S^{*}\left(I-I^{*}\right)\right]\hfill \\ & +\left(E-{{\rm E}}^{*}\right)\left[-\left(a+\delta \right)\left(E-E^{*}\right)+\sigma \beta \left(V-V^{*}\right)I+\beta \left(S-S^{*}\right)I-\beta \left(S^{*}+\sigma V^{*}\right)I^{*}\right]\\ & +\left(I-I^{*}\right)\left[\gamma \left(E-E^{*}\right)-\left(\delta +\mu \right)\left(I-I^{*}\right)\right]\hfill \\ & +\left(Q-Q^{*}\right)\left[\delta \left({\rm I}-{{\rm I}}^{*}\right)-\left(\left(1-\kappa \right)\lambda +\kappa \rho +\mu \right)\left(Q-Q^{*}\right)\right]\hfill \\ & +\left(R-R^{*}\right)\left[\left(1-\kappa \right)\lambda \left(Q-Q^{*}\right)-\mu \left(R-R^{*}\right)\right]\hfill \end{array} \\ +\left(V-V^{*}\right)\left[\alpha \left(S-S^{*}\right)-\mu \left(V-V^{*}\right)-\sigma \beta \left(V-V^{*}\right)I-\sigma \beta V^{*}\left(I-I^{*}\right)\right], \end{gathered}$$

(22)

since according to [1] it holds that $E^{*}=\frac{\delta +\mu }{\gamma }{{\rm I}}^{*}$, $Q^{*}=\frac{\delta }{\left(1-\kappa \right)\lambda +\kappa \rho +\mu }{{\rm I}}^{*},$ $R^{*}=\frac{\left(1-\kappa \right)\lambda \delta }{\mu \left(\left(1-\kappa \right)\lambda +\kappa \rho +\mu \right)}{{\rm I}}^{*},$ and $S^{*}=\frac{\Lambda }{\beta {{\rm I}}^{*}+\mu +\alpha } .$ Additionally, we notice that $\Lambda -\left(\alpha +\mu \right)S^{*}+\beta S^{*}{I}^{*}=0,$ $\beta \left(S^{*}+\sigma V^{*}\right)I^{*}-\left(\gamma +\mu \right){{\rm E}}^{*}=0$, and $\alpha S^{*}-\sigma \beta V^{*}{{\rm I}}^{*}-\mu V^{*}=0$, which derive from Equation (22) of [1] (Section 3.4), leading to our Equation (22).

Now, Expression (22) can be represented in matrix form as

$$\frac{dL}{dt}={\mathit{X}}^T\mathit{B}\mathit{X}={\mathit{X}}^T{\mathit{B}}_1\mathit{X}+{\mathit{X}}^T{\mathit{B}}_2\mathit{X},$$

(23)

where

$${\mathit{B}}_1=\left(\begin{array}{cccccc}-\left(\alpha +\mu \right)& 0& -\beta S^{*}& 0& 0& 0\\ 0& -\left(\gamma +\mu \right)& 0& 0& 0& 0\\ 0& \gamma & -\left(\delta +\mu \right)& 0& 0& 0\\ 0& 0& \delta & -\left(\left(1-\kappa \right)\lambda +\kappa \rho +\mu \right)& 0& 0\\ 0& 0& 0& \lambda \left(1-\kappa \right)& -\mu & 0\\ \alpha & 0& -\sigma \beta V^{*}& 0& 0& -\mu \end{array}\right),$$

and

$${\mathit{B}}_2=\left(\begin{array}{ccccccc}-\beta {\rm I}& 0& 0& 0& 0& 0& 0\\ \beta {\rm I}& 0& \frac{\left(\gamma +\mu \right)\left(\delta +\mu \right)}{\gamma }& 0& 0& 0& \sigma \beta {\rm I}\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& -\sigma \beta {\rm I}\end{array}\right) .$$

Hence, the model’s endemic equilibrium is globally asymptotically stable when $R_0>1$. The eigenvalues of ${\mathit{B}}_2$ are all non-positive due to $I\left(t\right)$ being non-negative, while all eigenvalues of ${\mathit{B}}_1$ lie on the left complex plain since they represent negative real numbers. Consequently, we obtain ${\mathit{X}}^T{\mathit{B}}_1\mathit{X}<0,$ and ${\mathit{X}}^T{\mathit{B}}_2\mathit{X}\le 0$. Thus, $\frac{dL}{dt}={\mathit{X}}^T\mathit{B}\mathit{X}<0,$ when $R_0>1$, proving the global stability of the endemic equilibrium. □

Theorem 3. 

When $R_0<1$, the extended SEIR model converges exponentially to the DFE according to the maximum eigenvalue of matrix $\mathit{A}$. On the other hand, in case $R_0>1$, the system converges to the endemic equilibrium based on the maximum eigenvalue of matrix $\mathit{B}$.

Proof of Theorem 3. 

In order to determine the convergence rate of the suggested epidemiological model to the DFE, it is necessary to find a positive value for the parameter $k$ that fulfills the inequality

$$\frac{dL\left(t\right)}{dt}\le -kL\left(t\right),$$

(24)

where $L\left(t\right)$ still represents the Lyapunov function. Based on the above, we can lead to the epidemic system’s convergence rate, which is determined by $\frac{k}{2}$. The distinction of the two cases, $R_0<1$ and $R_0>1$, is included to ensure that the existence of the two examined endemic equilibria is satisfied, before we proceed to the investigation of their convergence rates. We require the most appropriate $k$ value that satisfies Expression (24). After substituting the formulas for $L\left(t\right)$ and its derivate with respect to time, we obtain

$${\mathit{X}}^T\mathit{A}\mathit{X}\le -k\frac{1}{2}{\mathit{X}}^T\mathit{X},$$

(25)

$${\mathit{X}}^T\left(\mathit{A}+\frac{k}{2}{\rm I}\right)\mathit{X}\le 0,$$

(26)

leading to the conclusion that matrix $\mathit{A}+\frac{k}{2}\mathit{I}$ must be negative semidefinite. For $R_0<1$, the eigenvalues $\left({\lambda }_i\right)$ of $\mathit{A}$ are all negative. Our goal is to determine the value of $k$, for which the eigenvalues $\left({\lambda }_i+\frac{k}{2}\right)$ of matrix $\mathit{A}+\frac{k}{2}\mathit{I}$ are also negative.

As a result, we culminate in ${\lambda }_i+\frac{k}{2}\le 0$, when $k\le -2{\lambda }_i$ for $i=1, \dots , 6,$ which leads to the selection of $k=-2\max \left\{{\lambda }_i, i=1,\dots ,6\right\}>0$. This validates that the convergence rate to DFE is equal to $-\max \left\{{\lambda }_i, i=1,\dots ,6\right\}>0$. Similarly, in case $R_0>1$ the convergence rate of the epidemiological model to the endemic equilibrium is based on the positive equivalent of the maximum eigenvalue of matrix $\mathit{B}.$ □

6. Conclusions

In this paper, we have identified several issues regarding the theoretical results that are presented in [1] and accounted for the non-negativity, boundedness, existence, and local stability of epidemic equilibria. Moreover, special emphasis is placed on examining the global stability analysis of the produced equilibria based on the LaSalle’s invariance principle, extending the theoretical investigation of the aforementioned paper.

It is important to underline that the global stability analysis can provide insights into the entire state space’s stability, not just a neighborhood around an equilibrium point. Global stability analysis is often more robust to uncertainties and parameter variations. It can reveal whether a system remains stable under a wide range of conditions, making it particularly valuable in fields like control theory and engineering, where parameter variations are common. More importantly, it reveals the long-term behavior of the system regardless of the initial condition.

At this point, we should emphasize that despite the aforementioned issues, the statistical methodology proposed in [1], has an important role in the field of mathematical modelling in epidemiology. Kalman filtering provides the best linear unbiased estimate of a system’s states in the presence of noise and uncertainty, while it optimally combines measurements and a priori system predictions [12,27,28]. It can adapt to changing system dynamics by adjusting the filter’s parameters. This makes it suitable for systems with time-varying characteristics. Additionally, Kalman filters are computationally efficient, making them applicable in real-time systems. Like the traditional Kalman filter, Ensemble Kalman filtering provides estimates of state uncertainty and consistency, aiding in decision-making processes. Also, by sampling from the state space it accomplishes the capturing of complex nonlinear dynamics and avoids filter divergence.

Global stability analysis of COVID-19 models provides crucial insights that are immensely valuable for practical applications. It helps in predicting the long-term behavior of the disease spread. Understanding the stability of the model equilibria allows for projections about the disease trajectory, aiding in preparedness and resource allocation. Also, by analyzing the stability of different equilibria within the models, researchers can assess the effectiveness of various intervention strategies. This insight guides policymakers in implementing control measures such as vaccination drives, social distancing, or lockdowns. It assists in resource allocation by estimating the potential severity and duration of the outbreak. Hospitals, medical supplies, and personnel can be strategically deployed based on the projected stability of the disease dynamics. Finally, analyzing the stability of the model against real-world data allows for model validation. Insights gained from the analysis can also contribute to refining the model by identifying areas where the model might deviate from observed patterns.

There are several analyses in the literature that propose the utilization of statistical methodologies like Kalman filters, aiming to provide estimations about the future state of the COVID-19 pandemic [1,12,29,30,31]. The employed statistical methodology holds promise in steering decisions concerning the short-term trajectory of the pandemic. Conversely, the stability analysis provided in our study furnishes insights into the extended patterns of the phenomenon, augmenting awareness around this public health concern in the long run. Therefore, we argue that both analyses present valuable insights into the pandemic, each offering unique viewpoints.

In future work, we find it intriguing to explore a hybrid epidemiological particle filter. This approach handles the uncertainty inherent in pandemic phenomena by integrating particle filtering, which offers an alternative way to address the uncertainties present in both the equations defining the state and the observations of such phenomena. Moreover, delving into the disease’s evolution using various stochastic methods like discrete or continuous time Markov chains holds significant promise aiming to examine interesting stochastic descriptors [32]. Finally, numerical methods for the computationally efficient solving of the ODE system can be investigated [33], as the establishment of methodologies of low complexity is always of interest in mathematical modelling [34,35,36].

Finally, in the case of COVID-19, given the ongoing circumstances, it remains difficult to curtail the transmission of the virus in the foreseeable future. The $R_0$ decreases during periods of lockdown, although it rises right after the easing of restrictions to values which are far higher than unity [37,38]. Also, even after the initialization of the vaccination campaigns, variants like omicron continue to spread rapidly [39,40].

Several variants have emerged even after the onset of the vaccination period, with the most widely known being the alpha, delta, and omicron variants, while the corresponding values for the delta variant ranged between 3.2 and 8 with a mean of 5.08 [41,42]. Moreover, according to the review of Liu and Rocklöv [43], the basic reproduction number for the omicron variant is 2.5 times greater than the respective reproduction number of the delta variant. Hence, according to the above comments, it becomes evident that the transmissibility of the virus will persist for quite a long-time interval. This perspective was strongly supported by many researchers even during the early stages of the pandemic [44,45]. Neither the establishment of lockdowns nor the vaccination campaigns, reduced the reproduction values less than unity for sufficiently long periods. Therefore, the eradication of the disease seems almost impossible.

As a result, public authorities may emphasize the reduction of severe infections, hospitalizations, and deaths as these are the main issues of concern for the entire population. Until now, this policy has shown a significant improvement of the confrontation against the pandemic’s drawbacks. Without a doubt, the systematic and timely vaccination of the population plays a pivotal role in realizing this objective.