A Kind of Time-Delayed COVID-19 Dynamical Model with Vaccination

Cheng’ao Li, Junliang Lu^*

School of Statistics and Mathematics, Yunnan University of Finance and Economics, Kunming, China.

DOI: 10.4236/am.2022.134025                

Abstract

In the paper, we study a kind of time-delayed novel coronavirus pneumonia dynamical model with vaccination. This model considers that people are vaccinated, and the human immune system has a series of processes, which need a certain time. We first obtain the disease-free equilibrium and the basic reproduction number R₀, and the system has a unique endemic equilibrium when R₀ > 1. Then we discuss the stability of the disease-free equilibrium and the endemic equilibrium with different delays τ. For τ = 0, using the Lyapunov function approach, we obtained the stability of disease-free equilibrium and the endemic equilibrium, respectively. For any delay τ ≠ 0, using the Routh-Hurwitz Criteria, we obtained that the disease-free equilibrium is locally asymptotically stable. We also find the critical value τ₀ at the endemic equilibrium, and obtain the condition that the system has a Hopf bifurcation at the endemic equilibrium. Finally, with the suitable choices of the parameters, some numerical simulations are presented in order to verify the effectiveness of the obtained theoretical results.

Keywords

SVIQR Epidemic Model, the Time-Delayed COVID-19 Dynamical Model, Lyapunov Functional, Hopf Bifurcation, Stability

1. Introduction

In December 2019, a kind of new atypical viral pneumonia cases, which is caused by a novel coronavirus (2019-nCoV), appeared in Wuhan City, Hubei Province of China. Since then, the novel coronavirus pneumonia (COVID-2019) has rapidly spread to other provinces of China and other countries with the Spring Festival passenger flew. The number of suspected/confirmed cases has increased rapidly [1]. The outbreak of COVID-2019 has attracted the attention of the countries all over the world. According to the cases of COVID-19 in Wuhan and Hubei, the number of COVID-19 is a significant upward trend [2] [3]. There is already abundant information on 2019-nCoV and its spread worldwide [4] [5], as well as the risks it poses to the population. Globally, as of 5:31 pm CET, 10 January 2022, there have been 305,914,601 confirmed cases of COVID-19, including 5,486,304 deaths, reported to World Health Organization (WHO) [6] [7]. Due to the sudden outbreak of COVID-2019, it is a new type of infectious disease, people have not found a new measure to analyze, prevent and treat it. Every outbreak of the coronavirus has caused huge economic losses to the affected countries. The epidemic of infectious diseases in history has seriously affected people’s health and led to huge losses to the global economy [8] [9]. At the same time, the high mortality rate of the epidemic has also caused varying degrees panic. COVID-2019 spreads mainly through human contact. Therefore, governments throughout China and other countries have formulated policies for controlling populations to cut off the spread of the virus and to reduce the infection rate. At the same time, scientific medical treatment can effectively alleviate the epidemic situation and eventually terminate the epidemic situation after a period. However, the long-term population control has brought many inconveniences to the human daily life, and it has brought a great impact on the economy. Therefore, it is very important that the government makes the timely and effective treatment measures to eliminate COVID-2019. It is necessary to use traditional infectious disease dynamics models to analyze and predict the development of the epidemic. These can be achieved by establishing mathematical models that simulate the spread of the epidemic. Through the study of related mathematical models, the relevant parameters obtained through historical data inversion can measure the current situation of the epidemic in many ways, and numerical calculations from the model. Simulation can also predict the development trend of the epidemic, all of which provide effective reference value for proposing better and more scientific epidemic prevention and control schemes.

The classical mathematical models describing the spread of infectious diseases include SI, SIS, SIR, SIRS, SEIR, and SEIRS model, etc [10] [11] [12]. Its main idea is to divide the population into different compartments: susceptible(S), latent (E), infectious (I), and removed (R), and to establish the ordinary differential equation systems through the communication mechanism of one group transferring to another group, so as to reveal the law of the epidemic spread. These models (and their variants) are used to study the transmission of various infectious diseases such as measles, smallpox, rabies, Ebola virus, and population dynamics. After the SARS epidemic in 2003, the mathematical models of SARS, MERS and other coronavirus transmission laws were gradually enriched [13] [14].

Because most traditional models do not consider the effect of latency on propagation delay, even considering into the latent (E) in the SEIR models of COVID-19, it is only assumed that the latent E is weakly infectious, which cannot depict the characteristic that 2019-nCoV can propagate in the latent period. In addition, the above model hardly considers the data delay caused by the time required for diagnosis. To sum up, a model of COVID-19 with time delay effect is needed to describe the novel crown pneumonia epidemic more appropriately. At the same time, with the development of the scientific medical treatment, there are many immune drugs. As of 9 January 2022, 9,126,987,353 vaccine doses have been administered [7]. The infectious disease dynamics model based on time-delayed dynamics system with vaccination proposed in this paper inverse the model parameters through public historical data. The parameter analysis shows the effectiveness of prevention and control measures of governments at all levels. Based on these parameters, we well simulated the development of the current epidemic and accurately predicted the future trend of the epidemic. The simulation results show that the epidemic situation can be controlled and gradually ended in a short time.

The arrangement of the paper is as follows. Firstly, using the SVIR model in [15] as our baseline model, we formulate the time-delayed COVID-19 dynamical model with vaccination, and show the model biologically reasonable and mathematically defined well in Section 2. We obtain the disease-free equilibrium, the endemic equilibrium, and a formula for the reproduction number, $R_0$, of the time-delayed COVID-19 dynamical model with vaccination by the next generation matrix [16] [17] in Section 3. In Section 4, when $R_0>1$, we prove that the model (5) has a unique positive equilibrium. Then we discuss the eigenvalues of the model at the equilibrium through Routh-Hurwitz criteria [18] for $R_0>1$ and $R_0<1$, respectively, so as to judge the stability of the disease-free equilibrium and the endemic equilibrium. When the delay $\tau =0$, we construct Lyapunov function [19] to judge the global stability of the equilibrium. According to the conditions for the existence of Hopf bifurcation of delay differential equations [20] [21], we find that the system has Hopf bifurcation when $R_0>1$ and $\tau ={\tau }_0$. Finally, we carry out numerical simulation and get different simulation results according to different data to verify the correctness of the theoretical results in Section 5. We finally give brief discussions of our findings and our future research plan in Section 6.

2. The Model Formulation

2.1. The Building of the Delay Dynamical Model with Vaccination

From the study of infectious diseases, prevention is more important than treatment. The researchers have established mathematical models to show the dynamic behavior of infectious diseases based on the occurrence, transmission and development law of infectious diseases and related social factors. Theoretical analysis and numerical simulation are also carried out through the model [15] [22] [23] [24] [25]. It can show the development process of infectious diseases by simulating dynamic behavior, analyzing the causes and key factors of infectious diseases by studying the impact of prevention and control methods on the epidemic of infectious diseases, predict the epidemic law and development trend, seek for the development law of various infectious diseases, and find the optimal solution of disease prevention and control. These researchers also provide theoretical and quantitative basis for prevention as well as decision-making. What's more, vaccination is the most effective way to control the spread of the epidemic [26] [27] [28]. If people are vaccinated, the human immune system has a series of processes, which need a certain time [29] [30]. And if people can only have limited immune protection after vaccination [31], the vaccinators are not susceptible (S) and recovered (R) individuals, but belong to a new group. According to this situation, we build a time-delayed SVIQR model with incomplete vaccination based on the SVIR model established by Liu et al. [15].

$\{\begin{array}{l}\frac{\text{d}S\left(t\right)}{\text{d}t}=\Lambda -\beta S\left(t\right)I\left(t\right)-\mu S\left(t\right)-\rho S\left(t-\tau \right),\hfill \\ \frac{\text{d}V\left(t\right)}{\text{d}t}=\rho S\left(t-\tau \right)-\sigma \beta V\left(t\right)I\left(t\right)-\mu V\left(t\right)-\alpha V\left(t\right),\hfill \\ \frac{\text{d}I\left(t\right)}{\text{d}t}=\beta \left(S\left(t\right)+\sigma V\left(t\right)\right)I\left(t\right)-\left(\mu +\nu +\delta +\kappa \right)I\left(t\right),\hfill \\ \frac{\text{d}Q\left(t\right)}{\text{d}t}=\kappa I\left(t\right)-\left(\mu +\gamma \right)Q\left(t\right),\hfill \\ \frac{\text{d}R\left(t\right)}{\text{d}t}=\alpha V\left(t\right)+\nu I\left(t\right)+\gamma Q\left(t\right)-\mu R\left(t\right),\hfill \end{array}$ (1)

with the initial conditions

$\begin{array}{l}S\left(\theta \right)={\phi }_1\left(\theta \right),V\left(\theta \right)={\phi }_2\left(\theta \right),I\left(\theta \right)={\phi }_3\left(\theta \right),Q\left(\theta \right)={\phi }_4\left(\theta \right),R\left(\theta \right)={\phi }_5\left(\theta \right),\\ {\phi }_i\left(\theta \right)\in C\left[-\tau ,0\right),{\phi }_i\left(\theta \right)\ge 0,{\phi }_i\left(0\right)>0,i=1,2,3,4,5.\end{array}$

where $S\left(t\right),V\left(t\right),I\left(t\right),Q\left(t\right)$ and $R\left(t\right)$ are the density of the susceptible, vaccinated, infected, quarantined and recovered population at time t, respectively. $\tau$ is the time required for the vaccine to produce protection. $\Lambda$ is the constant entering flux into the S-class, $\beta$ is the transmission rates for individuals in S-class, $\mu$ is the constant natural death rate, $\alpha$ is complete vaccination rate, $\rho$ is the vaccination rate, $\kappa$ is isolation rate, $\delta$ is the death rate because of CODIV-19. $\frac{1}{\nu }$, $\frac{1}{\gamma }$ are the average time of the infectious period and quarantine, respectively. $\sigma \beta$ is the transmission rates for individuals in V-class, $0\le \sigma \le 1$. When $\sigma =0$ and $\sigma =1$, it means that vaccination is completely effective and ineffective, respectively.

2.2. Positive Invariant Sets of the Model System

Below, we show that the system (1) is biologically reasonable and mathematically defined well. For the simplicity, we give the following notations: ${ℝ}_{+}^5=\left\{S\left(t\right)\ge 0\right\}\cup \left\{V\left(t\right)\ge 0\right\}\cup \left\{I\left(t\right)\ge 0\right\}\cup \left\{Q\left(t\right)\ge 0\right\}\cup \left\{R\left(t\right)\ge 0\right\}$ for the nonnegative orthant of space ${ℝ}^5$, and $\mathrm{int}\left({ℝ}_{+}^5\right)={ℝ}_{+}^5\backslash \left\{S\left(t\right)=V\left(t\right)=I\left(t\right)=Q\left(t\right)=R\left(t\right)=0\right\}$ for the interior of the ${ℝ}_{+}^5$. From the system (1), we have

$\begin{array}{l}S\left(t\right)=\exp \left(-{\displaystyle {\int }_0^t}\left(\mu +\beta I\left(\eta \right)\right)\text{d}\eta \right)\\ \times \left(S\left(0\right)+{\displaystyle {\int }_0^t}\left(\Lambda -\rho S\left(\eta -\tau \right)\right)\exp \left({\displaystyle {\int }_{\eta }^t}\left(\mu +\beta I\left(\xi \right)\right)\text{d}\xi \right)\text{d}\eta \right)\mathrm{,}\\ V\left(t\right)=\exp \left(-{\displaystyle {\int }_0^t}\left(\mu +\alpha +\sigma \beta I\left(\eta \right)\right)\text{d}\eta \right)\\ \times \left(V\left(0\right)+{\displaystyle {\int }_0^t}\rho S\left(\eta -\tau \right)\exp \left({\displaystyle {\int }_{\eta }^t}\left(\mu +\alpha +\sigma \beta I\left(\xi \right)\right)\text{d}\xi \right)\text{d}\eta \right)\mathrm{,}\end{array}$

$\begin{array}{l}I\left(t\right)=I\left(0\right)\exp \left(-{\displaystyle {\int }_0^t}\left(\beta S\left(\eta \right)+\sigma V\left(\eta \right)-\left(\mu +\upsilon +\delta +\kappa \right)\right)\text{d}\eta \right)\mathrm{,}\\ Q\left(t\right)=\exp \left(-\left(\mu +\gamma \right)t\right)\left(Q\left(0\right)+{\displaystyle {\int }_0^t}\kappa I\left(\eta \right)\exp \left(\left(\mu +\gamma \right)\left(t-\eta \right)\right)\text{d}\eta \right)\mathrm{,}\\ R\left(t\right)=\exp \left(-\mu t\right)\left(R\left(0\right)+{\displaystyle {\int }_0^t}\left(\alpha V\left(\eta \right)+\nu I\left(\eta \right)+\gamma Q\left(\eta \right)\right)\exp \left(\mu \left(t-\eta \right)\right)\text{d}\eta \right)\mathrm{.}\end{array}$ (2)

If we suppose that $\left(S\left(0\right)\mathrm{,}V\left(0\right)\mathrm{,}I\left(0\right)\mathrm{,}Q\left(0\right)\mathrm{,}R\left(0\right)\right)\in \mathrm{int}\left({ℝ}_{+}^5\right)$, according to the third equation of (2), for all $t\ge 0$, $I\left(t\right)>0$. Because for any $t\ge 0$, $\Lambda -\rho S\left(t-\tau \right)>0$, hence, from other equation of (2), we have $S\left(t\right)>0$, $V\left(t\right)>0$, $Q\left(t\right)>0$, $R\left(t\right)>0$.

From the initial conditions, for all $t\in \left[-\tau \mathrm{,0}\right)$, $S\left(t\right)>0$, $V\left(t\right)>0$, $I\left(t\right)>0$, $Q\left(t\right)>0$, $R\left(t\right)>0$, thus, $S\left(t\right)>0$, $V\left(t\right)>0$, $I\left(t\right)>0$, $Q\left(t\right)>0$, $R\left(t\right)>0$ for any $t\in \left[-\tau \mathrm{,}+\infty \right)$.

In addition, we suppose that the total population $N\left(t\right)=S\left(t\right)+V\left(t\right)+I\left(t\right)+Q\left(t\right)+R\left(t\right)$. Thus, adding all the equations for $S\left(t\right),V\left(t\right),I\left(t\right),Q\left(t\right)$, and $R\left(t\right)$ in (1) together, we get

$\frac{\text{d}N\left(t\right)}{\text{d}t}=\Lambda -\mu N\left(t\right)-\delta I\left(t\right)\mathrm{.}$

Because of $N\left(t\right)>I\left(t\right)>0$, we obtain

$\Lambda -\left(\mu +\delta \right)N\left(t\right)\le \frac{\text{d}N\left(t\right)}{\text{d}t}\le \Lambda -\mu N\left(t\right)\mathrm{.}$

Consequently, for any $t\in \left[\mathrm{0,}+\infty \right)$,

$\begin{array}{l}\frac{\Lambda }{\mu +\delta }+\left(N\left(0\right)-\frac{\Lambda }{\mu +\delta }\right)\exp \left(-\left(\mu +\delta \right)t\right)\le N\left(t\right)\\ \le N^0+\left(N\left(0\right)-N^0\right)\exp \left(-\mu t\right)\mathrm{,}\end{array}$ (3)

where $N^0$ is the initial population, $N\left(0\right)$ is the population at time $t=0$, which yields that $N\left(0\right)\exp \left(-\left(\mu +\delta \right)t\right)\le N\left(t\right)\le N^0+\left|N\left(0\right)-N^0\right|$, for $t\in \left[\mathrm{0,}{T}^{+}\right)$, where $T^{+}>0$. Thus, we obtain that $N\left(t\right)$ is bounded on $t\in \left[\mathrm{0,}{T}^{+}\right)$, which implies that $T^{+}=+\infty$. Let $t\to +\infty$ in (3), we obtain

$\frac{\Lambda }{\mu +\delta }\le \underset{t\to +\infty }{\lim \inf }N\left(t\right)\le \underset{t\to +\infty }{\lim \sup }N\left(t\right)\le N^0.$

For showing that the system (1) is biologically reasonable and mathematically defined well, we need to prove that the set

$D=\left\{\left(S\left(t\right),V\left(t\right),I\left(t\right),Q\left(t\right),R\left(t\right)\right)\in \mathrm{int}\left({ℝ}_{+}^5\right):\frac{\Lambda }{\mu +\delta }\le N\left(t\right)\le N^0\right\}$ (4)

is positively invariant for the system (1). Let $\left(S\left(0\right)\mathrm{,}V\left(0\right)\mathrm{,}I\left(0\right)\mathrm{,}Q\left(0\right)\mathrm{,}R\left(0\right)\right)\in D$, then $S\left(t\right),V\left(t\right),I\left(t\right),Q\left(t\right)$, and $R\left(t\right)$ are defined on $\left[\mathrm{0,}+\infty \right)$, with (3) holding for $N\left(t\right)$. Since $\left(S\left(0\right)\mathrm{,}V\left(0\right)\mathrm{,}I\left(0\right)\mathrm{,}Q\left(0\right)\mathrm{,}R\left(0\right)\right)\in D$, it follows from (3) that $\frac{\Lambda }{\mu +\delta }\le N\left(t\right)\le N^0$ for all $t\in \left[\mathrm{0,}+\infty \right)$. This means that D is a positively invariant set. Therefore, we have

Theorem 2.1. The set D is a positively invariant set for the model system (1).

Because the total population $N\left(t\right)=S\left(t\right)+V\left(t\right)+I\left(t\right)+Q\left(t\right)+R\left(t\right)$, the model system (1) is equivalent to the following model system,

$\{\begin{array}{l}\frac{\text{d}S\left(t\right)}{\text{d}t}=\Lambda -\beta S\left(t\right)I\left(t\right)-\mu S\left(t\right)-\rho S\left(t-\tau \right),\hfill \\ \frac{\text{d}V\left(t\right)}{\text{d}t}=\rho S\left(t-\tau \right)-\sigma \beta V\left(t\right)I\left(t\right)-\mu V(t)-\alpha V\left(t\right),\hfill \\ \frac{\text{d}I\left(t\right)}{\text{d}t}=\beta \left(S\left(t\right)+\sigma V\left(t\right)\right)I\left(t\right)-\left(\mu +\nu +\delta +\kappa \right)I\left(t\right),\hfill \\ \frac{\text{d}Q\left(t\right)}{\text{d}t}=\kappa I\left(t\right)-\left(\mu +\gamma \right)Q\left(t\right),\hfill \end{array}$ (5)

with the initial conditions and the positively invariant region $\Omega$ as follows,

$\begin{array}{l}S\left(\theta \right)={\phi }_1\left(\theta \right),V\left(\theta \right)={\phi }_2\left(\theta \right),I\left(\theta \right)={\phi }_3\left(\theta \right),Q\left(\theta \right)={\phi }_4\left(\theta \right),\\ {\phi }_i\left(\theta \right)\in C\left[-\tau ,0\right),{\phi }_i\left(\theta \right)\ge 0,{\phi }_i\left(0\right)>0,i=1,2,3,4\\ \Omega =\left\{\left(S\left(t\right),V\left(t\right),I\left(t\right),Q\left(t\right)\right)\in \mathrm{int}\left({ℝ}_{+}^4\right):\frac{\Lambda }{\mu +\delta }\le N\left(t\right)\le N^0\right\}.\end{array}$ (6)

3. The Equilibria of the Model System

When $t\to +\infty$, $S\left(t-\tau \right)=S\left(t\right)=S$, $V\left(t\right)=V$, $I\left(t\right)=I$, $Q\left(t\right)=Q$. Equilibrium points of system (5) are given by solutions $\left(S\mathrm{,}V\mathrm{,}I\mathrm{,}Q\right)\in {ℝ}_{+}^4$ of

$\{\begin{array}{l}\Lambda -\beta SI-\mu S-\rho S=0,\hfill \\ \rho S-\sigma \beta VI-\mu V-\alpha V=0,\hfill \\ \beta \left(S+\sigma V\right)I-\left(\mu +\nu +\delta +\kappa \right)I=0,\hfill \\ \kappa I-\left(\mu +\gamma \right)Q=0.\hfill \end{array}$ (7)

We easily obtain the disease-free equilibrium point $E^0=\left(S_0\mathrm{,}{V}_0\mathrm{,}{I}_0\mathrm{,}{Q}_0\right)$ of system (5) in the situation, where

$S_0=\frac{\Lambda }{\rho +\mu },V_0=\frac{\rho \Lambda }{\left(\rho +\mu \right)\left(\mu +\alpha \right)},I_0=0,Q_0=0.$

The basic reproduction number $R_0$ is known as the expected number of secondary cases produced by a typical infectious individual during its entire period of infectiousness in a completely susceptible population [17] [32] [33]. For system (5), we can find the new infected person matrix $\Phi$ and the removed infected person matrix $\Psi$, respectively:

$\Phi =\beta \left(S\left(t\right)+\sigma V\left(t\right)\right)I\left(t\right),\Psi =\left(\mu +\nu +\delta +\kappa \right)I\left(t\right),$

and hence, the basic reproduction number $R_0$ is

$R_0={\rho \left(\Phi {\Psi }^{-1}\right)|}_{E^0}=\frac{\beta \left(S_0+\sigma V_0\right)}{\mu +\nu +\delta +\kappa }=\frac{\Lambda \beta \left(\mu +\alpha +\sigma \rho \right)}{\left(\rho +\mu \right)\left(\mu +\alpha \right)\left(\mu +\nu +\delta +\kappa \right)}.$

Let $c=\mu +\nu +\delta +\kappa$, the basic reproduction number can be written as

$R_0=\frac{\Lambda \beta \left(\mu +\alpha +\sigma \rho \right)}{c\left(\rho +\mu \right)\left(\mu +\alpha \right)}\mathrm{.}$

Theorem 3.1. Let $t\to +\infty$, then if $R_0<1$, the system (5) has no positive equilibria; if $R_0>1$, the system (5) has a unique positive equilibrium point $E^{*}=\left(S^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{V}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{I}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{Q}^{\mathrm{<mo>\; *\; </mo>}}\right)$.

Proof. By the first equation of (7), we get $S=\frac{\Lambda }{\rho +\mu +\beta I}$. By the second equation of (7), we get $V=\frac{\rho S}{\mu +\alpha +\sigma \beta I}=\frac{\rho \Lambda }{\left(\rho +\mu +\beta I\right)\left(\mu +\alpha +\sigma \beta I\right)}$. By the fourth equation of (7), we get $Q=\frac{\lambda I}{\mu +\gamma }$. Then we substitute $S\mathrm{,}V$ and $Q$ into the third equation of (7) and simplify it:

$H\left(I\right)=h_1{I}^2+h_{2}I+h_3=0,$ (8)

where

$\begin{array}{l}{h}_1=-\sigma {\beta }^{2}c<0,\\ h_2=\beta \left(\Lambda \beta \sigma -c\left(\mu +\alpha +\sigma \left(\rho +\mu \right)\right)\right),\\ h_3=c\left(\mu +\alpha \right)\left(\rho +\mu \right)\left(R_0-1\right).\end{array}$ (9)

If $R_0<1$, it is ease to prove that $h_2<0,h_3<0$. Hence the equation $H\left(I\right)=0$ does not have positive roots, so the (5) has no positive equilibrium.

If $R_0>1$, then $h_3>0$. Hence the equation $H\left(I\right)=0$ has a positive root and a negative root. We take the positive root as $I^{\mathrm{<mo>\; *\; </mo>}}$ and obtain the positive equilibrium point $E^{*}=\left(S^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{V}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{I}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{Q}^{\mathrm{<mo>\; *\; </mo>}}\right)$ of the system (5), where

$\begin{array}{l}{S}^{*}=\frac{\Lambda }{\rho +\mu +\beta I^{*}},\\ V^{*}=\frac{\rho \Lambda }{\left(\rho +\mu +\beta I^{*}\right)\left(\mu +\alpha +\sigma \beta I^{*}\right)},\\ Q^{*}=\frac{\lambda I^{*}}{\mu +\gamma },\end{array}$ (10)

and $I^{\mathrm{<mo>\; *\; </mo>}}$ is the unique positive root of the quadratic equation for $R_0>1$. □

4. The Stability of the Equilibria of the System (5)

Now, we research the stability of the equilibria of the system (5).

4.1. Global Stability of the Equilibria When $\tau =0$

At first, we give the definition of the Lyapunov-candidate-function to research the stability of the equilibria. Let $L\mathrm{<mo>\; :\; </mo>}{ℝ}^n\to ℝ$ be a continuous scalar function, if L satisfies

$L\left(0\right)=\mathrm{0,}L\left(x\right)>\mathrm{0,}\forall x\in U\backslash \left\{0\right\}\mathrm{,}$

with U being a neighborhood region around $x=0$, then L is called a Lyapunov-candidate-function.

Next, we state the general Lyapunov stability theorem [34].

Lemma 4.1. Consider the following vector field

$\stackrel{\dot{}}{x}=f\left(x\right)\mathrm{,}x\in {ℝ}^n\mathrm{.}$ (11)

Let $x_0$ be a equilibrium point of (11) and let $L\mathrm{<mo>\; :\; </mo>}U\to ℝ$ be a $C^1$ function defined on some neighborhood U of $x_0$ such that

1) $L\left(x_0\right)=0$ and $L\left(x\right)>$ if $x\ne x_0$,

2) $\stackrel{\dot{}}{L}\left(x\right)\le 0$ in $U\backslash \left\{x_0\right\}$,

then $x_0$ is stable. Moreover, if

3) $\stackrel{\dot{}}{L}\left(x\right)<0$ in $U\backslash \left\{x_0\right\}$,

then $x_0$ is asymptotically stable.

The proof and the geometric interpretation of the Lemma 4.1 can be found in the second chapter in reference [34].

Now, we prove the stability of the equilibrium points of the system (5).

Theorem 4.1. If $R_0<1$, $\tau =0$, then the disease-free equilibrium $E^0$ of system (5) is globally asymptotically stable in $\Omega$.

Proof. We transform the disease-free equilibrium $E^0$ into the origin by transforming $X=S-S_0$, $Y=V-V_0$, $I=I$, $Q=Q$, then we get a new system (12) by the transformation.

$\{\begin{array}{l}\frac{\text{d}X}{\text{d}t}=-\left(\rho +\mu \right)X\left(t\right)-\beta S_{0}I\left(t\right)-\beta X\left(t\right)I\left(t\right),\hfill \\ \frac{\text{d}Y}{\text{d}t}=\rho X\left(t\right)-\sigma \beta V_{0}I\left(t\right)-\mu Y\left(t\right)-\alpha Y\left(t\right)-\sigma \beta Y\left(t\right)I\left(t\right),\hfill \\ \frac{\text{d}I}{\text{d}t}=\beta \left(S_0+\sigma V_0\right)I\left(t\right)-cI\left(t\right)+\beta X\left(t\right)I\left(t\right)+\sigma \beta Y\left(t\right)I\left(t\right),\hfill \\ \frac{\text{d}Q}{\text{d}t}=\kappa I\left(t\right)-\left(\mu +\gamma \right)Q\left(t\right).\hfill \end{array}$ (12)

The feasible region of system (12) is

${\Omega }_1=\left\{X>-S_0,Y>-V_0,I>0,Q>0,X+Y+I+Q\le \frac{\Lambda }{\mu }-S_0-V_0\right\}.$

We define a Lyapunov function $L_1$ as

$L_1=X^2+\frac{\mu }{\rho }{Y}^2+2S_{0}I\mathrm{.}$

The derivative of the Lyapunov function $L_1$ along the trajectories of system (12) is

$\begin{array}{c}{\frac{\text{d}{L}_1}{\text{d}t}|}_{(12)}=-2\left(\rho +\mu \right)X^2-2\beta S_{0}XI-2\beta X^{2}I+2S_0\left[\beta \left(S_0+\sigma V_0\right)-c\right]I\\ +2\beta S_{0}XI+2\sigma \beta S_{0}YI+2\mu XY-2\sigma \beta S_{0}YI-\frac{2{\mu }^2}{\rho }{Y}^2\\ -\frac{2\mu \alpha }{\rho }{Y}^2-2\frac{\mu }{\rho }\sigma \beta Y^{2}I\\ =-2\left[\left(\rho +\mu \right)X^2-\mu XY+\frac{{\mu }^2}{\rho }{Y}^2\right]-2S_0\left[c-\beta \left(S_0+\sigma V_0\right)\right]I\\ -2\beta X^{2}I-2\frac{\mu }{\rho }\sigma \beta Y^{2}I-\frac{2\mu \alpha }{\rho }{Y}^2\end{array}$ (13)

If $R_0<1$, then $c-\beta \left(S_0+\sigma V_0\right)>0$. Because the arithmetic mean is larger than or equals to the geometric mean, we have

$\left(\rho +\mu \right)X^2+\frac{{\mu }^2}{\rho }{Y}^2\ge 2\sqrt{\frac{\rho +\mu }{\rho }}\mu XY\ge \mu XY$. For $\forall X\mathrm{,}Y\mathrm{,}I\mathrm{,}Q\in {\Omega }_1$, there is ${\frac{\text{d}{L}_1}{\text{d}t}|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>12\; <mo\; stretchy="false">\; )\; </mo>}}<0$, so ${\frac{\text{d}{L}_1}{\text{d}t}|}_{\text{Eq}\mathrm{.\; <mo\; stretchy="false">\; (\; </mo>12\; <mo\; stretchy="false">\; )\; </mo>}}$ is negative in ${\Omega }_1$. Hence, the disease-free equilibrium $E^0$ of system (5) is globally asymptotically stable. The proof is complete. □

Theorem 4.2. If $R_0>1$, $\tau =0$, then the endemic equilibrium $E^{*}$ is globally asymptotically stable.

Proof. At first, we prove that the endemic equilibrium $E^{*}$ is locally asymptotically stable.

If $R_0>1$, we obtain that the Jacobian matrix of the system (5) at the endemic equilibrium point $E^{*}$ is

$\begin{array}{c}J\left(E^{*}\right)=\left(\begin{array}{cccc}-\mu -\rho {\text{e}}^{-\lambda \tau }-\beta I^{\mathrm{<mo>\; *\; </mo>}}& 0& -\beta S^{\mathrm{<mo>\; *\; </mo>}}& 0\\ \rho {\text{e}}^{-\lambda \tau }& -\mu -\alpha -\sigma \beta I^{\mathrm{<mo>\; *\; </mo>}}& -\sigma \beta V^{\mathrm{<mo>\; *\; </mo>}}& 0\\ \beta I^{\mathrm{<mo>\; *\; </mo>}}& \sigma \beta I^{\mathrm{<mo>\; *\; </mo>}}& 0& 0\\ 0& 0& \kappa & -\mu -\gamma \end{array}\right)\\ =\left(\begin{array}{cccc}-\frac{\Lambda }{S^{\mathrm{<mo>\; *\; </mo>}}}+\rho \left(1-{\text{e}}^{-\lambda \tau }\right)& 0& -\beta S^{\mathrm{<mo>\; *\; </mo>}}& 0\\ \rho {\text{e}}^{-\lambda \tau }& \frac{-\rho S^{\mathrm{<mo>\; *\; </mo>}}}{V^{\mathrm{<mo>\; *\; </mo>}}}& -\sigma \beta V^{\mathrm{<mo>\; *\; </mo>}}& 0\\ \beta I^{\mathrm{<mo>\; *\; </mo>}}& \sigma \beta I^{\mathrm{<mo>\; *\; </mo>}}& 0& 0\\ 0& 0& \kappa & -\mu -\gamma \end{array}\right)\end{array}$ (14)

When $\tau =0$, we get the characteristic equation of the matrix $J\left(E^{*}\right)$ is

$\begin{array}{l}\left({\lambda }^3+a_1{\lambda }^2+a_2\lambda +a_3\right)\left(\lambda +\mu +\gamma \right)=\mathrm{0,}\hfill \end{array}$ (15)

where

$\begin{array}{l}{a}_1=\frac{\Lambda }{S^{*}}+\frac{\rho S^{*}}{V^{*}}>0\\ a_2={\sigma }^2{\beta }^2{V}^{*}{I}^{*}+{\beta }^2{S}^{*}{I}^{*}+\frac{\rho \Lambda }{V^{*}}>0\\ a_3=\sigma {\beta }^2\rho S^{*}{I}^{*}+\frac{\rho \beta S^{{*}^2}{I}^{*}}{V^{*}}+\frac{{\sigma }^2{\beta }^2\Lambda V^{*}{I}^{*}}{S^{*}}>0\end{array}$ (16)

Obviously, the characteristic Equation (15) always has a negative real root $\lambda =-\mu -\gamma$. And the remaining roots are determined by the following equation,

$\begin{array}{l}{\lambda }^3+a_1{\lambda }^2+a_2\lambda +a_3=0\hfill \end{array}$ (17)

By Routh-Hurwitz Criteria, the coefficients of the Equation (17) satisfy

$\begin{array}{l}{H}_1=\left|\begin{array}{c}{a}_1\end{array}\right|>0\hfill \end{array}$ (18)

$\begin{array}{c}{H}_2=\left|\begin{array}{cc}{a}_1& a_3\\ 1& a_2\end{array}\right|=a_1{a}_2-a_3\\ =\frac{\rho {\Lambda }^2}{S^{*}{V}^{*}}+{\beta }^2{S}^{*}{I}^{*}\left(\mu +\beta I^{*}\right)+\frac{{\rho }^2\mu S^{*}}{V^{*}}\\ +\rho {\beta }^2{\left(1-\sigma \right)}^2{S}^{*}{I}^{*}+\rho \sigma {\beta }^2{S}^{*}{I}^{*}\\ >0\end{array}$ (19)

$\begin{array}{l}{H}_3=\left|\begin{array}{ccc}{a}_1& a_3& 0\\ 1& a_2& 0\\ 0& a_1& a_3\end{array}\right|=a_3{H}_2>0\hfill \end{array}$ (20)

therefore, the endemic equilibrium $E^{*}$ is locally asymptotically stable.

Below, we prove that the endemic equilibrium $E^{*}$ is globally asymptotically stable. When $\tau \mathrm{<mo>\; =\; </mo>0}$, delay differential equation system (5) becomes the ordinary differential equation system as follows,

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\Lambda -\beta S\left(t\right)I\left(t\right)-\mu S\left(t\right)-\rho S\left(t\right)\hfill \\ \frac{\text{d}V}{\text{d}t}=\rho S\left(t\right)-\sigma \beta V\left(t\right)I\left(t\right)-\mu V\left(t\right)-\alpha V\left(t\right)\hfill \\ \frac{\text{d}I}{\text{d}t}=\beta \left(S\left(t\right)+\sigma V\left(t\right)\right)I\left(t\right)-cI\left(t\right)\hfill \\ \frac{\text{d}Q}{\text{d}t}=\kappa I\left(t\right)-\left(\mu +\gamma \right)Q\left(t\right).\hfill \end{array}$ (21)

We define a Lyapunov function $L_2$ as

$L_2=S-S^{*}-S^{*}\ln \frac{S}{S^{*}}+V-V^{*}-V^{*}\ln \frac{V}{V^{*}}+I-I^{*}-I^{*}\ln \frac{I}{I^{*}}.$

The derivative of the Lyapunov function $L_2$ along the trajectories of system (21) is

$\begin{array}{c}{\frac{\text{d}{L}_2}{\text{d}t}|}_{(21)}=\Lambda -\mu S-\mu V-\alpha V-cI-\frac{\Lambda S^{*}}{S}+\left(\rho +\mu \right)S^{*}+\beta S^{*}I\\ -\beta I^{*}\left(S+\sigma V\right)+c{I}^{*}-\rho S\frac{V^{*}}{V}+\mu V^{*}+\alpha V^{*}+\sigma \beta V^{*}I\\ =2\Lambda +\rho S^{*}-\mu S-\rho S^{*}\frac{V}{V^{*}}-\left(\rho +\mu \right)\frac{S^{{*}^2}}{S}-\beta I^{*}\frac{S^{{*}^2}}{S}-\rho S\frac{V^{*}}{V}-\beta S{I}^{*}\\ =2\left[\left(\rho +\mu \right)S^{*}+\beta S^{*}{I}^{*}\right]+\rho S^{*}-\mu S-\rho S^{*}\frac{V}{V^{*}}-\beta I^{*}\frac{S^{{*}^2}}{S}\\ -\left(\rho +\mu \right)\frac{S^{{*}^2}}{S}-\rho S\frac{V^{*}}{V}-\beta S{I}^{*}\\ =-\left(\mu S^{*}+\beta S^{*}{I}^{*}\right)\left(\frac{S^{*}}{S}+\frac{S}{S^{*}}-2\right)-\rho S^{*}\left(\frac{S^{*}}{S}+\frac{V}{V^{*}}+\frac{S{V}^{*}}{S^{*}V}-3\right).\end{array}$ (22)

Because the arithmetic mean is larger than or equals to the geometric mean, we have

$\frac{S^{\mathrm{<mo>\; *\; </mo>}}}{S}+\frac{S}{S^{\mathrm{<mo>\; *\; </mo>}}}-2\ge \mathrm{0,}$

$\frac{S^{\mathrm{<mo>\; *\; </mo>}}}{S}+\frac{V}{V^{\mathrm{<mo>\; *\; </mo>}}}+\frac{S{V}^{\mathrm{<mo>\; *\; </mo>}}}{S^{\mathrm{<mo>\; *\; </mo>}}V}-3\ge 0.$

Hence for $\forall S\mathrm{,}V\mathrm{,}I\mathrm{,}Q\in {\Omega }_1$, there is ${\frac{\text{d}{L}_2}{\text{d}t}|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>21\; <mo\; stretchy="false">\; )\; </mo>}}\le 0$. So the endemic equilibrium

$E^{\mathrm{<mo>\; *\; </mo>}}$ of system (5) is stable.

In addition, if $S^{\mathrm{<mo>\; *\; </mo>}}\ne S\mathrm{,}{V}^{\mathrm{<mo>\; *\; </mo>}}\ne V$, then, ${\frac{\text{d}{L}_2}{\text{d}t}|}_{(21)}<0$, hence, the endemic

equilibrium $E^{*}$ is globally asymptotically stable. The proof is complete. □

4.2. Local Stability of the Equilibria When $\tau \ne 0$

Theorem 4.3. For any time delay $\tau \ne 0$, the disease-free equilibrium $E^0$ is locally asymptotically stable if $R_0<1$, and unstable if $R_0>1$.

Proof. We obtain that the Jacobian matrix (23) of the system (5) at the disease-free equilibrium $E^0$ is

$J\left(E_0\right)=\left(\begin{array}{cccc}-\mu -\rho {\text{e}}^{-\lambda \tau }& 0& -\beta S_0& 0\\ \rho {\text{e}}^{-\lambda \tau }& -\mu -\alpha & -\sigma \beta V_0& 0\\ 0& 0& \beta \left(S_0+\sigma V_0\right)-c& 0\\ 0& 0& \kappa & -\mu -\gamma \end{array}\right)$ (23)

Hence, we get the characteristic equation for the above matrix $J\left(E_0\right)$ is

$\left(\lambda +\mu +\rho {\text{e}}^{-\lambda \tau }\right)\left(\lambda -\beta \left(S_0+\sigma V_0\right)+c\right)\left(\lambda +\mu +\alpha \right)\left(\lambda +\mu +\gamma \right)=0.$ (24)

If $R_0<1$, then The characteristic Equation (24) always has three negative real roots

${\lambda }_1=\beta \left(S_0+\sigma V_0-c\right),{\lambda }_2=-\mu -\alpha ,{\lambda }_3=-\mu -\gamma .$

The other characteristic roots of (24) is determined by the following equation

$\lambda +\mu +\rho {\text{e}}^{-\lambda \tau }=0.$ (25)

We suppose that a root of the Equation (24) is imaginary number ${\lambda }_0$. Through calculation, we get that its conjugate ${\bar{\lambda }}_0$ also satisfies the Equation (24). So the root of the Equation (24) must be a real number. Hence we obtain

$\lambda =-\mu -\rho {\text{e}}^{-\lambda \tau }\le -\mu .$

Therefore, the disease-free equilibrium $E^0$ is locally asymptotically stable.

If $R_0>1$, then $\lambda =\beta \left(S_0+\sigma V_0\right)-c>0$. Hence the characteristic Equation (24) has at least one positive root. So the disease-free equilibrium $E^0$ is unstable. The proof is complete. □

When $\tau \ne 0$, we get the characteristic equation about the matrix (14) is

$\left({\lambda }^3+b_1{\lambda }^2+b_2\lambda +b_3+\left(b_4{\lambda }^2+b_5\lambda +b_6\right){\text{e}}^{-\lambda \tau }\right)\left(\lambda +\mu +\gamma \right)=\mathrm{0,}$ (26)

where

$\begin{array}{l}{b}_1=\frac{\Lambda }{S^{*}}+\frac{\rho S^{*}}{V^{*}}-\rho ,\\ b_2={\sigma }^2{\beta }^2{V}^{*}{I}^{*}+{\beta }^2{S}^{*}{I}^{*}+\frac{\rho \Lambda }{V^{*}}-\frac{{\rho }^2{S}^{*}}{V^{*}},\\ b_3=\frac{\rho \beta S^{{*}^2}{I}^{*}}{V^{*}}+\frac{{\sigma }^2{\beta }^2\Lambda V^{*}{I}^{*}}{S^{*}}-\rho {\sigma }^2{\beta }^2{V}^{*}{I}^{*},\\ b_4=\rho ,\\ b_5=\frac{{\rho }^2{S}^{*}}{V^{*}},\\ b_6=\rho \sigma {\beta }^2{S}^{*}{I}^{*}+\rho {\sigma }^2{\beta }^2{V}^{*}{I}^{*}.\end{array}$ (27)

The characteristic Equation (26) always has a negative real root $\lambda =-\mu -\gamma$. And the remaining roots are determined by the following equation,

${\lambda }^3+b_1{\lambda }^2+b_2\lambda +b_3+\left(b_4{\lambda }^2+b_5\lambda +b_6\right){\text{e}}^{-\lambda \tau }=0.$ (28)

We assume that the pure imaginary number $\lambda =i\omega$ is the root of the Equation (28), and substitute $\lambda =i\omega$ into the equation,

$\begin{array}{l}-{\omega }^{3}i-\left(b_4\left(\cos \left(\lambda \tau \right)-\sin \left(\lambda \tau \right)i\right)+b_1\right){\omega }^2+\left(b_5\left(\cos \left(\lambda \tau \right)-\sin \left(\lambda \tau \right)i\right)+b_2\right)\omega i\\ +b_6\left(\cos \left(\lambda \tau \right)-\sin \left(\lambda \tau \right)i\right)+b_3=0\end{array}$ (29)

We separate the real and imaginary parts of the Equation (29) and obtain the following system,

$\{\begin{array}{l}-b_1{\omega }^2+b_3=\left(b_4{\omega }^2-b_6\right)\cos \left(\lambda \tau \right)-b_5\omega \sin \left(\lambda \tau \right),\\ -{\omega }^3+b_2\omega =-\left(b_4{\omega }^2-b_6\right)\sin \left(\lambda \tau \right)-b_5\omega \cos \left(\lambda \tau \right).\end{array}$ (30)

This leads to

${\left(-b_1{\omega }^2+b_3\right)}^2+{\left(-{\omega }^3+b_2\omega \right)}^2={\left(b_4{\omega }^2-b_6\right)}^2+b_5^2{\omega }^2\mathrm{.}$ (31)

We make $\phi ={\omega }^2$, and then simplify the Equation (31)

$G\left(\phi \right)={\phi }^3+g_1{\phi }^2+g_2\phi +g_3=0,$ (32)

where

$\begin{array}{l}{g}_1=b_1^2-b_4^2-2b_2,\\ g_2=b_2^2+2b_4{b}_6-2b_1{b}_3-b_5^2,\\ g_3=b_3^2-b_6^2.\end{array}$ (33)

Theorem 4.4. For any time delay $\tau \ne 0$, if $R_0>1$ and the coefficients in the $G\left(\phi \right)$ meet the Routh-Hurwitz Criteria, then the endemic equilibrium $E^{*}$ is locally asymptotically stable.

Proof. If the coefficients in the $G\left(\phi \right)$ meet the Routh-Hurwitz Criteria, then there is no positive real root in the Equation (32). Hence, there is no pure imaginary root in the Equation (28). It is known from Theorem 4.2. that when $\tau =0$, all roots of the Equation (28) are negative real parts. According to the continuity of roots to parameters, for any time delay $\tau$, all roots of the Equation (28) have strict negative real parts. Therefore the endemic equilibrium $E^{*}$ is locally asymptotically stable. The proof is complete. □

Theorem 4.5. For any time delay $\tau \ne 0$, if $R_0>1$, then the endemic equilibrium $E^{*}$ is locally asymptotically stable when $b_3-b_6<0$ and $\tau \in \left[\mathrm{0,}{\tau }_0\right)$ ; The system (5) undergoes a Hopf bifurcation at $E^{*}$ when $\tau ={\tau }_k$ ; the endemic equilibrium $E^{*}$ is unstable when $\tau >{\tau }_0$.

Proof. Because $b_3+b_6=\sigma {\beta }^2\rho S^{*}{I}^{*}+\frac{\rho \beta S^{{*}^2}{I}^{*}}{V^{*}}+\frac{{\sigma }^2{\beta }^2\Lambda V^{*}{I}^{*}}{S^{*}}>0$, $g_3<0$ when $b_3-b_6<0$. Hence $G\left(0\right)<0$, $G\left(+\infty \right)>0$. Because $G\left(\phi \right)$ is a continuous function, at least one positive real number $\phi$ satisfies $G\left(\phi \right)=0$. Therefore the Equation (28) has a pair of pure imaginary roots $\pm i{\omega }_0$.

We get ${\tau }_k>0$ from Equation (18), and the Equation (28) has a pair of pure imaginary roots $\pm i{\omega }_0$ when $\tau ={\tau }_k$, where

$\begin{array}{l}{\tau }_k=\frac{1}{{\omega }_0}\arccos \frac{\left(-b_1{\omega }_0^2+b_3\right)\left(b_4{\omega }_0^2-b_6\right)-b_5{\omega }_0\left(-{\omega }_0^3+b_2{\omega }_0\right)}{{\left(b_4{\omega }_0^2-b_6\right)}^2+{\left(b_5{\omega }_0\right)}^2}\\ +\frac{2k\pi }{{\omega }_0},\left(k=0,1,2,\cdots \right).\end{array}$ (34)

Taking the derivative of the Equation (28) with respect to $\tau$,

$\begin{array}{l}\left(3{\lambda }^2+2b_1\lambda +b_2\right)\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)+\left(2b_4\lambda +b_5\right){\text{e}}^{-\lambda \tau }\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)\\ +\left(b_4{\lambda }^2+b_5\lambda +b_6\right){\text{e}}^{-\lambda \tau }\left(-\lambda -\tau \left(\frac{\text{d}\lambda }{\text{d}\tau }\right)\right)=0,\end{array}$ (35)

we obtain

$\begin{array}{c}{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}=\frac{\left(3{\lambda }^2+2b_1\lambda +b_2\right)+\left(2b_4\lambda +b_5\right){\text{e}}^{-\lambda \tau }}{\lambda \left(b_4{\lambda }^2+b_5\lambda +b_6\right){\text{e}}^{-\lambda \tau }}-\frac{\tau }{\lambda }\\ =\frac{3{\lambda }^2+2b_1\lambda +b_2}{-\lambda \left({\lambda }^3+b_1{\lambda }^2+b_2\lambda +b_3\right)}+\frac{2b_4\lambda +b_5}{\lambda \left(b_4{\lambda }^2+b_5\lambda +b_6\right)}-\frac{\tau }{\lambda }.\end{array}$ (36)

Thus,

$\begin{array}{l}\text{sign}\left\{{\frac{\text{d}\left(Re\lambda \right)}{\text{d}\tau }|}_{\tau ={\tau }_k}\right\}=\text{sign}\left\{{Re{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}|}_{\lambda =i{\omega }_0}\right\}\\ =\text{sign}\left\{{Re\left[\frac{3{\lambda }^2+2b_1\lambda +b_2}{-\lambda \left({\lambda }^3+b_1{\lambda }^2+b_2\lambda +b_3\right)}+\frac{2b_4\lambda +b_5}{\lambda \left(b_4{\lambda }^2+b_5\lambda +b_6\right)}\right]|}_{\lambda \mathrm{<mo>\; =\; </mo>}i{\omega }_0}\right\}\\ =\text{sign}\left\{\frac{3{\omega }_0^6+2\left(b_1^2-2b_2\right){\omega }_0^4+\left(b_2^2-2b_1{b}_3\right){\omega }_0^2}{{\left(-b_1{\omega }_0^3+b_3\right)}^2+{\left(-{\omega }_0^4+b_2{\omega }_0^2\right)}^2}+\frac{-2b_4^2{\omega }_0^4+\left(2b_4{b}_6-b_5^2\right){\omega }_0^2}{{\left(b_4{\omega }_0^3-b_6{\omega }_0\right)}^2+b_5^2{\omega }_0^4}\right\}\end{array}$ (37)

From (31), we get ${\left(-b_1{\omega }_0^2+b_3\right)}^2+{\left(-{\omega }_0^3+b_2{\omega }_0\right)}^2={\left(b_4{\omega }_0^2-b_6\right)}^2+b_5^2{\omega }_0^2$, hence

$\text{sign}\left\{{\frac{\text{d}\left(Re\lambda \right)}{\text{d}\tau }|}_{\tau ={\tau }_k}\right\}=\text{sign}\left\{{Re{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}|}_{\lambda =i{\omega }_0}\right\}=\text{sign}\left\{\frac{G^{\prime }\left({\omega }_0^2\right)}{{\left(b_4{\omega }_0^2-b_6\right)}^2+b_5^2{\omega }_0^2}\right\}\mathrm{.}$ (38)

When $b_3-b_6<0$, $G^{\prime }\left({\omega }_0^2\right)>0$, so $\text{sign}\left\{\frac{G^{\prime }\left({\omega }_0^2\right)}{{\left(b_4{\omega }_0^2-b_6\right)}^2+b_5^2{\omega }_0^2}\right\}=1$.

Hence ${\frac{\text{d}\left(Re\lambda \right)}{\text{d}\tau }|}_{\tau ={\tau }_k}>0$.

When $\tau ={\tau }_0$, the Equation (28) has a pair of pure imaginary roots $\pm i{\omega }_0$, and the other roots are not integral multiples of $i{\omega }_0$. And ${\frac{\text{d}\left(Re\lambda \right)}{\text{d}\tau }|}_{\tau =0}\ne 0$. Therefore System (5) undergoes a Hopf bifurcation at $E^{*}$.

The Equation (28) does not have a pair of pure virtual roots when $\tau \in \left[\mathrm{0,}{\tau }_0\right)$. According to the theorem 4.2 and the theorem 4.4, we obtain that the endemic equilibrium $E^{*}$ is locally asymptotically stable.

Because ${\frac{\text{d}\left(Re\lambda \right)}{\text{d}\tau }|}_{\tau ={\tau }_k}>0$, the Equation (28) has positive real number roots when $\tau >{\tau }_0$. Therefore, the endemic equilibrium $E^{*}$ is unstable. The proof is complete. □

5. Numerical Simulation

In this section, we give some examples to imitate the theoretical results and to verify the correctness of the theoretical results.

Example 1. We choose a set parameters $\Lambda =1$, $\beta =0.2$, $\mu =0.05$, $\rho =0.8$, $\sigma =0.1$, $\nu =0.45$, $\kappa =0.3$, $\gamma =0.75$, $\alpha =0.28$, $\delta =0.01$, from which one can obtain $R_0=0.3609<1$. Therefore the disease-free equilibrium $E^0$ of system (5) is globally asymptotically stable (The figures of the system (5) are like as Figure 1).

Figure 1. For $\tau =0$, the disease-free equilibrium $E^0$ of system (5) is globally asymptotically stable.

Figure 2. For $\tau =0$, The endemic equilibrium $E^{*}$ of system (5) is globally asymptotically stable.

Figure 3. For $\tau =2.30$, The endemic equilibrium $E^{*}$ of system (5) is locally asymptotically stable.

Figure 4. For $\tau =2.33674$, (a) The endemic equilibrium $E^{*}$ is periodic solution. (b) The system (5) undergoes a Hopf bifurcation at $E^{*}$.

Example 2. We choose another set of parameters $\Lambda =1$, $\beta =0.4$, $\mu =0.05$, $\rho =0.8$, $\sigma =0.4$, $\nu =0.25$, $\kappa =0.3$, $\gamma =0.35$, $\alpha =0.28$, $\delta =0.01$, from which one can obtain $R_0=1.5195>1$, ${\tau }_0=2.33674$ and $b_3-b_6=-0.0003291<0$.

The endemic equilibrium $E^{*}$ of System (5) is globally asymptotically stable when $\tau =0$ (The figures of the System (5) are like as Figure 2).

The endemic equilibrium $E^{*}$ is locally asymptotically stable when $b_3-b_6<0$ and $\tau =2.30\in \left[0,{\tau }_0\right)$ (The figures of the System (5) are like as Figure 3).

The system (5) undergoes a Hopf bifurcation at $E^{*}$ when $\tau =2.52\approx {\tau }_0$ (The figures of the system (5) are like as Figure 4).

6. Conclusion

In this paper, we construct a time-delayed COVID-19 system with vaccination to study the effect of time-delay $\tau$ on the stability of the equilibrium points. We get the basic reproduction number, $R_0$, by the next generation matrix, and the system has a unique endemic equilibrium when $R_0>1$. Due to $\frac{\partial R_0}{\partial \rho }=\frac{\Lambda \beta \mu \left(\sigma -1\right)}{c\left(\mu +\alpha \right){\left(\rho +\mu \right)}^2}<0$, so we can reduce the basic reproduction number by increasing the vaccination rate. That is, if the vaccination rate is increased and COVID-19 is actively and effectively prevented, it will help to control the transmission of COVID-19.

For delay $\tau =0$, using the Lyapunov function approach, we obtained that the disease-free equilibrium $E^0$ and the endemic equilibrium $E^{*}$ are globally asymptotically stable when $R_0<1$ and $R_0>1$, respectively. For any time delay $\tau \ne 0$, using the Routh-Hurwitz Criteria, we obtained that the disease-free equilibrium $E^0$ is locally asymptotically stable when $R_0<1$. We also find the critical value ${\tau }_0$ at the endemic equilibrium. The endemic equilibrium $E^{*}$ is locally asymptotically stable when $\tau \in \left[\mathrm{0,}{\tau }_0\right)$. The Hopf bifurcation appears in the system, and the endemic equilibrium forms periodic solutions when $\tau ={\tau }_0$. From a biological point of view, the time delay affects the spread of infectious diseases. Therefore, in the epidemic prevention and control measures, we should strengthen the development of vaccines to reduce the time to generate antibodies after vaccination, so that we can better control the spread of the epidemic.

Our analysis is still not complete. We have not been able to show the global stability of the infection-free equilibrium and the endemic equilibrium for any delay $\tau \ne 0$, which is planned in our future research.

Acknowledgements

The authors sincerely thank to the editors and the anonymous reviewers for their valuable comments and suggestions, which improved the quality of the paper. At the same time, we would like to thank the National Natural Science Foundation of China for its support to the paper.

Funding

The research is supported by the National Natural Science Foundation of China (Grant number 11961075).

Ethics Approval

Not unethical.

Consent for Publication

Consent for publication.

Availability of Data and Material

Data transparency.

Code Availability

Software application or custom code.

Authors’ Contributions

The main idea of this paper was proposed by the authors and the authors completed the final manuscript together.

Conflicts of Interest

The authors declare they have no financial interests.

References

[1]	Wuhan Municipal Health Commission (2020) Wuhan Municipal Commission of Health and Health on Pneumonia of New Coronavirus Infection. http://www.wuhan.gov.cn/gsgg/202012/t20201229_1572409.shtml
[2]	Health Commission of Hubei Province (2020) Epidemic Situation of New Coronavirus Infection in Hubei Province. https://www.hubei.gov.cn/
[3]	National Health Commission of the P.R. China (2020) Update on Pneumonia Outbreak of New Coronavirus Infection. http://www.nhc.gov.cn/xcs/yqtb/list_gzbd.shtml
[4]	World Health Organization (WHO) (2020) Novel Coronavirus (2019-nCoV). Situation Report 1, 21 January. https://www.who.int/docs/default-source/coronaviruse/situation-reports/20200121-sitrep-1-2019-ncov.pdf?sfvrsn=20a99c10_4
[5]	Worldometer. https://www.worldometers.info/coronavirus
[6]	Roser, M., Ritchie, H. and Ortiz-Ospina, E. (2020) Coronavirus Disease (COVID-19)-Research and Statistics. https://ourworldindata.org/coronavirus
[7]	Lopez, A.D., Mathers, C.D., Ezzati, M., Jamison, D.T., Murray, C.J., Lamptey, P.R., Johnson, J.L., Khan, M., Nayeri, F. and Amini, E. (2006) Infectious Diseases. Changes in Individual Behavior Could Limit the Spread of Infectious Diseases. Population Bulletin, 61.
[8]	Mameli, C. and Zuccotti, G.V. (2013) The Impact of Viral Infections in Children with Community-Acquired Pneumonia. Current Infectious Disease Reports, 15, 197-202. https://doi.org/10.1007/s11908-013-0339-z
[9]	Diekmann, O., Heesterbeek, H. and Britton, T. (2013) Mathematical Tools for Understanding Infectious Disease Dynamics. Princeton Series in Theoretical and Computational Biology. Princeton University Press, Princeton.
[10]	Hethcote, H.W. (2000) The Mathematics of Infectious Diseases. SIAM Review, 42, 599-653. https://doi.org/10.1137/S0036144500371907
[11]	Keeling, M.J. and Rohani, P. (2007) Modeling Infectious Diseases in Humans and Animals. Princeton University Press, Princeton. https://doi.org/10.1515/9781400841035
[12]	Donnelly, C.A., Ghani, A., Leung, G., et al. (2003) Epidemiological Determinants of Spread of Causal Agent of Severe Acute Respiratory Syndrome in Hong Kong. The Lancet, 361, 1761-1766. https://doi.org/10.1016/S0140-6736(03)13410-1
[13]	https://www.yicai.com/news/101285104.html
[14]	Lipsitch, M., Cohen, T., Cooper, B., Robins, J., et al. (2003) Transmission Dynamics and Control of Severe Acute Respiratory Syndrome. Science (New York, N.Y.), 300, 1966-1970. https://doi.org/10.1126/science.1086616
[15]	Sun, C., Wei, Y., Arino, J. and Khan, K. (2011) Effect of Media-Induced Social Distancing on Disease Transmission in a Two Patch Setting. Mathematical Biosciences, 230, 87-95. https://doi.org/10.1016/j.mbs.2011.01.005
[16]	Pereda, A., Vielva, L.A., Vegas, A. and Prieto, A. (2001) Analyzing the Stability of the FDTD Technique by Combining the von Neumann Method with the Routh-Hurwitz Criterion. IEEE Transactions on Microwave Theory Techniques, 49, 377-381. https://doi.org/10.1109/22.903100
[17]	Kim, A.V. (1999) The Lyapunov Function Method. Springer, Berlin. https://doi.org/10.1007/978-94-017-1630-7_8
[18]	Hassard, D.D., Kazarinoff, N.D. and Wan, Y.-H. (2006) Theory and Applications of Hopf Bifurcation. SIAM Review, 24, 498-499. https://doi.org/10.1137/1024123
[19]	Jiang, Z.C., Ma, W.B. and Wei, J.J. (2016) Global Hopf Bifurcation and Permanence of a Delayed SEIRS Epidemic Model. Mathematics Computers in Simulation, 122, 35-54. https://doi.org/10.1016/j.matcom.2015.11.002
[20]	Kermack, W.O. and McKendrick, A.G. (1991) Contributions to the Mathematical Theory of Epidemics. II. The Problem of Endemicity. Bulletin of Mathematical Biology, 53, 57-87. https://doi.org/10.1016/S0092-8240(05)80041-2
[21]	Grossman, Z. (1980) Oscillatory Phenomena in a Model of Infectious Diseases. Theoretical Population Biology, 18, 204-243. https://doi.org/10.1016/0040-5809(80)90050-7
[22]	Rihan, F.A. and Anwar, M.N. (2012) Qualitative Analysis of Delayed SIR Epidemic Model with a Saturated Incidence Rate. International Journal of Differential Equations, 2012, Article ID: 408637. https://doi.org/10.1155/2012/408637
[23]	Elazzouzi, A., Alaoui, A.L., Tilioua, M. and Tridane, A. (2019) Global Stability Analysis for a Generalized Delayed SIR Model with Vaccination and Treatment. Advances in Difference Equations, 2019, Article No. 532. https://doi.org/10.1186/s13662-019-2447-z
[24]	Zhang, X.H., Jiang, D.Q., Hayat, T. and Ahmad, B. (2017) Dynamical Behavior of a Stochastic SVIR Epidemic Model with Vaccination. Physica A: Statistical Mechanics and Its Applications, 483, 94-108. https://doi.org/10.1016/j.physa.2017.04.173
[25]	Li, X.Z., Wang, J. and Ghosh, M. (2010) Stability and Bifurcation of an SIVS Epidemic Model with Treatment and Age of Vaccination. Applied Mathematical Modelling, 34, 437-450. https://doi.org/10.1016/j.apm.2009.06.002
[26]	Kumar, A., Srivastava, P.K. and Gupta, R.P. (2019) Nonlinear Dynamics of Infectious Diseases via Information-Induced Vaccination and Saturated Treatment. Mathematics and Computers in Simulation, 157, 77-99. https://doi.org/10.1016/j.matcom.2018.09.024
[27]	Buri, N., Mudrinic, M. and Vasovi, N. (2001) Time Delay in a Basic Model of the Immune Response. Chaos Solitons Fractals, 12, 483-489. https://doi.org/10.1016/S0960-0779(99)00205-2
[28]	Wang, K., Wang, W., Pang, H. and Liu, X. (2007) Complex Dynamic Behavior in a Viral Model with Delayed Immune Response. Physica D: Nonlinear Phenomena, 226, 197-208. https://doi.org/10.1016/j.physd.2006.12.001
[29]	Liu, X., Takeuchi, Y. and Iwami, S. (2008) SVIR Epidemic Models with Vaccination Strategies. Journal of Theoretical Biology, 253, 1-11. https://doi.org/10.1016/j.jtbi.2007.10.014
[30]	Sda, B. and Sba, C. (2021) Global Dynamics of SVIR Epidemic Model with Distributed Delay and Imperfect Vaccine. Results in Physics, 25, Article ID: 104245.
[31]	Diekmann, O., Heesterbeek, J. and Metz, J. (1990) On the Definition and the Computation of the Basic Reproduction Ratio R0 in Models for Infectious Diseases in Heterogeneous Populations. Journal of Mathematical Biology, 28, 365-382. https://doi.org/10.1007/BF00178324
[32]	Driessche, P. and Watmough, J. (2002) Reproduction Numbers and Sub-Threshold Endemic Equilibria for Compartmental Models of Disease Transmission. Mathematical Biosciences, 180, 29-48. https://doi.org/10.1016/S0025-5564(02)00108-6
[33]	Kuniya, T. (2013) Global Stability of a Multi-Group SVIR Epidemic Model. Nonlinear Analysis: Real World Applications, 14, 1135-1143. https://doi.org/10.1016/j.nonrwa.2012.09.004
[34]	Wiggins, S. (2003) Introduction to Applied Nonlinear Dynamical Systems and Chaos. 2nd Edition, Springer-Verlag, New York.