- Research
- Open access
- Published:
Dynamics analysis and optimal control of delayed SEIR model in COVID-19 epidemic
Journal of Inequalities and Applications volume 2024, Article number: 66 (2024)
Abstract
The coronavirus disease 2019 (COVID-19) remains serious around the world and causes huge deaths and economic losses. Understanding the transmission dynamics of diseases and providing effective control strategies play important roles in the prevention of epidemic diseases. In this paper, to investigate the effect of delays on the transmission of COVID-19, we propose a delayed SEIR model to describe COVID-19 virus transmission, where two delays indicating the incubation and recovery periods are introduced. For this system, we prove its solutions are nonnegative and ultimately bounded with the nonnegative initial conditions. Furthermore, we calculate the disease-free and endemic equilibrium points and analyze the asymptotical stability and the existence of Hopf bifurcations at these equilibrium points. Then, by taking the weighted sum of the opposite number of recovered individuals at the terminal time, the number of exposed and infected individuals during the time horizon, and the system cost of control measures as the cost function, we present a delay optimal control problem, where two controls represent the social contact and the pharmaceutical intervention. Necessary optimality conditions of this optimal control problem are exploited to characterize the optimal control strategies. Finally, numerical simulations are performed to verify the theoretical analysis of the stability and Hopf bifurcations at the equilibrium points and to illustrate the effectiveness of the obtained optimal strategies in controlling the COVID-19 epidemic.
1 Introduction
The coronavirus disease 2019 (COVID-19) is an epidemic caused by the SARS-CoV-2 virus, which affects mostly the human respiratory system. Due to the strong infectivity, fast transmission speed, and long incubation period, it has brought significant losses to human life and the global economy [1, 2]. Although many prevention and control measures have been taken to mitigate the disease spread, the world is still suffering from the infection and death cases of COVID-19. According to the World Health Organization report, as of 2 November 2023, the total number of infections caused by the COVID-19 epidemic has reached 697,318,367, with more than 6,934,066 deaths [3]. Thus, there is an urgent need to study the COVID-19 virus transmission behavior and provide effective control measures to prevent the spread of the disease.
Mathematical models have been proven effective means in describing and understanding the complex transmission dynamics of epidemics [4–6]. In this regard, several mathematical models have been investigated to formulate and control COVID-19 spread [7–9]. SIR (Susceptible, Infectious, Removed) models have been proposed to simulate and predict the COVID-19 pandemic in [10, 11]. The stability of the SEIR (Susceptible, Exposed, Infectious, Removed) model is analyzed in [12]. Stability analysis of a spatial extension of the SEIR model is carried out in [13]. Optimal control of both non-pharmacological and pharmacological interventions for a modified SEIR epidemiological model is studied in [14]. Optimal control of COVID-19 spread considering quarantine effect on people with diabetes is investigated in [15]. Stability analysis and optimal control of a COVID-19 system, including susceptible, exposed, symptomatically infected, asymptomatically infected, hospitalized, and recovered individuals, are discussed in [16]. Bifurcation analysis and optimal control of a discrete SIR system are studied in [17]. Optimal control of the Omicron and Delta strains in COVID-19 is investigated in [18]. Although the aforementioned results are certainly valid and interesting, time delays are ignored in the dynamics analysis and optimal control problems. In fact, time delays exist in COVID-19 virus transmission since there is an incubation period for the exposed individual to manifest COVID-19 symptoms and signs, and a recovery period is required for the infected individual to become a recovered individual [19, 20].
Motivated by the above issue, in the current paper, we propose a delayed SEIR model to describe COVID-19 virus transmission in which two delays representing the incubation period and the recovery period are introduced. For this delayed SEIR system, we prove that its solutions with the nonnegative initial conditions are nonnegative and ultimately bounded. Furthermore, we calculate the disease-free and endemic equilibrium points and analyze the locally asymptotical stability and Hopf bifurcations of these two equilibrium points. Then, by taking the weighted sum of the opposite number of recovered individuals at the terminal time, the number of exposed and infected individuals during the time horizon, and the system cost of control measures as the cost function, we present a delay optimal control problem, where two controls denote the social contact and the pharmaceutical intervention. We also derive the necessary optimality conditions to characterize the optimal controls. Based on the necessary optimality conditions, we develop a numerical approach for solving the delay optimal control problem, which is different from the reported methods in [21–26]. Finally, numerical simulations are performed to verify the theoretical analysis of the stability and the existence of Hopf bifurcations at equilibrium points and demonstrate the effectiveness of computed optimal strategies in controlling the COVID-19 epidemic. Compared with the existing literature, the main contributions and innovations of this paper include: (i) a novel SEIR system with incubation and recovery delays is proposed for COVID-19 virus transmission; (ii) the stability and Hopf bifurcation of the equilibria in the delayed SEIR model are thoroughly analyzed; and (iii) an effective solution approach is developed for solving the delayed optimal control problem with both the terminal and integral costs.
The paper is organized as follows: In Sect. 2, the delayed SEIR system is proposed, and some important properties are proved. In Sect. 3, the stability and existence of Hopf bifurcations at the disease-free and endemic equilibrium points are analyzed. In Sect. 4, we present the delay optimal control problem and derive the corresponding necessary optimality conditions. In Sect. 5, we illustrate the numerical simulation results. In Sect. 6, we provide conclusions.
2 Delayed SEIR model
In the COVID-19 virus transmission process, the primary source of transmission is social contact between susceptible and exposed individuals who are usually asymptomatic. Another transmission source is direct contact between frontline workers and infected individuals. Some exposed individuals usually recover due to their strong natural immune systems or use over-the-counter medications, while others require hospitalization after a certain incubation period. Many infected individuals recover from the disease after extensive treatment, but unfortunately, for some, the infection costs their lives. As is well-known, a susceptible individual does not show the corresponding symptoms quickly after infection with the virus; that is, there is an incubation period for the exposed individual to be an infective individual. Besides, it is required for the infected individual to spend some time to be a recovered individual. Thus, time delays must be incorporated into the COVID-19 virus transmission process. Based on SEIR system [14], we propose the following delayed SEIR system to describe the COVID-19 virus transmission:
where t is the time of process; \(S(t)\), \(E(t)\), \(I(t)\), and \(R(t)\) are, respectively, the susceptible, exposed, infected, and recovered individuals; \(N(t)=S(t)+E(t)+I(t)+R(t)\) is the total number of individuals affected by the outbreak; Λ is the recruitment in the human population; d is the natural mortality rate; \({\beta _{1}}\) is the transmission rate due to social contact; \({\beta _{2}}\) is the transmission rate due to frontline contact; ε is the infection rate; \({\gamma _{1}}\) is the recovery rate of infectious individuals; \({\gamma _{2}}\) is the immune recovery rate; μ is the probability of death due to COVID-19; \(\tau _{1}\) and \(\tau _{2}\) are two delay arguments that indicate the incubation and recovery periods, respectively. A flow chart of the delayed SEIR system (1) is shown in Fig. 1.
For system (1), the initial conditions are given by
where \(t \in [-\tau ,0]\) with \(\tau =\max \{\tau _{1},\tau _{2}\}\); and \(\varphi _{i}\), \(i=1,2,3,4\), are given functions. For system (1) with the initial conditions (2), we have the following theorem:
Theorem 1
If the initial conditions are as follows \(\varphi _{i} \geq 0\), \(i=1,2,3,4\), then the solutions of system (1) are nonnegative and ultimately bounded.
Proof
It is easy to show by Theorem 5.2.1 [27] that \(S(t)\geq 0\), \(E(t)\geq 0\), \(I(t)\geq 0\) and \(R(t)\geq 0\) for \(t\geq 0\). Hence, the nonnegative cone \(\mathbb{R}^{4}_{+}\) is invariant for system (1) with nonnegative initial conditions \(\varphi _{i} \geq 0\), \(i=1,2,3,4\). Furthermore, from (1), we have
which gives
Then, the delayed SEIR system (1) has a biologically feasible range as indicated below:
which implies the ultimate boundedness of the solutions. □
3 Equilibrium points, stability and Hopf bifurcation
In this section, we will analyze the local stability of system (1) at the equilibrium points and establish the existence of Hopf bifurcations at the equilibrium points.
3.1 Equilibrium points
For system (1), the disease-free and the endemic equilibrium points are defined as
where
Now, we consider the basic reproduction number \(R_{0}\) [28], which can be used to determine whether a disease is an epidemic. To begin with, we rewrite system (1) as \(\dot{x} = f - v\), where
Then, by taking E and I as the infection compartments, we obtain the following Jacobian matrices F and V at the disease-free equilibrium point \(E_{0}\):
Thus, \(R_{0}\) can be defined as
Furthermore, the endemic equilibrium point \(E_{1}\) defined in (3) can be rewritten as
Therefore, system (1) has a positive endemic equilibrium point when \(R_{0} > 1\).
3.2 Stability and Hopf bifurcation of disease-free equilibrium point
The Jacobian matrix of system (1) at \({E_{0}}\) is given by
For (6), two negative eigenvalues are given by \(\lambda _{1}=\lambda _{2}=-d\). The other two eigenvalues are determined by the following characteristic equation:
where \(A_{1}^{1} =\gamma _{2} + d - \beta _{1} + \mu + d\); \(A_{2}^{1} =(\gamma _{2} + d - \beta _{1})(\mu + d)\); \(B_{0}^{1} =\varepsilon \); \(B_{1}^{1} =\varepsilon (\mu + d - \beta _{2})\); \(C_{0}^{1} =\gamma _{1}\); \(C_{1}^{1} = \gamma _{1}(\gamma _{2} + d - \beta _{1})\); and \(D_{0}^{1} =\varepsilon \gamma _{1}\).
Now, we study the following cases for delay arguments \(\tau _{1}\) and \(\tau _{2}\) at \(E_{0}\).
Case 3.2.1 If \(\tau _{1}=\tau _{2}=0\), then equation (7) becomes
It follows from (5) that
where \(R_{1}= \frac{\beta _{1} }{\varepsilon + \gamma _{2} + d}\); and \(R_{2}= \frac{\beta _{2} \varepsilon }{(\varepsilon + \gamma _{2} + d)(\mu + \gamma _{1} +d)}\). Obviously, \(R_{1}>0\), and \(R_{2}>0\). Let
Then, we have
and
By the Routh-Hurwitz theorem [29], we obtain the following theorem:
Theorem 2
Suppose that \({\tau _{1}=\tau _{2}=0}\).
(i) If \({R_{0}<1}\), then \({E_{0}}\) is locally asymptotically stable;
(ii) If \({R_{0}>1}\), then \(E_{0}\) is unstable.
Case 3.2.2 If \(\tau _{1}=0\), and \(\tau _{2}>0\), then equation (7) turns to
Next, we will discuss the Hopf bifurcation of (10). Obviously, \(\lambda =i\omega\) (\(\omega >0\)) is a pure imaginary root of the equation (10) if and only if
Separating the real part from the imaginary part in (11) gives
By adding the squares of two equations in (12), we obtain
Let \(z=\omega ^{2}\), \(h_{01}=(A_{2}^{1} + B_{1}^{1})^{2} - (C_{1}^{1} + D_{0}^{1})^{2}\) and \(h_{11}=(A_{1}^{1} + B_{0}^{1})^{2} - 2(A_{2}^{1} + B_{1}^{1}) - (C_{0}^{1})^{2}\). Then, equation (13) can be rewritten as
Without loss of generality, we assume that equation (14) has two positive roots, denoted by \(z_{1}\) and \(z_{2}\). Then, two positive roots of equation (13) are
Furthermore, substituting (15) into (12) and solving the resulting equations, we have
where \(M_{1} =C_{1}^{1} + D_{0}^{1}\); \(M_{2} =C_{0}^{1}\); \(M_{3}(\omega _{2}^{j}) = ({\omega _{2}^{j}} )^{2} - (A_{2}^{1} + B_{1}^{1})\); and \(M_{4}(\omega _{2}^{j}) =-(A_{1}^{1} + B_{0}^{1})\omega _{2}^{j}\). Clearly,
Thus, we can define
Now, we consider the roots of equation (10) as a function of \(\tau _{2}\) and the root \(\lambda _{1}(\tau _{2})=\xi _{1}(\tau _{2})+i \omega _{1}(\tau _{2})\) of equation (10) such that
where \(\tau _{2}^{0}\) and \(\omega _{2}^{0}\) are as defined in (16). Substituting \(\lambda _{1}(\tau _{2})\) into equation (10) and taking the derivative with respect to \(\tau _{2}\), we have
Then,
where \(M_{1}^{1}= A_{1}^{1} + B_{0}^{1}\); \(S_{1}^{1}= 2\omega _{2}^{0}\); \(N_{1}^{1}= (A_{1}^{1} + B_{0}^{1}) ({\omega _{2}^{0}} )^{2}\); \(T_{1}^{1}= ({\omega _{2}^{0}} )^{3} - (A_{2}^{1} + B_{1}^{1})\omega _{2}^{0}\); \(Q_{1}^{1}= C_{0}^{1}\); \(P_{1}^{1} = -C_{0}^{1} ({\omega _{2}^{0}} )^{2}\); and \(O_{1}^{1} = (C_{1}^{1}+D_{0}^{1})\omega _{2}^{0} \). It follows that if the condition
is satisfied, then
Therefore, by Lemma 2.2 in [30] and Hopf bifurcation theory [31], we obtain the following theorem:
Theorem 3
Suppose that \(\tau _{1} = 0\), and \(\tau _{2} > 0\). If equation (14) has positive roots and condition (17) is met, then
(i) when \(\tau _{2} \in (0,\tau _{2}^{0})\), \(E_{0}\) is locally asymptotically stable;
(ii) when \(\tau _{2}>\tau _{2}^{0}\), \(E_{0}\) is unstable;
(iii) when \(\tau _{2}=\tau _{2}^{0}\), system (1) has a Hopf bifurcation at \(E_{0}\).
Case 3.2.3 If \(\tau _{1}>0\), and \(\tau _{2}=0\), then equations (7) and (13) become
and
respectively, where \(A_{1}^{1}\), \(A_{2}^{1}\), \(B_{0}^{1}\), \(B_{1}^{1}\), \(C_{0}^{1}\), \(C_{1}^{1}\), and \(D_{0}^{1}\) are as defined in (7). Using similar argument as for Case 3.2.2, we have
where \(\omega _{1}^{1}\) and \(\omega _{1}^{2}\) are two positive roots of (19); \(\bar{M}_{1} =B_{1}^{1} + D_{0}^{1}\); \(\bar{M}_{2} =B_{0}^{1}\); \(\bar{M}_{3}(\omega _{1}^{j}) = ({\omega _{1}^{j}} )^{2} - (A_{2}^{1} + C_{1}^{1})\); and \(\bar{M}_{4}(\omega _{1}^{j}) =-(A_{1}^{1} + C_{0}^{1})\omega _{1}^{j}\). Moreover, we consider the root of equation (18) as \(\lambda _{2}(\tau _{1})=\xi _{2}(\tau _{1})+i\omega _{2}(\tau _{1})\) with \(\xi _{2}(\tau ^{0}_{1})=0\) and \(\omega _{2}(\tau _{1}^{0})=\omega _{1}^{0}\). Then, we have
where \(\bar{M}_{1}^{1}= A_{1}^{1} + C_{0}^{1}\); \(\bar{S}_{1}^{1}= 2\omega _{1}^{0}\); \(\bar{N}_{1}^{1}= (A_{1}^{1} + C_{0}^{1}) ({\omega _{1}^{0}} )^{2}\); \(\bar{T}_{1}^{1}= ({\omega _{1}^{0}} )^{3} - (A_{2}^{1} + C_{1}^{1}) \omega _{1}^{0}\); \(\bar{Q}_{1}^{1}= B_{0}^{1}\); \(\bar{P}_{1}^{1} = -B_{0}^{1} ({\omega _{1}^{0}} )^{2}\); and \(\bar{O}_{1}^{1} = (B_{1}^{1}+D_{0}^{1})\omega _{1}^{0}\). Obviously, if condition
is satisfied, then
Thus, we have the following theorem:
Theorem 4
Suppose that \(\tau _{2} = 0\), and \(\tau _{1} > 0\). If equation (19) has positive roots and condition (20) is met, then
(i) when \(\tau _{1} \in (0,\tau _{1}^{0})\), \(E_{0}\) is locally asymptotically stable;
(ii) when \(\tau _{1}>\tau _{1}^{0}\), \(E_{0}\) is unstable;
(iii) when \(\tau _{1}=\tau _{1}^{0}\), system (1) has a Hopf bifurcation at \(E_{0}\).
Case 3.2.4 If \(\tau _{1}>0\), and \(\tau _{2}>0\), then as shown in [32], it is difficult to directly analyze the Hopf bifurcation with two delays. Thus, further, we only discuss the Hopf bifurcation of equilibrium point for system (1) with one delay and the other delay lying in its stability region. Without loss of generality, we investigate the Hopf bifurcation at \(E_{0}\) for system (1) with \(\tau _{1}>0\) and \(\tau _{2} \in (0, \tau _{2}^{0})\), where \(\tau _{2}^{0}\) is as given in Case 3.2.2. Let \(\tau _{2}\in (0, \tau _{2}^{0})\) be arbitrary but fixed. Then, \(\lambda =i\omega\) (\(\omega >0\)) is a pure imaginary root of equation (7) if and only if
Separating the real part from the imaginary part in (21) gives
By adding the squares of two equations in (22), we obtain
where \(l_{1}(\omega )=-2C_{0}^{1}\sin ({\omega \tau _{2}})\); \(l_{2}(\omega )=(A_{1}^{1})^{2} - 2A_{2}^{1} + (C_{0}^{1})^{2} - (B_{0}^{1})^{2} + 2(-C_{1}^{1} + A_{1}^{1}C_{0}^{1})\cos{\omega \tau _{2}}\); \(l_{3}(\omega )=2(A_{2}^{1}C_{0}^{1} + B_{0}^{1}D_{0}^{1} - A_{1}^{1}C_{1}^{1}) \sin ({\omega \tau _{2}})\); and \(l_{4}(\omega )=(A_{2}^{1})^{2} + (C_{1}^{1})^{2} - (B_{1}^{1})^{2} - (D_{0}^{1})^{2} + 2(A_{2}^{1}C_{1}^{1} - B_{1}^{1}D_{0}^{1})\cos{\omega \tau _{2}} \). Let
Clearly, if \(l_{4}(0)<0\), then \(G_{1}(0)<0\). Besides, \(\lim_{\omega \to +\infty} G_{1}(\omega )=+\infty \). Hence, equation (23) has at least one positive root denoted by \(\tilde{\omega}_{1}\). Furthermore, substituting \(\tilde{\omega}_{1}\) into (22) and solving the resulting equations, we have
where \(M_{5}(\tilde{\omega}_{1})=(B_{1}^{1} + D_{0}^{1} \cos ({ \tilde{\omega}_{1} \tau _{2}}))\); \(M_{6}(\tilde{\omega}_{1})=(B_{0}^{1}\tilde{\omega}_{1} - D_{0}^{1} \sin ({\tilde{\omega}_{1} \tau _{2}}))\); \(M_{7}(\tilde{\omega}_{1}) = -[(-{\tilde{\omega}_{1}}^{2} + A_{2}^{1}) + C_{0}^{1} \tilde{\omega}_{1} \sin ({\tilde{\omega}_{1} \tau _{2}}) + C_{1}^{1} \cos ({\tilde{\omega}_{1} \tau _{2}})]\); and \(M_{8}(\tilde{\omega}_{1})=-(A_{1}^{1}\tilde{\omega}_{1} + C_{0}^{1} \tilde{\omega}_{1} \cos ({\tilde{\omega}_{1} \tau _{2}}) - C_{1}^{1} \sin ({ \tilde{\omega}_{1} \tau _{2}})) \). Obviously,
Thus, we can define
As discussed in Case 3.2.3, we consider the roots of equation (7) as a function of \(\tau _{1}\) and the root \(\lambda _{3}(\tau _{1})= \xi _{3}(\tau _{1})+i\omega _{3}(\tau _{1})\) of equation (7) such that
Substituting \(\lambda _{3}(\tau _{1})\) into equation (7) and taking the derivative with respect to \(\tau _{1}\), we have
Then,
where \(M_{3}^{1}= A_{1}^{1} + C_{0}^{1} \cos (\tilde{\omega}^{0}_{1} \tau _{2}) - \tau _{2} C_{0}^{1} \tilde{\omega}^{0}_{1} \sin (\tilde{\omega}^{0}_{1} \tau _{2}) - \tau _{2} C_{1}^{1} \cos (\tilde{\omega}^{0}_{1} \tau _{2})\); \(S_{3}^{1}= 2\tilde{\omega}^{0}_{1} - C_{0}^{1} \sin (\tilde{\omega}^{0}_{1} \tau _{2}) - \tau _{2} C_{0}^{1} \tilde{\omega}^{0}_{1} \cos ( \tilde{\omega}^{0}_{1} \tau _{2}) + \tau _{2} C_{1}^{1} \sin ( \tilde{\omega}^{0}_{1} \tau _{2})\); \(N_{3}^{1}= -A_{1}^{1} (\tilde{\omega}^{0}_{1})^{2} - C_{0}^{1}( \tilde{\omega}^{0}_{1})^{2}\cos (\tilde{\omega}^{0}_{1} \tau _{2}) + C_{1}^{1} \tilde{\omega}^{0}_{1}\sin (\tilde{\omega}^{0}_{1} \tau _{2})\); \(T_{3}^{1}= - (\tilde{\omega}^{0}_{1})^{3} + A_{2}^{1} \tilde{\omega}^{0}_{1} + C_{0}^{1}(\tilde{\omega}^{0}_{1})^{2}\sin (\tilde{\omega}^{0}_{1} \tau _{2}) + C_{1}^{1}\tilde{\omega}^{0}_{1}\cos (\tilde{\omega}^{0}_{1} \tau _{2})\); \(Q_{3}^{1}= B_{0}^{1} - \tau _{2}^{0} D_{0}^{1} \cos (\tilde{\omega}^{0}_{1} \tau _{2}^{0})\); \(R_{3}^{1}= \tau _{2}^{0} D_{0}^{1} \sin (\tilde{\omega}^{0}_{1} \tau _{2}^{0})\); \(P_{3}^{1} = -B_{0}^{1} (\tilde{\omega}^{0}_{1})^{2}+ D_{0}^{1} \cos (\tilde{\omega}^{0}_{1} \tau _{2}^{0})\); and \(O_{3}^{1} = B_{1}^{1} \tilde{\omega}^{0}_{1} - D_{0}^{1} \sin ( \tilde{\omega}^{0}_{1} \tau _{2}^{0})\). It follows that if the condition
is satisfied, then
As a result, by the Hopf bifurcation theory [31], we obtain the following theorem:
Theorem 5
Suppose that \(\tau _{1} > 0\), and \(\tau _{2} \in (0,\tau _{2}^{0})\) with \(\tau ^{0}_{2}\) as defined in Case 3.2.2. If equation (23) has positive roots and condition (26) is met, then
(i) when \(\tau _{1} \in (0,\tilde{\tau}^{0}_{1})\), \(E_{0}\) is locally asymptotically stable;
(ii) when \(\tau _{1}>\tilde{\tau}^{0}_{1}\), \(E_{0}\) is unstable;
(iii) when \(\tau _{1}=\tilde{\tau}^{0}_{1}\), system (1) has a Hopf bifurcation at \(E_{0}\).
3.3 Stability and Hopf bifurcation of endemic equilibrium point
The Jacobian matrix of the system (1) at \({E_{1}}\) is given by
where \(N^{*}=S^{*}+E^{*}+I^{*}+R^{*}\); \(Y= \frac{\beta _{1} (N^{*}-S^{*})E^{*}}{(N^{*})^{2}} + \frac{\beta _{2} (N^{*} -S^{*})I^{*}}{(N^{*})^{2}}\); \(Y_{1}=\frac{\beta _{1} S^{*}N^{*} - \beta _{1} S^{*}E^{*} - \beta _{2} S^{*}I^{*}}{(N^{*})^{2}}\); \(Y_{2}=\frac{\beta _{2} S^{*}N^{*} - \beta _{2} S^{*}I^{*} - \beta _{1} S^{*}E^{*}}{(N^{*})^{2}}\); \(Y_{3}= \varepsilon e^{-\lambda \tau _{1}} +\gamma _{2} + d\); and \(Y_{4} =-(\mu + d + \gamma _{1} e^{-\lambda \tau _{2}})\). For (27), one negative eigenvalue is given by \(\lambda _{1}=-d \). The other eigenvalues are determined by the following characteristic equation:
where \(A_{1}^{2}=(Y + d)+(\gamma _{2} + d)+(\mu + d) - Y_{1}\); \(A_{2}^{2}=(Y + d)(\gamma _{2} + d) +(Y + d)(\mu + d) + (\gamma _{2} + d)(\mu + d) - (\mu + 2d)Y_{1} + \gamma _{2} {\frac{\beta _{1} S^{*}E^{*} + \beta _{2} S^{*}I^{*}}{(N^{*})^{2}}}\); \(A_{3}^{2}=(Y + d)(\gamma _{2} + d) (\mu + d) -Y_{1}(\mu + d)d + { \frac{\beta _{1} S^{*}E^{*} + \beta _{2} S^{*}I^{*}}{(N^{*})^{2}}} (\mu + d) \gamma _{2} \); \(B_{0}^{2}=\varepsilon \); \(B_{1}^{2}=(Y + d)\varepsilon + (\mu + d)\varepsilon - Y_{2} \varepsilon \); \(B_{2}^{2}=(Y + d)(\mu + d)\varepsilon - \varepsilon d Y_{2}\); \(C_{0}^{2}=\gamma _{1}\); \(C_{1}^{2}=(Y + d)\gamma _{1} + (\gamma _{2} + d)\gamma _{1} - Y_{1} \gamma _{1}\); \(C_{2}^{2}=(Y + d)(\gamma _{2} + d)\gamma _{1} + \gamma _{1} \gamma _{2} {\frac{\beta _{1} S^{*}E^{*} + \beta _{2} S^{*}I^{*}}{(N^{*})^{2}}} - \gamma _{1}\,dY_{1}\); \(D_{0}^{2}=\varepsilon \gamma _{1}\); and \(D_{1}^{2}={\frac{\beta _{1} S^{*}E^{*} + \beta _{2} S^{*}I^{*}}{(N^{*})^{2}}} \varepsilon \gamma _{1} + (Y + d)\varepsilon \gamma _{1}\).
Now, we investigate the following cases for \(\tau _{1}\) and \(\tau _{2}\) at the endemic equilibrium point \(E_{1}\).
Case 3.3.1 If \(\tau _{1}= \tau _{2}= 0\), then equation (28) turns to
Let
Obviously,
and
Then, we have
Recall that the positive endemic equilibrium point \(E_{1}\) exists if and only if \(R_{0}> 1\). Thus, by the Routh-Hurwitz theorem [29], we have the following theorem:
Theorem 6
Suppose that \(\tau _{1} =\tau _{2}= 0\). If \(R_{0}>1\), then \(E_{1}\) is locally asymptotically stable.
Case 3.3.2 If \(\tau _{1}=0\), and \(\tau _{2} > 0\), then equation (28) turns to
To discuss the Hopf bifurcation of (30), we observe that \(\lambda =i\omega\) (\(\omega >0\)) is a pure imaginary root of the equation (30) if and only if
Separating the real part from the imaginary part in (31) gives
By adding the squares of two equations in (32), we obtain
Let \(z = \omega ^{2}\), \(h_{03}=(A_{3}^{2} + B_{2}^{2})^{2} - (C_{2}^{2} + D_{1}^{2})^{2}\), \(h_{13}=(A_{2}^{2} + B_{1}^{2})^{2} - 2(A_{1}^{2} + B_{0}^{2})(A_{3}^{2} + B_{2}^{2}) - (C_{1}^{2} + D_{0}^{2})^{2} + 2C_{0}^{2}(C_{2}^{2} + D_{1}^{2})\) and \(h_{23}=(A_{1}^{2} + B_{0}^{2})^{2} - 2(A_{2}^{2} + B_{1}^{2}) - (C_{0}^{2})^{2}\). Then, (33) can be rewritten as
Without loss of generality, we assume that equation (34) has three positive roots denoted by \(\bar{z}_{j}\), \(j = 1,2,3\). Then, three positive roots of equation (33) are
Furthermore, substituting (35) into (32) and solving the resulting equations, we have
where \(\hat{F}={A_{1}^{2} + B_{0}^{2}}\); \(\hat{H}={A_{2}^{2} + B_{1}^{2}}\); \(\hat{I}={A_{3}^{2} + B_{2}^{2}}\); \(\hat{J}={C_{1}^{2} + D_{0}^{2}}\); and \(\hat{K}={C_{2}^{2} + D_{1}^{2}}\). Clearly,
Thus, we can define
Now, we consider the roots of equation (30) as a function of \(\tau _{2}\) and the root \(\tilde{\lambda}_{1}(\tau _{2})=\alpha _{1}( \tau _{2})+i\tilde{\omega}_{1}(\tau _{2})\) of equation (30) such that
where \(\bar{\tau}_{2}^{0}\) and \(\bar{\omega}_{2}^{0}\) are as defined in (36). Substituting \(\tilde{\lambda}_{1}(\tau _{2})\) into equation (30) and taking the derivative with respect to \(\tau _{2}\), we have
Then,
where \(M_{1}^{2}= -3(\bar{\omega}_{2}^{0})^{2} + (A_{2}^{2} + B_{1}^{2})\), \(S_{1}^{2}= 2(A_{1}^{2} + B_{0}^{2})\bar{\omega}_{2}^{0}\), \(N_{1}^{2}= -(\bar{\omega}_{2}^{0})^{4} + (A_{2}^{2} + B_{1}^{2})(\bar{\omega}_{2}^{0})^{2}\), \(T_{1}^{2}= (A_{1}^{2} + B_{0}^{2})(\bar{\omega}_{2}^{0})^{3} - (A_{3}^{2} + B_{2}^{2}) \bar{\omega}_{2}^{0}\), \(Q_{1}^{2}= C_{1}^{2} + D_{0}^{2}\), \(R_{1}^{2}= 2C_{0}^{2} \bar{\omega}_{2}^{0}\), \(P_{1}^{2} = -(C_{1}^{2} + D_{0}^{2})(\bar{\omega}_{2}^{0})^{2}\), and \(O_{1}^{2} = -C_{0}^{2}(\bar{\omega}_{2}^{0})^{3} + (C_{2}^{2} + D_{1}^{2})\bar{\omega}_{2}^{0} \). It follows that if the condition
is satisfied, then
Hence, by the Hopf bifurcation theory [31] and Lemma 2.1 in [33], we can obtain the following theorem:
Theorem 7
Suppose that \(\tau _{1} = 0\), and \(\tau _{2} > 0\). If equation (34) has positive roots and condition (37) is met, then
(i) when \(\tau _{2} \in (0,\bar{\tau}_{2}^{0})\), \(E_{1}\) is locally asymptotically stable;
(ii) when \(\tau _{2}>\bar{\tau}_{2}^{0}\), \(E_{1}\) is unstable;
(iii) when \(\tau _{2}=\bar{\tau}_{2}^{0}\), system (1) has a Hopf bifurcation at \(E_{1}\).
Case 3.3.3 If \(\tau _{1}>0\), and \(\tau _{2}=0\), then equations (28) and (33) become
and
respectively, where \(A_{1}^{2}\), \(A_{2}^{2}\), \(A_{3}^{2}\), \(B_{0}^{2}\), \(B_{1}^{2}\), \(B_{2}^{2}\), \(C_{0}^{2}\), \(C_{1}^{2}\), \(C_{2}^{2}\), \(D_{0}^{2}\), and \(D_{1}^{2}\) are as defined in (28). Using similar argument as given for Case 2.2.3, we have
where \(\bar{\omega}_{1}^{1}\), \(\bar{\omega}_{1}^{2}\), and \(\bar{\omega}_{1}^{3}\) are three positive roots of (39); \(\bar{F}={A_{1}^{2} + C_{0}^{2}}\); \(\bar{H}={A_{2}^{2} + C_{1}^{2}}\); \(\bar{I}={A_{3}^{2} + C_{2}^{2}}\); \(\bar{J}={B_{1}^{2} + D_{0}^{2}}\); and \(\bar{K}={B_{2}^{2} + D_{1}^{2}}\). Furthermore, we consider the root of equation (38) as \(\tilde{\lambda}_{2}(\tau _{1})=\alpha _{2}(\tau _{1})+i \tilde{\omega}_{2}(\tau _{1})\) with \(\alpha _{2}(\bar{\tau}_{1}^{0})=0\) and \(\tilde{\omega}_{2}(\bar{\tau}_{1}^{0})=\bar{\omega}_{1}^{0}\). Then, we have
where \(\bar{M}_{1}^{2}= -3(\bar{\omega}_{1}^{0})^{2} + (A_{2}^{2} + C_{1}^{2})\), \(\bar{S}_{1}^{2}= 2(A_{1}^{2} + C_{0}^{2})\bar{\omega}_{1}^{0}\), \(\bar{N}_{1}^{2}= -(\bar{\omega}_{1}^{0})^{4} + (A_{2}^{2} + C_{1}^{2})( \bar{\omega}_{1}^{0})^{2}\), \(\bar{T}_{1}^{2}= (A_{1}^{2} + C_{0}^{2})(\bar{\omega}_{1}^{0})^{3} - (A_{3}^{2} + C_{2}^{2}) \bar{\omega}_{1}^{0}\), \(\bar{Q}_{1}^{2}= B_{1}^{2} + D_{0}^{2}\), \(\bar{R}_{1}^{2}= 2B_{0}^{2} \bar{\omega}_{1}^{0}\), \(\bar{P}_{1}^{2} = -(B_{1}^{2} + D_{0}^{2})(\bar{\omega}_{1}^{0})^{2}\), and \(\bar{O}_{1}^{2} = -B_{0}^{2}(\bar{\omega}_{1}^{0})^{3} + (B_{2}^{2} + D_{1}^{2}) \bar{\omega}_{1}^{0} \). Clearly, if condition
is satisfied, then
Thus, we have the following theorem:
Theorem 8
Suppose that \(\tau _{1} > 0\), and \(\tau _{2} = 0\). If equation (39) has positive roots and condition (40) is met, then
(i) when \(\tau _{1} \in (0,\bar{\tau}_{1}^{0})\), \(E_{1}\) is locally asymptotically stable;
(ii) when \(\tau _{1}>\bar{\tau}_{1}^{0}\), \(E_{1}\) is unstable;
(iii) when \(\tau _{1}=\bar{\tau}_{1}^{0}\), system (1) has a Hopf bifurcation at \(E_{1}\).
Case 3.3.4 If \(\tau _{1}>0\), and \(\tau _{2}>0\), then using the similar discussion as for Case 3.2.4, we consider the Hopf bifurcation of \(E_{1}\) for the system (1) with \(\tau _{1}>0\) and \(\tau _{2} \in (0, \bar{\tau}_{2}^{0})\), where \(\bar{\tau}_{2}^{0}\) as given in Case 3.3.2. Let \(\tau _{2} \in (0,\bar{\tau}_{2}^{0})\) be arbitrary but fixed. Then, \(\lambda =i\omega\) (\(\omega >0\)) is a pure imaginary root of the equation (28) if and only if
Separating the real part from the imaginary part in (41) gives
By adding the squares of two equations in (42), we obtain
where \(\bar{l}_{1}(\omega ) = -2C_{0}^{2}\sin ({\omega \tau _{2}}) \); \(\bar{l}_{2}(\omega ) = (A_{1}^{2})^{2} - 2A_{2}^{2} + (C_{0}^{2})^{2} - (B_{0}^{2})^{2} + 2(A_{1}^{2}C_{0}^{2} - C_{1}^{2})\cos ({\omega \tau _{2}}) \); \(\bar{l}_{3}(\omega ) = 2(C_{2}^{2} + A_{2}^{2}C_{0}^{2} -A_{1}^{2}C_{1}^{2} + B_{0}^{2}D_{0}^{2})\sin ({\omega \tau _{2}})\); \(\bar{l}_{4}(\omega ) = (A_{2}^{2})^{2} + (C_{1}^{2})^{2} + 2B_{0}^{2}B_{2}^{2} - 2A_{1}^{2}A_{3}^{2} - 2C_{0}^{2}C_{2}^{2} - (B_{1}^{2})^{2} - (D_{0}^{2})^{2} - 2(A_{1}^{2}C_{2}^{2} - A_{2}^{2}C_{1}^{2} +A_{3}^{2}C_{0}^{2} + B_{1}^{2}D_{0}^{2} - B_{0}^{2}D_{1}^{2})\cos ({\omega \tau _{2}})\); \(\bar{l}_{5}(\omega ) = 2(A_{3}^{2}C_{1}^{2} -A_{2}^{2}C_{2}^{2} -B_{2}^{2}D_{0}^{2} + B_{1}^{2}D_{1}^{2})\sin ({\omega \tau _{2}})\); and \(\bar{l}_{6}(\omega ) = (C_{2}^{2})^{2} + (A_{3}^{2})^{2} - (B_{2}^{2})^{2} - (D_{1}^{2})^{2} + ( 2A_{3}^{2}C_{2}^{2} -2B_{2}^{2}D_{1}^{2})\cos ({ \omega \tau _{2}}) \). Let
Obviously, if \(\bar{l}_{6}(0)<0\), then \(G_{2}(0)<0\). In addition, \(\lim_{\omega \to +\infty} G_{2}(\omega )=+\infty \). Hence, equation (43) has at least one positive root denoted by \(\hat{\omega}_{1}\). Furthermore, substituting \(\hat{\omega}_{1}\) into (43) and solving the resulting equations, we have
where \(M_{9}(\hat{\omega}_{1})=- B_{0}^{2}(\hat{\omega}_{1})^{2} + B_{2}^{2} + D_{1}^{2}\cos ({\hat{\omega}_{1} \tau _{2}}) + D_{0}^{2} \hat{\omega}_{1} \sin ({\hat{\omega}_{1} \tau _{2}})\); \(M_{10}( \hat{\omega}_{1})= B_{1}^{2} \hat{\omega}_{1} - D_{1}^{2}\sin ({ \hat{\omega}_{1} \tau _{2}}) + D_{0}^{2} \hat{\omega}_{1} \cos ({ \hat{\omega}_{1} \tau _{2}})\); \(M_{11}(\hat{\omega}_{1})= A_{1}^{2} ( \hat{\omega}_{1})^{2} - A_{3}^{2} -(-C_{0}^{2}(\hat{\omega}_{1})^{2} + C_{2}^{2})\cos ({\hat{\omega}_{1} \tau _{2}}) - C_{1}^{2}\hat{\omega}_{1} \sin ({\hat{\omega}_{1} \tau _{2}})\); and \(M_{12}(\hat{\omega}_{1})=( \hat{\omega}_{1})^{3} - A_{2}^{2}\hat{\omega}_{1} +(-C_{0}^{2}( \hat{\omega}_{1})^{2} + C_{2}^{2})\sin ({\hat{\omega}_{1} \tau _{2}}) - C_{1}^{2}\hat{\omega}_{1} \cos ({\hat{\omega}_{1} \tau _{2}}) \). Clearly,
Thus, we can define
Next, we consider the roots of equation (28) as a function of \(\tau _{1}\) and the root \(\tilde{\lambda}_{3}(\tau _{1})=\alpha _{3}( \tau _{1})+i\tilde{\omega}_{3}(\tau _{1})\) of equation (28) such that
Substituting \(\tilde{\lambda}_{3}(\tau _{1})\) into equation (28) and taking the derivative with respect to \(\tau _{1}\), we have
Then,
where \(M_{3}^{2}= -3 (\hat{\omega}^{0}_{1})^{2} + A_{2}^{2} + 2C_{0}^{2} \hat{\omega}^{0}_{1}\sin (\hat{\omega}^{0}_{1} \tau _{2}) + C_{1}^{2} \cos (\hat{\omega}^{0}_{1} \tau _{2}) - \tau _{2}(-C_{0}^{2}( \hat{\omega}^{0}_{1})^{2} + C_{2}^{2})\cos (\hat{\omega}^{0}_{1} \tau _{2}) - \tau _{2} C_{1}^{2}\hat{\omega}^{0}_{1} \sin ( \hat{\omega}^{0}_{1} \tau _{2})\); \(S_{3}^{2}= 2 A_{1}^{2}\hat{\omega}^{0}_{1} + 2C_{0}^{2}\hat{\omega}^{0}_{1} \cos (\hat{\omega}^{0}_{1} \tau _{2}) - C_{1}^{2}\sin (\hat{\omega}^{0}_{1} \tau _{2}) + \tau _{2}(-C_{0}^{2}(\hat{\omega}^{0}_{1})^{2} + C_{2}^{2}) \sin (\hat{\omega}^{0}_{1} \tau _{2}) - \tau _{2} C_{1}^{2} \hat{\omega}^{0}_{1} \cos (\hat{\omega}^{0}_{1} \tau _{2})\); \(N_{3}^{2}= -(\hat{\omega}^{0}_{1})^{4} + A_{2}^{2}(\hat{\omega}^{0}_{1})^{2} + (C_{0}^{2}(\hat{\omega}^{0}_{1})^{3} - C_{2}^{2}\hat{\omega}^{0}_{1}) \sin (\hat{\omega}^{0}_{1} \tau _{2}) + C_{1}^{2}(\hat{\omega}^{0}_{1})^{2} \cos (\hat{\omega}^{0}_{1} \tau _{2})\); \(T_{3}^{2}= A_{1}^{2}(\hat{\omega}^{0}_{1})^{3} - A_{3}^{2} \hat{\omega}^{0}_{1} + (C_{0}^{2}(\hat{\omega}^{0}_{1})^{3} - C_{2}^{2}\hat{\omega}^{0}_{1})\cos (\hat{\omega}^{0}_{1} \tau _{2}) - C_{1}^{2}(\hat{\omega}^{0}_{1})^{2}\sin (\hat{\omega}^{0}_{1} \tau _{2})\); \(Q_{3}^{2}= B_{1}^{2} + D_{0}^{2}\cos (\hat{\omega}^{0}_{1} \tau _{2}) -\tau _{2}D_{0}^{2}\hat{\omega}^{0}_{1} \sin (\hat{\omega}^{0}_{1} \tau _{2})- \tau _{2}D_{1}^{2}\cos (\hat{\omega}^{0}_{1} \tau _{2})\); \(R_{3}^{2}= 2B_{0}^{2}\hat{\omega}^{0}_{1} - D_{0}^{2}\sin ( \hat{\omega}^{0}_{1} \tau _{2})- \tau _{2}D_{0}^{2}\hat{\omega}^{0}_{1} \cos (\hat{\omega}^{0}_{1} \tau _{2}) + \tau _{2}D_{1}^{2} \sin (\hat{\omega}^{0}_{1} \tau _{2})\); \(P_{3}^{2} = - B_{1}^{2}(\hat{\omega}^{0}_{1})^{2} - D_{0}^{2}( \hat{\omega}^{0}_{1})^{2} \cos (\hat{\omega}^{0}_{1} \tau _{2})+ D_{1}^{2} \hat{\omega}^{0}_{1}\sin (\hat{\omega}^{0}_{1} \tau _{2})\); and \(O_{3}^{2} = -B_{0}^{2}(\hat{\omega}^{0}_{1})^{3} + B_{2}^{2} \hat{\omega}^{0}_{1}+ D_{0}^{2}(\hat{\omega}^{0}_{1})^{2} \sin ( \hat{\omega}^{0}_{1} \tau _{2})+ D_{1}^{2}\hat{\omega}^{0}_{1}\cos ( \hat{\omega}^{0}_{1} \tau _{2}) \). Furthermore, if condition
is satisfied, then
Thus, by the Hopf bifurcation theory [31], we can obtain the following theorem:
Theorem 9
Suppose that \(\tau _{1} > 0\), and \(\tau _{2} \in (0,\bar{\tau}_{2}^{0})\) with \(\bar{\tau}^{0}_{2}\) as defined in Case 3.3.2. If equation (43) has positive roots and condition (46) is met, then
(i) when \(\tau _{1} \in (0,\hat{\tau}_{1}^{0})\), \(E_{1}\) is locally asymptotically stable;
(ii) when \(\tau _{1}>\hat{\tau}_{1}^{0}\), \(E_{1}\) is unstable;
(iii) when \(\tau _{1}=\hat{\tau}_{1}^{0}\), system (1) has a Hopf bifurcation at \(E_{1}\).
4 Delay optimal control problem
In this section, we will consider the delay optimal control problem in the COVID-19 epidemic and will derive the corresponding necessary optimality conditions.
The outbreak of the COVID-19 epidemic can be controlled by two interventions: (i) reduction and suppression of social contact through masking, home isolation, etc. and (ii) pharmacological interventions, leading to the use of novel treatments that minimize mortality and shorten hospital stays. Let \(u_{1}(t)\) and \(u_{2}(t)\) be the social contact and the pharmaceutical intervention, respectively. Then, the controlled SEIR model in the COVID-19 virus transmission can be formulated as:
where \(T>0\) is a given terminal time. It should be noted that the suppression of social contacts is generally costly, and it will result in economy recession. Pharmacological intervention is also expensive because it requires extensive clinical trials prior to experimental treatment. Therefore, we assume that
where \(u_{i}^{\max}\) (≤1), \(i=1,2\), are the upper bounds for controls \(u_{i}(t)\). Let \(\mathcal{U}\) be the set of all Borel measurable functions \(u_{i}: [0,T]\rightarrow \mathbb{R}\), \(i=1,2\), satisfying constraint (48).
Let \(x(t)= (S(t), E(t), I(t), R(t))^{\top}\), \(u(t)=(u_{1}(t), u_{2}(t))^{\top}\), and the right-hand side of system (47) be \(\tilde{f}(x(t),x(t-\tau _{1}),x(t-\tau _{2}),u(t))\). Then, system (47) with the condition (2) in a finite time horizon \([-\tau ,T]\) can be rewritten as
where \(\varphi (t)=(\varphi _{1}(t),\varphi _{2}(t),\varphi _{3}(t), \varphi _{4}(t))^{\top}\) is the initial condition.
During the COVID-19 spread, social contacts can be mitigated and suppressed by non-pharmacological interventions, such as wearing masks, maintaining social distancing, testing and isolation, and closing businesses, while the length of hospital stay can be minimized by pharmacological interventions and using new treatment modalities. Thus, using both non-pharmacological and pharmacological interventions as the control strategies, our goal is to maximize the number of recovered individuals at the terminal time, as well as to minimize the number of exposed and infected individuals during the time horizon, and the system cost of control measures. Therefore, the cost function in controlling the COVID-19 epidemic can be stated as
where \(q_{i}>0\), \(i = 1,2,3\), \(r_{1}>0\) and \(r_{2}>0\) are weighting coefficients. As a result, we propose the following delay optimal control problem:
To solve Problem (DOCP), we establish the following necessary optimality conditions.
Theorem 10
Let \(u^{*} \in \mathcal{U}\) be the optimal control of Problem (DOCP), and let \(x^{*}=x(\cdot |u^{*})\) be the corresponding solution of system (49). Then, there exists a costate vector \(\lambda ^{*}(t)=(\lambda _{1}^{*}(t),\lambda _{2}^{*}(t),\ldots , \lambda _{4}^{*}(t))^{\top}\) satisfying
with the conditions
where \(\bar{x}^{*}(t)=x_{1}^{*}(t) + x_{2}^{*}(t) + x_{3}^{*}(t) + x_{4}^{*}(t)\); \(W_{1} =\frac{\beta _{1} x_{2}^{*}(t)(1-{u}_{1}^{*}(t))(\bar{x}^{*}(t)-x_{1}^{*}(t))}{(\bar{x}^{*}(t))^{2}}\); \(W_{2} = \frac{\beta _{2} x_{3}^{*}(t)(1-{u}_{1}^{*}(t))(\bar{x}^{*}(t) - x_{1}^{*}(t))}{(\bar{x}^{*}(t))^{2}}\); \(W_{3} =\frac{\beta _{1} x_{1}^{*}(t) (1-{u}_{1}^{*}(t))(\bar{x}^{*}(t) - x_{2}^{*}(t))}{(\bar{x}^{*}(t))^{2}}\); \(W_{4} =\frac{\beta _{2} x_{1}^{*}(t)(1-{u}_{1}^{*}(t))(\bar{x}^{*}(t)-x_{3}^{*}(t))}{(\bar{x}^{*}(t))^{2}}\); \(W_{5} =\frac{\beta _{2} (1-{u}_{1}^{*}(t)) x_{1}^{*}(t) x_{3}^{*}(t)}{(\bar{x}^{*}(t))^{2}}\); and \(W_{6} =\frac{\beta _{1} (1-{u}_{1}^{*}(t)) x_{1}^{*}(t) x_{2}^{*}(t)}{(\bar{x}^{*}(t))^{2}}\). Furthermore, the optimal control \(u^{*}\) can be expressed as
Proof
Using proof similar to Theorem 1 in [34], we can complete the proof. □
5 Numerical simulations
In this section, we will carry out some numerical simulations to verify the stability and the Hopf bifurcation of equilibrium points in Sect. 3. Moreover, Problem (DOCP) will be solved based on the derived necessary optimality conditions in Sect. 4.
5.1 Numerical simulations of stability and Hopf bifurcation
In numerical simulations, two sets of parameter values for system (1) are chosen (Table 1). In addition, the initial condition is chosen as \(\varphi (t)=(99999, 1, 0, 0)^{\top}\) for \(t \in [-\tau ,0]\). For the first set of parameter values (i.e., the second row in Table 1), we obtain \(R_{0}=0.7037(<1)\). In this case, system (1) has a unique disease-free equilibrium point \(E_{0}=(30,0,0,0)\). For \(\tau _{1}=\tau _{2}=0\), we have \(a_{0}^{1}= 1\), \(a_{1}^{1}=1.65\), \(a_{2}^{1}=0.32\) and \(a_{3}^{1}=0\). Obviously, \(a_{1}^{1} a_{2}^{1} > a_{0}^{1} a_{3}^{1}\) and \(E_{0}\) is locally stable from Theorem 2 (Fig. 2). For \(\tau _{1}=0\) and \(\tau _{2}>0\), we obtain \(\omega _{2}^{0}=0.3820\) and \(\tau _{2}^{0}=3.8658\). As shown in Theorem 3, \(E_{0}\) is asymptotically stable when \(\tau _{2} \in (0,3.8658)\), and the Hopf bifurcation occurs when \(\tau _{2}=3.8658\), which are illustrated in Figs. 3 and 4. For \(\tau _{1}>0\) and \(\tau _{2}=0\), we obtain \(\omega _{1}^{0}=0.9840\) and \(\tau _{1}^{0}=1.4683\). The corresponding asymptotic stability and the Hopf bifurcation of \(E_{0}\) are shown in Figs. 5 and 6. For \(\tau _{1}>0\) and \(\tau _{2}=2 \in (0,3.8658)\), we get \(\tilde{\omega}_{1}^{0}=1.2902\) and \(\tilde{\tau}_{1}^{0}=1.2340\). The corresponding asymptotic stability and the Hopf bifurcation of \(E_{0}\) are plotted in Figs. 7 and 8.
For the second set of parameter values (i.e., the third row in Table 1), we obtain \(R_{0}=1.1746(>1)\). In this case, system (1) has a unique endemic equilibrium point \(E_{1}=(8.5135, 0.0787, 0.0708, 0.7236)\). For \(\tau _{1}=\tau _{2}=0\), we have \(a_{0}^{2}=1\), \(a_{1}^{2}=1.16365\), \(a_{2}^{2}=0.4634\) and \(a_{3}^{2}=0.0047\). Obviously, \(a_{1}^{2} a_{2}^{2} > a_{0}^{2} a_{3}^{2}\) and \(E_{1}\) is asymptotically stable from Theorem 6 (Fig. 9). For \(\tau _{1}=0\), \(\tau _{2}>0\), we get \(\bar{\omega}_{2}^{0}=0.2881\) and \(\bar{\tau}_{2}^{0}=4.1563\). Using Theorem 7, we conclude that \(E_{1}\) is asymptotically stable when \(\tau _{2} \in (0,4.1563)\) and that the Hopf bifurcation occurs when \(\tau _{2}=4.1563\) (Figs. 10 and 11). For \(\tau _{1}>0\) and \(\tau _{2}=0\), we obtain \(\bar{\omega}_{1}^{0}=1.0086\) and \(\bar{\tau}_{1}^{0}=1.2326\). The corresponding asymptotic stability and the Hopf bifurcation are plotted in Figs. 12 and 13. For \(\tau _{1}>0\) and \(\tau _{2}=3.5 \in (0,4.1563)\), we obtain \(\hat{\omega}_{1}^{0}=1.8961\) and \(\hat{\tau}_{1}^{0}=0.5437\). The corresponding asymptotic stability and the Hopf bifurcation are shown in Figs. 14 and 15.
5.2 Numerical results of optimal control problem
In this section, we will present numerical results for solving Problem (DOCP). Here, the delayed SEIR system (49) was solved using software package DDE23 in Matlab with delay arguments \(\tau _{1}=0.5\) and \(\tau _{2}=0.7\) and the initial condition \(\varphi (t) = (99999, 1, 0, 0)^{\top}\) [14]. The weight coefficients in the cost functional (50) are chosen as \(q_{1} = 10^{-8}\), \(q_{2} = q_{3} = 10^{-4}\), \(r_{1} = 50\) and \(r_{2} = 1\) [14]. In addition, we choose two sets of parameter values listed in Table 1 for calculating the optimal control strategies. As for the control boundaries, we assume that \(u^{\max}_{1} = 0.07\) and \(u^{\max}_{2} = 0.1\). Note that Problem (DOCP) was solved using forward-backward sweep method [35] in conjunction with the necessary optimality conditions in Theorem 10.
To explore the effects of each control means, we set up the following control scheme.
(i) Strategy A: Social contact only (\(u_{1}\)).
(ii) Strategy B: Pharmaceutical intervention only (\(u_{2}\)).
(iii) Strategy C: Social contact + Pharmaceutical intervention (\(u_{1}\), \(u_{2}\)).
For the first set of parameters (i.e., the second row in Table 1), we solve Problem (DOCP) and obtain the optimal control strategies \(u_{1}^{*}\) and \(u^{*}_{2}\) shown in Fig. 16. The state trajectories corresponding to different control strategies, i.e., \(u_{1}=u_{2}=0\), \(u_{1}=u_{1}^{*}\) and \(u_{2}=0\), \(u_{1}=0\) and \(u_{2}=u_{2}^{*}\), and \(u_{1}=u_{1}^{*}\) and \(u_{2}=u_{2}^{*}\), are plotted in Fig. 17. Similarly, for the second set of parameters (i.e., the third row in Table 1), we also solve Problem (DOCP) and obtain the optimal control strategies \(\bar{u}^{*}_{1}\) and \(\bar{u}_{2}^{*}\) depicted in Fig. 18. The state trajectories corresponding to different control strategies, i.e., \(u_{1}=u_{2}=0\), \(u_{1}=\bar{u}_{1}^{*}\) and \(u_{2}=0\), \(u_{1}=0\) and \(u_{2}=\bar{u}_{2}^{*}\), and \(u_{1}=\bar{u}_{1}^{*}\) and \(u_{2}=\bar{u}_{2}^{*}\), are plotted in Fig. 19.
From Figs. 16 and 18, we can see that, for the social contact \(u_{1}\), it begins with the maximal value of 0.07, keeps the maximal value for a period of time, and then decreases to zero. This is mainly because, in the case of a sudden outbreak of COVID-19, quarantine measures are effective in stopping the spread of the disease by cutting off the route of transmission in the real world. As for the pharmacological intervention \(u_{2}\), since the pharmacological intervention for the infected individuals is not immediately administered, it starts from the minimal value of zero, rapidly increases to the maximal value of 0.1, maintains at this maximal value for a period of time, and ultimately reduces to zero.
From Figs. 17(b), 17(c), 19(b), and 19(c), it can be seen that the number of exposed and infected individuals under no control is the highest, while the lowest number of exposed and infected individuals is under strategy C. Moreover, the implementation of \(u_{1}\) can decrease the number of both exposed and infectious individuals, while the implementation of \(u_{2}\) is particularly effective in reducing the number of infected individuals. Nevertheless, relying solely on one control measure (strategies A and B) or not implementing any control measures at all is less effective than the optimal control strategy C.
For strategy C, we also solved the optimal control problem without delay. The computed optimal control strategies under two sets of parameter values are also illustrated in Figs. 16 and 18, respectively. Under the optimal control strategies in Figs. 16 and 18, the corresponding optimal state trajectories are depicted in Figs. 20 and 21. From Figs. 20 and 21, we observe that the peaks of infected individuals for the cases with \(\tau _{1}=0.5\) and \(\tau _{2}=0.7\) are higher than those without delay. This implies that time delays could aggravate the transmission of COVID-19. As a result, to minimize the number of infections, effective control strategies should be implemented as soon as possible.
6 Conclusion
In this paper, we have studied the dynamics analysis and optimal control of the delayed SEIR model in the COVID-19 epidemic. Two delays representing the incubation and recovery periods in COVID-19 virus transmission have been introduced. A key issue with delays describing the incubation and recovery periods is whether they cause sustained oscillations. This was investigated through Hopf bifurcation analysis. By using the characteristic equations of delay differential equations, the existence of Hopf bifurcations at the disease-free and endemic equilibrium points was established. It has been shown that under some conditions, delays representing the incubation and recovery periods may destabilize the disease-free and endemic equilibrium points and cause the population to fluctuate. From Theorems 3, 4, 5, 7, 8, and 9, we see that thresholds for two delays were identified to characterize the existence of Hopf bifurcations at the disease-free and endemic equilibrium points. In addition, two controls representing the social contact and pharmaceutical intervention have also been introduced. The delay optimal control problem was formulated, and its necessary optimality conditions were exploited to solve the optimal controls. Numerical simulation results indicate that the pharmacological intervention was more effective for hospitalized patients, whereas suppression of social contact reduced the number of exposed individuals. More importantly, the obtained optimal strategies are effective in preventing and controlling the spread of COVID-19 infection.
We note that the effect of vaccination is ignored in systems (1) and (47). Vaccination has played a major role in preventing the spread of COVID-19 [36–38]. Moreover, the fractional derivative can be regarded as the generalization of the traditional integer derivative, which shows many important properties that the integer derivative does not have [39–44]. Many scholars have applied fractional derivative differential equations to study the spread of COVID-19 [45–47]. As a result, it is worthwhile to investigate the fractional optimal control of the COVID-19 epidemic by incorporating vaccination. This will be left for future research work.
Data Availability
No datasets were generated or analysed during the current study.
References
Kow, C.S., Hasan, S.S.: Do sleep quality and sleep duration before or after COVID-19 vaccination affect antibody response? Chronobiol. Int. 38, 941–943 (2021)
Avadhani, A., Cardinale, M., Akintade, B.: COVID-19 pneumonia: what APRNs should know. Nurse Pract. 46, 22–28 (2021)
Coronavirus Cases. https://www.worldometers.info/coronavirus/
Kermack, W.O., McKendrick, A.G.: A contribution to the mathematical theory of epidemics. Proc. R. Soc. Lond. 115, 700–721 (1927)
Krämer, A., Kretzchmar, M., Kricheberg, K.: Modern Infectious Disease Epidemiology. Springer, New York (2010)
Li, M.Y.: An Introduction to Mathematical Modeling of Infectious Diseases. Springer, Cham (2018)
Julien, A., Portet, S.: A simple model for COVID-19. Infect. Dis. Model. 5, 309–315 (2020)
Awasthi, A.: A mathematical model for transmission dynamics of COVID-19 infection. Eur. Phys. J. Plus 138, 285 (2023)
Xu, C., Yu, Y., Ren, G., Sun, Y., Si, X.: Stability analysis and optimal control of a fractional-order generalized SEIR model for the COVID-19 pandemic. Appl. Math. Comput. 457, 128210 (2023)
Nesteruk, I.: Simulations and predictions of COVID-19 pandemic with the use of SIR model. Innov. Biosyst. Bioeng. 4, 110–121 (2020)
Cooper, I., Mondal, A., Antonopoulos, C.G.: A SIR model assumption for the spread of COVID-19 in different communities. Chaos Solitons Fractals 139, 110057 (2020)
Annas, S., Pratama, M.I., Rifandi, M., Sanusi, W., Side, S.: Stability analysis and numerical simulation of SEIR model for pandemic COVID-19 spread in Indonesia. Chaos Solitons Fractals 139, 110072 (2020)
Ahmed, N., Elsonbaty, A., Raza, A., Rfiq, M., Adel, W.: Numerical simulation and stability analysis of a novel reaction-diffusion COVID-19 model. Nonlinear Dyn. 106, 1293–1310 (2021)
Biswas, S.K., Ahmed, N.U.: Mathematical modeling and optimal intervention of COVID-19 outbreak. Quant. Biol. 1, 84–92 (2021)
Kouidere, A., EL Youssoufi, L., Ferjouchia, H., Balatif, O., Rachik, M.: Optimal control of mathematical modeling of the spread of the COVID pandemic with highlighting the negative impact of quarantine on diabetics people with cost-effectiveness. Chaos Solitons Fractals 145, 110777 (2021)
Deressa, C.T., Duressa, G.F.: Modeling and optimal control analysis of transmission dynamics of COVID-19: the case of Ethiopia. Alex. Eng. J. 60, 719–732 (2021)
Ahmed, M., Masud, M., Sarker, M.: Bifurcation analysis and optimal control of discrete SIR model for COVID-19. Chaos Solitons Fractals 174, 113899 (2023)
Guo, Y., Li, T.: Modelling the competitive transmission of the Omecron strain and Delta strain of COVID-19. J. Math. Anal. Appl. 526, 127283 (2023)
Chen, Y., Cheng, J., Jiang, Y., Liu, K.: A time delay dynamic system with external source for the local outbreak of 2019-nCoV. Appl. Anal. 101, 146–157 (2022)
Paul, S., Lorin, E.: Estimation of COVID-19 recovery and decease periods in Canada using delay model. Sci. Rep. 11, 23763 (2021)
Liu, C., Loxton, R., Teo, K.L.: A computational method for solving time-delay optimal control problems with free terminal time. Syst. Control Lett. 72, 53–60 (2014)
Liu, C., Loxton, R., Teo, K.L.: Optimal parameter selection for nonlinear multistage systems with time-delays. Comput. Optim. Appl. 59, 285–306 (2014)
Liu, C., Loxton, R., Teo, K.L.: Switching time and parameter optimization in nonlinear switched systems with multiple time-delays. J. Optim. Theory Appl. 163, 957–988 (2014)
Liu, C., Gong, Z., Teo, K.L., Sun, J., Caccetta, L.: Robust multi-objective optimal switching control arising in 1, 3-propanediol microbial fed-batch process. Nonlinear Anal. Hybrid Syst. 25, 1–20 (2017)
Liu, C., Loxton, R., Lin, Q., Teo, K.L.: Dynamic optimization for switched time-delay systems with state-dependent switching conditions. SIAM J. Control Optim. 56, 3499–3523 (2018)
Liu, C., Loxton, R., Teo, K.L., Wang, S.: Optimal state-delay control in nonlinear dynamic systems. Automatica 135, 109981 (2022)
Smith, H.L.: An Introduction to the Theory of Competitive and Cooperative Systems. Am. Math. Soc., Rhode Island (1995)
Van den Driessche, P., Watmough, J.: Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci. 180, 29–48 (2002)
Cesari, L.: Asymptotic Behavior and Stability Problems in Ordinary Differential Equations. Springer, Berlin (2012)
Bianca, C., Ferrara, M., Guerrini, L.: The Cai model with time delay: existence of periodic solutions and asymptotic analysis. Appl. Math. Inf. Sci. 7, 21–27 (2013)
Hassard, B.D., Kazarinoff, N.D., Wan, Y.H.: Theory and Applications of Hopf Bifurcation. Cambridge University Press, Cambridge (1981)
Li, H., Liu, X., Yan, R., Liu, C.: Hopf bifurcation analysis of a tumor virotherapy model with two time delays. Phys. A, Stat. Mech. Appl. 553, 124266 (2020)
Ruan, S., Wei, J.: On the zeros of a third degree exponential polynomial with applications to a delayed model for the control of testosterone secretion. Math. Med. Biol. 18, 41–52 (2001)
Ray, W.H., Soliman, M.A.: The optimal control of processes containing pure time delays – I necessary conditions for an optimum. Chem. Eng. Sci. 25, 1911–1925 (1970)
Lenhart, S., Workman, J.T.: Optimal Control Applied to Biological Models. CRC Press, London (2007)
Libotte, G.B., Lobato, F.S., Platt, G.M., Neto, A.J.S.: Determination of an optimal control strategy for vaccine administration in COVID-19 pandemic treatment. Comput. Methods Programs Biomed. 196, 105664 (2020)
Li, T., Guo, Y.: Modeling and optimal control of mutated COVID-19 (Delta strain) with imperfect vaccination. Chaos Solitons Fractals 156, 111825 (2022)
ElHassan, A., AbuHour, Y., Ahmad, A.: An optimal control model for Covid-19 spread with impacts of vaccination and facemask. Heliyon 9, e19848 (2023)
Liu, C.Y., Gong, Z., Yu, C., Wang, S., Teo, K.L.: Optimal control computation for nonlinear fractional time-delay systems with state inequality constraints. J. Optim. Theory Appl. 191, 83–117 (2021)
Gong, Z., Liu, C., Teo, K.L., Wang, S., Wu, Y.H.: Numerical solution of free final time fractional optimal control problems. Appl. Math. Comput. 405, 126270 (2021)
Liu, C., Gong, Z.H., Teo, K.L., Wang, S.: Optimal control of nonlinear fractional-order systems with multiple time-varying delays. J. Optim. Theory Appl. 193, 856–876 (2022)
Liu, C., Gong, Z., Wang, S., Teo, K.L.: Numerical solution of delay fractional optimal control problems with free terminal time. Optim. Lett. 17, 1359–1378 (2023)
Liu, C., Yu, C., Gong, Z., Cheong, H., Teo, K.L.: Numerical computation of optimal control problems with Atangana-Baleanu fractional derivatives. J. Optim. Theory Appl. 197, 798–816 (2023)
Liu, C., Zhou, T., Gong, Z., Yi, X., Teo, K.L., Wang, S.: Robust optimal control of nonlinear fractional systems. Chaos Solitons Fractals 175, 113964 (2023)
Panwar, V.S., Sheik Uduman, P.S., Gómez-Aguilar, J.F.: Mathematiacal modeling of coronavirus disease COVID-19 dynamics using CF, and ABC non-singular fractional derivatives. Chaos Solitons Fractals 145, 110757 (2021)
Pandey, P., Chu, Y.M., Gómez-Aguilar, J.F., Jahanshahi, H., Aly, A.A.: A novel fractional mathematicatical model of COVID-19 epidemic considering quarantine and latent time. Results Phys. 26, 104286 (2021)
Guo, Y., Li, T.: Fractional-order modeling and optimal control of a new online game addition model based on real data. Commun. Nonlinear Sci. Numer. 121, 107221 (2023)
Acknowledgements
The supports of National Natural Science Foundation of China (No. 12271307) and Shandong Province Natural Science Foundation of China (No. ZR2023MA054) are gratefully acknowledged.
Author information
Authors and Affiliations
Contributions
C.L. and J.G. wrote the main manuscript and C.L. and J.K. revised the manuscript. All authors reviewed the manuscript
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Liu, C., Gao, J. & Kanesan, J. Dynamics analysis and optimal control of delayed SEIR model in COVID-19 epidemic. J Inequal Appl 2024, 66 (2024). https://doi.org/10.1186/s13660-024-03140-2
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s13660-024-03140-2