Fuhong Zhang¹, Xiaolan Liu^2
1College of Mathematics and Statistics, Sichuan University of Science and Engineering, Zigong, China.
²Key Laboratory of Higher Education of Sichuan Province for Enterprise Informationalization and Internet of Things, Zigong, China.
DOI: 10.4236/jamp.2026.149169   PDF    HTML   XML   2 Downloads   39 Views   Citations

Abstract

If a natural disaster like the epidemic breaks out again, we all know that the consequences will be unimaginable. In order to describe the transmission process of infectious diseases more vividly, this paper constructs a new fractional-order SEIR model, and introduces two time-delay parameters: incubation period and recovery period. For this system, we calculate two equilibrium points, and focus on the stability of the local equilibrium point affected by time delay and hopf bifurcation. Finally, we calculate the critical time delay through the numerical calculator maple, and verify the local asymptotic stability at the equilibrium point and the result of hopf bifurcation through numerical simulation.

Keywords

Hopf Bifurcation, Fractional-Order Delayed SEIR, Stability

Share and Cite:

1. Introduction

Since the outbreak of coronavirus disease 2019 (COVID-19), the global public health system has faced very serious challenges for a long time. Moreover, the widespread spread of the virus will not only cause serious problems of human health and life safety, but also have an irreparable impact on global economic development and social stability. According to the mathematical modeling method, it is the core direction with important practical value in the study of infectious disease dynamics to deeply explore the dynamic law of the spread of novel coronavirus, accurately predict the evolution trend of the epidemic situation, and design scientific and effective prevention and control intervention programs accordingly. This paper systematically sorts out and summarizes the research results in related fields at home and abroad. The existing academic research provides a solid theoretical basis and research reference for the study of Hopf bifurcation of fractional-order multi-delay COVID-19 propagation model from two aspects: model system construction and dynamic analysis method. In the related research of biodynamic models, especially infectious disease transmission models, many scholars have achieved rich results. Among them, Reference [1] constructed an integer-order SEIR dynamic model with dual delay parameters of incubation period and recovery period based on the real transmission characteristics of the novel coronavirus epidemic, and focused on the mechanism of time delay on the stability of the equilibrium point and the Hopf bifurcation dynamic behavior of the system. This study further integrates control strategies and drug intervention-related variables, constructs an optimal control analysis framework with time-delay characteristics, and establishes the time-delay disturbance problem as the core research context, which provides an important reference for the research ideas of this paper. Based on the traditional SEIR model, Reference [2] extended the model into a spatial reaction-diffusion system. At the same time, it also introduced many realistic factors such as the randomness of population spatial migration and nonlinear infection rate, which enriched the research ideas of infectious disease modeling coupled with multi-dimensional realistic variables, and provided a new reference direction for the optimization and improvement of the model in this paper.

Fractional differential equations are widely used in epidemic modeling because they can describe the memory effect and historical cumulative dependence of the system. In the literature [3]-[9], etc., combined with the actual epidemic characteristics of asymptomatic infection, multiple physiological delays, saturated incidence, multi-group warehouse division, etc., different forms of fractional-order delay SEIR-like dynamic systems are constructed. It is fully confirmed that fractional-order operators combined with discrete delays can be more suitable for the complex evolution of infectious diseases. In addition, Reference [10] explored the bifurcation characteristics of the double-delay SEIR model under the dual propagation path. Reference [11] carried out the quantitative calculation of Hopf bifurcation based on the combination of saturation incidence and time delay. Reference [12] completed the threshold determination and global stability proof of the classical SEIR model. Reference [13] introduced the complex network topology to describe the non-uniform contact mode of the crowd, and improved the extended form of the traditional SEIR model from multiple angles. In [14]-[16], the optimal control problem of vaccination with state constraints is added to the SEIR model, and the optimal intervention scheme is solved by distinguishing different constraints, which provides an important reference for this paper to build a time-delay optimal control model with dual control variables of social contact control and drug treatment.

In terms of theoretical analysis and numerical simulation methodology, the existing research has formed a complete and mature technical system. In many literatures, the fixed point theory is widely used to prove the existence and uniqueness of the solution of the fractional order system. The next generation matrix method is used to solve the basic reproduction number, and the Lyapunov direct method is used to complete the local and global stability of the equilibrium point. The stability switching condition is determined by the root distribution of the characteristic equation. Based on the center manifold theorem and the normal form method, the direction, period and stability of Hopf bifurcation are derived. In the numerical simulation stage, Adams-Bashforth-Moulton predictor-corrector algorithm, finite difference method, network topology simulation and other tools are widely used to verify the theoretical derivation results. The above complete qualitative analysis theory and numerical calculation scheme lay a methodological foundation for the equilibrium point analysis, Hopf bifurcation condition derivation and numerical verification of the optimal control strategy of the fractional-order multi-delay new crown propagation model.

2. Statement of the Problem

We propose the following delayed SEIR model to describe the spread of the COVID-19 virus.

$\{\begin{array}{l}{}_0^C{D}_t^{q}S\left(t\right)=\Lambda -\frac{{\beta }_{1}S\left(t\right)E\left(t\right)+{\beta }_{2}S\left(t\right)I\left(t\right)}{N\left(t\right)}-dS\left(t\right)+\delta R\left(t\right),\\ {}_0^C{D}_t^{q}E\left(t\right)=\frac{{\beta }_{1}S\left(t\right)E\left(t\right)+{\beta }_{2}S\left(t\right)I\left(t\right)}{N\left(t\right)}-\epsilon E\left(t-{\tau }_1\right)-\left({\gamma }_2+d\right)E\left(t\right),\\ {}_0^C{D}_t^{q}I\left(t\right)=\epsilon E\left(t-{\tau }_1\right)-{\gamma }_{1}I\left(t-{\tau }_2\right)-\left(\mu +d\right)I\left(t\right),\\ {}_0^C{D}_t^{q}R\left(t\right)={\gamma }_{1}I\left(t-{\tau }_2\right)+{\gamma }_{2}E\left(t\right)-dR\left(t\right)-\delta R\left(t\right).\end{array}$ (1)

Here, the fractional order parameter $q\in \left(0,1\right]$ , $S\left(t\right)$ , $E\left(t\right)$ , $I\left(t\right)$ , and $R\left(t\right)$ denote the susceptible, exposed, infected, and recovered at time $t$ , respectively, and $N\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)R\left(t\right)$ denotes the total population tree. ${\tau }_1$ , ${\tau }_2$ represents the time delay parameters of two different physiological processes, where ${\tau }_1$ represents the latent period delay after infection, and ${\tau }_2$ represents the recovery period delay of infected individuals. The biological interpretation of the remaining parameters in the system Equation (1) is summarized in Table 1.

Table 1. Parameter descriptions of the fractional-order delayed SEIR epidemic model.

Symbol	Description
$\Lambda$	Constant population recruitment rate
$d$	Natural mortality rate
$\mu$	Disease-induced excess mortality rate of infected individuals
${\beta }_1$	Effective transmission rate of exposed individuals $E$
${\beta }_2$	Effective transmission rate of symptomatic infected individuals $I$
$\epsilon$	Conversion rate from exposed compartment $E$ to infected compartment $I$
${\gamma }_1$	Recovery rate of infected individuals $I$
${\gamma }_2$	Direct self-healing rate of exposed individuals $E$
$\delta$	Immunity waning rate of recovered individuals

In order to facilitate the subsequent stability analysis and Hopf bifurcation study of system (1), the definition of Caputo fractional derivative and several basic lemmas are given.

Definition 2.1 (Caputo Fractional derivative [17]-[19]). The Caputo definition of fractional derivatives can be written as:

${}_a^C{D}_t^{\alpha }f\left(t\right)=\frac{1}{\Gamma \left(n-\alpha \right)}{\displaystyle {\int }_a^t\frac{f^{\left(n\right)}\left(\tau \right)}{{(t-\tau )}^{\alpha -n+1}}\text{d}\tau },n-1<\alpha \le n,$ (2)

where $f\left(t\right)\in C^m\left(\left[a,\infty \right),ℝ\right)$ . In particular, if $0<\alpha \le 1$ , $a=0$ , Equation (2) can be written as:

${}_0^C{D}_t^{\alpha }f\left(t\right)=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^t\frac{f^{\prime }\left(\tau \right)}{{(t-\tau )}^{\alpha }}\text{d}\tau },0<\alpha \le 1,t>0.$ (3)

Definition 2.2 (Laplace Transform [20] [21] of Fractional Derivative). The Laplace transform of Caputo fractional derivative of order $\alpha$ ($n-1<\alpha \le n$ ) for the function $f\left(t\right)\in C^n\left(\left[a,\infty \right),ℝ\right)$ is

$ℒ\left\{{}_a^C{D}_t^{\alpha }f\left(t\right);s\right\}=s^{\alpha }F\left(s\right)-{\displaystyle \sum }_{k=0}^{n-1}{s}^{\alpha -k-1}{f}^{\left(k\right)}\left(a\right),$ (4)

where $F\left(s\right)$ is the Laplace transform of $f\left(t\right)$ , and $f^{\left(k\right)}\left(a\right)$ $\left(k=0,1,\cdots ,n-1\right)$ are the initial conditions. Obviously, if $f^{\left(k\right)}\left(a\right)=0$ for $k=0,1,\cdots ,n-1$ , Equation (4) can be written as

$ℒ\left\{{}_a^C{D}_t^{\alpha }f\left(t\right);s\right\}=s^{\alpha }F\left(s\right).$

Definition 2.3. [22] Consider the following $n$ -dimensional fractional-order system with time delay:

${}_a^C{D}_t^{\alpha }{u}_i\left(t\right)=f_i\left(u_1\left(t\right),\cdots ,u_n\left(t\right);\tau \right),i=1,2,\cdots ,n,$ (5)

where $0<\alpha \le 1$ and the time delay $\tau \ge 0$ . System (5) undergoes Hopf bifurcation at the equilibrium $u^{\text{*}}=\left(u_1^{\text{*}},u_2^{\text{*}},\cdots ,u_n^{\text{*}}\right)$ when $\tau ={\tau }_0$ if the following three conditions are satisfied:

i) All the eigenvalues ${\lambda }_j$ ($j=1,2,\cdots ,n$ ) of the coefficient matrix $J$ of the linearized system of (5) with $\tau =0$ satisfy

$\left|\text{arg}\left({\lambda }_j\right)\right|>\frac{\alpha \pi }{2}.$

ii) The characteristic equation of the linearized system of (5) has a pair of purely imaginary roots $\pm i{\omega }_0$ when $\tau ={\tau }_0$ .

iii)

$\text{Re}\left[{\frac{\text{d}s\left(\tau \right)}{\text{d}\tau }|}_{\tau ={\tau }_0,\omega ={\omega }_0}\right]\ne 0,$

where $\text{Re}[\cdot ]$ denotes the real part of the complex number.

Using the next-generation matrix method proposed by van den Driessche and Watmough [23], the basic reproduction number of the system (1) can be expressed as:

$R_0=\frac{{\beta }_1\left({\gamma }_1+\mu +d\right)+{\beta }_2\epsilon }{\left(\epsilon +{\gamma }_2+d\right)\left({\gamma }_1+\mu +d\right)}$

For system (1), the disease-free and the endemic equilibrium points are defined as

$E_0=\left(\Lambda /d,0,0,0\right)$ and $E_1=\left(S^{\text{*}},E^{\text{*}},I^{\text{*}},R^{\text{*}}\right)$

where

$\{\begin{array}{l}{S}^{*}=\frac{\Lambda \left[{\epsilon }^2{\gamma }_1+\epsilon K\left(d+{\gamma }_1\right)+K\left(d+{\gamma }_2\right)M+\epsilon {\gamma }_{2}M\right]}{K\left(d+\epsilon \right)\left[{\beta }_{1}d+{\beta }_2\epsilon +{\beta }_1{\gamma }_1+{\beta }_1\mu -\epsilon \mu \right]},\\ E^{*}=\frac{\Lambda \left[\left({\beta }_{1}M+{\beta }_2\epsilon -CM\right)+\epsilon \left({\beta }_1-C\right)\right]}{CK\left[{\beta }_{1}M+\epsilon \left({\beta }_2-\mu \right)\right]},\\ I^{*}=\epsilon \frac{\Lambda \left[\left({\beta }_{1}M+{\beta }_2\epsilon -CM\right)+\epsilon \left({\beta }_1-C\right)\right]}{CKM\left[{\beta }_{1}M+\epsilon \left({\beta }_2-\mu \right)\right]},\\ R^{*}=\frac{\left({\gamma }_1\epsilon +{\gamma }_{2}M\right)\Lambda \left[C{\beta }_{1}M+{\beta }_2\epsilon -CM\right]+\epsilon \left(C{\beta }_1-C\right)}{K^{2}CM\left[{\beta }_{1}M+\epsilon \left({\beta }_1-\mu \right)\right]},\end{array}$ (6)

where $M={\gamma }_1+\mu +d$ , $K=d+\delta$ , $C=\epsilon +{\gamma }_2+d$

We know that when the basic reproduction number $R_0<1$ , only the disease-free equilibrium exists and the disease dies out; $R_0=1$ corresponds to the critical bifurcation case; when $R_0>1$ , an endemic equilibrium exists and the disease persists.

For Caputo fractional-order systems with $0<q<1$ , the equilibria and the Jacobian matrix evaluated at each equilibrium are completely determined by the right-hand-side vector field of the system and are independent of the fractional order $q$ . Accordingly, the eigenvalues of the Jacobian matrix can be solved by exactly the same algebraic procedures as those used for integer-order systems. Nevertheless, local asymptotic stability of an equilibrium for fractional-order systems is no longer equivalent to the condition that all eigenvalues have negative real parts. Instead, the fractional-order eigenvalue criterion should be adopted: an equilibrium is locally asymptotically stable if and only if every eigenvalue $\lambda$ of the Jacobian matrix satisfies

$\left|\text{arg}\left(\lambda \right)\right|>\frac{q\pi }{2}.$

For the delay-free counterpart of our system, all eigenvalues of the Jacobian at the disease-free equilibrium satisfy $\left|\text{arg}\left(\lambda \right)\right|>\frac{\pi }{2}$ whenever $R_0<1$ . Combined with $0<q<1$ , we have $\frac{q\pi }{2}<\frac{\pi }{2}$ , so the fractional-order stability condition

$\left|\text{arg}\left(\lambda \right)\right|>\frac{q\pi }{2}$ is automatically fulfilled. Hence, $R_0<1$ guarantees the local

asymptotic stability of the disease-free equilibrium in the delay-free case. It is worth noting that this property holds only in the absence of time delays. When time delays ${\tau }_1,{\tau }_2$ are involved, the characteristic equation becomes transcendental, and the eigenvalue analysis based on the constant Jacobian matrix is invalid. We have to investigate the stability switches and Hopf bifurcation induced by time delays.

At the equilibrium point $E_1=\left(S^{\text{*}},E^{\text{*}},I^{\text{*}},R^{\text{*}}\right)$ , by applying the transformations

$\widehat{S}\left(t\right)=S\left(t\right)-S^{\text{*}},\widehat{E}\left(t\right)=E\left(t\right)-E^{\text{*}},\widehat{I}\left(t\right)=I\left(t\right)-I^{\text{*}},\widehat{R}\left(t\right)=R\left(t\right)-R^{\text{*}},$

system (1) can be converted into

$\{\begin{array}{l}{D}_q\widehat{S}\left(t\right)=l_{11}\widehat{S}+l_{12}\widehat{E}+l_{13}\widehat{I}+l_{14}\widehat{R},\\ D_q\widehat{E}\left(t\right)=l_{21}\widehat{S}+l_{22}\widehat{E}+l_{23}\widehat{I}+l_{24}\widehat{R}+l_{25}\widehat{E}\left(t-{\tau }_1\right),\\ D_q\widehat{I}\left(t\right)=l_{31}\widehat{S}+l_{32}\widehat{E}+l_{33}\widehat{I}+l_{34}\widehat{R}+l_{35}\widehat{E}\left(t-{\tau }_1\right)+l_{36}\widehat{I}\left(t-{\tau }_2\right),\\ D_q\widehat{R}\left(t\right)=l_{41}\widehat{S}+l_{42}\widehat{E}+l_{43}\widehat{I}+l_{44}\widehat{R}+l_{45}\widehat{I}\left(t-{\tau }_2\right),\end{array}$ (7)

where

$l_{11}=-\left(\frac{{\beta }_1{E}^{\text{*}}+{\beta }_2{I}^{\text{*}}}{N^{\text{*}}}+d\right),l_{12}=-\frac{{\beta }_1{S}^{\text{*}}}{N^{\text{*}}},l_{13}=-\frac{{\beta }_2{S}^{\text{*}}}{N^{\text{*}}}$

$l_{14}=\delta ,l_{21}=\frac{{\beta }_1{E}^{\text{*}}+{\beta }_2{I}^{\text{*}}}{N^{\text{*}}},l_{22}=-\frac{{\beta }_1{S}^{\text{*}}}{N^{\text{*}}}-{\gamma }_2-d$

$l_{23}=\frac{{\beta }_2{S}^{\text{*}}}{N^{\text{*}}},l_{24}=0,l_{25}=-\epsilon$

$l_{31}=0,l_{32}=0,l_{33}=-\left(\omega +d\right)$

$l_{34}=0,l_{35}=\epsilon ,l_{36}=-{\gamma }_1$

$l_{41}=0,l_{42}={\gamma }_2,l_{43}=0$

$l_{44}=-\left(d+\delta \right),l_{45}={\gamma }_1$

3. Stability and Hopf Bifurcation

In this section, we thoroughly investigate the local stability of the system (1) at its equilibrium points. Taking the time delay as the bifurcation parameter, we further demonstrate the existence of Hopf bifurcation at the equilibrium points of the system through case-by-case analysis.

Firstly, we need to discuss the stability of system (1) at ${\tau }_1={\tau }_2=0$ .

Case 3.1.1 If ${\tau }_1={\tau }_2=0$ , We can get the characteristic equation of system (1) is

${\lambda }^4+a_3{\lambda }^3+a_2{\lambda }^2+a_1\lambda +a_0=0,$ (8)

where

$a_3=-\left(l_{11}+l_{22}+l_{35}+l_{33}+l_{36}+l_{44}\right),$

$\begin{array}{l}{a}_2=l_{11}{l}_{22}-l_{12}{l}_{21}+l_{11}{l}_{25}+l_{11}{l}_{33}-l_{13}{l}_{31}+l_{11}{l}_{36}+l_{11}{l}_{44}-l_{14}{l}_{44}\\ +l_{22}{l}_{33}-l_{23}{l}_{30}+l_{22}{l}_{36}-l_{23}{l}_{35}+l_{25}{l}_{33}+l_{25}{l}_{34}+l_{22}{l}_{44}\\ -l_{24}{l}_{42}+l_{25}{l}_{44}+l_{33}{l}_{44}-l_{34}{l}_{43}-l_{34}{l}_{45}+l_{36}{l}_{44},\end{array}$

$\begin{array}{l}{a}_1=l_{11}{l}_{23}{l}_{32}-l_{11}{l}_{22}{l}_{33}+l_{12}{l}_{21}{l}_{33}-l_{12}{l}_{23}{l}_{31}-l_{13}{l}_{21}{l}_{32}+l_{13}{l}_{22}{l}_{31}\\ -l_{11}{l}_{22}{l}_{36}+l_{11}{l}_{23}{l}_{35}-l_{11}{l}_{25}{l}_{33}+l_{12}{l}_{21}{l}_{36}-l_{13}{l}_{21}{l}_{35}+l_{13}{l}_{25}{l}_{31}\\ -l_{11}{l}_{25}{l}_{36}-l_{11}{l}_{22}{l}_{44}+l_{11}{l}_{24}{l}_{42}+l_{12}{l}_{21}{l}_{44}-l_{12}{l}_{24}{l}_{41}-l_{14}{l}_{21}{l}_{42}\\ +l_{14}{l}_{22}{l}_{41}-l_{11}{l}_{25}{l}_{44}+l_{14}{l}_{25}{l}_{41}-l_{11}{l}_{33}{l}_{44}+l_{11}{l}_{34}{l}_{43}+l_{13}{l}_{31}{l}_{44}\\ -l_{13}{l}_{34}{l}_{41}-l_{14}{l}_{31}{l}_{43}+l_{14}{l}_{33}{l}_{41}+l_{11}{l}_{34}{l}_{45}-l_{14}{l}_{31}{l}_{45}-l_{11}{l}_{36}{l}_{44}\\ +l_{14}{l}_{36}{l}_{41}-l_{25}{l}_{33}{l}_{44}+l_{25}{l}_{34}{l}_{43}-l_{24}{l}_{35}{l}_{45}+l_{25}{l}_{34}{l}_{45}-l_{25}{l}_{36}{l}_{44},\end{array}$

$\begin{array}{c}{a}_0=l_{11}{l}_{22}{l}_{33}{l}_{44}-l_{11}{l}_{22}{l}_{34}{l}_{43}-l_{11}{l}_{23}{l}_{32}{l}_{44}+l_{11}{l}_{23}{l}_{34}{l}_{42}+l_{11}{l}_{24}{l}_{32}{l}_{43}-l_{11}{l}_{24}{l}_{33}{l}_{42}\\ -l_{12}{l}_{21}{l}_{33}{l}_{44}+l_{12}{l}_{21}{l}_{34}{l}_{43}+l_{12}{l}_{23}{l}_{31}{l}_{44}-l_{12}{l}_{23}{l}_{34}{l}_{41}-l_{12}{l}_{24}{l}_{31}{l}_{43}+l_{12}{l}_{24}{l}_{33}{l}_{41}\\ +l_{13}{l}_{21}{l}_{32}{l}_{44}-l_{13}{l}_{21}{l}_{34}{l}_{42}-l_{13}{l}_{22}{l}_{31}{l}_{44}+l_{13}{l}_{22}{l}_{34}{l}_{41}+l_{13}{l}_{24}{l}_{31}{l}_{42}-l_{13}{l}_{24}{l}_{32}{l}_{41}\\ -l_{14}{l}_{21}{l}_{32}{l}_{43}+l_{14}{l}_{21}{l}_{33}{l}_{42}+l_{14}{l}_{21}{l}_{32}{l}_{43}-l_{14}{l}_{22}{l}_{33}{l}_{41}-l_{14}{l}_{23}{l}_{31}{l}_{42}+l_{14}{l}_{23}{l}_{32}{l}_{41}\\ -l_{11}{l}_{22}{l}_{34}{l}_{45}+l_{11}{l}_{24}{l}_{32}{l}_{45}+l_{12}{l}_{21}{l}_{34}{l}_{45}-l_{12}{l}_{24}{l}_{31}{l}_{45}-l_{14}{l}_{21}{l}_{32}{l}_{45}-l_{14}{l}_{21}{l}_{32}{l}_{45}\\ +l_{11}{l}_{22}{l}_{36}{l}_{44}-l_{11}{l}_{23}{l}_{35}{l}_{44}+l_{11}{l}_{24}{l}_{35}{l}_{43}-l_{11}{l}_{24}{l}_{36}{l}_{42}+l_{11}{l}_{25}{l}_{33}{l}_{44}-l_{11}{l}_{25}{l}_{34}{l}_{43}\end{array}$

According to Definition condition i), to guarantee the local asymptotic stability of the equilibrium point $E_1=\left(S^{\text{*}},E^{\text{*}},I^{\text{*}},R^{\text{*}}\right)$ of system (1) at ${\tau }_1={\tau }_2=0$ , the following Assumption (A1) needs to be satisfied.

(A1) The coefficients $a_i$ of Equation (8) satisfy the following conditions

i) $a_i>0,i=0,1,2,3$ ;

ii) $a_3{a}_2>a_1$ ;

iii) $a_3{a}_2{a}_1>a_1^2+a_3{a}_0$ .

Claim 1. According to the Routh-Hurwitz criterion, if Assumption (A1) holds, all roots of Equation (8) have negative real parts, which implies that $\left|\arg \left(\lambda \right)\right|>\frac{\pi }{2}$

is satisfied for all characteristic roots. Therefore, the equilibrium $E_1\left(S^{*},E^{*},I^{*},R^{*}\right)$ is locally asymptotically stable under the delay-free condition.

To analyze the stability and Hopf bifurcation behavior of system (1), this paper performs the Laplace transform on both sides of the linearized equation (7), and further derives the corresponding characteristic matrix of the system

$\Delta \left(s\right)=\left(\begin{array}{cccc}{s}^q-l_{11}& -l_{12}& -l_{13}& -l_{14}\\ -l_{21}& s^q-l_{22}-l_{25}{\text{e}}^{-s{\tau }_1}& -l_{23}& 0\\ 0& -l_{35}{\text{e}}^{-s{\tau }_1}& s^q-l_{33}-l_{36}{\text{e}}^{-s{\tau }_2}& 0\\ 0& -l_{42}& -l_{45}{\text{e}}^{-s{\tau }_2}& s^q-l_{44}\end{array}\right)$ .

It follows from $\det \left(\Delta \left(s\right)\right)=0$ that the characteristic equation of the fractional-order system with two delays is given by

$P_0\left(s,q\right)+P_1\left(s,q\right){\text{e}}^{-s{\tau }_1}+P_2\left(s,q\right){\text{e}}^{-s{\tau }_2}+P_3\left(s,q\right){\text{e}}^{-s\left({\tau }_1+{\tau }_2\right)}=0,$ (9)

where the polynomial functions $P_0\left(s,q\right),P_1\left(s,q\right),P_2\left(s,q\right),P_3\left(s,q\right)$ are defined as

$\begin{array}{l}{P}_0\left(s,q\right)=s^{4q}+E_{30}{s}^{3q}+E_{20}{s}^{2q}+E_{10}{s}^q+E_{00},\\ P_1\left(s,q\right)=E_{31}{s}^{3q}+E_{21}{s}^{2q}+E_{11}{s}^q+E_{01},\\ P_2\left(s,q\right)=E_{32}{s}^{3q}+E_{22}{s}^{2q}+E_{12}{s}^q+E_{02},\\ P_3\left(s,q\right)=E_{23}{s}^{2q}+E_{13}{s}^q+E_{03}.\end{array}$

where

$\begin{array}{l}{E}_{30}=-\left(l_{11}+l_{44}+l_{22}+l_{33}\right),\\ E_{20}=\left(l_{22}+l_{33}\right)\left(l_{11}+l_{44}\right)+l_{22}{l}_{33}-l_{12}{l}_{21},\\ E_{10}=-l_{44}{l}_{24}{l}_{42}-C\left(l_{11}+l_{44}\right)l_{22}{l}_{33}+l_{12}{l}_{24}\left(l_{33}+l_{44}\right),\\ E_{00}=l_{44}{l}_{21}{l}_{53}+\left(l_{11}{l}_{22}-l_{12}{l}_{21}\right)l_{33}{l}_{44},\\ E_{31}=-l_{25},\\ E_{21}=l_{25}{l}_{33}+l_{11}+l_{33}{l}_{44}-l_{23}{l}_{35},\\ E_{11}=-l_{25}{l}_{11}{l}_{33}+l_{33}{l}_{44}+l_{11}{l}_{44}+\left(l_{11}+l_{44}\right)l_{23}{l}_{35},\\ E_{01}=l_{11}{l}_{44}{l}_{25}{l}_{33}-l_{23}{l}_{35},\\ E_{32}=-l_{36},\\ E_{22}=l_{36}\left(l_{22}+l_{11}+l_{44}\right),\\ E_{12}=-\left(l_{11}+l_{44}\right)l_{22}{l}_{36}-l_{11}{l}_{36}{l}_{44}-l_{12}{l}_{21}{l}_{33},\\ E_{02}=l_{14}{l}_{21}{l}_{36}+l_{11}{l}_{22}{l}_{31}{l}_{44}-l_{12}{l}_{21}{l}_{31}{l}_{44},\\ E_{23}=l_{25}{l}_{36},\\ E_{13}=-\left(l_{11}+l_{44}\right)l_{25}{l}_{36},\\ E_{03}=l_{14}{l}_{21}{l}_{35}{l}_{45}.\end{array}$

Next, we investigate the following cases for the delay parameters ${\tau }_1$ and ${\tau }_2$ at the equilibrium point $E_1$ .

Case 3.1.2 If ${\tau }_1=0,{\tau }_2>0$ , then Equation (9) becomes

$P_0\left(s,q\right)+P_1\left(s,q\right)+\left(P_2\left(s,q\right)+P_3\left(s,q\right)\right){\text{e}}^{-s{\tau }_2}=0$ (10)

Assume that $s=\text{i}\omega$ ($\omega >0$ ) is a purely imaginary root of Equation (6). Substitute $s=\left(\cos \frac{\pi }{2}+\text{i}\sin \frac{\pi }{2}\right)\omega$ into $P_k\left(s,q\right)$ for $k=0,1,2,3$ . Let

$P_k\left(\text{i}\omega ,q\right)=A_k+\text{i}{B}_k$ . By separating the real and imaginary parts of the characteristic Equation (10), we derive the following system:

$\begin{array}{l}\mathrm{Re}:\left(A_2+A_3\right)\cos \omega {\tau }_2+\left(B_2+B_3\right)\sin \omega {\tau }_2=-\left(A_0+A_1\right),\\ \mathrm{Im}:-\left(A_2+A_3\right)\sin \omega {\tau }_2+\left(B_2+B_3\right)\cos \omega {\tau }_2=-\left(B_0+B_1\right).\end{array}$

where

$A_0=w^{4q}+E_{30}{w}^{3q}\cos \frac{3q\pi }{2}+E_{20}{w}^{2q}\cos q\pi +E_{10}{w}^q\cos \frac{q\pi }{2}+E_{00},$

$B_0=w^{4q}+E_{30}{w}^{3q}\sin \frac{3q\pi }{2}+E_{20}{w}^{2q}\sin q\pi +E_{10}{w}^q\sin \frac{q\pi }{2},$

$A_1=E_{31}{w}^{3q}\cos \frac{3q\pi }{2}+E_{21}{w}^{2q}\cos q\pi +E_{11}{w}^q\cos \frac{q\pi }{2}+E_{01},$

$B_1=E_{31}{w}^{3q}\sin \frac{3q\pi }{2}+E_{21}{w}^{2q}\sin q\pi +E_{11}{w}^q\sin \frac{q\pi }{2},$

$A_2=E_{32}{w}^{3q}\cos \frac{3q\pi }{2}+E_{22}{w}^{2q}\cos q\pi +E_{12}{w}^q\cos \frac{q\pi }{2}+E_{02},$

$B_2=E_{32}{w}^{3q}\sin \frac{3q\pi }{2}+E_{22}{w}^{2q}\sin q\pi +E_{12}{w}^q\sin \frac{q\pi }{2},$

$A_3=E_{23}{w}^{2q}\cos q\pi +E_{13}{w}^q\cos \frac{q\pi }{2}+E_{03},$

$B_3=E_{23}{w}^{2q}\sin q\pi +E_{13}{w}^q\sin \frac{q\pi }{2}.$

For the convenience of subsequent calculations, we perform equivalent transformation on Equation (10) and obtain the following system of equations

$\begin{array}{l}\mathrm{Re}:C_1\cos \omega {\tau }_2+C_2\sin \omega {\tau }_2=-C_3,\\ \mathrm{Im}:C_2\cos \omega {\tau }_2-C_1\sin \omega {\tau }_2=-C_4,\end{array}$

where

$C_1=A_2+A_3,C_2=B_2+B_3,C_3=A_0+A_1,C_4=B_0+B_1.$

We can readily derive

$\cos \omega {\tau }_2=\frac{C_1{C}_3+C_2{C}_4}{-C_1^2-C_2^2}=H_1\left(\omega \right),$

$\sin \omega {\tau }_2=\frac{-C_1{C}_4+C_2{C}_3}{-C_1^2-C_2^2}=H_2\left(\omega \right).$

By the identity ${\cos }^2\omega {\tau }_2+{\sin }^2\omega {\tau }_2=1$ , we have

$H_1^2\left(\omega \right)+H_2^2\left(\omega \right)=1.$ (11)

Assume that Equation (11) has at least one positive real root $\omega$ , which can be computed by the symbolic computation software Maple. Accordingly, the critical delay ${\tau }_{20}$ corresponding to the bifurcation point is defined as

${\tau }_{20}=\frac{1}{\omega }\underset{k=0,1,2,\cdots }{\min }\left(\arccos H_1\left(\omega \right)+2k\pi \right),$

where $\omega$ is the root of Equation (11).

Finally, we consider the transversality condition, namely Condition iii in Definition 2.3. Let $s\left(\tau \right)=\gamma \left(\tau \right)+i\omega \left(\tau \right)$ denote a root of Equation (6) in the neighborhood of $\tau ={\tau }_{20}$ , which satisfies $\gamma \left({\tau }_{20}\right)=0$ and $\omega \left({\tau }_{20}\right)={\omega }_{20}>0$ . Differentiating both sides of Equation (10) with respect to ${\tau }_2$ , we obtain

$\begin{array}{l}{P^{\prime }}_0\left(s,q\right)\frac{\text{d}s}{\text{d}{\tau }_2}+{P^{\prime }}_1\left(s,q\right)\frac{\text{d}s}{\text{d}{\tau }_2}+\left({P^{\prime }}_2\left(s,q\right)+{P^{\prime }}_3\left(s,q\right)\right){\text{e}}^{-s{\tau }_2}\frac{\text{d}s}{\text{d}{\tau }_2}\\ -\left(P_2\left(s,q\right)+P_3\left(s,q\right)\right){\text{e}}^{-s{\tau }_2}\left(s+{\tau }_2\frac{\text{d}s}{\text{d}{\tau }_2}\right)=0\end{array}$

where ${P^{\prime }}_i\left(s,q\right)$ stands for the derivative of the function $P_i\left(s,q\right)$ ($i=0,1,2,3$ ). Hence,

$\begin{array}{c}\frac{\text{d}s}{\text{d}{\tau }_2}=\frac{s\left[P_2\left(s,q\right)+P_3\left(s,q\right)\right]{\text{e}}^{-s{\tau }_2}}{{P^{\prime }}_0\left(s,q\right)+{P^{\prime }}_1\left(s,q\right)+\left[{P^{\prime }}_2\left(s,q\right)+{P^{\prime }}_3\left(s,q\right)\right]{\text{e}}^{-s{\tau }_2}-{\tau }_2\left[P_2\left(s,q\right)+P_3\left(s,q\right)\right]{\text{e}}^{-s{\tau }_2}}\\ =\frac{N}{D}\end{array}$ (12)

Let $P_k\left(i{\omega }_0,q\right)={\mu }_k+i{\epsilon }_k$ for $k=1,2,3$ , and ${P^{\prime }}_l\left(i{\omega }_0,q\right)={\hslash }_l+i{\ell }_l$ for $l=0,1,2,3$ . Substituting these expressions into Equation (8), we obtain

${\mathrm{Re}\left(\frac{\text{d}s}{\text{d}{\tau }_2}\right)|}_{{\tau }_2={\tau }_{20},\omega ={\omega }_{20}}=\frac{N_1{D}_1+D_2{N}_2}{D_1^2+D_2^2}$

where

$N_1={\omega }_0\left[-\left({\epsilon }_2+{\epsilon }_3\right)\text{cos}\left({\omega }_0{\tau }_{20}\right)+\left({\mu }_2+{\mu }_3\right)\text{sin}\left({\omega }_0{\tau }_{20}\right)\right],$

$N_2={\omega }_0\left[\left({\mu }_2+{\mu }_3\right)\text{cos}\left({\omega }_0{\tau }_{20}\right)+\left({\epsilon }_2+{\epsilon }_3\right)\text{sin}\left({\omega }_0{\tau }_{20}\right)\right],$

$\begin{array}{l}{D}_1={\hslash }_0+{\hslash }_1+\left({\hslash }_2+{\hslash }_3\right)\text{cos}\left({\omega }_{20}{\tau }_{20}\right)+\left({\ell }_2+{\ell }_3\right)\text{sin}\left({\omega }_{20}{\tau }_{20}\right)\\ -{\tau }_{20}\left({\mu }_2+{\mu }_3\right)\text{cos}\left({\omega }_{20}{\tau }_{20}\right)+{\tau }_{20}\left({\epsilon }_1+{\epsilon }_2\right)\text{sin}\left({\omega }_{20}{\tau }_{20}\right),\end{array}$

$\begin{array}{l}{D}_2={\hslash }_0+{\hslash }_1+\left({\hslash }_2+{\hslash }_3\right)\text{cos}\left({\omega }_{20}{\tau }_{20}\right)-\left({\hslash }_2+{\hslash }_3\right)\text{sin}\left({\omega }_{20}{\tau }_{20}\right)\\ -{\tau }_{20}\left({\epsilon }_2+{\epsilon }_3\right)\text{cos}\left({\omega }_{20}{\tau }_{20}\right)+{\tau }_{20}\left({\mu }_2+{\mu }_3\right)\text{sin}\left({\omega }_{20}{\tau }_{20}\right).\end{array}$

In summary, we propose the second assumption and obtain the following theorem1.

(A2):

${\mathrm{Re}\left(\frac{\text{d}s}{\text{d}{\tau }_2}\right)|}_{{\tau }_2={\tau }_{20},\omega ={\omega }_{20}}=\frac{N_1{D}_1+D_2{N}_2}{D_1^2+D_2^2}\ne 0$

Theorem 1. Assume that conditions (A1) and (A2) hold, then the equilibrium point $\left(S^{*},E^{*},I^{*},R^{*}\right)$ of the system (1) is asymptotically stable when ${\tau }_2\in \left[0,{\tau }_{20}\right)$ ; when ${\tau }_2={\tau }_{20}$ , Hopf bifurcation occurs near the equilibrium point.

Case 3.1.3 If ${\tau }_1>0,{\tau }_2=0$ , then Equation (9) becomes

$P_0\left(s,q\right)+P_2\left(s,q\right)+\left(P_1\left(s,q\right)+P_3\left(s,q\right)\right){\text{e}}^{-s{\tau }_1}=0.$ (13)

Separating the real and imaginary parts, we obtain the following system of trigonometric equations:

$\begin{array}{l}\mathrm{Re}:\left(A_1+A_3\right)\cos \left(\omega {\tau }_1\right)+\left(B_1+B_3\right)\sin \left(\omega {\tau }_1\right)=-\left(A_0+A_2\right),\\ \mathrm{Im}:\left(B_1+B_3\right)\cos \left(\omega {\tau }_1\right)-\left(A_1+A_3\right)\sin \left(\omega {\tau }_1\right)=-\left(B_0+B_2\right).\end{array}$

For simplicity, we denote

$G_1=A_1+A_3,G_2=B_1+B_3,G_3=A_0+A_2,G_4=B_0+B_2.$

Then the equation is rewritten as

$\begin{array}{l}\mathrm{Re}:G_1\cos \left(\omega {\tau }_1\right)+G_2\sin \left(\omega {\tau }_1\right)=-G_3,\\ \mathrm{Im}:G_2\cos \left(\omega {\tau }_1\right)-G_1\sin \left(\omega {\tau }_1\right)=-G_4.\end{array}$

By Cramer’s rule, we solve for $\text{cos}\left(\omega {\tau }_1\right)$ and $\text{sin}\left(\omega {\tau }_1\right)$ :

$\begin{array}{l}\cos \left(\omega {\tau }_1\right)=\frac{G_1{G}_3+G_2{G}_4}{-G_1^2-G_2^2}=k_1\left(\omega \right),\\ \sin \left(\omega {\tau }_1\right)=\frac{-G_1{G}_4+G_2{G}_3}{-G_1^2-G_2^2}=k_2\left(\omega \right).\end{array}$

Using the fundamental trigonometric identity

${\text{sin}}^2\left(\omega {\tau }_1\right)+{\text{cos}}^2\left(\omega {\tau }_1\right)=1,$

we arrive at the constraint

$K_1^2\left(\omega \right)+K_2^2\left(\omega \right)=1.$ (14)

Equation (14) guarantees the existence of a positive real root $\omega$ . The critical delay ${\tau }_{10}$ can be computed numerically with Maple 13 via the formula

${\tau }_{10}=\frac{1}{\omega }\underset{k=0,1,2,\cdots }{\min }\left(\arccos K_1\left(\omega \right)+2k\pi \right),$

The parameters are defined the same as those in Subsection case 3.1.2. Taking the partial derivative on both sides of Equation (13), we finally obtain:

${\mathrm{Re}\left(\frac{\text{d}s}{\text{d}{\tau }_1}\right)|}_{{\tau }_1={\tau }_{10},\omega ={\omega }_{10}}=\frac{E_1{F}_1+E_2{F}_2}{F_1^2+F_2^2}$

where

$E_1={\omega }_{10}\left[-\left({\epsilon }_1+{\epsilon }_3\right)\text{cos}\left({\omega }_{10}{\tau }_{10}\right)+\left({\mu }_1+{\mu }_3\right)\text{sin}\left({\omega }_{10}{\tau }_{10}\right)\right],$

$E_2={\omega }_{10}\left[\left({\mu }_1+{\mu }_3\right)\text{cos}\left({\omega }_{10}{\tau }_{10}\right)+\left({\epsilon }_1+{\epsilon }_3\right)\text{sin}\left({\omega }_{10}{\tau }_{10}\right)\right],$

$\begin{array}{l}{F}_1={\hslash }_0+{\hslash }_2+\left({\hslash }_1+{\hslash }_3\right)\text{cos}\left({\omega }_{10}{\tau }_{10}\right)+\left({\ell }_1+{\ell }_3\right)\text{sin}\left({\omega }_{10}{\tau }_{10}\right)\\ -{\tau }_{10}\left({\epsilon }_1+{\epsilon }_3\right)\text{cos}\left({\omega }_{10}{\tau }_{10}\right)-{\tau }_{10}\left({\mu }_1+{\mu }_3\right)\text{sin}\left({\omega }_{10}{\tau }_{10}\right),\end{array}$

$\begin{array}{l}{F}_2=\left({\ell }_0+{\ell }_2\right)+\left({\ell }_1+{\ell }_3\right)\text{cos}\left({\omega }_{10}{\tau }_{10}\right)-\left({\hslash }_1+{\hslash }_3\right)\text{sin}\left({\omega }_{10}{\tau }_{10}\right)\\ -{\tau }_{10}\left({\mu }_1+{\mu }_3\right)\text{cos}\left({\omega }_{10}{\tau }_{10}\right)+{\tau }_{10}\left({\epsilon }_1+{\epsilon }_3\right)\text{sin}\left({\omega }_{10}{\tau }_{10}\right).\end{array}$

Similar to the above case, we will get similar assumptions and theorems.

(A3):

${\mathrm{Re}\left(\frac{\text{d}s}{\text{d}{\tau }_1}\right)|}_{{\tau }_1={\tau }_{10},\omega ={\omega }_{10}}=\frac{E_1{F}_1+E_2{F}_2}{F_1^2+F_2^2}\ne 0$

Theorem 2. Assume that conditions (A1) and (A3) hold, then the equilibrium point $\left(S^{*},E^{*},I^{*},R^{*}\right)$ of the system (1) is asymptotically stable when ${\tau }_1\in \left[0,{\tau }_{10}\right)$ , when ${\tau }_1={\tau }_{10}$ , Hopf bifurcation occurs near the equilibrium point.

Case 3.1.4 If ${\tau }_1={\tau }_2=\tau >0$ , then equation (9) becomes

$P_0\left(s,q\right)+\left[P_1\left(s,q\right)+P_2\left(s,q\right)\right]{\text{e}}^{-s\tau }+P_3\left(s,q\right){\text{e}}^{-2s\tau }=0.$ (15)

Multiply both sides of the equation (15) by ${\text{e}}^{s\tau }$ , we have

$P_0\left(s,q\right){\text{e}}^{s\tau }+\left[P_1\left(s,q\right)+P_2\left(s,q\right)\right]+P_3\left(s,q\right){\text{e}}^{-s\tau }=0.$ (16)

Separating the real and imaginary parts, we obtain the following system of trigonometric equations:

$\begin{array}{l}\mathrm{Re}:\left(A_0+A_3\right)\cos \left(\omega \tau \right)+\left(B_3-B_0\right)\sin \left(\omega \tau \right)=-\left(A_1+A_2\right),\\ \mathrm{Im}:\left(B_0+B_3\right)\cos \left(\omega \tau \right)+\left(A_0-A_3\right)\sin \left(\omega \tau \right)=-\left(B_1+B_2\right).\end{array}$

For simplicity, we denote

$L_1=A_0+A_3,L_2=B_3-B_0,L_3=A_1+A_2,L_4=B_0+B_3,L_5=A_0-A_3,L_6=B_1+B_2.$

Then the system is rewritten as

$\begin{array}{l}\mathrm{Re}:L_1\cos \left(\omega \tau \right)+L_2\sin \left(\omega \tau \right)=-L_3,\\ \mathrm{Im}:L_4\cos \left(\omega \tau \right)-L_5\sin \left(\omega \tau \right)=-L_6.\end{array}$

By Cramer’s rule, we solve for $\text{cos}\left(\omega \tau \right)$ and $\text{sin}\left(\omega \tau \right)$ :

$\begin{array}{l}\text{cos}\left(\omega \tau \right)=\frac{L_2{L}_6-L_3{L}_5}{L_1{L}_5-L_2{L}_4}=M_1\left(\omega \right),\\ \text{sin}\left(\omega \tau \right)=\frac{L_3{L}_4-L_1{L}_6}{L_1{L}_5-L_2{L}_4}=M_2\left(\omega \right).\end{array}$

Using the fundamental trigonometric identity

${\text{sin}}^2\left(\omega \tau \right)+{\text{cos}}^2\left(\omega \tau \right)=1,$

we arrive at the constraint

$M_1^2\left(\omega \right)+M_2^2\left(\omega \right)=1.$ (17)

Equation (17) guarantees the existence of a positive real root $\omega$ . The critical delay ${\tau }_0$ can be computed numerically with Maple 13 via the formula

${\tau }_0=\frac{1}{\omega }\underset{k=0,1,2,\cdots }{\min }\left(\arccos M_1\left(\omega \right)+2k\pi \right).$

The parameters are defined the same as those in Subsection case 3.1.2. Taking the partial derivative on both sides of Equation (16), we finally obtain:

${\mathrm{Re}\left(\frac{\text{d}s}{\text{d}\tau }\right)|}_{\tau ={\tau }_0,\omega ={\omega }_0}=\frac{X_1{Y}_1+X_2{Y}_2}{Y_1^2+Y_2^2}$

where

$X_1={\omega }_0\left[-\left({\epsilon }_3+{\epsilon }_0\right)\text{cos}\left({\omega }_0{\tau }_0\right)+\left({\mu }_3-{\mu }_0\right)\text{sin}\left({\omega }_0{\tau }_0\right)\right],$

$X_2={\omega }_0\left[\left({\mu }_3-{\mu }_0\right)\text{cos}\left({\omega }_0{\tau }_0\right)+\left({\epsilon }_3+{\epsilon }_0\right)\text{sin}\left({\omega }_0{\tau }_0\right)\right],$

$\begin{array}{l}{Y}_1={\hslash }_1+{\hslash }_2+\left({\hslash }_0+{\hslash }_3\right)\text{cos}\left({\omega }_0{\tau }_0\right)+\left({\ell }_3-{\ell }_0\right)\text{sin}\left({\omega }_0{\tau }_0\right)\\ -{\tau }_0\left(\left({\epsilon }_0+{\epsilon }_3\right)\text{sin}\left({\omega }_0{\tau }_0\right)\right)-{\tau }_0\left(\left({\mu }_3-{\mu }_0\right)\text{cos}\left({\omega }_0{\tau }_0\right)\right),\end{array}$

$\begin{array}{l}{Y}_2=\left({\ell }_0+{\ell }_2\right)+\left({\ell }_0+{\ell }_3\right)\text{cos}\left({\omega }_0{\tau }_0\right)+\left({\hslash }_0-{\hslash }_3\right)\text{sin}\left({\omega }_0{\tau }_0\right)\\ +{\tau }_0\left({\epsilon }_0-{\epsilon }_3\right)\text{cos}\left({\omega }_0{\tau }_0\right)+{\tau }_0\left({\mu }_1+{\mu }_3\right)\text{sin}\left({\omega }_0{\tau }_0\right).\end{array}$

Similar to the above case, we will get similar assumptions and theorems.

(A4):

${\mathrm{Re}\left(\frac{\text{d}s}{\text{d}\tau }\right)|}_{\tau ={\tau }_0,\omega ={\omega }_0}=\frac{X_1{Y}_1+X_2{Y}_2}{Y_1^2+Y_2^2}\ne 0$

Theorem 3. Assume that conditions (A1) and (A4) hold, then the equilibrium point $\left(S^{*},E^{*},I^{*},R^{*}\right)$ of the controlled system (1) is asymptotically stable when $\tau \in \left[0,{\tau }_0\right)$ , when $\tau ={\tau }_0$ , Hopf bifurcation occurs near the equilibrium point.

4. Numerical Results

In this section, the influence of time delay as a bifurcation parameter on the stability of the equilibrium point of the system is intuitively felt through numerical examples. In this paper, two groups of parameters are selected as shown in Table 2. The first set of values makes the basic reproduction number $R_0<1$ , so there is only one disease-free equilibrium. However, this paper mainly studies the local equilibrium point, so we use the second set of data, which makes $R_0>1$ , there is a local equilibrium point, so we start numerical simulation. The fractional order parameter $q=0.9$ is selected, and the critical frequency and critical bifurcation point are calculated by numerical calculation maple. The following five pictures show the stability changes of the five equilibrium points due to different time delays.

Table 2. Parameter values in the first and second groups.

Parameters	First	Second
$\lambda$	1.50	0.50
${\beta }_1$	0.50	1.00
${\beta }_2$	0.30	0.30
$d$	0.05	0.05
$\delta$	0.01	0.01
${\gamma }_1$	0.55	0.50
${\gamma }_2$	0.10	0.01
$\mu$	0.20	0.20
$\epsilon$	1.20	1.20

Figure 1. ${\tau }_1={\tau }_2=0$ $E_1$ is asymptotically stable.

Figure 2. When ${\tau }_1=0$ , $E_1$ is asymptotically stable for ${\tau }_2=5.16<{\tau }_{20}$ .

Figure 3. When ${\tau }_1=0$ , $E_1$ is asymptotically stable for ${\tau }_2=6>{\tau }_{20}$ .

Figure 4. When ${\tau }_2=0$ , $E_1$ is asymptotically stable for ${\tau }_1=1.1<{\tau }_{10}$ .

Figure 5. When ${\tau }_2=0$ , $E_1$ is asymptotically stable for ${\tau }_1=1.4>{\tau }_{10}$ .

From Figures 1-5, we observe that small time delays preserve stability of the endemic equilibrium $E_1$ , while crossing the critical delay threshold triggers Hopf bifurcation and periodic oscillations. Figure 1 shows that when τ1 = τ2 = 0, the time-series of each compartment gradually converges, and the phase trajectories settle at a fixed point, so $E_1$ is asymptotically stable. For τ1 = 0, Figure 2 illustrates that $E_1$ remains asymptotically stable for τ2 = 5.16 < τ20. In Figure 3, as τ2 = 6 > τ20, periodic oscillations appear in all population curves and limit cycles form in phase portraits, which implies the onset of Hopf bifurcation. When τ2 = 0, Figure 4 presents stable dynamics for τ1 = 1.1 < τ10. By contrast, Figure 5 shows sustained periodic fluctuations once τ1 = 1.4 > τ10, confirming that a supercritical Hopf bifurcation takes place. Collectively, these simulations verify that each delay can destabilize E1 and induce periodic solutions when exceeding its corresponding critical value. From an epidemiological perspective, this warns public-health authorities that slow roll-out of quarantine policies will produce oscillatory waves of infections. To mitigate such periodic resurgence, intervention measures should be activated as early as possible to keep the effective response delay under the computed critical threshold.

5. Conclusions and Suggestions

For the fractional-order SEIR epidemic dynamics model constructed in this paper, the existence and uniqueness of the positive equilibrium point of the system are first strictly demonstrated. Then, the time delay term is selected as the core bifurcation parameter. By classifying and discussing the different situations of time delay, the local asymptotic stability of the positive equilibrium point is deeply analyzed, and the occurrence conditions of Hopf bifurcation of the system are determined. Based on the classical theory of Hopf bifurcation of fractional-order delay differential systems, a set of new sufficient criteria for determining the birth of Hopf bifurcation are derived. At the same time, the critical delay threshold of bifurcation of the system is accurately solved by Maple symbolic computing software.

Our mathematical and numerical results carry clear epidemiological implications. When the intervention response delay is kept below the critical threshold (${\tau }_2<{\tau }_{20}$ ), the endemic equilibrium remains asymptotically stable, meaning that infection cases will gradually converge to a steady level. Once intervention delay exceeds ${\tau }_{20}$ , Hopf bifurcation occurs and sustained periodic oscillations emerge in susceptible, exposed and infected populations, which corresponds to recurrent epidemic waves in reality. This indicates that delayed implementation of quarantine measures may trigger repeated resurgence of infections, even if other control parameters remain unchanged. Therefore, shortening the response lag of isolation interventions is crucial to suppressing periodic epidemic oscillations. Although full calibration against real-world case data is not the primary goal of this theoretical study, our parameter setup refers to epidemiological literature on latent periods and intervention timelines, so that our numerical observations can offer qualitative reference for practical public-health decision-making.

It should be noted that the parameter group used for numerical simulation in this paper is taken from typical ranges in existing epidemic-modelling literature for the purpose of qualitative verification of theoretical conclusions. The exact numerical values of critical delay thresholds ${\tau }_{10}$ and ${\tau }_{20}$ are parameter-dependent. Direct quantitative application of these threshold figures for real-world disease calls for further calibration against real surveillance data, which will be considered in our future research.

Funding

This work is supported by the National Natural Science Foundation of China (Grant No. 11872043), Natural Science Foundation of Sichuan Province (Grant No. 2023NSFSC1299), the Scientific Research and Innovation Team Program of Sichuan University of Science and Engineering (Grant No. SUSE652B002), Opening Fund of Key Laboratory of Higher Education of Sichuan Province for Enterprise Informationalization and Internet of Things (Grant No. 2023WZJ02).

Author Contributions

Conceptualization, F.H.Z. and X.L.L.; methodology, F.H.Z.; software, F.H.Z.; validation, F.H.Z., X.L.L., and Z.Z.; formal analysis, F.H.Z.; investigation, F.H.Z.; resources, X.L.L.; data curation, F.H.Z.; writing—original draft preparation, F.H.Z.; writing—review and editing, X.L.L.; visualization, F.H.Z.; supervision, X.L.L.; project administration, X.L.L. All authors have read and agreed to the published version of the manuscript.