Stability analysis of time-delayed SAIR model for duration of vaccine in the context of temporary immunity for COVID-19 situation

1.   Introduction

In early 2020, COVID-19 (Corona Virus Disease 2019) outbroke among the population and became a worldwide pandemic. The COVID-19 is a kind of infectious disease caused by SARS-CoV-2 pathogen (severe acute respiratory syndrome coronavirus 2) that is highly contagious, universally susceptible and widely transmitted. In some serious areas, people must be quarantined at home or in concentrated isolation once exposed to the source of infection, which is detrimental to the development and progress of society. More importantly, the COVID-19 not only affects the social production and life, but also threatens our health. According to the statistics from World Health Organization, more than 200 million people have been diagnosed with COVID-2019 and more than 6 million people have died from this until November 1, 2022. Thereinto, about 20% of the patients require hospitalization and the death rate is nearly up to 2% ^[1]. In particular, in the early stage of the outbreak, the treatment process was difficult because people had never been exposed to the disease. Although we have better research on treatments with the continuous understanding of the virus, the test for human beings is far more than that. The COVID-19 continues mutating as the transmission progress among humans. In the first half of 2021, the first virus strain SARS-CoV-2 Delta variant of concern (B.1.617.2) originating from India, was discovered ^[2]. Compared to the previous strain, the Delta strain showed the higher transmissibility and drug resistance. In November of the same year, the new virus strain Omicron was found with more genetic mutations and easier transmission noted by The World Health Organization. As a result, vaccination is considered the most effective measure to control the epidemic situation.

For different infectious diseases, scholars have developed different models of infectious diseases to study their dynamical properties and behaviors ^[3,4,5,6]. Accordingly, there are always corresponding vaccines to enhance herd immunity, so as to control the epidemic at the root. In recent studies, many researchers have added vaccine effects into their infectious disease models ^{[7,8,9,10,11]}. In the study of Arino et al. ^[7], an $SLIRA$ model, representing susceptible-latent-infected-recovered-asymptomatic, was established for the avian influenza (H5N1) with vaccination and antiviral therapy. They analyzed how to control the epidemic situation if it evolved into a strain with transmission among humans. In ^[8], Greenhalgh et al. established an $SIRS$ model and analyzed recurrent epidemic cycles of infected diseases when the vaccine-induced immunity was only temporary. Ho and Chao ^[9] took vaccine capacity, strain mismatch and priority group into considerations and constructed an $SEIR$ compartmental model aiming at studying the influenza vaccination policy. An $SVIR$ model with distributed delay and imperfect vaccine was constructed, which included the susceptible, vaccinated, infected and recovered populations ^[10]. In ^[11], an $SVIS$ cyclic epidemic vaccination model was established by Kabir et al.

Since the outbreak of COVID-19, the whole society has attached great importance to it. Many scholars have attempted to study the effect of COVID-19 vaccines on controlling the epidemic situation ^[12,13,14]. Vishaal et al. ^[12] presented a modified age-structured $SIR$ model to analyze the effect of vaccination and age-targeted distributions. They found that population age-distribution with vaccination had a significant effect on transmission and mortality rates. In ^[13], Shen et al. developed a dynamic model with ten compartments of COVID-19 transmission in the four most severely affected states. They evaluated the vaccine effectiveness and coverage required to suppress the COVID-19 outbreaks under certain conditions. There is also a generalized epidemiological $SEIR$ model for COVID-19 vaccine established by Rajaeia et al. ^[14]. It included the susceptible ($S$), exposed ($E$), infected ($I$), quarantined ($Q$), recovered ($R$), deceased ($D$) and insusceptible ($P$) populations to manage the COVID-19 by comparing the differences before and after its vaccine invention.

Previously, many models used ordinary differential equation (ODE) systems to describe the dynamic behavior of diseases. However, there were always some considerable differences between the realistic behavior and the response of their mathematical models. In 1979, Cooke first put forward "time-delay" to study the spread of epidemic diseases which made the models more realistic ^[15]. Since then, many researchers have attempted to consider time-delay in their models ^{[16,17,18,19,20]}. In ^[16], an $SVEIR$ model with time-delay for imperfect vaccine was developed. They considered the time from the susceptible or the vaccinated to the exposed. Yang and Wang ^[17] evaluated the effect of the latency when they established the $SEIRS$ epidemic model. And in ^[18], a delayed $SEIS$ epidemic model was constructed regarding the time used to cure the infectious people. An $SVEIR$ model was constructed ^[19] considering the time used to cure the infectious people as well. Din et al. ^[20] developed an infectious disease model for hepatitis B virus and introduced a time delay to represent the incubation period of the virus in the diffusion term.

Analysis of stability and Hopf bifurcation has always been the focus of epidemic models ^{[21,22,23,24,25,26,27,28]}. In ^[21], an age-structured $SIR$ epidemic model was constructed by Kuniya et al. and was studied the sustained periodic solution through the Hopf bifurcation. Duan et al. ^[22] formulated and studied an age-structured epidemic model with age of recovery. They considered the existence and stability of equilibria and the bifurcation of periodic solutions. Cai and Wang ^[23] demonstrated that the stability of solutions only changed at every turning points of the bounded bifurcation. ^{[24,25,26,27,28]} discussed the dynamical properties of Hopf bifurcation in different systems.

The motivation of our paper is as follows. Firstly, it is reported that there are some cases of infection that have been vaccinated, suggesting that vaccine-induced immunity is not permanent, so it is necessary to strengthen the immunization. But how long our original vaccination will last and when to schedule a mass booster vaccination? That's where we're going to study. Secondly, as the virus continues to mutate, the infection rate gradually increases. What effect will this have on the time of booster vaccination? That's also what we're looking at here. Thirdly, we know that not all countries can accomplish booster vaccination in a short time under the constraints of human, material and financial resources. In order to ensure that the epidemic is within the controllable range, we want to figure out the suitable time of booster vaccination with above limitation.

Through the study of this paper, we have a more accurate grasp of when to carry out booster vaccination. By simulating the current epidemic state, a reasonable time range for booster vaccination is proposed, which can ensure that the epidemic can reach a stable state instead of being out of control. That is meaningful to the formulation of epidemic prevention and control policies.

The structure of this paper is as follows. In Section 2, considering the characteristics of COVID-19, we construct a delayed Susceptible-Asymptomatic-Infected-Recovered (hereinafter referred to as $SAIR$) model associated with temporary immunity. In Section 3, we analyse the existence of equilibria and discuss the existence of Hopf bifurcation. In Section 4, we calculate the normal form of Hopf bifurcation by using multiple time scales method. And in Section 5, we show some numerical simulations and obtain the conclusion of booster vaccination time by analysis. Finally we show the conclusions and propose some suggestions in Section 6.

2.   Mathematical modelling

The susceptible people ($S$) can be infected by COVID-19 through contact with the source of infection. They first become the asymptomatic people ($A$: one group of people with the virus latent in their bodies but with negative results of nucleic acid detection, and another group of people who have positive results but do not show any symptoms). Some of them become the recovered people ($R$) through the immune response of their autoimmune system, while others become the infected people ($I$) with symptoms such as fever and weakness since they are less fit and cannot fully fight off the virus. Such people is either cured, also called the recovered people ($R$), or die from the failure of resuscitation. For those who are vaccinated, since the vaccine is the inactivated virus, similar to the mechanism of self-healing or cure after infection, it also causes the body to produce an immune response and thus produce antibodies. Therefore, we classify vaccinated people as the recovered people ($R$) as well, which can be understood as recovering from the risk of infection.

We consider the two effects of vaccination rates. In general, if there are no diagnosed cases of COVID-19 in a district, the vaccination can also be organized to enhance herd immunity and reduce the risk of infection. So we denote this part as the natural vaccination rate. On the other hand, when there are infection people in an area, it will be in a state of emergency. In addition to the necessary protection, isolation and treatment, the vaccination process will be accelerated as well. However, we believe that the rate of vaccination is not simply linear with respect to the number of infections, but grows in a much faster way. Since all nonlinear functions can be Taylor expanded, we approximate the expansion to its simplest second-order form, denoted as $\alpha S^2$, where $\alpha$ is the vaccine distribution rate which can be interpreted as vaccine distribution from other areas to speed up the vaccination process.

Taking all these considerations into account, we obtain the specific conversions between the four cabins that are given in Figure 1. The specific definitions of variables and paraments are given in the Table 1. In this table, all parameters and variables are positive.

For the COVID-19 vaccine, antibodies will be produced in susceptible people after vaccination, but this immunity is not permanent. After a period of time, the antibody efficiency will gradually decrease until it disappears. So there is a time delay from the recovered people to the susceptible people. Meanwhile, although there is also a period of time when an asymptomatic person turn into infected one, it can be ignored comparing with the time delay from the recovered people to the susceptible people. Combining the above considerations, we use $\tau$ for the time from the recovered people ($R$) to the suspected people ($S$), which is the time for vaccine failure. Thus, we construct the following delayed differential equation:

where the meanings of variables and parameters are given in Table 1, and $\tau$ is the time of vaccine failed.

When the vaccine fails, the antibodies disappear and the person move from $R$ to $S$. But under normal circumstances, in order to prevent the failure of the vaccine to lead to a large area of infection, which leads to the epidemic out of control, we often improve the population immunization rates through booster vaccination. Therefore, it is particularly important to select the right time for booster vaccination to control the epidemic.

Generally speaking, the booster vaccination time is often earlier than the expiration time of the vaccine. But when the booster shot is injected, the amount of antibody in the body is significantly more than that when not vaccinated, and the antiviral efficiency will be greatly improved. So we can believe that the original antibody has no immunity relative to the existing antibody. In addition, due to the constant mutation of the virus, even if the antibody did not disappear, the original antibody has been unable to resist the attack of the mutant strain, it will also have the possibility of infection. Considering the above two points, we can think that the vaccine failure time is consistent with the booster vaccination time.

Let $C([ - \tau, 0], \mathbb{R}_{ + 0}^4)$ be the Banach space of continuous functions mapping interval $[- \tau, 0]$ into $\mathbb{R}_{ + 0}^4$. For any $\phi \in C([ - \tau, 0], \mathbb{R}_{ + 0}^4)$, we define the norm $\Vert \phi \Vert$ = $\sup\limits_{\theta\in [ - \tau, 0]}\left|(\theta)\right|$ and the initial condition of the system (2.1) is

The system (2.1) has a unique solution under the above initial conditions. Let $(S(t), A(t), I(t), R(t))$ be an arbitrary solution of the system under this initial condition, we give the discussion about the positivity and boundedness of solution of the system (2.1).

1) The positivity of solution

When $t\in[0, \tau]$, we solve each equation of system (2.1) and obtain:

In this paper, since we are mainly considering the issues related to booster vaccination, we assume that most people have received their first vaccination. The state is mainly located in the middle and late stages of the epidemic with a small number of infected people relative to the beginning, that is, the number of the recovered people $R$ (the vaccinated people and the cured people) is so far greater than the infected people $I$, so that $B + \gamma R(t-\tau) - mI^2(t) > 0$ is reasonable. By the expression of $S(t)$, we can obtain that $S(t) > 0$ holds when $t\in[0, \tau]$.

Now we prove $I(t) > 0$. We assume that $I(t)$ is not always positive for $t\in[0, \tau]$ and make $t_1$ be the first time that $I(t_1) = 0$, $I'(t_1) < 0$. According to the third equation of the system (2.1), we can obtain $I'(t_1) = \rho_1A(t_1) > 0$. In addition,

Thus, $A(t_1) > 0$ and $I'(t_1) > 0$. Since $t_1$ is the first time that pass through $t$-axis and makes $I(t) < 0$, so $I'(t_1) < 0$. The two conclusion we obtain are contradictory and $I(t) > 0$. In the same way, we can also obtain $A(t)\ge 0$ when $t\in[0, \tau]$.

Finally, we prove $R(t)\ge 0$. Similarly, using the proof by contradiction, we assume that $t_2$ be the first time that $R(t_2) = 0$, $R'(t_2) < 0$. That is to say, for any $t\in[0, t_2]$, we have $R(t) > 0$, so $R(t_2-\tau) > 0$ holds due to $t_2-\tau \in [-\tau, 0]$. By the forth equation of system (2.1), we have $R'(t_2) = uI(t_2)+\rho _2A(t_2)+\alpha S(t_2)+mI(t_2)^2-\gamma R\left(t_2-\tau \right).$ If $\tau = 0$, we can obtain $R'(t_2) = uI(t_2)+\rho _2A(t_)+\alpha S(t_2)+mI(t_2)^2 > 0$ under the condition of $S(t) > 0, A(t) > 0, I(t) > 0$, which is contradictory with assumption. So $R(t)\ge 0$ when $\tau = 0$. If $\tau > 0$, we focus on the expression of $R(t)$ in Eq (2.2). In Subsection 5.1, we find that the values of parameters $u, \rho_2, m$ are larger than $\gamma$ through the analysis of the actual, so it is reasonable to assume that $uI(t)+\rho _2A(t)+\alpha S(t)+mI(t)^2-\gamma R\left(t-\tau \right) > 0$, which can ensures that $R(t) > 0$ according to the expression of $R(t)$ in Eq (2.2).

Overall, when $t\in[0, \tau]$, the all solutions of system (2.1) are positive after we consider the realistic situation and give some assumptions. By the similar approach, it can be shown that all solutions of the system (2.1) are also positive on $[\tau, 2\tau]$ under some restrictions. This conclusion can be extended to the interval $[n\tau, (n+1)\tau]\; (n\in N)$. Therefore, for any $t > 0$, all solutions of system (2.1) are positive with the restrictions.

2) The boundedness of solution

Let $N(t) = S(t)+A(t)+I(t)+R(t)$, and $N(t)$ represents the total number of people at time $t$. Adding four equations of system (2.1), we obtain $N'(t) = B-dN(t)-cI(t)$. According to the discussion of the positivity of solution, we can obtain that $I(t) > 0$ under certain condition. Thus, $N'(t) = B-dN(t)-cI(t)\leq B-dN(t)$. Based on the comparison theorem, we can get:

Therefore, $\mathop {\lim }\limits_{t \to \infty } {\text{sup}}\; N(t) \leq \frac{B }{d }$ and the solution $S(t)$, $A(t)$, $I(t)$, $R(t)$ of the system (2.1) is bounded when $t > 0$.

Under the conditions and assumptions that the positivity of solution exists we discussed above, we mark the feasible region of the system (2.1) is,

Although the positivity and boundedness of our model are conditional, our analysis of the epidemic shows that the present state is satisfying the conditions for the existence of positivity and boundedness, and therefore our model is reasonable and meaningful.

3.   Stability of equilibria and existence of Hopf bifurcation

In order to study the existence of equilibrium, Driessche introduced the next generation matrix method in ^[29]. By this way, they defined the basic reproduction number $\mathcal{R}_0$, that is, the number of secondary infections resulting from a single primary infection into an otherwise susceptible population ^[30]. Since then, many models have discussed the existence conditions of equilibria by the basic reproduction number $\mathcal{R}_0$. According to this method, we can obtain the basic reproduction number $\mathcal{R}_0$ of our model:

The basic reproduction value $\mathcal{R}_0$ is a very important tool when we discuss the spread of a disease, and the threshold condition $\mathcal{R}_0 = 1$ completely determines whether an infectious disease can invade the population or not. When the value is less than 1, it predicts that the disease is relaxed, that is, there is no epidemic prevalence in the population. On the contrary, when the value is greater than 1, the disease will be transformed into an epidemic ^[31,32].

System (2.1) has two equilibria

where

with

The equilibrium $E_1 = (S_{1}^{*}, A_{1}^{*}, I_{1}^{*}, R_{1}^{*})$ always exists. For the equilibrium $E_2 = (S_{2}^{*}, A_{2}^{*}, I_{2}^{*}, R_{2}^{*})$, if $\mathcal{R}_0 > 1$ holds, $q$ is negative, so $S_{2}^*$, $A_{2}^*$, $I_{2}^*$, $R_{2}^*$ are all positive, which means that equilibrium $E_2$ is meaningful under this condition. On the contrary, if $\mathcal{R}_0 < 1$, we can deduce that $I_{2}^* < 0$, so it is nonsignificant in reality. Therefore, we have the following conclusion.

Lemma 3.1. If $\mathcal{R}_0 < 1$, the system (2.1) has unique equilibrium: the disease-free equilibrium $E_1 = (S_{1}^{*}, A_{1}^{*}, I_{1}^{*}, R_{1}^{*})$. If $\mathcal{R}_0 > 1$, the system not only has the equilibrium $E_1$, but also has an endemic equilibrium $E_2 = (S_{2}^{*}, A_{2}^{*}, I_{2}^{*}, R_{2}^{*})$.

Then we make linear transformation: $\tilde{S} = S-S_{k}^{*}$, $\tilde{A} = A-A_{k}^{*}$, $\tilde{I} = I-I_{k}^{*}$, $\tilde{R} = R-R_{k}^{*}$, $k = 1, 2.$ Substitute $\tilde{S}$, $\tilde{A}$, $\tilde{I}$, $\tilde{R}$ into system (2.1) and drop the tilde for convenience, the system (2.1) becomes:

Since we are more concerned about the impact of antibody failure time on the epidemic, we take time delay $\tau$ as bifurcation parameter when discussing the Hopf bifurcation.

3.1.   Analysis for the disease-free equilibrium $E_1$

Firstly, we consider $\mathcal{R}_0 < 1$. The characteristic equation of the linearized system (3.3) about $E_1$ is as follows:

where $\varPi _1 = -\beta _1S_{1}^{*}+\rho _1+\rho _2+u+2d+c, \varPi _2 = \left(-\beta _1S_{1}^{*}+\rho _1+\rho _2+d \right) \left(u+d+c \right) -\rho _1\beta _2S_{1}^{*}.$

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

where $\varPi _1$ and $\varPi _2$ are given by Eq (3.4). So we can obtain the eigenvalues: $\lambda_{11} = -\alpha < 0$, $\lambda_{12} = -\alpha-d-\gamma < 0$ and $\lambda_{13}$, $\lambda_{14}$ satisfying equation $\lambda ^2+\varPi _1\lambda +\varPi _2 = 0$. Due to $\mathcal{R}_0 < 1$, we have:

Thus, all the roots of Eq (3.5) have negative real parts. We can conclude the disease-free equilibrium $E_1$ is locally asymptotically stable when $\tau = 0$.

When $\tau > 0$, $\lambda_{11} < 0$, $\lambda_{13} < 0$ and $\lambda_{14} < 0$ is available with $\mathcal{R}_0 < 1$. So we just need to consider the following equation.

Substitute $\lambda = {\text{i}}\omega \left(\omega > 0 \right)$ into Eq (3.6) and separate the real and imaginary parts, we can obtain:

Thus,

Adding the square of the two equations in Eq (3.8), we obtain:

Therefore, we show the following assumption:

Under (H1), Eq (3.9) has unique positive root: $\omega _0 = \sqrt{\gamma ^2-\left(\alpha +d \right) ^2}$. Then we substitute $\omega _0$ into Eq (3.8). Since $\omega_0 > 0$, $r > 0$ we can get $\sin(\omega\tau) = \frac{\omega}{r} > 0$ and

Let $\lambda\left(\tau \right) = a\left(\tau \right)+{\mathrm{i}}\omega\left(\tau \right)$ be the root of Eq (3.6). If $\tau = \tau_{0}^{\left(j\right)}$, Eq (3.6) has a pair of imaginary roots: $\pm {\mathrm{i}}\omega_0$ with $a(\tau_{0}^{(j)}) = 0$, $\omega(\tau_0^{(j)}) = \omega_0 = \sqrt{\gamma ^2-\left(\alpha +d \right) ^2}$ $\left(j = 0, 1, 2, \cdots\right)$, where $\tau_0^{(j)}$ is given by Eq (3.10).

Lemma 3.2. If (H1) holds, $z_k = \omega_k^2$ and $h'(z_k)\neq0$, we have the transversality condition:

When $\mathcal{R}_0 > 1$, we can get $\lambda_{13}\lambda_{14} = \varPi_2 < 0$ according to Eq (3.4), the characteristic equation has roots whose real parts are greater than zero. Thus, the disease-free equilibrium $E_1$ is unstable.

Theorem 3.1. For system (2.1), there always exists equilibrium $E_1$. And we have the following conclusions.

1) If $\mathcal{R}_0 < 1$, the equilibrium $E_1$ is locally asymptotically stable when $\tau = 0$. Further, if (H1) holds, the system undergoes Hopf bifurcation when $\tau = \tau _0^{(j)} \left(j = 0, 1, 2, \cdots\right)$. $E_1$ is locally asymptotically stable for $0 < \tau < \tau _0^{(0)}$ and is unstable for $\tau > \tau _0^{(0)}$, where $\tau_0^{(j)}$ is given by Eq (3.10).

2) If $\mathcal{R}_0 > 1$, the disease-free equilibrium $E_1$ is unstable.

3.2.   Analysis for the endemic equilibrium $E_2$

For the endemic equilibrium $E_2$, the characteristic equation of system (3.3) becomes:

where

When $\tau = 0$, the Eq (3.11) can be transformed into the following form:

where $q_1$, $q_2$, $q_3$, $q_4$, $p_1$, $p_2$, $p_3$ are shown in Eq (3.11), $a_1 = q_1+\gamma$, $a_2 = q_2+\gamma p_1$, $a_3 = q_3+\gamma p_2$, $a_4 = q_4+\gamma p_3$. According to the Routh-Hurwitz criterion, we show the following hypothesis:

If (H2) is satisfied, all eigenvalues of Eq (3.12) have negative real parts, the endemic equilibrium $E_2$ of system (2.1) is locally asymptotically stable when $\tau = 0$.

When $\tau > 0$, we try to discuss the existence of Hopf bifurcation. We assume that $\lambda = {\mathrm{i}}\omega \left(\omega > 0\right)$ is a pure imaginary root of Eq (3.11). Substitute it into Eq (3.11) and separate the real and imaginary parts, we can obtain:

Equation (3.13) derives to:

Adding the square of the two equations in Eq (3.13) and letting $z = \omega^2$, we get:

where

Discussion about the roots of Eq (3.15) is similar to that in ^[33]. According to this, we give the following lemma.

Lemma 3.3. Assume that $\delta_1 = \frac{c_2}{2}-\frac{3}{16}c_1^2$, $\delta_2 = \frac{c_1^3}{32}-\frac{c_1c_2}{8}+c_3$, $\Phi = (\frac{\delta_2}{2})^2+(\frac{\delta_1}{3})^3$, we can get the following results.

1) If $c_4 < 0$, Eq (3.15) has at least one positive root.

2) If $c_4\geq 0$, Eq (3.15) has positive roots if and only if $z_1 > 0$ and $h(z_1) < 0$ when $\Phi\geq0$. On the contrary, when $\Phi < 0$, Eq (3.15) has positive root if and only if there exists at least one $z^*\in \{z_1, z_2, z_3\}$, such that $z^* > 0$ and $h(z*)\leq0$, where $z_1$, $z_2$ and $z_3$ are the roots of $h'(z)$.

If $h(z)$ satisfies the existence condition for positive roots in Lemma 3.3, we hypothesize the Eq (3.15) has $n \left(n = 1, 2, 3, 4 \right)$ positive roots marked as $z_1 < z_2 < \dots < z_{n}$ $\left(n = 1, 2, 3, 4 \right)$. Substituting $\omega_{n} = \sqrt{z_{n}}$ into Eq (3.14), we can get the expression of $\tau$:

where $Q_{n} = {\sin(\omega\tau)}|_{\omega = {\omega _{n}}, \tau = \tau _{n}^{(j)}}, {P_{n}} = \cos(\omega\tau) |_{\omega = {\omega _{n}}, \tau = \tau _{n}^{(j)}},$ $\sin(\omega\tau)$, $\cos(\omega\tau)$ are given in Eq (3.14).

Let $\lambda\left(\tau \right) = a\left(\tau \right)+ {\mathrm{i}}\omega\left(\tau \right)$ be the root of Eq (3.11). If $\tau = \tau_{n}^{\left(j\right)}$, characteristic Eq (3.11) has a pair of imaginary roots: $\pm {\mathrm{i}} \omega_{n}$ with $a(\tau_{n}^{j}) = 0, \omega(\tau_{n}^{j}) = \omega_{n} \left(j = 0, 1, 2, \cdots, n = 1, 2, 3, 4\right)$.

Lemma 3.4. If (H2) holds and $z_{n} = \omega_{n}^{2}$, $h'\left(z_n \right)\neq0$, we have the transversality condition:

Theorem 3.2. When $\mathcal{R}_0 > 1$ and (H2) hold, the equilibrium $E_2$ is locally asymptotically stable for $\tau = 0$. When $\tau > 0$, we have the following conclusions.

1) If $h(z)$ does not satisfy the condition for the existence of positive root in Lemma 3.3, the equilibrium $E_2$ is locally asymptotically stable when $\tau > 0$.

2) If $h(z)$ satisfies the condition of the existence of positive root in Lemma 3.3 and only has one positive root $z_1$, the system (2.1) undergoes Hopf bifurcation at $E_2$ when $\tau = \tau_{1}^{(j)}$ and $h'(z_1) > 0$. The equilibrium $E_2$ is locally asymptotically stable for $\forall\; 0 < \tau < \tau_{1}^{(0)}$ and is unstable for $\forall$ $\tau > \tau_{1}^{(0)}$.

3) If $h(z)$ satisfies the condition of the existence of positive root in Lemma 3.3 and has two positive roots $z_1$, $z_2$, the system (2.1) undergoes Hopf bifurcation at $E_2$ when $\tau = \tau_{1}^{(j)}$ and $\tau = \tau_{2}^{(j)}\; (j = 0, 1, 2, \cdots)$.Assuming that $z_2 < z_1$, we can get $h'(z_1) > 0, h'(z_2) < 0$. Thus, if $\tau_{1}^{(0)} < \tau_{2}^{(0)}$, there exists $k$ that makes: $0 < \tau_1^{(0)} < \tau_{2}^{(0)} < \tau_{1}^{(1)} < \tau_{2}^{(1)} < \dots < \tau_{1}^{(k)} < \tau_{1}^{(k+1)}$. When $\tau\in\left[ 0, \left.\tau _{1}^{(0)} \right) \right.\cup \bigcup_{i = 1}^k{(\tau _{2}^{\left(i-1 \right)}, \tau _{1}^{\left(i \right)})}$, the equilibrium $E_2$ is locally asymptotically stable and when $\tau\in\bigcup_{i = 0}^{k-1}{(\tau_1^{(i)}, \tau_2^{(i)})}\cup (\tau_1^{(k)}, +\infty)$, the equilibrium $E_2$ is unstable.

4) If $h(z)$ satisfies the condition of the existence of positive root in Lemma 3.3 and has three positive roots $z_1$, $z_2$, $z_3$, the system (2.1) undergoes Hopf bifurcation at $E_2$ when $\tau = \tau_{1}^{(j)} (n = 1, 2, 3)$. Assuming that $z_3 < z_2 < z_1$, we have $h'(z_1) > 0, h'(z_2) < 0, h'(z_3) > 0$. Similar to the analysis of (3), the stability of equilibrium $E_2$ switches with the increase of $\tau$. And finally, the equilibrium $E_2$ is unstable.

5) If $h(z)$ satisfies the condition of the existence of positive root in Lemma 3.3 and has four positive roots $z_1$, $z_2$, $z_3$ and $z_4$, the system (2.1) undergoes Hopf bifurcation at $E_2$ when $\tau = \tau_{1}^{(j)} (n = 1, 2, 3, 4)$. Assuming that $z_4 < z_3 < z_2 < z_1$, we can obtain $h'(z_1) > 0, h'(z_2) < 0, h'(z_3) > 0, h'(z_4) < 0$. Similar to the analysis of (3), the stability of equilibrium $E_2$ switches with the increase of $\tau$. Finally, the equilibrium $E_2$ is unstable.

4.   Normal form and bifurcation analysis

In this section, we will deduce the normal form of Hopf bifurcation for system (3.3) with the help of multiple time scales method. Assuming that Eq (3.4) or (3.11) has a pair of pure imaginary roots $\lambda = \pm{\mathrm{i}}\omega_k$ at equilibrium $E_k = (S_k^*, A_k^*, I_k^*, R_k^*)\; (k = 1, 2)$, we can obtain the normal form as follows. The calculating details can be referred to Appendix.

where $\chi_k, M_k, H_k$ are given in Eqs (A11) and (A16) in Appendix.

Then, we let $G = \gamma {\text{e}}^{{\text{i}}\theta_k}$ and substitute it into Eq (A17). The normal form of Hopf bifurcation in polar coordinate is as follows:

The nontrivial equilibrium of Eq (4.2) corresponds to the periodic solution of system (2.1). So we can discuss the stability of periodic solution in system (2.1) by studying the stability of nontrivial equilibrium in Eq (4.2).

Theorem 4.1. If $\frac{\mathit{{\text{Re}}}\left(\chi _k M_k \right) \tau _{\varepsilon}}{\mathit{{\text{Re}}}\left(\chi _k H_k \right) \tau _c} < 0 \; (k = 1, 2)$ holds, the Eq (4.2) has nontrivial fixed point $r^* = \sqrt{-\frac{\mathit{{\text{Re}}}\left(\chi _k M_k \right) \tau _{\varepsilon}}{\mathit{{\text{Re}}}\left(\chi _k H_k \right) \tau _c}}$, so the system (2.1) has periodic solution around the equilibrium $E_k$.

1) If $\operatorname{Re}(\chi _k M_k){{\tau }_{\varepsilon }} < 0$, the periodic solution of system (2.1) is unstable. When $\tau_{\varepsilon} > 0\; (\tau_{\varepsilon} < 0)$, the unstable periodic solution of Hopf bifucation is forward (backward).

2) If $\operatorname{Re}(\chi _k M_k){{\tau }_{\varepsilon }} > 0$, the periodic solution of system (2.1) is stable. When $\tau_{\varepsilon} > 0\; (\tau_{\varepsilon} < 0)$, the stable periodic solution of Hopf bifurcation is forward (backward).

Remark 1:

In this paper, in order to determine the time for booster vaccination, we propose an $SAIR$ epidemic model with time delay for exploring the impact of booster vaccination time on the epidemic. We obtain the basic reproduction number $\mathcal{R}_0$ and the following conclusions.

1) $\mathcal{R}_0 < 1$: when $\tau = 0$, the disease-free equilibrium $E_1$ is locally asymptotically stable, which means the epidemic can be effectively controlled; when $\tau = \tau _{0}^{\left(j \right)}\; (j = 1, 2, \cdots)$, equilibrium $E_1$ undergoes Hopf bifurcation, where $\tau _{0}^{\left(j \right)}$ is critical value of Hopf bifurcation near $E_1$ given by Eq (3.10). In order to effectively control the epidemic, we need to control $\tau$ within $\tau_0^{(0)}$. That is to say, we should vaccinate the booster shot within the critical time while the original antibody is invalid in a sense. The booster vaccination can generate antibodies, which can improving herd immunity and further control the outbreak of the epidemic. $E_2$ is unstable in this condition.

2) $\mathcal{R}_0 > 1$: the disease-free equilibrium $E_1$ is unstable and when $\tau = 0$, $E_2$ is locally asymptotically stable under the preconditions of Theorem 3.2, the epidemic can be controlled effectively; when $\tau = \tau_n^{(j)}\; (j = 1, 2, \cdots, \; n = 1, 2, 3, 4)$, the equilibrium $E_2$ undergoes Hopf bifurcation if Eq (3.15) has positive root, where $\tau_n^{(j)}$ is critical value of Hopf bifurcation near $E_2$ given by Eq (3.16). In order to control the epidemic, we also need to control the booster vaccination time within the time delay that can make epidemic stability.

With analysis of equilibria, we can also get the direction of Hopf bifurcation by calculating the normal form according to Appendix. We can also obtain the stability of periodic solutions and the direction of Hopf bifurcation in Theorem 4.1.

5.   Numerical simulations

In this section, we will carry out numerical simulations to verify our theoretical analysis. Then, we will simulate the effect of virus mutation on the critical time of the booster vaccination by controlling the infection rates $\beta_1$ and $\beta_2$. Finally, we will estimate the time required to complete a round of vaccination worldwide through fitting and prediction, and will propose reasonable suggestions about the booster vaccination time by analysis of the results.

5.1.   Determination for parameter values

Based on official statistics (https://ourworldindata.org/covid-vaccinations; http://2019ncov.chinacdc.cn/2019-nCoV/global.html), we can obtain the data of cure rates and vaccination rates for 211 countries. In order to ensure that the data can reflect the average, we take countries with more than 400,000 infections and eliminate outliers and missing values. Then, we screen cure rates for 52 countries and vaccination rates for 116 countries. Since most COVID-19 vaccines are inactivated vaccines with two doses, we considered the vaccination rates of all two doses. As a result, we make bar chart and line chart respectively, as shown in Figure 2.

For Figure 2(a), we can clearly see that the cure rates of 52 countries are almost at the same level, so we choose the mean value: 0.9383 as the value of $u$. As for vaccination rate, it's not difficult to find that the vaccination rates of these countries are mostly in the range of 0.2 to 0.6 see in Figure 2(b), so we choose $c = 0.3089, c = 0.4126\in(0.3, 0.6)$ respectively.

Based on the above consideration and the actual situation, we take the following two groups of parameters:

1) $B = 10$, $d = 0.00714$, $c = 0.0484$, $\alpha = 0.3098$, $\gamma = 0.50$, $u = 0.9383$, $\beta_1 = 0.00046$, $\beta_2 = 0.00092$, $\rho_1 = 0.0025$, $\rho_2 = 0.74$, $m = 0.6$;

2) $B = 100$, $d = 0.00714$, $c = 0.0484$, $\alpha = 0.4126$, $\gamma = 0.35$, $u = 0.9383$, $\beta_1 = 0.00011$, $\beta_2 = 0.00089$, $\rho_1 = 0.1528$, $\rho_2 = 0.84$, $m = 0.6$.

5.2.   Simulation and verification

For the group of parameters (1): $B = 10$, $d = 0.00714$, $c = 0.0484$, $\alpha = 0.3098$, $\gamma = 0.50$, $u = 0.9383$, $\beta_1 = 0.00046$, $\beta_2 = 0.00092, \rho_1 = 0.0025$, $\rho_2 = 0.74$, $m = 0.6$, we calculate the non-negative equilibrium $E_1 = [S_1^*, A_1^*, I_1^*, R_1^*] = [869, 0, 0,531]$ according to Eq (3.2). This group of parameters meets (H1), so $\mathcal{R}_0 < 1$, the equilibrium $E_1$ is locally asymptotically stable when $\tau = 0$. Substitute parameters (1) into Eqs (3.8)–(3.10), we get $\omega_0 = 0.3867, \tau_0^{(0)} = 5.8372, \sin(\omega_0\tau) = 0.7734, \cos(\omega_0\tau) = -0.6339$. According to Theorem 3.1, we know that $E_1$ is locally asymptotically stable for any $0 < \tau < \tau _0^{(0)} = 5.8372$ and it is unstable for any $\tau > \tau _0^{(0)} = 5.8372$. Then we calculate the normal form of Hopf bifurcation and obtain $\operatorname{Re}(\chi _k M_k) > 0, \operatorname{Re}(\chi _k H_k) > 0$ from Eqs (A11) and (A16). According to Theorem 4.1, we can get when $\tau_\varepsilon < 0$, $\operatorname{Re}(\chi _k M_k){{\tau }_{\varepsilon }} < 0$, the periodic solution is unstable and the direction of Hopf bifurcation is backward.

When $\tau = 0$, it means people consistently vaccinate to avoid infection. In this case, the virus constantly mutates, making the original antibody unable to defend against the mutated strains, so the immune system continues to fail. We choose the initial values [1000,100, 50,500], the equilibrium is locally asymptotically stable shown in Figure 3. Although there are fluctuations in 0–4 years, it tends to be stable after the 4th year. There are no infected people and the disease is eliminated. This means that in this case, it is efficient for people to vaccinate against disease consistently, but this situation is under the perfect condition. Considering that it takes a certain time to complete a round of universal vaccination and the antibody are not instantaneous failure, this situation is basically impossible.

When $\tau = 1.5\in(0, \tau _0^{(0)})$ with $\tau_0^{(0)} = 5.8372$, which means that we will get the booster vaccination 1.5 years later, we still choose the initial values [1000,100, 50,500] and the figure of numerical simulation is shown in Figure 4.

From Figure 4, we find that the disease-free equilibrium $E_1$ is also locally asymptotically stable. The fluctuation in the first eight years is obvious, but it gradually stabilizes after 6th years, and the disease can also be eliminated. Although its regional stability speed is less than that $\tau = 0$, it can still ensure stability in the short time. That is to say, when we do the booster vaccination 1.5 years later, the immune effectiveness of the original antibody is lost in a sense, and people can vaccinate again to deal with the constantly changeable virus. This will enhance the immunity of the population and control the epidemic effectively.

Through the analysis of Hopf bifurcation direction, we obtain: when $\tau = \tau_0^{(0)} = 5.8372$, $\operatorname{Re}(\chi _k M_k) > 0, \operatorname{Re}(\chi _k H_k) > 0, \tau_\varepsilon < 0$, so there is backward Hopf bifurcation near $\tau_0^{(0)}$, and the bifurcating periodic solution is unstable. In order to control the epidemic effectively, we hope to avoid this situation.

Remark 2: According to the numerical simulations of parameters (1), we can get: when $\tau < \tau_0^{(0)}$, the shorter time of booster vaccination is, the better the epidemic can be controlled. However, if $\tau > \tau_0^{(0)}$, the equilibrium $E_2$ is unstable. In other words, if people are vaccinated after more than $\tau_0^{(0)}$ years, the epidemic will be difficult to control.

For the group of parameters (2): $B = 100$, $d = 0.00714$, $c = 0.0484$, $\alpha = 0.4126$, $\gamma = 0.35$, $u = 0.9383$, $\beta_1 = 0.00011$, $\beta_2 = 0.00089$, $\rho_1 = 0.1528$, $\rho_2 = 0.84$, $m = 0.6$, we can also calculate the equilibria $E_1 = [S_1^*, A_1^*, I_1^*, R_1^*] = [6498, 0, 0, 7507],$ $E_2 = [S_2^*, A_2^*, I_2^*, R_2^*]$ $= [4051,307, 47, 9279]$ by Eq (3.2). We find that (H1) is not satisfied and (H2) is available. So the equilibrium $E_1$ is unstable and $E_2$ is locally asymptotically stable when $\tau = 0$. Substitute this group of parameters into Eqs (3.11)–(3.15), we obtain that $h(z)$ has two different positive real roots $z_1 \approx 0.5658 > z_2 \approx 0.2869$. According to Eq (3.16), we have $\tau_1^{(0)} = 3.1114 < \tau_2^{(0)} = 7.2011 < \tau_1^{(1)} = 11.4645 < \tau_1^{(2)} = 18.9323 < \tau_2^{(1)} = 19.8176.$

In accordance with Theorem 3.2, the system (2.1) has the conclusion: if $\tau\in(0, \tau_1^{(0)})\bigcup(\tau_2^{(0)}, \tau_1^{(1)})$, the equilibrium $E_2$ is locally asymptotically stable, and if $\tau\in(\tau_1^{(0)}, \tau_2^{(0)})\bigcup(\tau_1^{(1)},$ $+\infty)$, $E_2$ is unstable. Besides, the system undergoes Hopf bifurcation near $E_2$. Through Eqs (A10)–(A17) and Theorem 3.2, we get: when $\tau = \tau_1^{(0)}$, $\operatorname{Re}(\chi_2M_2) > 0$, $\operatorname{Re}(\chi_2H_2) < 0$, $\tau_\varepsilon > 0$, the periodic solution is stable; when $\tau = \tau_2^{(0)}$, $\operatorname{Re}(\chi_2M_2) < 0$, $\operatorname{Re}(\chi_2H_2) < 0$, $\tau_\varepsilon < 0$, the periodic solution is stable; and when $\tau = \tau_1^{(1)}$, $\operatorname{Re}(\chi_2M_2) > 0$, $\operatorname{Re}(\chi_2H_2) < 0$, $\tau_\varepsilon > 0$, the periodic solution is also stable.

When $\tau = 0$ which means people consistently vaccinate to avoid infection with instantaneous failure of immunity, we can see the equilibrium $E_2$ is asymptotically stable with the initial value [5000,300,100, 10000] shown in Figure 5. In this figure, there are fluctuations in 0–6 years, but it tends to be stable after the 6th year. Similar to the situation of parameters (1), although the epidemic can be effectively controlled in a short time, this situation does not exist.

When $\tau = 1\in(0, \tau_1^{(0)})$ where $\tau_1^{(0)} = 3.1114$, which means that people get booster vaccination after one year along with failure of original immunity, we can get the equilibrium $E_2$ is also asymptotically stable shown in Figure 6. Although there are fluctuations in the epidemic in the previous 8 years, it tends to be stable after 10th year.

When $\tau\in(\tau_1^{(0)}, \tau_2^{(0)}) = (3.1114, 7.2011)$, which means that people get booster vaccination between 3.1114 and 7.2011 years, the equilibrium $E_2$ is unstable. Since $\operatorname{Re}(\chi_2M_2) > 0, \operatorname{Re}(\chi_2H_2) < 0, \tau_\varepsilon > 0$ at $\tau = \tau_1^{(0)}$, system (2.1) has forward Hopf bifurcation near $\tau_1^{(0)}$ and the periodic solution is stable. This is not ideal for controlling the outbreak.

When $\tau = 8 \in(\tau_2^{(0)}, \tau_1^{(1)})$, which means that people get booster vaccination after 8 years and the original immunity fails at this time, the results of the numerical simulation are presented in Figure 7 with the initial value [5000,300,100, 10000]. Although the equilibrium is asymptotically stable in the long term, it continues to fluctuate in the next two decades and the epidemic can not be effectively controlled in a short time. Therefore, it will not be adopted.

Remark 3: Through numerical simulations, it can be found that as time delay $\tau$ changes, the stability of the equilibrium $E_2$ constantly switches. However, in the short term, the epidemic can be controlled effectively only when $\tau < \tau_1^{(0)}$. And when $\tau > \tau_1^{(0)}$, the outbreak of epidemic will continue. Because the small variation of time delay will lead to great changes near $\tau_1^{(0)}$, therefore, much emphasis should be put on controlling the booster vaccination time.

Next, we will simulate the virus mutation based on the second group of parameters. To keep the epidemic under control, the booster vaccination time must be kept within $\tau_1^{(0)}$. Otherwise, antibody will fail, the epidemic will be difficult to control.

As the virus mutates and selective expression of genes, the infection rate and the drug resistance increases, so we consider that the infection rates $\beta_1$ and $\beta_2$ change within a certain range to simulate the virus mutation. And we use the critical time delay $\tau_1^{(0)}$ to reflect the efficiency of vaccine. We find that when $\beta_1\in(0.00009, 0.00014)$, $\beta_2\in(0.00053, 0.00185)$ and other parameters remain unchanged, the stability and the direction of Hopf bifurcation is the same as situation of parameter (2). Therefore, we choose $\beta_1 = 0.00009, 0.00010, 0.00011, 0.00012$ and $\beta_2$ is changed within (0.00053, 0.00185), and we obtain a graph about $\tau_1^{(0)}$ with respect to the infection rates.

According to Figure 8, we can clearly see that $\tau_1^{(0)}$ decreases as $\beta_2$ increases when $\beta_1$ is fixed, and the delay $\tau_1^{(0)}$ decreases for the four decreasing values of $\beta_1$ when $\beta_2$ is fixed. In fact, as the virus continues to mutate and the rate of infection increases, the greater the probability that the original antibody will fail, requiring us to take action for booster vaccination in a shorter period of time. This is consistent with the results of our simulations. Therefore, in order to keep the outbreak under control and rather than the sustained outbreak, we need to reduce the booster vaccination time when the virus mutates.

5.3.   Forecast and suggestion

It takes a certain amount of time for a region to complete a round of vaccination which cannot be ignored in the study of the booster vaccination time. So we should consider this in our model.

Now, we use the available data to predict the time to complete a round of vaccination. We obtain the time-varying data of the number of people vaccinated per 100 people in each country from the official data (https://ourworldindata.org/grapher/covid-vaccination-doses-per-capita?time=latest). And then we select the average as the level of the world through the data of 2020/12/1–2021/10/12. We choose the initial time as $t = 0$ with day as a unit, perform fitting with logistic function on the data and get the equation as:

The fitting and prediction result is shown in Figure 9. Through calculation, it is concluded the correlation coefficient between predicted value and actual value $R^2 = 0.9803$, with strong correlation and good fitting effect.

In fact, due to individual variability, not all people can be vaccinated, especially for elderly people with some underlying diseases. So it is not possible to achieve the vaccination rate of 100%, which is consistent with the results of our fitting. In the Figure 9, we find that the growth of vaccination is already slow when $t$ is at 450 days, so we consider that vaccination is complete at $t = 450$ days, which means that a round of vaccination can be completed in 1.23 years with the world average.

From the simulation results in Subsection 5.2, we can obtain that under the second group of parameters, the booster vaccination can eventually reach stability within 3.11 years. Based on the fitting result in this subsection, the time to complete a round of vaccination is predicted to 1.23 year. Therefore, when the vaccination process starts within 3.11-1.23 = 1.88 years, it can ensure that the round of vaccination can be completed and the final stabilization can be guaranteed.

Next we illustrate the more general relationship between fitting results and system (2.1). Firstly, we give the definitions of some function and region in Table 2.

The virus mutation affects the infection rate $\beta_2$ most predominantly, so we define the booster vaccination time function as $t_1 = \tau_1^{(0)}(\beta_2)$ expressed as the relationship between the critical booster vaccination time and the infection rate $\beta_2$. However, since vaccination is a large-scale behavior, the time to complete a round of vaccination is not negligible, assuming $t_0$. Taking this factor into account in the system, we define the restricted booster vaccination time function as $t_2 = \tau_1^{(0)}(\beta_2)-t_0$, expressed as the relationship between the critical booster vaccination time that allows vaccination to proceed smoothly and the infection rate $\beta_2$. With the above two functions, we can divide the $\beta_2- t$ plane into three parts: $D_1 = \{(\beta_2, t)|t > \tau_1^{(0)}(\beta_2)\}$, $D_2 = \{(\beta_2, t)|\tau_1^{(0)}(\beta_2)-t_0 < t < \tau_1^{(0)}(\beta_2)\}$ and $D_3 = \{(\beta_2, t)|0 < t < \tau_1^{(0)}(\beta_2)-t_0\}$, as shown in Figure 10(a). $D_1$ describes the region where the critical booster vaccination time beyond the critical booster vaccination time, which is unstable according to the previous simulations. $D_3$ describes the region where the booster vaccination time is less than the critical booster vaccination time after removing the time to complete a round of vaccination, which is clearly the stable state. The region between the two functions $D_1$, $D_2$ is not able to complete the vaccination before reaching instability with the current vaccination process, but if we intervene appropriately to speed up the vaccination process, we can complete the vaccination on time and eventually reach the stable state. So we call this region as buffer region.

Taking the second set of parameters and the results of the fit and prediction, we get $\beta_2 = 0.00089$, $t_0 = 1.23$, and we draw the function $t_1$, $t_2$ and the region $D_1$–$D_3$, as shown in Figure 10(b). Nowadays, most countries in the world require booster vaccination at 6 months, so we mark the point (0.00089, 0.5) in the figure as the current state. It is obvious that this point is in region $D_3$, so the current state is the stable state. However, as the virus mutates, the transmission rate $\beta_2$ will gradually increase, which we represent by the red line with arrow in Figure 10(b). We can clearly see that the state point reaches the region $D_2$ when $\beta_2 > 0.0073$, at which time we think we cannot complete a round of vaccination. However, we can change the size of the region $D_1$ and $D_2$ by shortening the vaccination process, so that the state point reaches again within the new stable region $D_3$.

Remark 4: For the numerical simulation results, we have the following conclusions.

1) Through the simulation of the two groups, we find that the epidemic would be effectively controlled when the booster vaccination time is controlled within the critical delay $\tau_0^{(0)}$ or $\tau_1^{(0)}$, while the epidemic would be difficult to be controlled if the booster vaccination time was larger than the critical delay.

2) For the mutant strain, the infection rate will increase with the continuous mutation of the virus, which will lead to the continuous reduction of the critical booster vaccination time. That is detrimental to the control of the epidemic.

3) Not all countries are able to complete a round of mass vaccination within a short period of time, taking human, material and financial factors into account. For the current world average, we should keep the vaccination time within 1.88 years. Therefore, it is reasonable to carry out booster vaccination at 6 months, because it can keep the epidemic at a stable state. But as the virus mutates and the infection rate increases, we may reach a state within the buffer zone, at which point the vaccination process needs to be accelerated so that the state is in the new stable zone.

6.   Conclusions

In this paper, considering of the characteristics of COVID-19 vaccine, we have constructed an $SAIR$ model with time delay to study the booster vaccination time. We have studied the stability of the equilibria and the existence of Hopf bifurcation with the help of the basic reproduction number $\mathcal{R}_0$. Then we have analysed the stability and direction of Hopf bifurcation periodic solution by calculating the normal form with the multiple time scales method.

In Section 5, we carried out numerical simulations. Firstly, the changes of epidemic situation are simulated through two groups of parameters to verify the theoretical analysis results. Secondly, we considered the impact of mutant strains on the critical time of booster vaccination. With the continuous variation of virus and the continuous increase of infection rate, the critical time of booster vaccination should gradually be decreased to ensure the safety and effectiveness of epidemic prevention and control. Then, using the official vaccination data, we fitted and predicted the time required to complete a round of universal vaccination. Through the analysis, we define three regions: stable region, unstable region, and buffer region. Based on this, we judge that the current 6-month vaccination is reasonable, but with the mutation of the virus, it is necessary to accelerate the booster vaccination process or advance the booster vaccination time appropriately.

In the future, as the virus continues to mutate, it is possible that the virus will have a high rate of transmission but will gradually become less aggressive to humans, and COVID-19 may become a disease like influenza. What kind of dynamic phenomenon will be generated at that time? This is an issue for future research to explore.

Acknowledgments

This study was funded by Fundamental Research Funds for the Central Universities of China (Grant No. 2572022DJ06).

Conflict of interest

The authors declare that they have no competing interests.

Appendix

System (3.3) can be written as:

where

Then, make transformation: $t\mapsto t/\tau$ for convenience, we have:

where $A$, $B$, $F$ are given in Eq (A1).

We linearize system (A2): $\widetilde X'(t) = \tau AX(t) + \tau BX(t - 1)$. Let $h_{k}\; (k = 1, 2)$ be the eigenvector of the linear operator corresponding to the eigenvalue ${\text{i}}\omega_{k}\tau$ and let $h_{k}^{*} \; (k = 1, 2)$ be the normalized eigenvector of the adjoint operator of the liner operator corresponding to the eigenvalues $-{\text{i}}\omega_{k}\tau$ satisfying the inner product $\left < h_{k}, h_{k}^{*}\right > = {\overline {h_{k}^{*}} ^{\text{T}}}h_{k} = 1$. By a simple calculation, we get:

where

We consider $\tau$ as a bifurcation parameter and take a time scale $\tau = \tau_c+\varepsilon \tau_\varepsilon$, where $\tau_c = \tau_n^{(j)} (n = 0, 1, 2, 3, 4)$ in Eq (3.11) or (3.16). $\varepsilon$ is the scale parameter without dimension and $\tau_\varepsilon$ is the disturbance parameter. Thus, $X(t)$ can be written as:

where

$X'(t)$ can be written as:

where $D_i = \frac{\partial}{\partial T_i} (i = 1, 2, 3, \cdots)$ is differential operator.

Using the Taylor expansion, we can obtain the formula of $X(t-1)$:

where $X_{i1} = X_i(t-1, \varepsilon t, \varepsilon^2t, \cdots)\; (i = 1, 2, \cdots).$

Then we substitute Eqs (A4)–(A6) into Eq (A1). Corresponding to the coefficient of $\varepsilon$ on both sides of the equation, we get the following expression:

Corresponding to the coefficient of $\varepsilon^2$ on both sides of the equation, we obtain:

For $\varepsilon^3-$order terms, we have:

The solution of Eq (A7) can be expressed as the following form:

where $h_k$ is given in (A3). Substitute Eq (A10) into the right hand side of Eq (A8), and mark the coefficient before ${\text{e}}^{{\text{i}}\omega\tau_cT_0}$ as vector $m_1$. Through solvability condition $< h_k^*, m_1 > = 0$, we can get the expression of $\frac{\partial G}{\partial T_1}$:

where $\chi_k = \frac{1}{1+\gamma\tau_c(h_4\bar h_1^*-h_4\bar h_4^*){\text{e}}^{-{\text{i}}\omega_k\tau_c}}, M_k = {\text{i}}\omega_k$.

Assume the solutions of (A8) are as follows:

Then substitute Eq (A12) into Eq (A8), we have:

where

We can solve $f_1, f_2, f_3, f_4, g_1, g_2, g_3, g_4, l_1, l_2, l_3, l_4$ from Eqs (A13)–(A15). Thus, we can get the expressions of $S_{2, k}, A_{2, k}, I_{2, k}, R_{2, k}$ according to Eq (A12). Then we substitute Eqs (A10)–(A12) into the right expression of Eq (A9), and note the coefficient of ${\text{e}}^{{\text{i}}\omega_k\tau_cT_0}$ as vector $m_2$. On the basis of solvability condition $< h_k^*, m_2 > = 0$, we obtain the expression of $\frac{\partial G}{\partial T_2}$. Because $\tau_\varepsilon^2$ has less impact on normal form, so we can ignore $\tau_\varepsilon^2G$ term. Therefore, we have:

where

Then, we let $G\rightarrow G/\varepsilon$. Therefore, we get the normal form of Hopf bifurcation for system (2.1):

where $\chi_k, M_k, H_k$ are given in Eqs (A11) and (A16).