Stability and bifurcation analysis of $SIQR$ for the COVID-19 epidemic model with time delay

1.   Introduction

COVID-19 is a respiratory infectious disease and was first reported in Wuhan, China. It is caused by a novel coronavirus named SARS-CoV-2. SARS-COV-2 appeared after severe acute respiratory syndrome coronavirus (SARS) and Middle East respiratory syndrome coronavirus (MERS) ^[1]. Coronaviruses are enveloped nonsegmented positive-sense RNA viruses and are broadly distributed in mammals, including humans. They can destroy the human respiratory systems and cause severe acute respiratory syndrome ^[2]. Infected individuals usually cough, have a fever, and have difficulty breathing. Some of them even die of COVID-19. Moreover, there are many ways to spread COVID-19, such as direct transmission, aerosol transmission and contact transmission. This caused COVID-19 to spread rapidly worldwide. As a global infectious disease, COVID-19 has caused more than one million deaths. We attempt to establish a model to analyze the spread of COVID-19 and provide some advice to make COVID-19 controllable.

In recent investigations, many researchers have studied different epidemic models applied to the spread of COVID-19. In ref. ^[3] and ref. ^[4], Varotsos et al. established an $SPRD$ model that included people that were susceptible, infected, recovered and dead. The model was based on COVID-19 data and was established to predict COVID-19 spread. We think there is no need to include people who died as a main part of the model because the proportion of people who died is very small. In ref. ^[5], León et al. developed an $SEIARD$ mathematical model to investigate the COVID-19 outbreak in Mexico. The $SEIARD$ model included infection with symptoms, infection without symptoms (asymptomatic), recovery from symptomatic infection and recovery from asymptomatic infection. If we consider only the infectivity of COVID-19, both infection with symptoms and infection without symptoms (asymptomatic) can infect the suspected cases. There is no need to divide them into infected with symptoms and infected without symptoms (asymptomatic). Both the recovery from symptomatic infection and the recovery from asymptomatic infection have very small possibilities for reinfection. There is no need to divide them into recovered from symptomatic infection and recovered from asymptomatic infection. In ref. ^[6], an $SEIR$ epidemic model was developed for COVID-19. In ref. ^[7] and ref. ^[8], two $SIR$ epidemic models were developed for COVID-19. For COVID-19, many countries have taken quarantine measures, which play an important part in the spread of COVID-19. These authors did not consider the influence of these quarantine measures. The quarantine measures taken by many countries make the proportion of quarantined people large. These quarantine measures can prevent infected individuals from infecting suspected individuals. In addition, exposure to air pollution can damage the heart and lungs and increase vulnerability to more serious coronavirus effects ^[9]. These quarantine measures can reduce the influence of air pollution. We think it is necessary to consider quarantined people because of the importance of these quarantine measures and the number of quarantined people. We attempt to construct an $SIQR$ epidemic model including quarantined people to describe the spread of COVID-19.

In the spreading process of COVID-19, there is a time delay from infection to isolation. There are some studies about other epidemics with time delays. Liu et al. ^[10] developed an $SEIRU$ epidemic model with a time delay before an infected person can transmit the infection to another person. They evaluated the effect of the latency period on the dynamics of COVID-19. Xu ^[11] and Yang et al. ^[12] also developed epidemic models with time delays before an infected person can transmit the infection to another person. Although this kind of time delay exists, human beings have difficulty changing this kind of time delay. In addition, Wang et al. ^[13] and Liu et al. ^[14] made time delays in their models the sojourn times in an infective state. Lu et al. ^[15] presented an $SIQR$ model with a time delay from infection to recovery. The influence of the time delay from infection to recovery was discussed in detail. According to their conclusions, the smaller the time delay from infection to isolation was, the better COVID-19 was controlled. We decide to focus on the time delay from infection to isolation and discuss the effect of the time delay from infection to isolation on the spread of COVID-19, which is the main difference between our work and the works of others.

Stability and bifurcation have great significance to epidemic models, and some studies have analyzed stability and bifurcation for some epidemic models. In ref. ^[16], Greenhalgh investigated the stability of some $SEIRS$ epidemiological models with vaccination and temporary immunity. Greenhalgh proposed that there was a threshold parameter $R_0$, and the disease could persist if and only if $R_0$ exceeded one. In addition, disease-free equilibrium always existed and was locally stable if $R_0 < 1$ and unstable if $R_0 > 1$. However, for $R_0 > 1$, the endemic equilibrium was unique and locally asymptotically stable. Based on this conclusion, we think it is necessary for us to calculate the threshold parameter $R_0$ in our model. In ref. ^[17], Xie et al. investigated the global stability of endemic equilibrium of an $SIS$ epidemic model in complex networks. They proved that the endemic equilibrium was globally asymptotically stable by using a combination of the Lyapunov function method and a monotone iterative technique. However, the Hopf bifurcation analysis was not mentioned in this paper. In refs. ^[18,19], two $SIS$ models were proposed. In ref. ^[20], an $SIQRS$ model was proposed. The stability of the models was analyzed, but bifurcation analysis was not mentioned. In refs. ^{[21,22,23,24]}, the stability and bifurcation of some $SIR$ models were analyzed by different authors. In refs. ^{[25,26,27,28]}, different $SIR$ models were developed. The stability and bifurcation of them were analyzed. In ref. ^[29], an $SIR$ epidemic model with time delay was proposed. In ref. ^[30], an $SIR$ model with information-dependent vaccination was proposed. We compare the methods they used to analyze the stability and bifurcation of these $SIR$ models. Then, we use the multiple time scales method, which was also used in ref. ^[21]. The multiple time scales method is systematic. It is a standard method for calculating the normal form of Hopf bifurcation. It can be directly applied to the original nonlinear dynamical system, which is described by ordinary differential equations and delayed differential equations, without application of the center manifold theory ^[31].

It is obvious that the isolation measures taken by many countries play an important role in controlling the spread of COVID-19. The time delay also has a considerable influence on the spread. Thus, we attempt to develop one mathematical model with a time delay to investigate the spread of COVID-19 and analyze the stability of this model. In some studies such as refs. ^[3,4,5], the authors forecasted how the epidemic situation changed with detailed data. However, we decide to focus on the stability of equilibria in our model, calculate the critical time delay from infection to isolation and analyze the dynamical property of our model to discuss the tendency of the spread of COVID-19.

The main objective of this paper includes constructing the DDE model for COVID-19 and analyzing the stability of the model, the existence of Hopf bifurcation and the stability of the bifurcating periodic solution. The rest of the paper is organized as follows. In Section 2, we construct the DDE model according to the characteristics of the spread of COVID-19. In Section 3, we analyze the existence and stability of the equilibria and the existence of Hopf bifurcation associated with the COVID-19 model with a time delay. We calculate the critical time delay from infection to isolation and analyze the dynamic properties of our model. In Section 4, we derive the normal form of the Hopf bifurcation for the above model and analyze the stability of the bifurcating periodic solution. In Section 5, we present the results of simulations to verify the correctness of our analysis. Finally, the conclusion is drawn in Section 6.

2.   Mathematical modeling

Suspected individuals can be infected with COVID-19 by infected individuals. After individuals are infected with COVID-19, some of them will be quarantined, and some of them will die of COVID-19. The quarantined individuals will not infect others. For COVID-19, we ignore the effect of reinfection because the probability of reinfection is very small. In addition, we ignore the probability of transforming the suspected individuals into recovered individuals directly since the probability is also small. We divide people into four kinds (suspected, infected, quarantined and recovered) based on their infectivity. The specific conversion between the four kinds of people is given in Figure 1.

The variables and parameters in Figure 1 are given in Table 1.

For COVID-19, there is a latent period after the suspected individuals are infected. In the latent period, the infected individuals have no symptoms. However, these asymptomatic infected individuals can also infect suspected individuals. If we consider only their infectivity, there is little difference between asymptomatic infected individuals and symptomatic infected individuals. Thus, both of them can be divided into infected. Additionally, some suspected individuals may be quarantined, and they will not be infected. This means that the number of suspected persons who can be infected will decrease. Considering the influence of the suspected individuals who are quarantined, we can reduce the transmission rate from $S$ to $I$, and there is no need to include the suspected individuals who are quarantined in the model.

Based on Figure 1, we construct the DDE model (2.1):

where the definitions of parameters and variables are presented in Figure 1 and Table 1, and $\tau > 0$ denotes the time delay from infection to isolation.

The initial condition of system (2.1) is ${\phi _S}(\theta) \ge 0, {\phi _I}(\theta) \ge 0, {\phi _Q}(\theta) \ge 0, {\phi _R}(\theta) \ge 0$, where $\Phi = ({\phi _S}(\theta), {\phi _I}(\theta), {\phi _Q}(\theta), {\phi _R}(\theta)) \in C([- \tau, 0], \mathbb{R} _{ + 0}^4)$ for $\theta \in [- \tau, 0]$, and $C([- \tau, 0], \mathbb{R}_{ + 0}^4)$ is the Banach space of continuous functions mapping interval $[- \tau, 0]$ into $\mathbb{R}_{ + 0}^4$. For the system (2.1), $\mathbb{R}_{+0}^4 = \{ (S, I, Q, R)\left| {S \ge 0, I \ge 0, Q \ge 0, R \ge 0} \right.\}$.

According to the initial condition of system (2.1), we present a theorem addressing the nonnegativity and the boundedness of the solution to system(2.1).

Theorem 2.1. If $S(0) \ge 0, I(0) \ge 0, Q(0) \ge 0, R(0) \ge 0$, the solution $S^*(t)$, $I^*(t)$, $Q^*(t)$, $R^*(t)$ to system (2.1) with $\tau = 0$ is nonnegative and bounded when $t > 0$.

Proof. First, we prove $S^*(t)\ge 0$ when $t\ge 0$ under the initial condition of system (2.1).

We assume that $S^*(t)$ is not always nonnegative for $t\ge 0$ and make $t_1$ the first time that $S^*(t_1) = 0$, $S'(t_1) < 0$. According to the first equation of system (2.1), we can obtain $S'(t_1) = \Lambda > 0$. The two conclusions we obtain are contradictory.

Therefore, $S^*(t)\ge 0$ when $t > 0$. In the same way, $I^*(t)\ge 0$, $Q^*(t)\ge 0$, $R^*(t)\ge 0$ when $t > 0$. The solution to system (2.1) is positive when $t > 0$.

Then, let $N(t) = S(t)+I(t)+Q(t)+R(t)$, and $N(t)$ represent the total number of people at time $t$.

Add four equations of Eq. (2.1), and we obtain $N'(t) = \Lambda-\mu(S(t)+I(t)+Q(t)+R(t))$. Then, we can obtain $\mathop {\lim }\limits_{t \to \infty } \rm{sup}N(t) = \frac{\Lambda }{\mu }$. Thus, the solution $S^*(t)$, $I^*(t)$, $Q^*(t)$, $R^*(t)$ to system (2.1) is bounded when $t > 0$.

Remark 1: We prove if $S(0) \ge 0, I(0) \ge 0, Q(0) \ge 0, R(0) \ge 0$, the solution $S^*(t)$, $I^*(t)$, $Q^*(t)$, $R^*(t)$ to system (2.1) with $\tau = 0$ is positive. It is not easy for us to prove that the solution to system (2.1) is positive when $\tau > 0$. However, according to our numerical simulation, we can find that when system (2.1) is stable, the solution to system (2.1) is always positive, which is not contradictory to the positivity of the solution to system (2.1).

Next, we consider the dynamic phenomena of system (2.1).

3.   Analysis of stability

3.1.   Existence and stability of equilibria

In this section, the system (2.1) is considered, and we first determine the equilibria of the above system. Obviously, the system (2.1) has two equilibria:

with

Transferring the equilibria $E_k (k = 1, 2)$ to the origin, we make $w = S+S_k^*$, $x = I+I_k^*$, $y = Q+Q_k^*$, $z = R+R_k^*$. We can obtain $S = w-S_k^*$, $I = x-I_k^*$, $Q = y-Q_k^*$ and $R = z-R_k^*$. Substituting these equations into the model (2.1), we use $S$, $I$, $Q$ and $R$ to represent $w$, $x$, $y$ and $z$. We obtain the following model:

In this subsection, we first demonstrate the stability of the equilibrium $E_1 = (S_{1}^{*}, I_{1}^{*}, Q_{1}^{*}, R_{1}^{*}) = (\frac{\Lambda }{\mu }, 0, 0, 0).$ We obtain the characteristics equation of the linearized system of Eq. (3.3) at the equilibrium $E_1$ as follows:

When $\tau = 0$, Eq. (3.4) becomes:

We show the following assumption:

(H1) $\beta S_{1}^{*}-(\varepsilon +\mu +{{\alpha }_{1}}+{{\gamma }_{1}}) < 0.$

Under (H1), the four roots of Eq. (3.5) have negative real parts due to $\mu > 0, \alpha_1 > 0, \gamma_1 > 0$, and the equilibrium $E_1$ is locally asymptotically stable when $\tau = 0$.

Similarly, for the stability of the other equilibrium $E_2 = (S_{2}^{*}, I_{2}^{*}, Q_{2}^{*}, R_{2}^{*}),$ we obtain the characteristic equation of the linearized system of Eq. (3.3) at the equilibrium $E_2$ as follows:

When $\tau = 0$, Eq. (3.6) becomes:

If (H1) is not satisfied, the four roots of Eq. (3.7) have negative real parts due to $\mu > 0, \alpha_2 > 0, \gamma_2 > 0$, and the equilibrium $E_2$ of system (2.1) is locally asymptotically stable when $\tau = 0$.

We can find that $E_1$ is the disease-free equilibrium, and $E_2$ is the endemic equilibrium. We calculate the basic reproduction number $R_0$, which is the number of suspected individuals who are infected by the same infectious individual and can estimate the infectiousness of an infectious disease. According to the system (2.1), we can obtain the new infections matrix $\mathcal{F}$ and the transition matrix $\mathcal{V}$.

Then, we make $F_0$ represent the derivative of $\mathcal{F}$ at $E_1$ and $V_0$ represent the derivative of $\mathcal{V}$ at $E_1$:

The inverse of $V_0$ is:

And we can obtain:

The maximum eigenvalue of ${F_0}V_0^{ - 1}$ is $R_0$:

We compare $R_0$ with (H1), and we find that if $R_0 < 1$, the inequality of (H1) exists constantly, and if $R_0 > 1$, (H1) is not satisfied. This represents if $R_0 < 1$, the disease-free equilibrium $E_1$ is locally asymptotically stable when $\tau = 0$ and if $R_0 > 1$, the endemic equilibrium $E_2$ is locally asymptotically stable when $\tau = 0$.

3.2.   Existence of Hopf bifurcation

For equilibrium $E_1$, when $\tau > 0$, let $\lambda = \mathrm{i}\omega_0\; (\omega_0 > 0)$ be a root of the following equation:

Substituting $\lambda = i\omega_0$ into Eq. (3.8) and separating the real and imaginary parts, we have

Thus, we obtain:

Adding the square of two equations (3.9), we can obtain

Thus, $\omega_0 = \sqrt{{{\varepsilon }^{2}}-{{[\beta S_{1}^{*}-(\mu +{{\alpha }_{1}}+{{\gamma }_{1}})]}^{2}}}$ when ${{{\varepsilon }^{2}}-{{[\beta S_{1}^{*}-(\mu +{{\alpha }_{1}}+{{\gamma }_{1}})]}^{2}}} \ge 0$, and substituting it into Eq. (3.10). Since $\omega_0 > 0$ and $\varepsilon \ge 0$, we can obtain $Q_1 > 0$ and obtain:

where $Q_1, P_1$ are given in Eq. (3.10), $j = 0, 1, 2\cdots.$

Let $\lambda(\tau) = \alpha(\tau)+i\omega(\tau)$ be the root of Eq. (3.8) satisfying $\alpha(\tau_1^{(j)}) = 0, \omega(\tau_1^{(j)}) = \omega_0, j = 0, 1, 2\cdots$. Then, we have the transversality conditions

Thus, the system (2.1) undergoes a Hopf bifurcation at equilibrium $E_1$ of system (2.1) when $\tau = \tau_1^{(j)}, j = 0, 1, 2, \cdots$.

For the other equilibrium $E_2$, when $\tau > 0$, let $\lambda = \mathrm{i}\omega\; (\omega > 0)$ be the root of the following equation:

Substituting it into the equation and separating the real and imaginary parts, we have

Equation (3.14) can be presented as:

Adding the square of two equations Eq. (3.14), let $z = \omega^2$, and

where ${{c}_{1}} = (\beta I_{2}^{*}+\mu)^2-2{{\beta }^{2}}S_{2}^{*}I_{2}^{*}, {{c}_{0}} = {{({{\beta }^{2}}S_{2}^{*}I_{2}^{*})}^{2}}-2{{\beta }^{2}}S_{2}^{*}I_{2}^{*}(\beta I_{2}^{*}+\mu)\varepsilon.$

Therefore, we have the following assumption:

(H2) $c_1^2-4c_0 > 0, c_1 < 0, c_0 > 0$.

(H3) $c_0 < 0$.

If (H2) holds, then Eq. (3.16) has two positive roots, $z_1$ and $z_2$ ($z_1 < z_2$). If (H3) holds, then Eq. (3.16) has only one positive root $z_3$. Without loss of generality, substituting $\omega_k = \sqrt{z_k}\; (k = 1, 2, 3)$ into Eq. (3.15), we obtain:

where $Q_2, P_2$ are given in Eq. (3.15) and ${Q_{2, k}} = {\left. {{Q_2}} \right|_{\omega = {\omega _k}, \tau = \tau _{2, k}^{(j)}}}, {P_{2, k}} = {\left. {{P_2}} \right|_{\omega = {\omega _k}, \tau = \tau _{2, k}^{(j)}}}, k = 1, 2, 3.$

Note that since $Q_1 = \frac{\omega_0}{\varepsilon}$ is always nonnegative, we can obtain $\tau_1^{(j)}$ directly. However, whether $Q_2$ is positive depends on the parameters of Eq. (3.15), and we need to obtain $\tau_{2, k}^{(j)}$, $(j = 1, 2, 3)$ according to the value of $Q_2$.

Furthermore, let $\lambda(\tau) = \alpha(\tau)+\mathrm{i}\omega_k(\tau)$ be the root of Eq. (3.17) satisfying $\alpha(\tau_{2, k}^{(j)}) = 0$, $\omega(\tau_{2, k}^{(j)}) = \omega_k\; (k = 1, 2, 3;\; j = 0, 1, 2\ldots)$.

If (H2) or (H3) holds and $z_k = \omega_k^{2} (k = 1, 2, 3)$, then we can deduce ${{\left. \operatorname{Re}(\frac{\mathrm{d}\tau }{\rm{d}\lambda }) \right|}_{\tau = \tau _{2, k}^{(j)}}}$:

There is a Hopf bifurcation at equilibrium $E_2$ of system (2.1) when $\tau = \tau_{2, k}^{(j)}, j = 0, 1, 2\cdots$.

Then, we present a stability theorem of system (2.1)'s equilibria $E_1$, $E_2$ and the existence of Hopf bifurcation.

Theorem 3.1. We consider the model (2.1):

(1) If (H1) holds, the equilibrium $E_1$ of the model (2.2) undergoes a Hopf bifurcation when $\tau = \tau_1^{(j)}, j = 0, 1, 2\cdots$, and we obtain the following. When $\tau\in [0, \tau_1^{(0)})$, the equilibrium $E_1$ is locally asymptotically stable. When $\tau\in[\tau_1^{(0)}, +\infty)$, the equilibrium $E_1$ is locally asymptotically unstable.

(2) Under the condition that (H1) is not satisfied, and if (H2) or (H3) holds, the equilibrium $E_2$ of the model (2.1) undergoes a Hopf bifurcation when $\tau = \tau_{2, k}^{(j)}(k = 1, 2, 3)$, and we obtain the following.

(a) If (H2) holds, $h(z)$ has two positive roots, $z_1$ and $z_2$; then we assume $z_1 < z_2$, and we obtain ${h}'({{z}_{1}}) < 0, {h}'({{z}_{2}}) > 0$. Then, $\exists \; m\in N$, which can make $0 < \tau _{2, 2}^{(0)} < \tau _{2, 1}^{(0)} < \tau _{2, 2}^{(1)} < \tau _{2, 1}^{(1)} < \cdots < \tau_{2, 1}^{(m-1)} < \tau _{2, 2}^{(m)} < \tau _{2, 2}^{(m+1)}$. When $\tau \in [0, \tau _{2, 2}^{(0)})\cup \bigcup\limits_{l = 1}^{m}{(\tau _{2, 1}^{(l-1)}}, \tau _{2, 2}^{(l)})$, the equilibrium of the model is locally asymptotically stable. When $\tau \in \bigcup\limits_{l = 0}^{m-1}{(\tau _{2, 2}^{(l)}}, \tau _{2, 1}^{(l)})\cup(\tau_{2, 2}^{(m)}, +\infty)$, the equilibrium is locally asymptotically unstable.

(b) If (H3) holds, $h(z)$ has only one positive root $z_3$, when $\tau\in [0, \tau_{2, 3}^{(0)})$, the equilibrium $E_2$ is locally asymptotically stable. When $\tau\in(\tau_{2, 3}^{(0)}, +\infty)$, the equilibrium $E_2$ is unstable.

4.   Normal form of Hopf bifurcation

Without loss of generality, we denote the critical value $\tau = \tau^*$ when the characteristic Eq. (3.4) and Eq. (3.6) have eigenvalue $\lambda = i\omega_k^{*}, k = 1, 2$, at which system (2.1) undergoes a Hopf bifurcation at equilibrium $E_k = (S_{k}^{*}, I_{k}^{*}, Q_{k}^{*}, R_{k}^{*}), k = 1, 2$.

First, we transform the equilibrium to the origin; then, we use the multiple time scales method, and system (2.1) can be written as

where $X(t) = (\hat{S_k}(t), \hat{I_k}(t), \hat{Q_k}(t), \hat{R_k}(t))^{\mathrm{T}}$, $X(t-\tau) = (\hat{S_k}(t-\tau), \hat{I_k}(t-\tau), \hat{Q_k}(t-\tau), \hat{R_k}(t-\tau))^\mathrm{T}$,

We make $\widetilde {\dot X}(t) = AX(t) + BX(t - \tau)$ the linear system of system (4.1). We assume that $h_k(k = 1, 2)$ is the eigenvector of the eigenvalue $\lambda \rm{ = }i\omega_k^{*} (k = 1, 2)$ of the linear system and that $h_k^{*}$ is the eigenvector of the eigenvalue $\lambda \rm{ = }-i\omega_k^{*}$ of the linear system. We make $\widetilde {\rm I}$ the unit matrix, and $h_k$ satisfies $(\lambda \widetilde {\rm I} - A - B{\rm{e}^{ - \lambda \tau }}){h_k} = 0$, where $\lambda = i\omega_k^*$, and $h_k^*$ satisfies ${(\lambda \widetilde {\rm I} - A - B{\rm{e}^{ - \lambda \tau }})^\rm{T}}{h_k}^* = 0$, where $\lambda = -i\omega_k^*$. Additionally, $\langle {{h}_k^{*}}, h_k \rangle = {\overline {h_k^*} ^\rm{T}}h_k = 1$, and we can obtain

where ${{d}_{k}} = {{(\frac{-{{\beta }^{2}}S_{k}^{*}I_{k}^{*}}{{{(\mu +\beta I_{k}^{*})}^{2}}+{{\omega_k^* }^{2}}}+1)}^{-1}}, k = 1, 2$.

We use the multiple time scale method to deduce the normal form of the Hopf bifurcation and assume the solution to Eq. (4.1) is

where $X({{T}_{0}}, {{T}_{1}}, {{T}_{2}}, \cdots) = {{[S({{T}_{0}}, {{T}_{1}}, {{T}_{2}}, \cdots), I({{T}_{0}}, {{T}_{1}}, {{T}_{2}}, \cdots), Q({{T}_{0}}, {{T}_{1}}, {{T}_{2}}, \cdots), R({{T}_{0}}, {{T}_{1}}, {{T}_{2}}, \cdots)]}^\mathrm{T}}$, ${{X}_{k}}({{T}_{0}}, {{T}_{1}}, {{T}_{2}}, \cdots) = \; {{[{{S}_{k}}({{T}_{0}}, {{T}_{1}}, {{T}_{2}}, \cdots), {{I}_{k}}({{T}_{0}}, {{T}_{1}}, {{T}_{2}}, \cdots), {{Q}_{k}}({{T}_{0}}, {{T}_{1}}, {{T}_{2}}, \cdots), {{R}_{k}}({{T}_{0}}, {{T}_{1}}, {{T}_{2}}, \cdots)]}^\mathrm{T}}$, $k = 1, 2$, $T_i = \epsilon^i t, i = 0, 1, 2, \cdots$ and $T_i$ is the scaling transform in the time direction.

The derivative with respect to $t$ is

where ${{D}_{i}} = \frac{\partial }{\partial {{T}_{i}}}, i = 0, 1, 2, \cdots.$

Note that ${{X}_{j}} = {{({{S}_{j}}, {{I}_{j}}, {{Q}_{j}}, {{R}_{j}})}^\mathrm{T}} = {{X}_{j}}({{T}_{0}}, {{T}_{1}}, {{T}_{2}}, \cdots), {{X}_{j, {{\tau }_{c}}}} = {{({{S}_{j, {{\tau }_{c}}}}, {{I}_{j, {{\tau }_{c}}}}, {{Q}_{j, {{\tau }_{c}}}}, {{R}_{j, {{\tau }_{c}}}})}^\mathrm{T}} = {{X}_{j}}({{T}_{0}}-{{\tau }_{c}}, {{T}_{1}}, {{T}_{2}}, \cdots), j = 1, 2, \cdots.$

Then, we can obtain

We consider that $\tau$ is the bifurcation parameter, and we set $\tau = \tau_c+{{\epsilon}} \tau_{\varepsilon}$, where $\tau_{c} = \tau_k^{(j)}(k = 1, 2)$ is the critical value of the Hopf bifurcation, $\tau_{\varepsilon}$ is the perturbation parameter, and $\epsilon$ is the dimensionless parameter.

Using Taylor series expansion of $X(t-\tau)$, we obtain that

where $X_{j, \tau_c} = X_j(T_0-\tau_c, T_1, T_2, \cdots), j = 1, 2, \cdots.$

Then, substituting the solution with the multiple scales Eqs. (4.3)-(4.5) into Eq. (4.1) and balancing the coefficients of $\epsilon^{n} (n = 1, 2, 3)$, we obtain a set-ordered linear differential equation:

Since $\pm i{{\omega }_{k}^{*}}(k = 1, 2)$ are the eigenvalues of the characteristic Eq. (4.1), the solution to Eq. (4.6) can be expressed in the form of Eq. (4.7):

where $G_1$ and ${\bar G}_1$ represent the coordinates on the center manifold, and the linear part of the center manifold can be presented as equation (4.7).

Next, for the $\epsilon^2$-order terms, we can obtain the following equations:

We substitute Eq. (4.7) into the right expression of Eq. (4.8), and the coefficients before $\mathrm{e}^{i \omega_k^* T_0}$ are denoted by the vector $m_1$. According to the solvability condition $\langle {{h}_k^{*}}, m_1 \rangle = 0$, we can solve $\frac{\partial G}{\partial {{T}_{1}}}$ as follows:

where ${{M}_{k}} = \varepsilon {\mathrm{e}^{-i\omega_k^{*} {{\tau }_{c}}}}{{[\frac{\beta I_{k}^{*}(i\omega_k^{*} +\mu +\beta I_{k}^{*})}{{{(\mu +\beta I_{k}^{*})}^{2}}+{{\omega_k^{*} }^{2}}}\frac{\beta S_{k}^{*}}{i\omega_k^{*} +\mu +\beta I_{k}^{*}}+1-\varepsilon {{\tau }_{c}}{\mathrm{e}^{-i\omega_k^{*} {{\tau }_{c}}}}]}^{-1}}, k = 1, 2$.

We assume:

Substituting Eq. (4.10) into Eq. (4.8), we obtain:

Next, balancing the $\epsilon^{3}$-order terms, we can obtain equations as follows:

Substituting Eq. (4.7) and Eq. (4.10) into the right expression of Eq. (4.12), and the coefficients before $\rm{e}^{i \omega_k T_0}$ are denoted by the vector $m_2$. According to the solvability condition $\langle {{h}_1^{*}}, m_2 \rangle = 0$ and note that $\tau_\varepsilon^2$ is small enough for small unfolding parameter $\tau_\varepsilon$, we ignore the term $\tau_\varepsilon^{2}G_1$. Then, we have:

where

where $g_k$ are given in Eq. (4.11) and $h_k$ are given in Eq. (4.2).

Let $G_1 \mapsto (G_1/\varepsilon)$, and we can obtain the normal form of Hopf bifurcation of system (2.1) as:

where $M_k$ is given in Eq. (4.9), and $H_k, \chi_k$ are given in Eq. (4.13).

Let $G = r{\rm{e}^{i\theta_k^{*} }}(k = 1, 2)$ and substitute it into Eq. (4.14), and we can obtain the normal form of Hopf bifurcation in polar coordinates:

where $M_k$ is given in Eq. (4.9), and $H_k, \chi_k$ are given in Eq. (4.13) $(k = 1, 2)$.

According to the normal form of bifurcation in polar coordinates, there is a theorem as follows:

Theorem 4.1. For system (4.15), if $\frac{\operatorname{Re}(M_k){{\tau }_{\varepsilon }}}{\operatorname{Re}(H_k\chi_k)} < 0 (k = 1, 2)$ holds, then system (2.1) exists periodic solutions near equilibrium $E_k$:

(1) If $\operatorname{Re}(M_k){{\tau }_{\varepsilon }} < 0$, the bifurcating periodic solutions are unstable.

(2) If $\operatorname{Re}(M_k){{\tau }_{\varepsilon }} > 0$, the bifurcating periodic solutions are locally asymptotically stable.

5.   Numerical simulations

In this section, according to the natural mortality in China and the data presented in ref. ^[32], we choose $\Lambda = 100, \mu = 0.00713, \varepsilon = 0.8, \gamma_1 = 0.1, \gamma_2 = 0.1, \alpha_1 = 0.025, \alpha_2 = 0.025, \beta = 0.00001$. According to Eq. (3.2), we obtain $S_{1}^{*} = \frac{\Lambda }{\mu } = 14025, I_{1}^{*} = 0, Q_{1}^{*} = 0, R_{1}^{*} = 0, S_2^{*} = 93213, I_2^{*} = -605.719, Q_2^{*} = -3667.4, R_2^{*} = -59932.$ Obviously, assumption (H1) holds and equilibrium $E_1$ is locally asymptotically stable when $\tau = 0$. However, the equilibrium $E_2$ is unstable for $I_2^{*} < 0$.

Remark 2: If $I_2^{*} < 0$, the number of infected individuals is less than zero, which is inconsistent with the facts. Additionally, according to our theoretical analysis, when $I_2^{*} < 0$, the equilibrium $E_2$ is unstable, which is consistent with the facts.

Substituting these parameter values into Eqs. (3.10)-(3.12), by using MATLAB, we can obtain $\omega_0 = 0.8000, Q_1 = 0.0102, P_1 = 0.9999, \tau_1^{(0)} = 1.9509$.

Thus, the equilibrium $E_1$ is locally asymptotically stable when $\tau\in[0, \tau_1^{(0)})$, and Hopf bifurcation occurs near the equilibrium $E_1$ when $\tau = \tau_1^{(0)}$. According to Eq. (4.9) and Eq. (4.13), we obtain $\rm{Re}(M_1) < 0$, $\rm{Re}({{H}_{1}}{{\chi }_{1}}) > 0$. The bifurcating periodic solution is unstable due to Theorem 4.1. We find that if the time delay from infection to isolation is over the critical value, the epidemic will not be controlled under a stable situation, and the epidemic situation will be more severe.

When $\tau = 0$, we choose the initial value $(1500, 1000,100,200)$, and the solution corresponding to a locally asymptotically stable equilibrium is shown in Figure 2. Thus, when $\tau = 0$, once the individuals are infected, they will be quarantined and they will not infect others, and the epidemic situation can be controlled into a stable situation.

When $\tau = 1.5\in (0, \tau_1^{(0)})$, we choose initial values $(1,500, 1,000,100,200)$, and the equilibrium $E_1$ of system (2.1) is locally asymptotically stable. See Figure 3.

When $\tau = 2\in (\tau_1 ^{(0)}, +\infty)$, we choose the initial values $(1,500, 1,000,100,200)$, and the equilibrium $E_1$ is unstable, as shown in Figure 4. This means that if $\tau$ is larger than the critical value $\tau_1^{(0)}$, the epidemic situation cannot be controlled into a stable state.

According to Figures 2-4, we find that the time delay from infection to isolation has a great influence on the epidemic situation in $E_1$, and with the increasing time delay, the epidemic situation will become increasingly severe.

According to the natural mortality in China and the data presented in ref. ^[32], we choose another group of parameters, namely, $\Lambda = 200, \mu = 0.00713, \varepsilon = 0.8, \gamma_1 = 0.05, \gamma_2 = 0.05, \alpha_1 = 0.025, \alpha_2 = 0.025, \beta = 0.0001$. According to Eq. (3.2), we obtain $S_1^{*} = 28,050, I_1^{*} = 0, Q_1^{*} = 0, R_1^{*} = 0, S_{2}^{*} = 8,821.300, I_{2}^{*} = 155.424, Q_{2}^{*} = 1,513.9, R_{2}^{*} = 11,703.$ We can find that assumption (H1) does not hold and the equilibrium $E_1$ is unstable. However, the equilibrium $E_2$ is locally asymptotically stable when $\tau = 0$.

Substituting these parameter values into Eq. (3.16), we can obtain $c_1 = -0.0269, c_2 = -0.00030938$, and Eq. (3.16) has only one positive root, $z_1 = 0.0356$. According to Eq. (3.15)and Eq. (3.17), we can obtain $\omega_2 = 0.1887, Q_2 = 0.9892, P_1 = 0.1463, \tau_{2, 3}^{(0)} = 0.7782$ by using MATLAB.

Thus, the equilibrium $E_2$ is locally asymptotically stable for $\tau\in[0, \tau_{2, 3}^{(0)})$, and Hopf bifurcation occurs near the equilibrium $E_2$ when $\tau = \tau_{2, 3}^{(0)}$. According to Eq. (4.9) and Eq. (4.13), we obtain $\rm{Re}(M_1) > 0, \rm{Re}({{H}_{1}}{{\chi }_{1}}) > 0$. Thus, $\rm{Re}(M_1)\tau_{\varepsilon } < 0$, the bifurcation solution is unstable. In $E_2$, the epidemic situation cannot be controlled under a stable situation.

When $\tau = 0$, we choose the initial values $(1,000, 50, 2,000, 1,000)$, and the solution corresponding to a locally asymptotically stable equilibrium is shown in Figure 5.

When $\tau = 0.5\in (0, \tau_{2, 3}^{(0)})$, we choose the initial value $(1,000, 50, 2,000, 1,000)$, and the equilibrium $E_2$ of system (2.1) is locally asymptotically stable. See Figure 6.

When $\tau\in (\tau_{2, 3}^{(0)}, +\infty), \tau_{2, 3}^{(0)} = 0.7782$, the equilibrium $E_2$ is unstable. We choose $\tau = 0.8 > \tau_{2, 3}^{(0)} = 0.7782$ and the initial value $(1,000, 50, 2,000, 1,000)$. See Figure 7.

According to numerical simulations, we can find that the smaller $\tau$ is, the better the control effect is. The results of the numerical simulation are consistent with this fact. Then, we provide some explanations for the above results.

Remark 3: For the two groups of parameter values that are $\Lambda = 100, \mu = 0.00713, \varepsilon = 0.8, \gamma_1 = 0.1, \gamma_2 = 0.1, \alpha_1 = 0.025, \alpha_2 = 0.025, \beta = 0.00001$ and $\Lambda = 200, \mu = 0.00713, \varepsilon = 0.8, \gamma_1 = 0.05, \gamma_2 = 0.05, \alpha_1 = 0.025, \alpha_2 = 0.025, \beta = 0.0001$, we have the following conclusions.

For the first group of parameter values, we can control $\tau < 1.9509$, and the epidemic situation can gradually become stable. However, if $\tau > 1.9509$, the epidemic cannot be controlled into a stable situation. This means that if we can make the time delay from infection to isolation within 1.9509, we can make the number of infected individuals gradually tend to zero, and the epidemic situation will be controllable. However, if we cannot make the time delay from infection to isolation within 1.9509, the epidemic situation will be uncontrollable.

For the second group of parameter values, we can control $\tau < 0.7782$, and the epidemic situation can be gradually stable. If $\tau > 0.7782$, the epidemic cannot be controlled into a stable situation. This means that if we can make the time delay from infection to isolation within 0.7782, we can make the number of infected individuals gradually tend to stabilize, and the epidemic situation will be controllable. However, if we cannot make the time delay from infection to isolation within 0.7782, the epidemic situation will be uncontrollable.

Remark 4: According to our theoretical analysis, with the increase in $\tau$, the epidemic situation will be increasingly severe, and COVID-19 will further threaten the health of humans. To control the epidemic situation, we need to take some measures to shorten the time from infection to isolation. In addition, for epidemic situations in different areas, we can deduce the corresponding parameter values and use the model we constructed to obtain the stability determination and the strategies to control epidemic situations.

Our work is reproducible. For different regions, we can calculate the different critical time delays by changing variables such as the transmission rate from $S$ to $I$, the transmission rate from $I$ to $Q$ and the transmission rate from $I$ to $R$. We can obtain the critical time delay based on Eqs. (3.10)-(3.12) and Eqs. (3.15)-(3.17). Then, according to Theorem 3.1, we know that the equilibria are locally asymptotically stable if the time delay is within the critical time delay. Finally, we can obtain the normal form of Hopf bifurcation by using the multiple time scales method and determine whether the epidemic situation can be controlled into a stable situation based on Theorem 4.1. The time delay from infection to isolation is the most important variable in our model. We calculate the critical time delay based on the parameter values. We choose parameter values according to the natural mortality in China and the data presented in ref. ^[32]. Then, we change the value of the time delay from infection to isolation, and we know that the equilibria are locally asymptotically stable if the time delay is within the critical time delay through numerical simulations.

6.   Conclusion

In this paper, we constructed a DDE according to the characteristics of the COVID-19 epidemic phenomenon based on the $SIQR$ model. We considered the existence and stability of equilibria in the above delayed $SIQR$ model. We also analyzed the existence and dynamic properties of Hopf bifurcation associated with both equilibria. We chose two groups of parameter values according to the natural mortality in China and the data presented in ref. ^[32]. Numerical simulations were carried out to verify the analytical results.

According to the results of numerical simulations, the epidemic solution would be more severe with increasing the time delay from infection to isolation. In addition, we predicted the stability of epidemic solutions in different areas and provided effective strategies to control the epidemic.

Acknowledgments

This study was funded by Fundamental Research Funds for the Central Universities of China. (Grant No. 2572019BC14), the Heilongjiang Provincial Natural Science Foundation of China (Grant No. LH2019A001) and College Students Innovations Special Project funded by Northeast Forestry University of China (No. 202010225035).

Conflict of interest

All authors declare no conflicts of interest in this paper.