Mathematical modeling and dynamic analysis of SIQR model with delay for pandemic COVID-19

1.   Introduction

Coronavirus disease 2019 (COVID-19) remains an on-going global pandemic at present. The World Health Organization declared COVID-19 as a Public Health Emergency of International Concern (PHEIC) on January 30, 2020. According to data released by Johns Hopkins University, there are 37, 213, 592 confirmed cases and 1, 072, 959 deaths in 188 countries and regions around the world by October 11, 2020 ^[1]. The disease is caused by a novel coronavirus named severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) ^[2]. Although most SARS-CoV-2-infected cases have asymptomatic or mild-to-moderate diseases, around 10$\%$ of those infected may develop severe pneumonia and other associated organ malfunctions ^[3].

In recent investigations, many scholars studied different epidemic models of COVID-19. In Zhang's study ^[4], a new mathematical model (SEIRD) was proposed, which is constructed with five classes including susceptible, exposed, infected, recovered and deaths to describe the possibility of transmission in a given general population. Bhadauria et al. ^[5] studied the SIQ model by using the stability theory of nonlinear ordinary differential equations. Li et al. ^[6] developed a numerical method preserving positivity for a stochastic SIQS epidemic model. In Ref. ^[7], Higazy proved the existence of a stable solution for the fractional order COVID-19 SIDARTHE model. A $\theta {\rm{ - }}$SEIHQRD model which is more relevant to COVID-19 was developed by Ramos et al. ^[8]. Nisar et al. ^[9] constructed the SIRD model and verified the correctness. Batistela et al. ^[10] proposed an SIRSi model of COVID-19. Especially, Paré et al. ^[11] proposed traditional group models, continuous-time and discrete-time versions of the models with non-trivial networks on simple SIR-based models, and supported the need for networked models through presenting a set of simulations.

In the propagation process of COVID-19, we consider that there is a time-delay from infection to recovery. At present, there are some researches on epidemic models with time-delay. Zhu et al. ^[12] constructed a time delay reaction-diffusion model that is closer to the actual spread of the COVID-19 epidemic with considering the time delay effect of infected persons during the spread of the epidemic. In Ref. ^[13], the mathematical modeling of COVID-19 fatality trends was constructed by Scheiner et al., which was based on infection-to-death delay rule. Wei et al. ^[14] considered the time delay from susceptible individuals to infected individuals, thus, proposed a new SVEIR epidemic disease model with time delay, and analyzed the dynamic behavior of the model under pulse vaccination. In Ref. ^[15], Mukherjee considered an S-I epidemic model with time delay, the time delay is the immune period and the incubation period, and gave an estimate on the length of delays for system which is stable in the absence of delays remains stable.

At the same time, many researchers analyzed the stability of the COVID-19 model. In Ref. ^[16], Araz dealt with a mathematical model about COVID-19 spread, and analyzed global and local stability for the considered model. Annas et al. ^[17] carried out the stability analysis and numerical simulation of the SEIR model on the spread of COVID-19 in the research. Besides, some scholars conducted mathematical analysis on the principle of COVID-19 infection. Almocera et al. ^[18] studied an in-host model, and the stability of a unique positive equilibrium point, with viral load ${{V}^{\rm{*}}}$, suggested that the virus may replicate fast enough to overcome T cell response and cause infection. In Ref. ^[19], Samanta analyzed the stability of the proposed model to control the epidemic.

And some scholars predicted the epidemic situation. León et al. ^[20] proposed an SEIARD mathematical model and attempted to forecast the evolution of the outbreak. Youssef et al. ^[21] carried out numerical verification and predictions of the proposed SEIR model, and compared the results with the real data due to the spreading of the COVID-19 in Saudi Arabia.

In addition, some scholars used the models with several compartments to put forward analysis and opinions on epidemic prevention and control. In Ref. ^[22], Carli et al. proposed a multi-region SIRQTHE model and an optimal control approach, which supported governments in defining the most effective strategies to be adopted during post-lockdown mitigation phases in a multi-region scenario. Giordano et al. ^[23] established a SIDARTHE model that predicted the course of the epidemic to help plan an effective control strategy. On the basis of the expected utility theory, Odagaki ^[24] carried out a theoretical framework to find out an optimum strategy for minimizing the maximum number of infeced and for controlling the outbreak of pandemic. In Ref. ^[25], Saldaña et al. developed a compartmental epidemic model of the COVID-19 epidemic outbreak to evaluate the theoretical impact of plausible control interventions.

In this paper, we establish the COVID-19 epidemic model based on the SIQR model. However, these large deviations between model predictions and the actually recorded numbers stem from the uncertainty of the underlying SIQR model parameters: they may not be sufficiently well known for the novel COVID-19 pandemic yet. Generally, considering the latent period gives rise to models with the incorporation of delays to solve problems. Therefore, time-delay has important biologic meaning in epidemic models. Differently, by describing the time-delay from the infection period to the recovery period, we propose a new delay SIQR epidemic model with horizontal transmission in our paper, and obtain the critical treatment time through calculation and stability analysis, which can contribute greatly to medicine.

The rest of the paper is organized as follows: In Section 2, considering the treatment time of epidemic, we present a delayed differential equation for COVID-19. In Section 3, we analyze the existence and stability of equilibria and the existence Hopf bifurcation of the system. In Section 4, we derive the normal form of Hopf bifurcation for above model, then determine the direction of Hopf bifurcation and the stability of bifurcating periodic solutions. In Section 5, we present the simulated results to verify the correctness of the theoretical analysis, and show the effect of treatment time on epidemic. Finally, some conclusions are shown in Section 6.

2.   Modeling

COVID-19 is similar to other infectious diseases. The susceptible population can be infected by COVID-19 carriers. Some of the infected people are quarantined, and some are killed by COVID-19. Both infected and quarantined people may be cured after treatment. A feature of COVID-19 can be obtained by combining a large number of cases, that is, cured infected people are difficult to be infected again. In Ref. ^[26], Liu et al. proposed an infectious disease model as follows:

where $S$, $I$, $Q$, $R$ denote the numbers of susceptible, infective, quarantined and removed, $N = S + I + Q + R$ is the number of total population individuals. The parameter $\Lambda$ is the recruitment rate of $S$ corresponding to births and immigration; $\beta$ denotes the average number of adequate contacts; $\mu$ is the natural death rate; $\gamma$ and $\epsilon$ denote the recover rates from group $I$, $Q$ to $R$, respectively; $\delta$ denotes the removal rate from $I$; $\alpha$ is the disease-caused death rate of $I$ and $Q$. The parameters involved in system (2.1) are all positive constants.

On the basis of reference ^[26], we define $\beta$ as the contact rate, which eliminates the calculation of the total population number $N$ and simplifies the form of the equation. It is a good way to reduce the error of the model results by reducing the parameters. Combined with practical reports, it is known that some quarantined patients will receive treatment, and the mortality of treated patients will be greatly reduced. Therefore, unlike the model in Ref.^[26], we believe that the COVID-19 mortality of infected and quarantined person is different. Meanwhile, the uncertainty of the model parameters is also the main reason for the difference between the predicted results and the actual. Here, we consider the time-delay ($\tau$) from infection to recovery process on the basis of the model in Ref.^[26], and the critical treatment time obtained by the calculation with time-delay is helpful to the treatment of COVID-19, which is rarely studied before and should be paid attention to in severe epidemic areas. In particular, the model can change the parameters according to the situation of different regions to get the corresponding critical treatment time, and take reasonable measures to control the epidemic effectively. Finally, we present a new COVID-19 epidemic model, and the relationships between the four populations (susceptible population($S$), infected population($I$), quarantined population($Q$), recovered population($R$)) are obtained, as shown in Figure 1.

Thus, we construct the following COVID-19 epidemic model:

where $\Lambda, \; \mu, \; \beta, \; \delta, \; \gamma, \; \epsilon, \; \alpha_1, \; \alpha _2$ are parameters; $S, \; I, \; Q, \; R$ are control variables and $\tau$ is the time-delay. The specific definitions are given in the Table 1.

3.   Stability analysis of equilibrium and existence of Hopf bifurcation

In this section, system (2.2) is considered. Obviously, system (2.2) has two equilibria:

where

$S_1^* = \frac{\Lambda}{\mu}$, $I_1^* = 0$, $Q_1^* = 0$, $R_1^* = 0$, $S_2^*$ = $\frac{{{\gamma}{\rm{ + }}{\alpha _1}{\rm{ + }}\delta {\rm{ + }}\mu }}{\beta }$, $I_2^*$ = $\frac{{\Lambda \beta - \mu ({\gamma} + {\alpha _1} + \delta + \mu)}}{{\beta ({\gamma} + {\alpha _1} + \delta + \mu)}}$,

$Q_2^*$ = $\frac{{\delta [\Lambda \beta - \mu ({\gamma} + {\alpha _1} + \delta + \mu)]}}{{\beta ({\gamma} + {\alpha _1} + \delta + \mu)({\alpha _2} + {\epsilon} + \mu)}}$, $R_2^*$ = $\frac{{[\Lambda \beta - \mu ({\gamma} + {\alpha _1} + \delta + \mu)][{\gamma}({\epsilon} + {\alpha _2} + \mu) + {\epsilon}\delta]}}{{\beta \mu ({\gamma} + {\alpha _1} + \delta + \mu)({\alpha _2} + {\epsilon} + \mu)}}$.

Firstly, we consider the first equilibrium $E_1 = (S_1^*, I_1^*, Q_1^*, R_1^*) = (\frac{\Lambda}{\mu}, 0, 0, 0)$. Transferring the equilibrium to the origin and linearizing the system around it, we obtain the characteristic equation of the linearized system as follows:

where $S_1^*$ = $\frac{\Lambda }{\mu }$.

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

Eq. (3.2) has four roots: ${\lambda _1}{\rm{ = }}{\lambda _2}{\rm{ = - }}\mu, \; {\lambda _3}{\rm{ = - }}{\epsilon} - {\alpha _2} - \mu, \; {\lambda _4} = \beta S_1^* - \delta - {\alpha _1} - \mu- {\gamma}$, due to $\mu > 0, {\epsilon} > 0, {\alpha _2} > 0$, so in the actual situation, ${\lambda _1} < 0$, ${\lambda _2} < 0$, ${\lambda _3} < 0$.

We consider the following assumption:

(H1) $\frac{{\beta \Lambda }}{\mu } < \mu + \delta + {\alpha _1} + {\gamma}$.

When (H1) holds, all the roots of Eq. (3.2) have negative real parts, and the equilibrium $E_1$ is locally asymptotically stable when $\tau = 0$.

When $\tau > 0$, let $\lambda = \mathrm{i}\omega\; (\omega > 0)$ be a root of Eq. (3.1). Due to ${\lambda _1} < 0$, ${\lambda _2} < 0$, ${\lambda _3} < 0$, actually, we only need to discuss the following equation:

Substituting $\lambda = \mathrm{i}\omega\; (\omega > 0)$ into Eq. (3.3) and separating the real and imaginary parts, we obtain

Eq. (3.4) leads to

Adding the square of two equations of Eq. (3.5), we have

Therefore, we give the following assumption:

(H2) ${(\beta S_1^* - \delta - {\alpha _1} - \mu)^2} - \gamma^2 < 0$.

If (H2) holds, then Eq. (3.6) has one positive root $\omega_0 = \sqrt {\gamma^2 - {{(\beta S_1^* - \delta - {\alpha _1} - \mu)}^2}}$. Substituting $\omega_0$ into Eq. (3.5), we get

where

Lemma 3.1. If (H2) holds, when $\tau = \tau_1^{(j)}\; (\; j = 0, 1, 2, \cdots)$, then Eq. (3.1) has a pair of pure imaginary roots $\pm \mathrm{i}\omega_0$, and all the other roots of Eq. (3.1) have nonzero real parts.

Furthermore, let $\lambda(\tau) = \alpha(\tau)+\mathrm{i}\omega(\tau)$ be the root of Eq. (3.1) satisfying $\alpha(\tau_1^{(j)}) = 0$, $\omega(\tau_1^{(j)}) = \omega_0\; (\; j = 0, 1, 2, \cdots)$.

Lemma 3.2. If (H2) holds, we have the following transversality conclusions:

${\rm{Re}}(\frac{{{\rm{d}}\tau }}{{{\rm{d}}\lambda }})\left| {_{\tau = \tau _1^{(j)}}} \right. = {\rm{Re}}{(\frac{{{\rm{d}}\lambda }}{{{\rm{d}}\tau }})^{ - 1}}\left| {_{\tau = \tau _1^{(j)}}} \right. = \frac{1}{{\gamma^2}} > 0$, where $j = 0, 1, 2, \cdots$.

Secondly, for the other equilibrium $E_2$ = $(S_2^*, I_2^*, Q_2^*, R_2^*)$ of the system (2.2), similarly, transferring the equilibrium to the origin and linearizing the system around it, we obtain the characteristic equation of the linearized system as follows:

where $S_2^*$ = $\frac{{{\gamma}{\rm{ + }}{\alpha _1}{\rm{ + }}\delta {\rm{ + }}\mu }}{\beta }$, $I_2^*$ = $\frac{{\Lambda \beta - \mu ({\gamma} + {\alpha _1} + \delta + \mu)}}{{\beta ({\gamma} + {\alpha _1} + \delta + \mu)}}$.

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

Eq. (3.9) has four roots: ${\lambda _1}{\rm{ = - }}\mu$, ${\lambda _2}{\rm{ = - }}{\epsilon} - {\alpha _2} - \mu$, ${\lambda _3} = \frac{{ - \mu - \beta S_2^* + \sqrt {{{(\mu + \beta S_2^*)}^2} - 4{\beta ^2}S_2^*I_2^*} }}{2}4$, ${\lambda _4} = \frac{{ - \mu - \beta S_2^* - \sqrt {{{(\mu + \beta S_2^*)}^2} - 4{\beta ^2}S_2^*I_2^*} }}{2}$. Due to $\mu > 0$, ${\epsilon} > 0$, ${\alpha _2} > 0$, in the actual situation, ${\lambda _1} < 0$, ${\lambda _2} < 0$.

Note that $\beta ^2S_2^*I_2^* > 0$, we consider the following assumption:

(H3) $\frac{{\beta \Lambda }}{\mu } > \mu + \delta + {\alpha _1} + {\gamma}$.

Therefore, under the assumption (H3), all the roots of Eq. (3.9) have negative real parts, and the equilibrium $E_2 = (S_2^*, I_2^*, Q_2^*, R_2^*)$ is locally asymptotically stable when $\tau = 0$.

When $\tau > 0$, let $\lambda = \mathrm{i}\omega\; (\omega > 0)$ be a root of Eq. (3.8). Due to $\mu > 0$, ${\epsilon} > 0$, ${\alpha _2} > 0$, actually, we only need to consider the equation:

Substituting $\lambda = \mathrm{i}\omega\; (\omega > 0)$ into Eq. (3.10) and separating the real and imaginary parts, we have

Eq. (3.11) leads to

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

where ${c_1} = {(\mu + \beta I_2^*)^2} - 2{\beta ^2}S_2^*I_2^*$, ${c_0} = {({\gamma}\mu + {\gamma}\beta I_2^* - {\beta ^2}S_2^*I_2^*)^2} - {({\gamma}\mu + {\gamma}\beta I_2^*)^2}$. Therefore, we give the following assumptions:

(H4) ${c_0} < 0.$

(H5) $c_1^2 - 4{c_0} > 0$, ${c_1} < 0, {c_0} > 0$.

If (H4) holds, then Eq. (3.13) has only one positive real root $z_1$. If (H5) holds, then Eq. (3.13) has two positive real roots $z_2$ and $z_3$. Substituting $\omega_k = \sqrt{z_k}\; (\; k = 1, 2, 3)$ into Eq. (3.12), we get

where

Lemma 3.3. If (H4) or (H5) holds, when $\tau = \tau_{2, k}^{(j)}\; (\; k = 1, 2, 3;\; j = 0, 1, 2, \cdots)$, then Eq. (3.8) has a pair of pure imaginary roots $\pm \mathrm{i}\omega_k$, and all the other roots of Eq. (3.8) have nonzero real parts.

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

Lemma 3.4. If (H4) or (H5) holds, and ${z_k} = \omega _k^2$, $h'({z_k}) \ne 0$, then we have the following transversality conclusions:

${\rm{Re}}(\frac{{{\rm d}\tau }}{{{\rm d}\lambda }}) \left| {_{\tau = \tau _{2, k}^{(j)}}} \right. = {\mathop{\rm Re}\nolimits} {(\frac{{{\rm d}\lambda }}{{{\rm d}\tau }})^{ - 1}}\left| {_{\tau = \tau _{2, k}^{(j)}}} \right. = \frac{{h'(z_k)}}{{\gamma^2[{{(\mu + \beta I_2^*)}^2} + {\omega _k^2}]}}\ne 0$, $k = 1, 2, 3, \; j = 0, 1, 2, \cdots$.

Theorem 3.1. We show the conclusion associated with two equilibria of the system (2.2).

(1) If the assumptions (H1) and (H2) hold, the equilibrium $E_1$ of the system (2.2) undergoes Hopf bifurcation at $\tau = \tau_1^{(j)}\; (\; j = 0, 1, 2, \cdots)$, where $\tau_1^{(j)}$ is given by Eq. (3.7), and we have: when $\tau \in [0, \tau _1^{(0)})$, the equilibrium $E_1$ is locally asymptotically stable, and the equilibrium $E_1$ is unstable when $\tau > \tau _1^{(0)}$.

(2) If the assumptions (H4) or (H5) holds, the equilibrium $E_2$ of the system (2.2) undergoes Hopf bifurcation at $\tau = \tau_{2, k}^{(j)}\; (\; k = 1, 2, 3;\; j = 0, 1, 2, \cdots)$, where $\tau_{2, k}^{(j)}$ is given by Eq. (3.14), and

(a) If the assumptions (H3) and (H4) hold, $h(z)$ has one positive root $z_1$, then when $\tau \in [0, \tau _{2, 1}^{(0)})$, the equilibrium $E_2$ is locally asymptotically stable, and the equilibrium $E_2$ is unstable when $\tau > \tau _{2, 1}^{(0)}$.

(b) If the assumptions (H3) and (H5) hold, $h(z)$ has two positive roots $z_2$ and $z_3$, we suppose $z_2 < z_3$, then $h'({z_2}) < 0, h'({z_3}) > 0$, note that $\tau _{2, 2}^{(0)} > \tau _{2, 3}^{(0)}$. Then $\exists$ $m \in N$ makes $0 < \tau _{2, 3}^{(0)} < \tau _{2, 2}^{(0)} < \tau _{2, 3}^{(1)} < \tau _{2, 2}^{(1)} < \cdots < \tau _{2, 2}^{(m-1)} < \tau _{2, 3}^{(m)} < \tau _{2, 3}^{(m + 1)}$. When $\tau \in [0, \tau _{2, 3}^{(0)})\cup \bigcup\limits_{l = 1}^m {} (\tau _{2, 2}^{(l-1)}, \tau _{2, 3}^{(l)})$, the equilibrium $E_2$ of the system (2.2) is locally asymptotically stable, and when $\tau \in \bigcup\limits_{l = 0}^{m-1} {}(\tau _{2, 3}^{(l)}, \tau _{2, 2}^{(l)}) \cup (\tau _{2, 3}^{(m)}, + \infty)$, the equilibrium $E_2$ of the system (2.2) is unstable.

4.   Direction of Hopf bifurcation and stability of periodic solution

In this section, we discuss the normal form of Hopf bifurcation for the system (2.2) by using the multiple time scales method. Combining with actual situation, we concern about the impact of treatment time on epidemic control. Therefore, we consider the time-delay $\tau$ as a bifurcation parameter, let $\tau {\rm{ = }}{\tau _c} + \varepsilon {\tau _\varepsilon }$, where $\tau_c$ is the critical value of Hopf bifurcation given in Eq. (3.7) or Eq. (3.14) respectively, $\tau _\varepsilon$ is the disturbance parameter, and $\varepsilon$ is the dimensionless scale parameter. We suppose the characteristic Eq. (3.1) and Eq. (3.8) have eigenvalue $\lambda = {\rm i}\omega^{(k)}\; (\; k = 1, 2)$, where $\omega^{(1)} = \omega_0$, $\omega^{(2)} = \omega_1$, $\omega_2$ or $\omega_3$, at which system (2.2) undergoes a Hopf bifurcation at equilibrium $E_k = (S_k^*, I_k^*, Q_k^*, R_k^*)$, $k = 1, 2$, respectively.

Then system (2.2) can be written as

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

We suppose $h_k$, $h_k^*$($k = 1, 2$) are the eigenvector of the corresponding eigenvalue $\lambda = {\rm i}\omega^{(k)}$, $\lambda = -{\rm i}\omega^{(k)}$ respectively of Eq. (4.1) for equilibrium $E_k$, and satisfies $\langle {h_k^*, h_k} \rangle = {\overline {h_k^*} ^{\rm T}}h_k = 1$. By simple calculation, we can get:

where ${d_k} = \frac{{{\beta ^2}S_k^*I_k^*}}{{{\beta ^2}S_k^*I_k^* - {{(-\lambda + \mu + \beta I_k^*)}^2}}}, \; k = 1, 2.$

We suppose the solution of Eq. (4.1) as follows:

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)]^{\rm 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)]^{\rm T}}.$

The derivative with respect to t is transformed into:

where ${D_i} = \frac{\partial }{{\partial {T_i}}}$, $i = 0, 1, 2, \cdots.$

We make ${X_j} = {({S_j}, {I_j}, {Q_j}, {R_j})^{\rm 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}})^{\rm T}} = {X_j}({T_0} - {\tau _c}, {T_1}, {T_2}, \cdots)$, $j = 1, 2, 3, \cdots.$

From Eq. (4.3) we can get:

The Taylor expansion of $X(t - \tau)$ is carried out:

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

Substituting Eqs. (4.3) $\sim$ (4.5) into Eq. (4.1), and balancing the coefficients before $\varepsilon$ on both sides of the equation, the following expression is obtained:

Thus Eq. (4.6) has the following solution form:

The expression of the coefficient before ${\varepsilon ^2}$ is as follows:

Substituting Eq. (4.7) into the right hand side of Eq. (4.8), and the coefficient vector of ${{\rm e}^{{\rm i}\omega^{(k)} {T_0}}}$ is denoted by $m_1$. According to the solvable condition $\langle {h_k^*, {m_1}} \rangle = 0$, the expression of $\frac{{\partial G}}{{\partial {T_1}}}$ can be obtained as follows:

where ${N_k} = \frac{{{{({\rm i}\omega^{(k)} + \mu + \beta I_k^*)}^2}{\gamma} \cdot {\rm i}\omega^{(k)} {{\rm e}^{ - {\rm i}\omega^{(k)} {\tau _c}}}}}{{{{({\rm i}\omega^{(k)} + \mu + \beta I_k^*)}^2} - {\beta ^2}S_k^*I_k^* - {{({\rm i}\omega^{(k)} + \mu + \beta I_k^*)}^2}{\gamma}{\tau _c}{{\rm e}^{ - {\rm i}\omega^{(k)} {\tau _c}}}}}, \; k = 1, 2.$

Since $\tau_\varepsilon$ is a disturbance parameter, we only consider its effect on the linear part. It has little effect on the high order, so it can be ignored. Therefore, we ignore the part containing $\tau_\varepsilon$ in the higher order. We suppose:

where

The expression of the coefficient before ${\varepsilon ^3}$ is:

Substituting Eq. (4.7), Eq. (4.10) and Eq. (4.11) into the right hand side of Eq. (4.12), and the coefficient vector of ${{\rm e}^{{\rm i}\omega^{(k)} {T_0}}}$ is denoted by $m_2$. According to the solvable condition $\langle {h_k^*, {m_2}} \rangle = 0$, the expression of $\frac{{\partial G}}{{\partial {T_2}}}$ can be obtained as follows:

where

Let $G \to G/{\varepsilon }$, then the deduce to third-order normal form of Hopf bifurcation of system (2.2) is:

where $N_k$ is given in Eq. (4.9), and $\xi _k$ is given in Eq. (4.13), $k = 1, 2$.

By substituting $G = r{{\rm e}^{{\rm i}\theta }}$ into Eq. (4.14), the following normal form of Hopf bifurcation in polar coordinates can be obtained:

where $N_k$ is given in Eq. (4.9), and $\xi _k$ is given in Eq. (4.13), $k = 1, 2$.

According to the normal form of Hopf bifurcation in polar coordinates, we just need to consider the first equation in system (4.15). Thus, there is the following theorem:

Theorem 4.1. For the system (4.15), when $\frac{{{\mathop{\rm Re}\nolimits} (N_k){\tau _\varepsilon }}}{{{\mathop{\rm Re}\nolimits} ({\xi _k})}} < 0\; (\; k = 1, 2)$, there is a nontrivial fixed point $r = \sqrt { - \frac{{{\mathop{\rm Re}\nolimits} ({N_k}){\tau _\varepsilon }}}{{{\mathop{\rm Re}\nolimits} ({\xi _k})}}}$, and system (2.2) has periodic solution:

(1) If ${\mathop{\rm Re}\nolimits} (N_k){\tau _\varepsilon } < 0$, then the periodic solution reduced on the center manifold is unstable.

(2) If ${\mathop{\rm Re}\nolimits} (N_k){\tau _\varepsilon } > 0$, then the periodic solution reduced on the center manifold is stable.

5.   Numerical simulations

In this section, according to the data presented in Refs. ^[27] and ^[28], we let $\Lambda = 10.48$, $\beta = 0.0004$, $\delta = 0.7$, ${\alpha _1} = 0.0484$, ${\alpha _2} = 0.022778$, ${\gamma} = 0.42386$, ${\epsilon} = 0.475$, $\mu = 0.00714$. Then, system (2.2) only has one nonnegative equilibrium $E_1 = (S_1^{\rm{*}}, I_1^{\rm{*}}, Q_1^{\rm{*}}, R_1^{\rm{*}}) \approx (1467.8, 0, 0, 0).$ Obviously, the assumption (H1) holds, thus the equilibrium $E_1$ is locally asymptotically stable when $\tau$ = 0.

Substituding these parameter values into the Eqs. (3.5)$\sim$(3.7), using MATLAB, we can obtain ${\omega _0}\approx0.3900$, $Q_0\approx 0.9201$, $P_0\approx-0.4284$, $\tau _1^{(0)}\approx5.1629$. According to the Theorem 3.1, the equilibrium $E_1$ is locally asymptotically stable at $\tau \in [0, \tau _1^{(0)})$, and when $\tau = \tau _1^{(0)}$, Hopf bifurcation occurs near the equilibrium $E_1$. Then, we obtain ${\mathop{\rm Re}\nolimits} (N_1) > 0$, ${\mathop{\rm Re}\nolimits} (\xi_1) > 0$ from Eqs. (4.9)$\sim$(4.13), thus according to the Theorem 4.1, the system (2.2) has backward periodic solution and the periodic solution is unstable when $\tau_\varepsilon < 0$.

When $\tau = 0$, we choose the initial value $(1473, 1, 2, 1)$ and the corresponding locally asymptotically stable equilibrium is shown in Figure 2.

When $\tau = 4.2 \in [0, \tau _1^{(0)})$, we choose the initial value $(1473, 1, 2, 1)$, the equilibrium $E_1$ is also locally asymptotically stable (see Figure 3). Thus, the epidemic will be controlled eventually.

When $\tau = 5.5 \in (\tau _1^{(0)}, +\infty)$, we choose the initial value $(1473, 1, 2, 1)$, and the equilibrium $E_1$ is unstable shown in Figure 4.

It can be seen from Figure 2$\sim$ Figure 4, when $\tau \in [0, \tau _1^{(0)})$, the equilibrium $E_1$ of system (2.2) is locally asymptotically stable (see Fig. 2 and Fig. 3). When $\tau \in (\tau _1^{(0)}, +\infty)$, the equilibrium $E_1$ of system (2.2) is unstable, and accompanied by the fluctuation with increased amplitude (see Fig. 4). There is no periodic solution, so we believe the theoretical analysis is correct.

Remark 1: Through numerical simulations, it can be found that when the treatment time $\tau < \tau_1^{(0)}$, the shorter treatment time is, the faster the epidemic tends to be stable. If $\tau > \tau_1^{(0)}$, the epidemic cannot be controlled effectively, and accompanied by the fluctuation with increased amplitude with time. Thus, we obtain the critical treatment time $\tau_1^{(0)}$ for the epidemic.

According to the data presented in Refs. ^[27] and ^[28], we choose another set of parameters, $\Lambda = 2.48$, $\beta = 0.86834$, $\delta = 0.3$, ${\alpha _1} = 0.0484$, ${\alpha _2} = 0.022778$, ${\gamma} = 0.42386$, ${\epsilon} = 0.475$, $\mu = 0.00714$. We obtain $E_1 = (S_1^{\rm{*}}, I_1^{\rm{*}}, Q_1^{\rm{*}}, R_1^{\rm{*}})\approx(347.3, 0, 0, 0)$, $E_2 = (S_2^{\rm{*}}, I_2^{\rm{*}}, Q_2^{\rm{*}}, R_2^{\rm{*}})\approx (0.8976, 3.1737, 1.8857, 313.8526)$. Apparently, the assumption (H1) does not hold, but the assumption (H3) holds, therefore, the equilibrium $E_1$ is unstable and the equilibrium $E_2$ is locally asymptotically stable when $\tau = 0$.

Substituding these parameter values into the Eq. (3.13), we get $c_1\approx3.3383$, $c_0\approx-0.4174$. It satisfies the assumption (H4), therefore, Eq. (3.13) has only one positive root $z_1\approx0.1204$. Substituting these parameter values into the Eqs. (3.12)$\sim$(3.14), using MATLAB, we can obtain ${\omega _1}\approx0.3474$, $Q_1\approx 0.5926$, $P_1\approx-0.8055$, $\tau _{2, 1}^{(0)}\approx7.2175$. According to the Theorem 3.1, the equilibrium $E_2$ is locally asymptotically stable at $\tau \in [0, \tau _{2, 1}^{(0)})$, and when $\tau = \tau _{2, 1}^{(0)}$, Hopf bifurcation occurs near the equilibrium $E_2$. Then we obtain ${\mathop{\rm Re}\nolimits} (N_2) > 0$, ${\mathop{\rm Re}\nolimits} (\xi_2) < 0$ by Eqs. (4.9)$\sim$(4.13), according to the Theorem 4.1, the system (2.2) has stable periodic solution near $E_2$, and the periodic solution reduced on the center manifold is stable.

When $\tau$ = 0, we choose the initial value $(0.85, 3, 1.9, 300)$ and the corresponding locally asymptotically stable equilibrium is shown in Figure 5.

When $\tau = 4 \in [0, \tau _{2, 1}^{(0)})$, we choose the initial value $(0.85, 3, 1.9, 300)$, and the periodic solution of this model is locally asymptotically stable shown in Figure 6.

When $\tau = 5 \in [0, \tau _{2, 1}^{(0)})$, we choose the initial value $(0.85, 3, 1.9, 300)$, the equilibrium $E_2$ has a stable equilibrium is shown in Figure 7.

It can be seen from Figure 5$\sim$Figure 7, when $\tau \in [0, \tau _{2, 1}^{(0)})$, the equilibrium $E_2$ of system (2.2) is locally asymptotically stable. Especially, if $\tau$ is larger, the time required for the system (2.2) to be controlled stable is longer.

When $\tau = 7.23 > \tau _{2, 1}^{(0)} = 7.2175$, we choose the initial value $(0.85, 3, 1.9, 300)$, the model has stable periodic solution shown in Figure 8.

When $\tau = 7.5 > 7.23$, we choose the initial value $(0.85, 3, 1.9, 300)$, and the equilibrium $E_2$ is unstable shown in Figure 9.

When $\tau = 7.9 > 7.5$, we choose the initial value $(0.85, 3, 1.9, 300)$, and the equilibrium $E_2$ is unstable shown in Figure 10.

It can be seen from the Figure 8$\sim$Figure 10, when $\tau \in (\tau _{2, 1}^{(0)}, +\infty)$, the equilibrium $E_2$ of system (2.2) is unstable. When $\tau$ approaches the critical value $\tau _{2, 1}^{(0)}$, the equilibrium $E_2$ of system (2.2) exhibits periodic fluctuation and bifurcates stable periodic solutions (see Fig. 8). As $\tau$ increasing, the fluctuation phenomenon lasts for a period of time, and shows increasing volatility trends (see Fig. 9) and disappears eventually (see Fig. 10). According to the Theorem 3.1 and Theorem 4.1, we know that the theoretical analysis is correct.

Remark 2: Through numerical simulations, it can be found that when the treatment time $\tau < \tau _{2, 1}^{(0)}$, the epidemic can be controlled effectively, and the smaller the time-delay $\tau$ is the faster the epidemic is controlled. When the treatment time $\tau \in (\tau _{2, 1}^{(0)}, +\infty)$ is close to $\tau _{2, 1}^{(0)}$, the epidemic will repeatedly occur periodically, otherwise, the epidemic cannot be controlled. Since the small variation in treatment time will lead to the epidemic polarization, therefore, whether to grasp the treatment time is essential to control the epidemic.

Remark 3: In Ref. ^[26], Liu et al. used the Milstein's Higher Order Method mentioned to illustrate their main results. And they obtained that the disease is persistent. Similarly, in our simulation, the epidemic is persistent and tends to stable with time.

In Ref. ^[29], the study carried out numerical simulation. First of all, they changed the value of parameters to get the "without control'' and "with control'' graphs, from which we can see that different populations tend to be stable in the end. At the same time, through comparison, it is concluded that the control technique is useful in controlling the disease. Our simulation results show that the system is locally asymptotically stable in the critical treatment time, and the disease will be effectively controlled. Then they fit real data with the infected class of their model and found that the number of cases grows exponentially with time, the disease would be serious without applying the proper optimal control strategies. And we conclude that if not treat within the critical time, that is to say, without timely and effective treatment, the epidemic will be uncontrollable.

In a word, although the model studied is different from Ref. ^[29], and we focus on the impact of time delay on the epidemic situation, but the overall idea and conclusion are consistent, so we can still verify the correctness of our simulation results. Moreover, we emphasize the time-delay of treatment, so we also obtain the critical treatment time, which can take more targeted measures to control the epidemic. To effectively control the epidemic situation, we need to take measures not only in the external environment, but also in the internal medical care, and now there is little research on the treatment direction. Therefore, the critical treatment time we get is helpful for the alleviation of the disease.

6.   Conclusion

In this paper, according to the propagation characteristics of COVID-19, we have constructed the SIRQ epidemic model with time delay for the COVID-19. We have also analyzed the stability of the equilibria and the existence of Hopf bifurcation associated with both equilibria. Then, we have used the multiple time scales method to derive the normal form of Hopf bifurcation for above COVID-19 epidemic model. Finally, We have chosen two groups of parameter values according to the data presented in Refs. ^[27] and ^[28] for numerical simulations to verify the correctness of the theoretical analysis. Compared with the numerical simulation results of Refs. ^[26] and ^[29], we have got the conclusion that our numerical simulation results were accuracy.

The numerical simulations showed that the small change of time-delay $\tau$ leads to the epidemic from stable to uncontrollable: the smaller $\tau$ is the better the epidemic controlled, but with $\tau$ increasing, the epidemic will occur repeatedly and outbreak eventually. Therefore, it is significant that find the critical treatment time to control the epidemic. In addition, we would predict the critical treatment time of epidemic in different regions according to the relevant parameters, so as to effectively control the epidemic in treatment, and realize the major breakthrough of medicine in the epidemic situation.

As for the part of numerical simulation, according to the opinions of reviewers, we will use more accurate data to get more valuable conclusions in further research.

Acknowledgments

The authors are extremely grateful to the anonymous referees and editor for their careful reading, valuable comments and helpful suggestions, which have helped us to improve the presentation of this work significantly. 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

The Authors declare that this work has no conflict of interest.