Mathematical modeling of intervention and low medical resource availability with delays: Applications to COVID-19 outbreaks in Spain and Italy

1.   Introduction

In the past decade, infectious diseases have been frequently threatening human lives on a large scale. After the first epidemic model of smallpox was proposed by a pioneering mathematician Daniel Bernoulli in the eighteenth century ^[1,2], a huge number of epidemic models have been proposed and studied to uncover the underlying spread mechanisms and to control the spread of infectious diseases. Kermack and Mckendrick ^[3] introduced a seminal SIR compartmental model in 1927 to study the plague disease in Mumbai and succeeded in revealing its epidemiological transition. Since then, mathematical modeling has been growing as a vital tool to suggest public health responses against the spread and transmission of infectious diseases. The SIR and SEIR type models have been performed a crucial role in modern mathematical epidemiology ^[4,5].

Numerous factors significantly affect the dynamical behavior of a particular epidemic model, such as the demographics, recovery or treatment rates, and incidence rates ^[6,7]. In particular, different types of incidence functions can ensure appropriate dynamical characteristics of the associated infectious diseases ^[8,9]. In classical SIR/SEIR models, the general incidence rate $\frac{\beta SI}{N},$ bilinear incidence rate $\beta SI$ and linear recovery rate $\gamma I$ are commonly used. The constant $\beta$ represents the average number of adequate contacts between susceptible and an infected individual per unit time, and $\gamma$ represents the per capita recovery rate of infected individuals. Since the infection risk increases with the increase in the number of infective individuals, therefore, the mobility of individuals probably influences these numbers. When a particular infectious disease appears and spreads in a region/community, intervention strategies play a crucial role in controlling the disease ^[10,11].

The impacts of intervention strategies, such as mask-wearing, border screening, isolation, quarantine, or communications through the mass media (to convey the information to the public about the epidemic/outbreak and probably risk-reducing behavior), play their significant roles in managing effectual interventions to control the spread of disease and expectantly remove epidemic diseases ^{[12,13,14,15]}. For instance, during the outbreak of the H1N1 influenza pandemic in 2009 ^[16,17] and the outbreak of SARS in 2003 ^[18,19], the intervention strategies (such as postponing conferences, closing restaurants/schools, isolating infectives, etc.) were chosen by the Chinese government. Recently, a new pandemic COVID-19 appeared and affected more than 200 countries and territories worldwide. To control the spread of COVID-19, governments of different countries strictly implemented the lockdown ^[20]. Besides lockdown, various governments are also adopting numerous steps and implementing different intervention strategies, such as social distancing, tracing close contacts, washing hands for 30 seconds, etc ^[21]. In many countries, these strategies impressively worked so that the contacts between susceptible and infected individuals per unit time are being decreased sufficiently, and therefore, the reduced incidence rate ^{[22,23,24,25]}. This showed the importance of considering the infection forces that comprises the adaptation of individuals to infection risks under intervention strategies.

In classical epidemic models, the recovery rate of infected individuals is commonly considered to be proportional to the number of infectives ^[26,27]. However, in general, this is unacceptable because the implicit assumption is that treatment resources are sufficient ^[28]. Every society should have a limited capacity of resources for treatment. If medical resources are available in large quantities, the community needs to pay unnecessary costs during non-epidemic periods, and if available resources are low, the community has the risk of a disease outbreak without needed treatment ^[26]. The recovery rate depends on the resources of the health system available to the community, particularly, the capacity of the hospital settings, efficiency, and effectiveness of treatment ^[29,30]. The prime factor is the number of health employees, including nurses, physicians, pharmacists, and other health care workers (HCW). The amenities of the hospital, such as equipment and medical apparatus, medicines, and the number of hospital beds are other noteworthy factors that are important for effective and safe prevention, treatment, and diagnosis of patients ^[28].

Shan et al. ^[28] studied an SIR model with the impact of the number of beds in hospitals. The authors considered the standard incidence rate and nonlinear recovery rate. They studied the complex dynamics of the model including numerous bifurcations such as backward bifurcation, saddle-node bifurcation, Hopf bifurcation, Bogdanov-Takens bifurcation. The study recommends that maintaining an adequate number of beds in the hospital is crucial to control the disease. Li et al. ^[29] investigated the dynamics of an SIR model using nonlinear incidence rate and nonlinear recovery rate. The model exhibited complicated dynamics and suggested that a sufficient number of beds is critical to control the disease. Mu et al. ^[30] proposed an SI-SIR model for avian influenza with a nonlinear recovery rate under available hospital resources. Zhao et al. ^[31] investigated the complex dynamics of Zika virus transmission via mathematical modeling with limited medical resources. In ^{[28,29,30,31]}, the authors considered a similar recovery rate that decreases with a decrease in the number of beds/resources in the hospitals. In the case of low availability of resources for treatment, the recovery rate primarily grows with an increase in the number of infected individuals and approaches the maximum and then starts decreasing. When deliveries of treatment (immunization, medicine, etc.) are depleted due to a large number of infected individuals, the available resources for treatment become very low.

Time delays are inevitable in almost all realistic systems. Introducing time delays in epidemic models can be important, and time delays include the latency period ^[32,33], the delay in recovery ^[34], the delay in media broadcast ^[35], the delay in awareness programs ^[36] etc. Different types of delay differential equations (DDEs) have fruitfully been utilized to study infectious diseases ^[37,38,39]. Greenhalgh et al. ^[40] proposed a mathematical model with multiple time delays, one delay due to the memory disappearance of aware individuals, and the other delay between the rising awareness and the disease appearance. Zhou et al. ^[41] analyzed an SIR model with media coverage incorporating time delay and studied delay induced local and global Hopf-bifurcation. The authors found that the delay in media coverage does not influence the stability of the disease-free equilibrium. In ^[42], the authors proposed a mathematical model with media coverage incorporating time delay. Recently, Liu et al. ^[43] studied the delay induced local Hopf bifurcation of an SEIR model with multiple time delays caused by the latent period and the recovery period.

A common assumption about intervention strategies in mathematical models ^[10,11] is that the impact of intervention strategies on the transmission dynamics is effective, i.e., the number of infected individuals instantly decrease. The recovery rate is assumed to be proportional to the number of infectives, but many researchers argued that this assumption is questionable in reality as we have discussed earlier. As an improvement, we propose an SEIR model with infection force under intervention strategies and a nonlinear recovery function (under low availability of treatment resources), where treatment declines (due to resources limitation) after attaining its maximum value. Further, we incorporate two time delays, the first delay represents the latent period of interventions and the second delay represents the period for curing the infected individuals. The two main objectives of the study are $(i)$ to investigate the impact of low availability of resources for treatment in the presence of intervention strategies, and $(ii)$ to assess the role of time delays in the transmission dynamics. Furthermore, we validate our model and estimate the parameter values based on the available COVID-19 data in Spain and Italy.

The remaining paper is organized as follows. In Section 2, we formulate an SEIR model (both non-delayed and delayed) with infection force under intervention strategies and recovery rate (under low treatment resources). In Section 3, we analyze the non-delayed system (2.2). In Section 4, we analyze the delayed system (2.3). Persistence has been shown in Section 5. In Section 6, numerical simulations have been presented for both non-delayed system (2.2) and delayed system (2.3). In Section 7, a case study has been discussed for COVID-19 outbreaks in Spain and Italy. The paper ends with conclusion and future developments in Section 8.

2.   Mathematical formulation of the model

We start with the classical SEIR model in which the total population is divided into four compartments: susceptible individuals (S), exposed individuals (E), infected individuals (I), and recovered individuals (R). Susceptible individuals are healthy but can be infected via contacting with infectives. Individuals who are exposed to disease pathogen and may be infected via contacts are called exposed individuals. Exposed individuals may or may not develop the disease, and these individuals are typically not infectious. Infected individuals are those who have been infected and are capable of transmitting the disease to susceptible individuals. Recovered individuals are those who were infected but have been healthy now, and they have immune protection for a long time. Under the above assumptions, the SEIR model with incidence rate under intervention strategies and recovery rate under low availability of medical resources is provided by

associated with the state space $\mathbb{R}^{4}_{+} = \left\lbrace(S, E, I, R):S > 0, E > 0, I > 0, R > 0 \right\rbrace$. All parameters are assumed to be positive. Here $A$ is the rate at which new individuals (including immigrants and newborns) enter into the susceptible population; $\mu$ is the natural mortality rate; $q$ is the rate at which exposed individuals become infectious; $\delta$ is the disease induced death rate. The infection force $F(I)$ and recovery function $G(I)$ in system (2.1) are the functions of infective individuals. These two functions play important roles in governing the disease transmission dynamics. There are various nonlinear transmission rates and nonlinear recovery rates proposed by researchers, for instance, ^{[10,11,13,14,28,44,45]} and references therein.

System (2.1) comprises the adaptation of human behavior under intervention strategies. In fact, when a new infectious disease emerges, both the infection probability and contact rate increase since people have less knowledge about the disease. Further, when the number of infective individuals get larger and the disease becomes more serious, psychological factors lead people to modify their behavior and implement suitable measures/intervention to reduce the opportunities of contact and the probability of infection. For instance, $F(I)$ may decrease when the number of infectives increases since the government may tend to contain the infectious disease by intervention strategies. It has also been interpreted as the psychological effect ^[14]. Mathematically, this idea could be modeled as follows: the infection force $F(I)$ increases when $I$ is small and decreases when $I$ is large (shown in Figure 1). For notational convenience, we assume that the infection force $F(I)$ could be factorized into $\frac{\beta I}{f(I)}$, where $\frac{1}{f(I)}$ signifies the impact of intervention strategies on the decrease of effective contact coefficient $\beta$ ^[27]. In the absence of intervention strategies, i.e., $f(I) = 1$, it is notable that the infection force takes the form of well-known bilinear transmission $\beta SI$. Here, we have the following assumptions for $f(I)$:

● $H(1): f(0) > 0$ and $f'(I) > 0$ for $I > 0$,

● $H(2):$ There is a $\xi$ such that $\left(\frac{I}{f(I)}\right)' > 0$ for $0 < I < \xi$ and $\left(\frac{I}{f(I)}\right)' < 0$ for $I > \xi$.

These assumptions explain the impacts of intervention strategies determined by a critical value $\xi :$ if $0 < I < \xi$, the infection force increases, while if $I > \xi$, the infection force decreases.

In case of low availability of medical resources for a large number of infectives, the recovery function primarily grows with increase of the number of infectives and approaches the maximum and then declines (shown in Figure 1). We have the following assumptions for $g(I):$

● $H(3):$ $g(0) > 0$ and $g'(I) > 0$ for $I > 0$,

● $H(4):$ There is a $\eta$ such that $\left(\frac{I}{g(I)}\right)' > 0$ for $0 < I < \eta$ and $\left(\frac{I}{g(I)}\right)' < 0$ for $I > \eta$.

These assumptions explain the impacts of treatment policies determined by a critical value $\eta :$ if $0 < I < \eta$, the recovery function increases, while if $I > \eta$, the recovery function decreases.

Under the above assumptions on the infection force and recovery function, the new SEIR model (2.1) takes the following form:

Furthermore, the functions $f(I)$ and $g(I)$ could be chosen from a combination of any of the following formats:

1. $f(I) \; {\rm{or}} \; g(I) = 1 +\alpha I^p, \quad p > 1,$

2. $f(I)\; {\rm{or}}\; g(I) = \exp(\alpha I),$

3. $f(I)\; {\rm{or}}\; g(I) = 1+I + \alpha I^p, \quad p > 1,$

where $\alpha$ and $p$ are positive constants. For the case of intervention strategies, the parameter $\alpha$ could be understood as the strength of interventions, and for the case of treatment, the parameter $\alpha$ could be understood as limitation on the treatment availability. In particular, if we take $g(I) = 1$, then model (2.2) becomes same as the model discussed in ^[11].

In Figure 2, the infection force has been plotted with respect to infected individuals and the strength of interventions. It depicts that when the strength of interventions is zero, then incidence function always increases. However, in presence of interventions, the incidence function first increases and then decreases. Moreover, the maximum of $F(I)$ becomes smaller with an increase in the strength of interventions. Similarly in Figure 3, the recovery function has been plotted with respect to infected individuals and limitation on treatment resources. We see that when the resource limitation is zero, then the recovery function increases. However, in presence of resource limitations, the recovery function first increases and then decreases. In addition, the maximum of $G(I)$ becomes smaller with an increase in resource limitations.

Delays naturally exist and can play an important role in disease dynamics. We include two time delays in system (2.2), the first delay $\tau_1$ represents the latent period of the intervention strategies and the second delay $\tau_2$ represents the period for curing the patients. Delays in intervention and curing processes are significant in infectious disease modeling. For example, the delay in intervention strategies, including the time duration for individuals' responses to the reported infection and the reporting delay, have been observed in H1N1-2009 ^[35]. The model (2.2) with these time delays is provided by

The initial conditions for the model (2.3) are defined as

where $\theta \in [-\tau, 0]$, $\tau = \max\left\lbrace \tau_{1}, \tau_{2} \right\rbrace$, $(\phi_{1}, \phi_{2}, \phi_{3}, \phi_{4})\in C([-\tau, 0], \mathbb{R}^{4})$, $\phi_{i}(0) > 0, i = 1, 2, 3, 4,$ and $C([-\tau, 0], \mathbb{R}^{4})$ is the Banach space of continuous functions from $[-\tau, 0]$ to $\mathbb{R}^{4}$ equipped with the sup-norm. It could easily be established from the elementary theory of functional differential equations ^[49] that the system (2.3) possesses a unique solution with initial conditions (2.4).

3.   Dynamics of non-delayed system (2.2)

3.1.   Positivity and boundedness

It is essential to show that all the population variables are nonnegative for all $t \geq 0$, which indicates that any trajectory that begins with positive initial condition will stay positive for $t \geq 0$.

Theorem 3.1. The closed region

is positively invariant for system (2.2).

Proof. See Appendix A.

3.2.   Equilibrium analysis

The system (2.2) has the following two equilibria:

1. The disease-free equilibrium $D_0 (\frac{A}{\mu}, 0, 0, 0),$ which always exists. Here we define the basic reproduction number $R_0$ for system (2.2) as

where $\frac{1}{\mu + \delta +\frac{\gamma}{g(0)}}$ is the life expectancy of infectious persons. $\frac{A}{\mu}$ signifies the number of susceptible persons at the starting of the infectious process and $\frac{\beta A}{\mu f(0)}$ represents the value when all the individuals are susceptible. $\frac{q}{\mu +q}$ represents the fraction of exposed individuals who survive in class $E$ and become infected. Hence the basic reproduction number $R_0$ is biologically well interpreted.

2. The endemic equilibrium $D^*(S^*, E^*, I^*, R^*),$ where

and $I^*$ is a unique positive root of the following equation

Let

then we can see that $Q(0) = \frac{A \beta}{\mu f(0)}-\frac{\gamma}{q}\frac{(\mu +q)}{g(0)}-\frac{(\mu +q)(\mu +\delta)}{q} > 0,$ if $R_0 > 1.$ Let $Q_1(I) = \frac{A\beta}{\beta I +\mu f(I)}$ and $Q_2(I) = \frac{(\mu +q)(\mu +\delta)}{q}+\frac{\gamma}{q}\frac{(\mu +q)}{g(I)}$, then we have $\lim\limits_{I \rightarrow \infty}Q_1(I) = 0$, $\lim\limits_{I \rightarrow \infty}Q_2(I) = \frac{(\mu +q)(\mu +\delta)}{q}$, $Q_1^{'}(I) < 0$ and $Q_2^{'}(I) < 0$. We also observe that $Q_1(0) > Q_2(0)$ when $R_0 > 1$. Hence $Q_1(I)$ and $Q_2(I)$ intersect each other only once when $R_0 > 1$, which proves the existence and uniqueness of $D^*$.

We observe that $Q_1(0) < Q_2(0)$ when $R_0 < 1$ and both $Q_1(I)$ and $Q_2(I)$ are monotonically decreasing. Thus, $Q_1(I)$ and $Q_2(I)$ do not intersect each other when $R_0 < 1$. Therefore, there is no endemic equilibrium when $R_0 < 1$.

On the other hand, by choosing $f(I) = 1+\alpha_1 I^2$ and $g(I) = 1+\alpha_2 I^2$, we see the existence of $D^*$ by plotting the two nullclines in SI-plane.

3.3.   Local stability behavior

Theorem 3.2. 1. The disease-free equilibrium $D_0$ is locally asymptotically stable when $R_0 < 1$ and unstable when $R_0 > 1$.

2. System (2.2) undergoes a transcritical bifurcation when $R_0 = 1$ and the endemic equilibrium $D^*$ is locally asymptotically stable when $R_0 > 1$.

Proof. For the proof of Theorem 3.2, interested readers may refer to Appendix A.

The threshold quantity $R_0$ denotes the average number of new disease infection generated by a single infection introduced into a community of susceptible individuals. The result of Theorem 3.2 shows that the disease can be eliminated in the community when the basic reproduction number $R_0 < 1$ and the initial population sizes are within the attraction basin of $D_0$.

3.4.   Global stability behavior

Theorem 3.3. The disease-free equilibrium $D_0$ is globally asymptotically stable when $R_0 < 1.$

Proof. For the proof of Theorem 3.3, one may refer to Appendix A.

Theorem 3.4. The endemic equilibrium $D^*$ is globally asymptotically stable when $R_0 > 1.$

Proof. Refer to Appendix A.

Remark 3.1. 1. One can observe that in the absence of intervention strategies, i.e., $f(I) = 1$, the basic reproduction number is given by $R_0 = \frac{\beta Aq}{\mu (\mu +q)(\mu +\delta +\frac{\gamma}{g(0)})}$. Thus if $f(0) = 1$, which means that if we perform intervention strategies at a suitable level of infection, $R_0$ remains same. However, by our assumption, $f$ is an increasing function, which means $R_0$ is greater than that in the presence of interventions. From Theorem 3.4, $D^*$ is globally asymptotically stable, thus intervention strategies decreases the endemic level of disease.

2. Further, if there is no limitation on resources of treatment for large infective individuals, i.e., $g(I) = 1$, the basic reproduction number is given by $R_0 = \frac{\beta Aq}{\mu f(0) (\mu +q)(\mu +\delta +\gamma)}$. Thus if $g(0) = 1$, which means that if there are low availability of resources for treatment at some level of infection, $R_0$ remains same. However, by our assumption, $g$ is an increasing function, which means $R_0$ would be less than that in case of low availability of resources for treatment. Since from Theorem 3.4, $D^*$ is globally asymptotically stable, therefore, low resources of treatment increase the endemic level of disease. Clearly, it could be concluded that the disease persists along with the low availability of resources for treatment when the basic reproduction number is greater than one.

4.   Dynamics of delayed system (2.3)

4.1.   Positivity and boundedness

Define $L(t) = \min\limits_{ t \geq 0}\left\lbrace S(t), E(t), I(t), R(t)\right\rbrace$. Clearly, $L(0) = \min \left\lbrace S(0), E(0), I(0), R(0) \right\rbrace > 0$. We need to show that $L(t) > 0$ for all $t\geq 0$. Suppose that there exists a $t_0 > 0$ such that $L(t_0) = 0$ and $L(t) > 0$ for all $t \in [0, t_0)$. Here, we need to discuss the following four cases: $(i)$ $L(t_0) = S(t_0);$ $(ii)$ $L(t_0) = E(t_0);$ $(iii)$ $L(t_0) = I(t_0);$ $(iv)$ $L(t_0) = R(t_0)$. We only give the proof of case $(iv)$. The remaining cases could be discussed similarly.

Let $L(t_0) = R(t_0)$. Since $L(t) > 0$ for all $t \in [0, t_0)$, $I(t-\tau_{2}) > 0$,

integrating from $0$ to $t_0$, we obtain

which leads a contradiction. Hence $S(t), E(t), I(t), R(t)$ are positive for all $t \geq 0$. For the boundedness of system (2.3), one can follow Section 3. The closed region $\mho$ defined in (3.1) is also invariant for system (2.3).

4.2.   Disease-free equilibrium and its stability

The first three equations of system (2.3) do not include the variable $R$. Therefore, for simplicity, we ignore the forth equation and analyze the following reduced system containing only first three equations:

Since delay does not affect the number of equilibria, system (4.1) always has a disease-free equilibrium $D_0$ and a unique endemic equilibrium $D^*$, which have already been defined in Section 3. The basic reproduction number $R_0$ has also been defined in Section 3.

Theorem 4.1. If $V_{11}^{2}-V_{13}^{2}-2V_{12} > 0$ and $V_{12}^{2}-V_{14}^{2} > 0$, $D_0$ is locally asymptotically stable when $\tau_{1} > 0, \tau_{2} > 0$ and $R_0 < 1$.

2. $D_0$ is unstable when $R_0 > 1$.

Proof. For the proof of Theorem 4.1, one may refer to Appendix A.

Remark 4.1. From the proof of Theorem 4.1 in Appendix A, we observe that if $\tau_{1} > 0, \tau_{2} = 0$, $D_0$ is locally asymptotically stable when $R_0 < 1$. Further if $\tau_1 = 0$ and $\tau_{2} > 0$, then $D_0$ is not locally asymptotically stable. Thus the result of Theorem 4.1 shows that time delay in the interventions ($\tau_{1}$) does not affect the stability of $D_0$, however, the result for the time delay in the recovery ($\tau_{2}$) suggests that only a restriction on the basic reproduction number $R_0$ will not be sufficient for the local stability of disease free equilibrium.

4.3.   Stability of $D^*$ and local Hopf bifurcation

In this subsection, we discuss the stability of $D^*$ and establish the local Hopf bifurcation of system (4.1). The characteristic equation of system (4.1) at $D^*$ is given by

where

Eq. (4.2) is equivalent to

where

For the stability and Hopf bifurcation, we have the following results:

Remark 4.2. If $\tau_{1} = \tau_{2} = 0$. In this case, system (4.1) becomes the non-delayed system (2.2). We have already discussed this case in Section 3.

Lemma 4.2. ^[54] For the transcendental equation

as $\tau_{1}, \tau_{2}, \dots \tau_{m}$ vary, the sum of the orders of the zeros of $x(\lambda, e^{-\lambda \tau_{1}}, \ldots , e^{-\lambda \tau_{m}})$ in the open right half plane can change only if a zero occurs on or crosses the imaginary axis.

Lemma 4.3. If $M_{11}M_{13}+M_{12}M_{14} > 0$, then $\left(\frac{d\lambda}{d\tau_{1}}\right)^{-1}_{\lambda = i\omega_{1}', \tau_{1} = \tau_{1_j}} > 0\; (j = 0, 1, 2, \dots)$ holds.

Similarly one can prove $\left(\frac{d\lambda}{d\tau_{1}}\right)^{-1}_{\lambda = i\omega_{1}^{\pm}, \tau_{1} = \tau_{1_k}} > 0\; (k = 0, 1, 2, \dots)$ and $\left(\frac{d\lambda}{d\tau_{1}}\right)^{-1}_{\lambda = i\omega_{1}^{1, 2, 3}, \tau_{1} = \tau_{1_l}} > 0\; (l = 0, 1, 2, \dots)$.

Theorem 4.4. For reduced system (4.1), the following conclusions hold:

1. If $c_{11} > 0, c_{12} > 0, c_{13} > 0$, then $D^*$ is locally asymptotically stable for $\tau_{1}\in (0, +\infty)$.

2. If $c_{13} < 0, c_{11} > 0, c_{12} > 0$ or $c_{13} < 0, c_{11} > 0, c_{12} < 0$ or $c_{13} < 0, c_{11} < 0, c_{12} < 0$, there exists a $\xi_{1_j}$ such that $0 < \xi_{1_j}\leq \tau_{1_j}$. Then $D^*$ is locally asymptotically stable for $\xi_{1_j} < \tau_{1} < \tau_{1_j}$ and unstable for $\xi_{1_{j+1}} > \tau_{1} > \tau_{1_j}$. Additionally, system (4.1) undergoes a supercritical Hopf bifurcation at $D^*$ when $\tau_{1} = \tau_{1_j}$.

3. If $c_{11} > 0$, $c_{12} < 0, c_{13} > 0$ or $c_{11} < 0$, $c_{12} > 0, c_{13} > 0$ or $c_{11} < 0$, $c_{12} < 0, c_{13} > 0$, there exists a $\sigma_{1_k}$ such that $0 < \sigma_{1_k}\leq \tau_{1_k}$. Then $D^*$ is locally asymptotically stable for $\sigma_{1_k} < \tau_{1} < \tau_{1_k}$ and unstable for $\sigma_{1_{k+1}} > \tau_{1} > \tau_{1_k}$. Additionally, system (4.1) undergoes a supercritical Hopf bifurcation at $D^*$ when $\tau_{1} = \tau_{1_k}$.

4. If $c_{11} < 0$, $c_{12} > 0, c_{13} < 0$, there exists a $\sigma_{1_l}$ such that $0 < \sigma_{1_l}\leq \tau_{1_l}$. Then $D^*$ is locally asymptotically stable for $\sigma_{1_l} < \tau_{1} < \tau_{1_l}$ and unstable for $\sigma_{1_{l+1}} > \tau_{1} > \tau_{1_l}$. Additionally, system (4.1) undergoes a supercritical Hopf bifurcation at $D^*$ when $\tau_{1} = \tau_{1_l}$.

Lemma 4.5. If $M_{21}M_{23}+M_{22}M_{24} > 0$, then $Re\left(\frac{d\lambda}{d\tau_{2}}\right)^{-1}_{\lambda = i\omega_{2}', \tau_{2} = \tau_{2_j}} > 0\; (j = 0, 1, 2, \dots)$ holds.

Similarly one can also prove $\left(\frac{d\lambda}{d\tau_{2}}\right)^{-1}_{\lambda = i\omega_{2}^{\pm}, \tau_{2} = \tau_{2_k}} > 0\; (k = 0, 1, 2, \dots)$ and $\left(\frac{d\lambda}{d\tau_{2}}\right)^{-1}_{\lambda = i\omega_{2}^{1, 2, 3}, \tau_{2} = \tau_{2_l}} > 0\; (l = 0, 1, 2, \dots)$.

Theorem 4.6. For reduced system (4.1), the following conclusions hold:

1. If $d_{11} > 0, d_{12} > 0, d_{13} > 0$, then $D^*$ is locally asymptotically stable for $\tau_{2}\in (0, +\infty)$.

2. If $d_{13} < 0, d_{11} > 0, d_{12} > 0$ or $d_{13} < 0, d_{11} > 0, d_{12} < 0$ or $d_{13} < 0, d_{11} < 0, d_{12} < 0$, there exists a $\xi_{2_j}$ such that $0 < \xi_{2_j}\leq \tau_{2_j}$. Then $D^*$ is locally asymptotically stable for $\xi_{2_j} < \tau_{2} < \tau_{2_j}$ and unstable for $\xi_{2_{j+1}} > \tau_{2} > \tau_{2_j}$. Additionally, system (4.1) undergoes a supercritical Hopf bifurcation at $D^*$ when $\tau_{2} = \tau_{2_j}$.

3. If $d_{11} > 0$, $d_{12} < 0, d_{13} > 0$ or $d_{11} < 0$, $d_{12} > 0, d_{13} > 0$ or $d_{11} < 0$, $d_{12} < 0, d_{13} > 0$, there exists a $\sigma_{2_k}$ such that $0 < \sigma_{2_k}\leq \tau_{2_k}$. Then $D^*$ is locally asymptotically stable for $\sigma_{2_k} < \tau_{2} < \tau_{2_k}$ and unstable for $\sigma_{2_{k+1}} > \tau_{2} > \tau_{2_k}$. Additionally, system (4.1) undergoes a supercritical Hopf bifurcation at $D^*$ when $\tau_{2} = \tau_{2_k}$.

4. If $d_{11} < 0$, $d_{12} > 0, d_{13} < 0$ there exists a $\sigma_{2_l}$ such that $0 < \sigma_{2_l}\leq \tau_{2_l}$. Then $D^*$ is locally asymptotically stable for $\sigma_{2_l} < \tau_{2} < \tau_{2_l}$ and unstable for $\sigma_{2_{l+1}} > \tau_{2} > \tau_{2_l}$. Additionally, system (4.1) undergoes a supercritical Hopf bifurcation at $D^*$ when $\tau_{2} = \tau_{2_l}$.

Lemma 4.7. If $M_{31}M_{33}+M_{32}M_{34} > 0$, then $Re\left(\frac{d\lambda}{d\tau}\right)^{-1}_{\lambda = i\omega_{3}', \tau = \tau_{j}} > 0\; (j = 0, 1, 2, \dots)$ holds.

Similarly we can also prove $\left(\frac{d\lambda}{d\tau}\right)^{-1}_{\lambda = i\omega_{3}^{\pm}, \tau = \tau_{k}} > 0\; (k = 0, 1, 2, \dots)$ and $\left(\frac{d\lambda}{d\tau}\right)^{-1}_{\lambda = i\omega_{3}^{1, 2, 3}, \tau = \tau_{l}} > 0\; (l = 0, 1, 2, \dots)$.

Theorem 4.8. For reduced system (4.1), the following conclusions hold:

1. If $e_{11} > 0, e_{12} > 0, e_{13} > 0$, then $D^*$ is locally asymptotically stable for $\tau \in (0, +\infty)$.

2. If $e_{13} < 0, e_{11} > 0, e_{12} > 0$ or $e_{13} < 0, e_{11} > 0, e_{12} < 0$ or $e_{13} < 0, e_{11} < 0, e_{12} < 0$, there exists a $\xi_{3_j}$ such that $0 < \xi_{3_j}\leq \tau_{j}$. Then $D^*$ is locally asymptotically stable for $\xi_{3_j} < \tau < \tau_{j}$ and unstable for $\xi_{3_{j+1}} > \tau > \tau_{j}$. Additionally, system (4.1) undergoes a supercritical Hopf bifurcation at $D^*$ when $\tau = \tau_{j}$.

3. If $e_{11} > 0$, $e_{12} < 0, e_{13} > 0$ or $e_{11} < 0$, $e_{12} > 0, e_{13} > 0$ or $e_{11} < 0$, $e_{12} < 0, e_{13} > 0$, there exists a $\sigma_{3_k}$ such that $0 < \sigma_{3_k}\leq \tau_{k}$. Then $D^*$ is locally asymptotically stable for $\sigma_{3_k} < \tau < \tau_{k}$ and unstable for $\sigma_{3_{k+1}} > \tau > \tau_{k}$. Additionally, system (4.1) undergoes a supercritical Hopf bifurcation at $D^*$ when $\tau = \tau_{k}$.

4. If $e_{11} < 0$, $e_{12} > 0, e_{13} < 0$, there exists a $\sigma_{3_l}$ such that $0 < \sigma_{3_l}\leq \tau_{l}$. Then $D^*$ is locally asymptotically stable for $\sigma_{3_l} < \tau < \tau_{l}$ and unstable for $\sigma_{3_{l+1}} > \tau > \tau_{l}$. Additionally, system (4.1) undergoes a supercritical Hopf bifurcation at $D^*$ when $\tau = \tau_{l}$.

Lemma 4.9. If $M_{41}M_{43}+M_{42}M_{44} > 0$, then $Re \left(\frac{d\lambda}{d\tau_{1}}\right) ^{-1}_{\lambda = i\omega_{4}^{*}, \tau_{1} = \tau_{1_j}^{*}} > 0\; (j = 0, 1, 2, \dots)$ holds.

Theorem 4.10. For reduced system (4.1), If $\textbf{H1}$ holds, there exists a $\xi_{4_j}$ such that $0 < \xi_{4_j}\leq \tau_{1_j}^{*}$. Then $D^*$ is locally asymptotically stable for $\xi_{4_j} < \tau_{1} < \tau_{1_j}^{*}$ and unstable for $\xi_{4_{j+1}} > \tau_{1} > \tau_{1_j}^{*}$. Additionally, system (4.1) undergoes a supercritical Hopf bifurcation at $D^*$ when $\tau_{1} = \tau_{1_j}^{*}$.

Remark For the proofs of Lemmas 4.3, 4.5, 4.7, 4.9, and Theorems 4.4, 4.6, 4.8, 4.10, interested readers may refer to Appendix A.

Remark 4.3. Without loss of generality, in Subsection 4.3 we only discussed the case $Re \left(\frac{d\lambda}{d\tau_{j}}\right) ^{-1}_{\lambda = i\omega, \tau = \tau_{j}} > 0\; (j = 1, 2)$. Similarly, for the case of $Re \left(\frac{d\lambda}{d\tau_{j}}\right) ^{-1}_{\lambda = i\omega, \tau = \tau_{j}} < 0\; (j = 1, 2)$, subcritical Hopf bifurcation will appear instead of supercritical Hopf bifurcation. We do not provide details here.

4.4.   Stability and direction of Hopf bifurcation

In Subsection 4.3, we have established that the reduced system (4.1) possesses a family of periodic solutions bifurcating from $D^*$ at the different critical values of delay parameters $\tau_{1}$ and $\tau_{2}$. In this subsection, by using the center manifold theorem and normal form theory ^[57], the properties of the Hopf bifurcation at the critical value $\tau_{1}^{*}$ are determined. Throughout the section, we consider that $\tau_{2}^{*} < \tau_{1}^{*}$, where $\tau_{2}^{*} \in (0, {\tau_{2}}_{j})$. For detailed discussion related to properties of the Hopf bifurcation (direction and stability of Hopf bifurcation), one may refer to Appendix B.

Remark 4.4. The disappearance or appearance of a periodic orbit through a local change in the stability properties of a steady point is called Hopf bifurcation. There are two types of Hopf bifurcation, the bifurcation is known as supercritical if the bifurcated periodic solutions are stable and subcritical if they are unstable. The parameter $\mu_2$ determines the directions of the Hopf bifurcation: if $\mu_2 > 0 \; (\mu_2 < 0)$, then the Hopf bifurcation is supercritical (subcritical).

Theorem 4.11. If $Re \left\lbrace c_1(0)\right\rbrace < 0\; (Re \left\lbrace c_1(0)\right\rbrace > 0)$, then the reduced system (4.1) undergoes a supercritical (subcritical) Hopf bifurcation at $D^*$ when $\tau_{1}$ crosses its critical value $\tau_{1}^{*}$. Furthermore, the bifurcated periodic solutions occurring through Hopf bifurcations are orbitally asymptotically stable on the center manifold if $Re \left\lbrace c_1(0)\right\rbrace < 0$ and unstable if $Re \left\lbrace c_1(0)\right\rbrace > 0$.

4.5.   Global continuation of local Hopf bifurcation

In this subsection, we investigate the global continuation of periodic solutions bifurcating from $D^*$ for system (2.3) for fixed $\tau_2$ in the interval $(0, {\tau_{2}}_{j})$. Throughout this subsection, we follow the notations in Wu ^[59]. For simplification, we denote $\tau = \tau_{1}$. Let $z(t) = (z_1(t), z_2(t), z_3(t)) = (S(t), E(t), I(t))$, then system (2.3) can be rewritten as

where $z_t(\theta) = (z_1(t+\theta), z_2(t+\theta), z_3(t+\theta))^T = (z_{1t}(\theta), z_{2t}(\theta), z_{3t}(\theta))^T \in C([-\tau, 0], \mathbb{R}^3).$ It is clear that if $R_0 > 1$, then system (2.3) has a disease-free equilibrium $(D_0)$ and an endemic equilibrium $(D^*)$. Following Wu ^[59], we define

Lemma 4.12. Assume that $(\bar{z}, \tau, p)$ is an isolated center satisfying (A1-A4) in ^[59]. Denote by $l_{(\bar{z}, \tau, p)}$ the connected component of $(\bar{z}, \tau, p)$ in $\mho$. Then either

$(i)$ $l_{(\bar{z}, \tau, p)}$ is unbounded, or

$(ii)$ $l_{(\bar{z}, \tau, p)}$ is bounded, $l_{(\bar{z}, \tau, p)} \cap M$ is finite, and $\sum\limits_{(\bar{z}, \tau, p) \in l_{(\bar{z}, \tau, p)} \cap M} \gamma_{m}(\bar{z}, \tau, p) = 0$ for all $m = 1, 2, 3, \ldots ,$ where $\gamma_{m}(\bar{z}, \tau, p)$ is the $m^{th}$ crossing number of $(\bar{z}, \tau, p)$, if $m \in J(\bar{z}, \tau, p),$ or it is zero otherwise.

It is well recognized that if the condition (ii) of above lemma is not true, then $l_{(\bar{z}, \tau, p)}$ is unbounded. Thus, if the projections of $l_{(\bar{z}, \tau, p)}$ onto $p$-space and onto $z$-space are bounded, then the projection onto $\tau$-space is unbounded. Further, if we can show that the projection of $l_{(\bar{z}, \tau, p)}$ onto $\tau$-space is away from zero, then the projection of $l_{(\bar{z}, \tau, p)}$ onto $\tau$-space must include interval $[\tau, +\infty)$. Following this concept, we can prove our consequences on the global continuation of local Hopf bifurcation.

Lemma 4.13. If $R_0 > 1$, then all non-constant periodic solutions of system (2.3) with initial conditions (2.4) are uniformly bounded.

Lemma 4.14. If $R_0 > 1,$ then there does not exist any non-constant periodic solution of system (2.3) with period $\tau$.

Proof. Suppose there exists a non-constant periodic solution of system (2.3) with period $\tau$. Then system (2.2) also has a non-constant periodic solution. Both systems (2.2) and (2.3) have same equilibria, i.e., $D_0$ and $D^*$. Note that E-axis and I-axis are the invariable manifolds of system (2.2) and the orbits of system (2.2) do not cross each other. Thus, there is no solution that intersects the coordinate axis.

On the other hand, if system (2.2) has a periodic solution, then there must be an equilibrium in its interior ^[41] and $D_0$ is located on the coordinate axis. Thus, it can be concluded that the periodic orbit of system (2.2) must lie in the first octant. From Theorem 3.4, the positive equilibrium is globally asymptotically stable in $\mathbb{R}_+^4$, thus, there is no periodic orbit in the first quadrant. This completes the proof.

Theorem 4.15. Let $\omega_{4}^{*}$ and ${\tau_{1}^{*}}_{j}\; (j = 0, 1, 2, \ldots )$ be defined in case 5 in subsection 4.3. Then for each $\tau > {\tau_{1}^{*}}_{j}\; (j\geq 1)$, system (2.3) has at least $j+1$ periodic solutions.

Proof. It is sufficient to prove that the projection of $l_{(D^*, {\tau_{1}^{*}}_{j}, 2\pi / \omega_{4}^{*})}$ onto $\tau$-space is $[\tau, +\infty)$ for each $j > 0,$ where $\bar{\tau} \leq {\tau_{1}^{*}}_{j}$. The characteristic matrix of (4.4) at positive equilibrium $D^*$ takes the following form

where $J_{1i}, i = 1, 2, 3, 4, 5, 6, 7, 8, 9$ and $J_{20}$ are defined in subsection 4.3. From the discussion of local Hopf bifurcation, it is easy to verify that $(D^*, {\tau_{1}^{*}}_{j}, 2\pi / \omega_{4}^{*}), \; j = 1, 2, ...$ are isolated centers. There exists $\epsilon > 0, \delta > 0$ and a smooth curve $\lambda : ({\tau_{1}^{*}}_{j}-\delta, {\tau_{1}^{*}}_{j}+\delta) \rightarrow C$, such that $det(\Delta (\lambda(\tau))) = 0, \vert \lambda(\tau) - \omega_4^* \vert < \epsilon$ for all $\tau \in [{\tau_{1}^{*}}_{j}-\delta, {\tau_{1}^{*}}_{j}+\delta]$ and

Let

It is easy to verify that on $[{\tau_{1}^{*}}_{j}-\delta, {\tau_{1}^{*}}_{j}+\delta]\times \partial \Omega_{\epsilon, 2\pi/\omega_{4}^*}$,

if and only if $\eta = 0, \tau = {\tau_{1}^{*}}_{j}, p = 2\pi/\omega_{4}^*$. Thus, the hypotheses (A1-A4) in ^[59] are satisfied. Moreover, if we define

then we have the crossing number of isolated center $\Big(D^*, {\tau_{1}^{*}}_{j}, \frac{2\pi}{\omega_{4}^*} \Big)$ as follows:

By Theorem 3.2 in Wu ^[59], we conclude that the connected $l_{(D^*, {\tau_{1}^{*}}_{j}, 2\pi/\omega_{4}^*)}$ through $(D^*, {\tau_{1}^{*}}_{j}, 2\pi/\omega_{4}^*)$ in $\mho$ is nonempty. We have

Hence $l_{(D^*, {\tau_{1}^{*}}_{j}, 2\pi/\omega_{4}^*)}$ is unbounded.

From Eq. (.34), we see that, for $j\geq 1, 2\pi / \omega_{4}^* < {\tau_{1}^{*}}_{j}$. Then, we are in a position to show that the projection of $l_{(D^*, {\tau_{1}^{*}}_{j}, 2\pi/\omega_{4}^*)}$ onto $\tau$-space is $[\bar{\tau}, +\infty),$ where $\bar{\tau} < {\tau_{1}^{*}}_{j}$. Clearly, it follows from the proof of Lemma 4.14 that system (2.3) with $\tau = 0$ has no non-constant periodic solution. Hence, the projection of $l_{(D^*, {\tau_{1}^{*}}_{j}, 2\pi/\omega_{4}^*)}$ onto $\tau$-space is away from zero.

For a contradiction, we suppose that the projection of $l_{(D^*, {\tau_{1}^{*}}_{j}, 2\pi/\omega_{4}^*)}$ onto $\tau$-space is bounded. This means that the projection of $l_{(D^*, {\tau_{1}^{*}}_{j}, 2\pi/\omega_{4}^*)}$ onto $\tau$-space is included in an interval $(0, \tau^{*}).$ Note that $2\pi/\omega_{4}^* < {\tau_{1}^{*}}_{j}$ and applying Lemma 4.14, we have $p < \tau^{*}$ for $(z, \tau, p)$ belonging to $l_{(D^*, {\tau_{1}^{*}}_{j}, 2\pi/\omega_{4}^*)}$. This implies that the projection of the connected component $l_{(D^*, {\tau_{1}^{*}}_{j}, 2\pi/\omega_{4}^*)}$ onto $p$-space is bounded. In addition, from Lemma 4.13, we obtain that the projection of $l_{(D^*, {\tau_{1}^{*}}_{j}, 2\pi/\omega_{4}^*)}$ onto $z$-space is bounded if the projection of $l_{(D^*, {\tau_{1}^{*}}_{j}, 2\pi/\omega_{4}^*)}$ onto $\tau$-space is bounded. Thus, the connected component $l_{(D^*, {\tau_{1}^{*}}_{j}, 2\pi/\omega_{4}^*)}$ crossing through $(D^*, {\tau_{1}^{*}}_{j}, 2\pi/\omega_{4}^*)$ is bounded, which is a contradiction. This implies that the projection of $l_{(D^*, {\tau_{1}^{*}}_{j}, 2\pi/\omega_{4}^*)}$ onto $\tau$-space is $[\bar{\tau}, +\infty)$ for each $j\geq 1$, where $\bar{\tau} < {\tau_{1}^{*}}_{j}$. This completes the proof.

4.6.   Estimation of the length of delay to preserve stability

In this subsection, we determine the length of delay to preserve the stability of system (4.1) by the Nyquist criterion ^[55]. We consider the system (4.1) and space of all real-valued continuous functions on $[-\tau, +\infty)$ satisfying the initial conditions on $[-\tau, 0]$. We use the transformations $S = S^*+u_{1}, E = E^*+u_2, I = I^*+u_{3}$ and obtain the following linearized system:

where

Now taking Laplace transform of system (4.6), we obtain

where $K_{1}(s) = \int_{-\tau_{1}}^{0}e^{-s\tau_{1}}u_{3}(t)dt$, $K_{2}(s) = \int_{-\tau_{2}}^{0}e^{-s\tau_{2}}u_{3}(t)dt$, and $\bar{u}_{i}(s)$ are the Laplace transforms of $u_{i}(t), i = 1, 2, 3,$ respectively. The characteristic equation of Eq. (4.7) is

where $P_{1i}, i = 1, 2, 3, 4, 5, 6, 7, 8$ are given in subsection 4.3. The conditions for local asymptotic stability of $D^*$ are given by

where $\eta_{0}$ is the smallest positive root of Eq. (4.9).

We have already revealed that $D^*$ is locally asymptotically stable when $\tau_{1} = \tau_{2} = 0$. Hence, by continuity, all eigenvalues will have negative real parts for sufficiently small $\tau_{1} > 0, \tau_{2} > 0$ provided one can assure that no eigenvalues with positive real parts bifurcates from infinity as $\tau_{1}$ and $\tau_{2}$ increase from zero. This can be proved by Butler's lemma ^[56]. Now, we consider the conditions of Case 5, i.e., $\tau_{1} > 0$ and $\tau_{2}$ is in stable interval. Eq. (4.9) and (4.10) gives

If Eq. (4.11) and (4.12) hold simultaneously, these are sufficient conditions to assure stability. We shall use them to obtain an estimation of the length of delay. Our goal is to discover an upper bound $\eta^{+}$ on $\eta_{0}$ free of $\tau_{1}$ and $\tau_{2}$, and then estimate $\tau_{1}$ so that Eq. (4.12) holds for all values of $\eta$, $0\leq \eta \leq \eta^{+}$ and hence in particular at $\eta = \eta_{0}.$ From (4.11), we obtain

Maximize $P_{13}+(P_{18}-P_{16}\eta_{0}^{2})\cos (\eta_{0}\tau_{2})+P_{17}\eta_{0} \sin (\eta_{0}\tau_{2})+ P_{15} \cos(\eta_{0}\tau_{1})+P_{14}\eta_{0}\sin (\eta_{0}\tau_{1})$ subject to $|\sin(\eta_{0} \tau_{1})|\leq 1, |\sin(\eta_{0} \tau_{2})|\leq 1, |\cos(\eta_{0} \tau_{1})|\leq 1, |\cos(\eta_{0} \tau_{2})|\leq 1$, we obtain

Assume that Eq. (4.14) has a positive root such $\eta^{+}$, then obviously from inequality (4.14), we have $\eta_{0}\leq \eta^{+}$.

Now, rearranging the inequality (4.12), using Eq. (4.11), and assuming $|\sin(\eta_{0} \tau_{2})|\leq 1,$ $|\cos(\eta_{0} \tau_{2})|\leq 1$, we obtain

Using the bounds, we obtain

and

Hence, from inequality (4.15), we obtain

where

Hence, if

then the stability is preserved for $0 \leq \tau_{1} < \tau_{1}^{+}$.

Similarly, one can follow the above calculation to determine the length of delay $\tau_{2}$ where the stability is preserved, taking $\tau_{1}$ in its stable interval.

4.7.   Stability switching curves and crossing directions

To study stability switching, we need to seek purely imaginary characteristic roots. For this, we assume $\lambda = i\omega (\omega > 0)$ and substituting this in Eq. (4.3), we obtain

with

Since $|e^{-i\omega \tau_{2}}| = 1$, we have

which is equivalent to

After simplification, we obtain

Thus,

where $K_1(\omega) = {\rm{Re}}(m_1 \bar{m}_2)$ and $K_2(\omega) = {\rm{Im}}(m_1 \bar{m}_2)$.

If there is some $\omega$ such that $K_1(\omega)^2 + K_2 (\omega)^2 = 0,$ then

The right hand side of Eq. (4.20) is zero with any $\tau_{1}$, and

Therefore, if there is an $\omega$ such that both (4.21) and (4.22) are satisfied, then all $\tau_{1} \in \mathbb{R}_+$ are solution of (4.19).

If $K_1(\omega)^2 + K_2 (\omega)^2 > 0,$ then there exists some continuous function $\chi_1(\omega)$ such that

Indeed,

Therefore, Eq. (4.20) becomes

Obviously, a sufficient and necessary condition for the existence of $\tau_{1} \in \mathbb{R}_+$ satisfying the above equation is

Denote $\Omega^1$ to be $\omega \in \mathbb{R}_+$ satisfying (4.24). We note that (4.24) includes the case $K_1^2 + K_2^2 = 0$ which leads to (4.21) and (4.22).

Let

We have

By substituting the value of $\tau_{1}(\omega)$ given by (4.25), into (4.18) and we obtain an explicit formula for $\tau_{2}(\omega)$ unconditionally with each $\omega \in \Omega^1,$ i.e.,

Hence, the stability crossing curves are

Another way to find $\tau_{2}$ is similarly to the analysis of $\tau_{1}$, which gives

where

with the condition on $\omega$:

which defines a region $\Omega^2$.

By squaring both sides of the conditions (4.24) and (4.29), one can observe that (4.24) and (4.29) are equivalent. Hence

and $\Omega$ is called the crossing set.

Lemma 4.16. $\Omega$ consists of a finite number of intervals of finite length.

Proof. The theoretical proof of this lemma can easily be established similarly as Lemma 3.2 of ^[64]. Numerically, we show that $F(\omega) = (|m_3|^2 - |m_1|^2 - |m_2|^2)^2 - 4 (K_1^2 + K_2^2)$ has finite number of roots on $\mathbb{R}_+$ in Figure 21. If $F(0) > 0,$ then $F(\omega)$ has roots $0 < a_1 < b_1 < a_2 < b_2 < ... < a_N < b_N < +\infty$ and

If $F(0) < 0,$ then $F(\omega)$ has roots $0 < b_1 < a_2 < b_2 < ... < a_N < b_N < +\infty$ and

For any $\Omega_n$, we have a restriction on the range of $\chi_i(\omega), i = 1, 2$. We require $\chi_i(\omega)$ to be the smallest continuous branch with the property that there exists an $\omega_i \in \Omega_n,$ such that

Therefore, $k_i$ has a lower bound, denoted by $M_{i, n}.$

By (4.26), one can obtain either $\tau_{2, k_2}^{+}$ or $\tau_{2, k_2}^{-}$ (but not both) for a given $\tau_{1, k_1}^{+}$, and similarly for $\tau_{1, k_1}^{-}$. By knotty computation, one can verify that when $\tau_{1} = \tau_{1, k_1}^{+}(\omega),$ we have $\tau_{2} = \tau_{2, k_2}^{-}(\omega),$ and when $\tau_{1} = \tau_{1, k_1}^{-}(\omega),$ we have $\tau_{2} = \tau_{2, k_2}^{+}(\omega).$ Therefore

$\Gamma$ defined in (4.30) is the set of all stability switching curves on the $(\tau_{1}, \tau_{2})$-plane and $\Gamma_{k_1, k_2}^{\pm n}$ is continuous in $\mathbb{R}^2$.

For the crossing directions, we directly use Theorem 4.1 of ^[64] which is shown in Figure 22.

5.   Persistence

In this section, we prove the persistence result for both systems (2.2) and (2.3). Biologically, persistence means the survival of all populations in the future time. Mathematically, persistence means positive solutions do not have omega limit points on the boundary of the non-negative set. A population $x(t)$ is called uniformly persistent if there is an $\delta > 0,$ free of $x(0) > 0$ such that $\lim\limits_{ t \rightarrow \infty} x(t) > \delta$. A system is called uniformly persistent if each component persists uniformly.

Theorem 5.1. Both systems (2.2) and (2.3) are uniformly persistent in the interior of $\mho$ whenever $R_0 > 1$.

Proof. Refer the Appendix C.

6.   Numerical simulations

In this section, we present some numerical evaluations to verify our theoretical results and to exhibit different dynamics of systems (2.2) and (2.3).

6.1.   Model system (2.2)

First, we numerically validate the theoretical results obtained for non-delayed system (2.2). We consider $f(I) = 1+\alpha_1 I^2$ and $g(I) = 1+\alpha_2 I^2$ and set the following parameter values (there may exist another set of parameter values):

In this way, the system (2.2) can be rewritten as

In this case, $R_0 = 0.7619 < 1$ and the system (6.2) has a disease-free equilibrium $D_0(10.00, 0, 0, 0)$, which is globally asymptotically stable as shown in Figure 5(a). We observe that curve $S$ starts up very quickly, attains its maximum value and stabilize. The curves $E, I, R$ go down to zero.

Further, we increase the incidence rate as $\beta = 0.5.$ By simple calculation, we obtain $R_0 = 1.9048 > 1$ and endemic equilibrium $D^*(5.313, 3.348, 0.896, 0.353)$, which is globally asymptotically stable as shown in Figure 5(b). We observe that curve $S$ starts up very quickly, attains its maximum value and then goes down and stabilize at a certain value. The curves $E, I, R$ go up to a certain value and stabilize. We plot a transcritical bifurcation diagram in Figure 6.

In Figure 7, we show that by increasing the value of recovery rate ($\gamma$), the level of infected population decreases and by increasing the value of contact rate ($\beta$), the level of infected population increases. From this observation, it could also be concluded that if we increase the recovery rate (by increasing the availability of resources for treatment), the level of infected population could be reduced. Similarly, if we decrease the contact rate (by increasing the intervention strategies), the level of infected population could be reduced. In Figure 8, we plot the infected population with respect to time for different values of $\alpha_1$ and $\alpha_2$. We observe that for greater values of $\alpha_1$, the trajectory of the infected population settles down at a lower level, which means by imposing more intervention, the infection level can be reduced. For higher values of $\alpha_2,$ the infected population settles down at a high level, which means the level of infection increases with the limitation of treatment. In Figure 9, we show that the value of $R_0$ with respect to $\beta$ and $\gamma$. Here we observe that $\beta$ increases the value of $R_0$ and $\gamma$ decreases the value of $R_0$, which means the disease can be extinct by reducing the value of $\beta$ (by reducing contact between susceptible and infected individuals) and by increasing the value of $\gamma$ (by increasing the treatment). From Figure 9, it can also be observed that $\beta$ has more influence than $\gamma$ on $R_0$.

In Figure 10(a)(b), 11(a)(b) and 12(a)(b), we show the combined effect of various parameters on endemic level. Figure 10(a) represents the variation of endemic level with respect to two parameters $\beta$ and $\gamma$. We observe that by decreasing the value of contact rate $\beta$ and by increasing the value of recovery rate ($\gamma$), endemic level could be reduced. Figure 10(b) represents the variation of endemic level with respect to two parameters $\beta$ and $\alpha_1$. We observe that the parameters $\beta$ and $\alpha_1$ have equal impact on endemic level. From this observation, it can be concluded that by increasing the interventions and reducing the contact rate between susceptible and infected individuals, endemic level can be reduced. Figure 11(a) shows that by increasing the value of recovery rate $\gamma$ and increasing the strength of interventions $(\alpha_1)$, number of infected individuals could be reduced. Figure 11(b) demonstrates that by decreasing the value of contact rate $\beta$ and reducing the limitation on resources availability $(\alpha_2)$, number of infected individuals could be reduced. Figure 12(a) depicts the variation of endemic level with respect to two parameters $\gamma$ and $\alpha_2$. We note that by increasing the recovery rate ($\gamma$) and by reducing the limitation on resources availability ($\alpha_2$), endemic level can be reduced. Figure 12(b) declares that number of infected individuals could be reduced by increasing the value of $\alpha_1$ (by increasing the strength of interventions) and decreasing the value of $\alpha_2$ (by decreasing the limitation on resources availability). These results are suitable in real life when the disease spread fast and the infection level is high, for example, COVID-19, H1N1 influenza, SARS. Our study suggests that treatment resources should be available with intervention strategies.

6.2.   Model system (2.3)

In this section, we perform some numerical simulations to verify our theoretical results of system (2.3) in Section 3. We set the parameter values as

and then the system (2.3) can be rewritten as

By simple calculation, we obtain that $R_0 = 64.9351$ and the endemic equilibrium $D^*(248,$ $362,$ $32.93,$ $68.12)$.

For $\tau_{2} = 0, \tau_{1} > 0,$ by simple calculations, we obtain the critical value of $\tau_{1} = 3.2.$ By Theorem 4.4, we can deduce that $D^*$ is locally asymptotically stable when $\tau_{1} = [0, 3.2)$ and a supercritical Hopf bifurcation occur at the critical value $\tau_{1} = 3.2$. Figure 13 also depicts that $D^*$ is locally asymptotically stable for $\tau_{1} = 2 < \tau_{1}^{*} = 3.2$ and bifurcated periodic solutions feasible for $\tau_{1} = 10 > \tau_{1}^{*} = 3.2$.

For $\tau_{1} = 0, \tau_{2} > 0,$ we obtain the critical value of $\tau_{2} = 3.74.$ By Theorem 4.6, we can deduce that $D^*$ is locally asymptotically stable when $\tau_{2} = [0, 3.74)$ and a supercritical Hopf bifurcation occur at the critical value $\tau_{2} = 3.74$. It has been presented in Figure 14 that $D^*$ is locally asymptotically stable for $\tau_{2} = 1 < \tau_{2}^{*} = 3.74$ and bifurcated periodic solutions feasible for $\tau_{2} = 10 > \tau_{1}^{*} = 3.74$.

For $\tau_{1} = \tau_{2} = \tau > 0,$ we obtain the critical value of $\tau = 0.63.$ By Theorem 4.8, we can deduce that $D^*$ is locally asymptotically stable when $\tau = [0, 0.63)$ and a supercritical Hopf bifurcation occur at the critical value $\tau = 0.63$. It has been shown in Figure 15 that $D^*$ is locally asymptotically stable for $\tau = 0.5 < \tau^{*} = 0.63$ and bifurcated periodic solutions feasible for $\tau = 0.65 > \tau^{*} = 0.63$.

For $\tau_{1} > 0, \tau_{2} \in (0, 3.74],$ by simple calculations we obtain the critical value of $\tau_{1} = 0.05.$ By Theorem 4.10, we can deduce that $D^*$ is locally asymptotically stable when $\tau_{1} = [0, 0.05)$ and a supercritical Hopf bifurcation occur at the critical value $\tau_{1} = 0.05$. Figure 16 illustrates that $D^*$ is locally asymptotically stable for $\tau_{1} = 0.04 < \tau_{1}^{*} = 0.05$ and bifurcated periodic solutions feasible for $\tau_{1} = 0.16 > \tau_{1}^{*} = 0.05$.

Furthermore, we plot bifurcation diagrams using the time delays as bifurcation parameters for different cases in Figure 17 and obtain that $D^*$ is locally asymptotically stable before the critical values of delays and bifurcated periodic solutions occur after the critical values of delays. We also plot global branches of Hopf bifurcation in Figure 18 and obtain the global existence of periodic solutions.

Figure 19 represents the changes in critical value of $\tau_{1}$ with respect to $\beta$ and $\alpha_1$. The critical value of $\tau_{1}$ decreases with an increase in $\beta$ while increases with an increase in $\alpha_1$. Figure 20 represents the changes in critical value of $\tau_{2}$ with respect to $\gamma$ and $\alpha_2$. The critical value of $\tau_{2}$ increases with an increase in $\gamma$ while decreases with an increase in $\alpha_2$. At the same time, we give the variation values of $\tau_{1}$ and $\tau_{2}$ with respect to parameters $\beta, \alpha_1, \gamma$ and $\alpha_2$ in Table 1 and 2. These results indicate that these parameters also significantly influence the dynamics of the system with time delays such as critical value of bifurcation parameters.

For stability switching curves, we obtain that $F(0) = - 158.96 < 0$ and $F(\omega)$ has three roots $b_1 = 0.01575, a_2 = 0.02216,$ and $b_2 = 0.0667$ (see Figure 21). Thus

We have

Since $(\delta_{1}^{a}, \delta_{2}^{a}) \neq (\delta_{1}^{b}, \delta_{2}^{b})$, therefore by Theorem 3.1 of ^[64], $\Gamma$ is a set of continuous curves with their two end points on the axises (class Ⅱ) (see Figure 22).

We know that for system (2.3), the interior equilibrium is stable when $\tau_{1} = \tau_{2} = 0,$ i.e., no characteristic roots have positive real parts. Hence from Figure 22, the interior equilibrium is stable if and only if $(\tau_{1}, \tau_{2})$ is on the small bottom-left region of the $(\tau_{1}, \tau_{2})$-plane. As $\tau_{1}$ and $\tau_{2}$ increase, there is a trend that there are more and more characteristic roots with positive real parts.

7.   Application of the study to COVID-19 in Spain and Italy

Recently, a novel coronavirus (namely COVID-19) spread worldwide. Several countries are in the list of top countries, including the USA, India, Brazil, Russia ^[61]. But early on the disaster, Spain and Italy emerged as a major global coronavirus hotspot. Both countries have implemented strict lockdown from March 2020. Restrictions started in Spain on March 14 and in Italy on March 9 ^[62]. Following this, the daily new cases are decreased in Spain and Italy after March 2020 ^[61,63]. In Spain, the maximum number of reported cases was observed on April 2, 2020, which implies the number of infections was on an upward trend till April 2, 2020. After that day, the number of daily cases was observed to decrease asymptotically after imposing interventions, which is the most significant interpretation of the mathematical model. A similar case can be seen in Italy, the maximum number of reported cases was observed on March 21, 2020, after which the number of daily cases was observed to decrease asymptotically. Without loss of generality, we have chosen the data of Spain and Italy during March 1, 2020 to May 31, 2020. We have plotted the bar graph of new daily cases in Spain and Italy from March 1, 2020 to May 31, 2020, i.e., the real cases for 92 days in Figure 23. We calibrate the system (2.3) to fit the new daily cases and cumulative cases of COVID-19 in Spain and Italy and estimate some parameters of the system (2.3) by the least square method using real data of Spain and Italy from March 1, 2020 to May 31, 2020, which is shown in Figure 25. The cumulative number of cases is represented by $C(t),$ which is given as

and the daily new cases are given by $C_{d}(t)$, where

The demographic effects are not considered in this discussion because of the short epidemic time scale in comparison to demographic time scale, i.e., we take $A = 0, \mu = 0$ here. Other estimated parameters and initial conditions for system (2.3) are given in Table 3 and Table 4, respectively.

The basic reproduction number $R_0$ can be explicitly calculated as $R_0 = \frac{\beta N}{\gamma + \delta}$. Based on the parameter values in Table 3, we calculated the value of $R_0$ for both countries Spain and Italy are $5.3391$ and $9.7921$, respectively.

7.1.   Impact of intervention strength and time delays on the early phase of COVID-19 in Spain and Italy

This section aims to compare the outcomes of cumulative cases over time in Spain and Italy for different values of different parameters. Here, using the parameter values given in Table 3 as the baseline values, we have shown that how a model could be used to observe the effect of changing strength of interventions and time delays on epidemic dynamics. In particular, we illustrate the impact of the strength of interventions ($\alpha_1$), time delay in intervention ($\tau_1$), and time delay in recovery ($\tau_2$) on the early phase of pandemic COVID-19. In mathematical terms, we explore how our model dynamics depends on the parameters $\alpha_{1}, \tau_{1},$ and $\tau_{2}$. We compare the total cases generated by the model for different values of strength of interventions $(\alpha_1)$ and time delays $\tau_{1},$ and $\tau_{2}.$ We observe that if the strength of interventions and time delays in the model alter, total cases generated via models are different.

From Figure 26(a), we observe that the strength of intervention $(\alpha_{1})$ could have reduced the total number of cases at the early stage of epidemic in Italy. We note that Italy had recorded $2.32 \times 10^5$ total number of cases by May 31, 2020. We analyze that if the strength of interventions ($\alpha_1$) had increased, the total cases of COVID-19 would have decreased than that the actual cases. Figure 26(a) shows that if the strength of interventions had increased by $50 \% (\alpha_{1} = 3\times 10^{-8})$ from its baseline value $(\alpha_{1} = 2\times 10^{-8})$, then the cumulative cases would have decreased by $12 \% (2.04 \times 10^5)$ and if the strength of interventions had increased by $100 \% (\alpha_{1} = 4\times 10^{-8})$ from its baseline value, then the cumulative cases would have decreased by $18 \% (1.89 \times 10^5)$. If the strength of interventions had decreased by $50 \% (\alpha_{1} = 1\times 10^{-8})$ from its baseline value, then the cumulative cases would have increased by $34 \% (3.12 \times 10^5).$ We observe that the strength of interventions would have been beneficial in the early phase of COVID-19 (measured in terms of total confirmed cases of COVID-19).

We also illustrate that how time delays in intervention and recovery influence the outcomes of early cumulative cases in Italy. The results in Figure 26(b) show that the sooner the interventions are imposed, the greater reduction in the cumulative number of cases. We analyzed that if time delay $\tau_1$ increased by 2 times $(\tau_{1} = 2)$ by its baseline value $(\tau_{1} = 1)$, then cumulative number of cases would have been increased by $8 \% (2.52 \times 10^5)$ and if $\tau_1$ increased by 3 times $(\tau_{1} = 3)$ by its baseline value, then cumulative number of cases would have been increased by $25 \% (2.91 \times 10^5).$ Thus the timing of implementing the interventions is crucial and analyzing the impact of time delay in the system is significant. Imposing the interventions sufficiently early, may have the potential to reduce the total number of infected cases. Figure 27 shows that if the time delay $(\tau_{2})$ had increased by 3/2 times $(\tau_{2} = 3)$ by its baseline value $(\tau_{2} = 2)$ then number of cumulative cases would have been increased by $4 \%$ and if the time delay $(\tau_{2})$ had decreased by 2 times $(\tau_{2} = 4)$ by its baseline value $(\tau_{2} = 2)$ then number of cumulative cases would have been decreased by $2 \%$.

We find the same qualitative result (Figures 28 and 29) when we parameterize the model using the data of COVID-19 in Spain. These results also demonstrate that model system can produce observed disease dynamics, but total cases are different when strength of interventions and time delays are altered. From Figure 28(a), we observe that the strength of intervention $(\alpha_{1})$ would have reduced the total number of cases at the early stage of epidemic in Spain. We note that Spain had recorded $2.39 \times 10^5$ total number of cases by May 31, 2020. We analyze that if the strength of interventions ($\alpha_1$) had increased, the total cases of COVID-19 would have decreased. Figure 28(a) shows that if the strength of interventions had increased by $100 \% (\alpha_{1} = 2\times 10^{-8})$ from its baseline value $(\alpha_{1} = 1\times 10^{-8})$, then the cumulative cases would have decreased by $25 \% (1.796 \times 10^5)$ and if the strength of interventions had increased by $200 \% (\alpha_{1} = 3\times 10^{-8})$ from its baseline value, then the cumulative cases would have decreased by $38 \% (1.48 \times 10^5)$. If the strength of interventions had increased by $300 \% (\alpha_{1} = 4\times 10^{-8})$ from its baseline value, then the cumulative cases would have decreased by $45 \% (1.30 \times 10^5).$ Figure 28(b) illustrates that how time delay in intervention affects the outcomes of early cumulative cases in Spain. We analyzed that if time delay $\tau_1$ increased by 1.28 times $(\tau_{1} = 4.5)$ by its baseline value $(\tau_{1} = 3.5)$, then cumulative number of cases would have been increased by $21 \% (2.91 \times 10^5)$ and if $\tau_1$ decreased by 1.4 times $(\tau_{1} = 2.5)$ by its baseline value, then cumulative number of cases would have been decreased by $7 \% (2.22 \times 10^5).$ Figure 29 shows that if the time delay $(\tau_{2})$ had increased by 3/2 times $(\tau_{2} = 3)$ by its baseline value $(\tau_{2} = 2)$ then number of cumulative cases would have been increased by $5 \%$ and if the time delay $(\tau_{2})$ had decreased by 2 times $(\tau_{2} = 1)$ by its baseline value $(\tau_{2} = 2)$ then number of cumulative cases would have been decreased by $2 \%$.

8.   Discussion

In the case of many infectious diseases, either we are unable to perform sufficient experiments, or a deeper understanding of disease dynamics is essentially required. Therefore, modeling of epidemics (mathematical and computational methods) is one of the important tools to understand the disease dynamics. Mathematical modeling along with precise computations gives new insights about different important issues of disease dynamics and their prevention, e.g., effects of different important parameters: delay, treatment, epidemic trends, intervention strategies, control, and spread of infections. Intervention strategies and treatment have great impression to control different contagious diseases. In the last few years, models accompanying the impact of intervention strategies on the spread of infectious diseases showed impressive popularity. But these models overlooked the impact of the availability of treatment resources on the dynamics of infectious disease. These models also do not involve different time delays, such as time delay in interventions and cure rate. However, in many biological systems, time delays are inevitably present, and several epidemic model systems have been explored incorporating the delays ^{[35,46,47,49]}. Here we put our efforts to understand the impact of intervention strategies and available resources for the treatment/control of infectious diseases, with special emphasis on the current pandemic COVID-19.

In this work, we have studied an SEIR model system (2.2) with infection force under intervention strategies and recovery rate under the low availability of resources for treatment. Furthermore, we included the two time delays in system (2.2), first-time delay $\tau_1$ represents the latent period of the intervention strategies, and second-time delay $\tau_2$ is due to the period that is used to cure the infected individuals. For systems (2.2) and (2.3), we obtained the various theoretical and epidemiological findings.

(ⅰ) We investigated that the system (2.2) has two equilibria, namely, disease-free equilibrium $(D_0)$ and endemic equilibrium $(D^*)$. The existence of equilibria depends on the basic reproduction number $(R_0).$ For the non-delayed model, stability analysis for both disease-free and endemic equilibriums are performed. In particular, it has been shown that for disease spread, the basic reproduction number $(R_0)$ acts as a threshold value. When $R_0 < 1,$ only disease-free equilibrium $(D_0)$ exists and locally asymptotically stable. When $R_0 > 1,$ there exists two equilibria, disease-free equilibrium $(D_0)$ and endemic equilibrium $(D^*)$.

(ⅱ) System (2.2) undergoes a transcritical bifurcation at $R_0 = 1,$ which means disease die out whenever $R_0 < 1$ and disease persists whenever $R_0 > 1.$ We have obtained that $D_0$ is globally stable in its feasible region $R_0 < 1$ and $D^*$ is globally stable if $R_0 > 1.$ In numerical analysis, by choosing a particular combination of $f(I)$ and $g(I)$, we have performed the impact of various parameters in the disease dynamics (Figure 7, 8). We have shown the relation of the different parameters with $R_0$ (Figure 9). We have also illustrated the combined effect of various parameters on the endemic level (Figure 10, 11, 12). Our observations revealed that increasing the interventions, disease level settle down at low level, but limitation on treatment sets the infected population at higher level. We additionally observed that $R_0$ can be reduced by increasing interventions and decreasing the limitation on resources availability.

(ⅲ) For delayed system (2.3), we have obtained that when $R_0 < 1,$ $D_0$ is locally asymptotically stable with certain conditions (Theorem 4.1) when $\tau_1 > 0$ and $\tau_2 > 0.$ Further, when $R_0 > 1$, $D_0$ remain unstable when $\tau_1 > 0$ and $\tau_2 > 0$. The system (2.3) seems sensitive to delay lengths. It is observed by both theoretical results and numerical simulations. When $R_0 > 1$, the effect of two-time delay has been investigated by taking different cases and the threshold values of both the delays have been obtained in section 4.3. Analytically and numerically, we have shown that system $(2.3)$ is stable when delay parameters are less than their threshold values (Figure 13, 14, 15, 16), which means the disease could be controlled easily. Whenever these parameters cross these threshold values, periodic solutions will occur through Hopf bifurcation. It can also be understood by the bifurcation diagram (Figure 17). In additional, the direction and stability of Hopf bifurcation are studied by normal form theory and center manifold theorem. We have also studied the global continuation of local Hopf bifurcation via both analytically and numerically (section 4.5 and Figure 18). We concluded that the time delays in the system (2.3) affect the structure of the solutions substantially. The time delays play the roles to destabilize the system. Indeed, it is evidenced by our theoretical results. The length of delay to preserve stability has also been estimated by employing the Nyquist criterion. More importantly, the stability crossing curves along with crossing directions ^[47,48,64] have been discussed in two-delay parameter space (refer to the stability switching curve plots in Figure 22). Our system (2.3) demonstrates rich dynamics, such as Hopf bifurcation, and is suitable for large population densities.

(ⅳ) The dynamics of system (2.3) is also affected by the parameters $\beta, \alpha_1, \gamma, \alpha_2$. Numerical simulations indicate that the critical value of $\tau_{1}$ decreases with an increase in $\beta$ (cf. Figure 19 and Table 1). Therefore, we conclude that the impact of an increase in transmission rate is to early the occurrence of Hopf bifurcation. Similarly, the critical value of $\tau_{1}$ increases with an increase in $\alpha_1$, which means the impact of increasing the value of the parameter $\alpha_1$ (imposing more interventions) is to postpone the occurrence of Hopf bifurcation. It can also be concluded that the epidemic cycle is more difficult to occur as the intervention strategies increases. On the other hand, the critical value of $\tau_{2}$ increases with an increase in $\gamma$ (cf. Figure 20 and Table 2). Thus, the impact of an increase in recovery rate is to postpone the occurrence of Hopf bifurcation while the impact of the increasing the value of the parameter $\alpha_2$ (increasing limitations on treatment) is to early the occurrence of Hopf bifurcation because the critical value of $\tau_{2}$ decreases with increasing $\alpha_2$.

(ⅴ) Furthermore, we have applied our model to a case study based on the data of COVID-19 in Spain and Italy. We estimated the parameters of system (2.3) by using the least-square method. It has also been depicted that if we could be able to increase the strength of interventions sufficiently, total number of cases in Spain and Italy would have reduced in early phase of COVID-19 (refer to Subsection 7.1). In both countries (Spain and Italy), the similar effects could also be achieved by reducing time delays in intervention and recovery. Further, it is expected that the similar scenarios could also be observed in many other European countries such as UK, Germany, and France. As such, the outcomes related to different significant values may differ from actual ones. Here we should note that many attributes of the system also depend on the published/collected data. The data have been collected for a short time span and do not have any future predictions. The quantitative results have been illustrated for the early phase of the pandemic COVID-19 in only two countries, Spain and Italy. Our model could also be applicable for other countries and also for long time periods.

8.1.   Epidemiological implications and future scope

Our proposed model systems generalize the study ^[10,11] and become more realistic via investigation of recovery rate and delay parameters on the system dynamics. Treatment resources play a significant role in protecting the individuals against a particular disease and to control a contagious disease, for example, COVID-19. Thus, our results are the complement to the study by Zhou et al. ^[11]. The results found here could be favorable for measuring the effect of the intervention strategies and treatment in the control of infectious diseases. The existence of periodic solutions in epidemic models is significant and interesting both in mathematics and applications since it provides a satisfactory explanation for the breakout and recurrence of a disease, which may have profound implications for the prediction, prevention, and control of disease transmission. If the delay due to the latent period of the intervention strategies is large enough, it will lead to repetitive episodes of endemic and in this situation, the disease would be out of control. Our study suggests that imposing intervention strategies timely may be helpful in controlling the outbreak of disease. Further, our study also suggests that to control and predict the disease spread, we should trim the time delays. When the delay is not too large, ultimately, the epidemic disease will become endemic. Similarly, the time delay due to the period used to cure the infected individuals also leads to periodic solutions. Appropriately, we should truncate the delays in the model as much as possible so that we could control and predict the disease spread. The outcomes of our results in Section 7 emphasized the importance of the intervention strategies and time delays for the case of COVID-19. Policy makers should take care of reducing the time delays in such biological processes. These policies can be used by policymakers and other professionals to make decisions on how certain interventions such as lockdown, quarantine, isolation, and social distancing can effectively be implemented.

By considering the incidence rate under intervention strategies and recovery rate under low resources availability treatment, the system (2.2) does not have complex dynamics. It will be intriguing to consider the incidence and recovery rate in which the infection and recovery function first increases and attain its maximum, and then eventually tends to a saturation level due to crowding effect and limitation of treatment, respectively. It will be insightful to compare the dynamics of these cases with our results as future work. In future studies, according to the actual situation of COVID-19 in Spain and Italy, the pre-stage exposed and post-stage exposed period ^[65] can be further considered in the model to analyze the impact of intervention and medical resources on COVID-19 prevention and control.

Acknowledgments

The research work of Sarita Bugalia is supported by the Council of Scientific & Industrial Research (CSIR), India [File No. 09/1131(0025)/2018-EMR-I]. The research work of Jai Prakash Tripathi is supported by Science and Engineering Research Board (SERB), India [File No. ECR/2017/002786] and UGC-BSR Research Start-Up-Grant, India [No. F.30-356/2017(BSR)]. The research work of Hao Wang is partially supported by The Natural Sciences and Engineering Research Council of Canada (NSERC) Individual Discovery Grant RGPIN-2020-03911 and Discovery Accelerator Supplement Award RGPAS-2020-00090. The authors are grateful to the handling editor and reviewers for their helpful comments and suggestions that have improved the quality of the manuscript.

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendices

Appendix A

Proof of Theorem 3.1: From the first equation of system (2.2), we have

and

Hence, if $S(0) > 0,$ then $S(t) > 0$. In the similar fashion, one can also check the positivity of $E(t), I(t),$ and $R(t)$.

Now, we shall prove the boundedness of all the solutions of system (2.2). For this, we define a function

Taking the time derivative of $N(t)$ and using (2.2), we obtain

which implies

and

This implies that all the solutions of system (2.2) are bounded. Thus we proved the Theorem 3.1.

Proof of Theorem 3.2: 1. The Jacobian matrix of right hand function in system (2.2) at $D_0$ is given by

which has two eigenvalues $-\mu, -\mu$ (which are negative) and the other eigenvalues are the eigenvalues of the following matrix

The trace and determinant of matrix $B$ are given by $\mathrm{tr}(B) = -(\mu + \delta + \frac{\gamma}{g(0)}+\mu +q) < 0$ and $\det (B) = (\mu + \delta + \frac{\gamma}{g(0)})(\mu +q)-q\frac{\beta A}{\mu f(0)}.$ We observe that $\det (B) > 0$ if $R_0 < 1$ and $\det (B) < 0$ if $R_0 > 1.$ Therefore, by Routh-Hurwitz criterion ^[50], it can be concluded that $D_0$ is locally asymptotically stable when $R_0 < 1$ and unstable when $R_0 > 1$.

2. When $R_0 = 1,$ i.e., $\beta = \beta^* = \frac{\mu f(0) (\mu +q) \left(\delta +\mu +\frac{\gamma }{g(0)} \right)}{A q}$, then

One can easily observe that $J_{(D_0, \beta^*)}$ has a simple zero eigenvalue and all other three eigenvalues are $-\mu, -\mu, -(2\mu +q +\delta +\frac{\gamma}{g(0)}),$ which are negative. Computing the right eigenvector $w$ and left eigenvector $v$ of $J_{D_0}$ with respect to zero eigenvalue, we have

Now from Theorem 4.1 of ^[51], we need to calculate the bifurcation constants $a$ and $b$. Denote $S = x_1, \, E = x_2, \, I = x_3,$ and $R = x_4.$ After some algebraic calculation, we obtain

Therefore, from Theorem 4.1 of ^[51], the disease-free equilibrium $(D_0)$ changes its stability from stable to unstable at $R_0 = 1.$ Moreover, there exists a positive equilibrium as $R_0$ crosses value 1, which is locally asymptotically stable. Therefore, system (2.2) undergoes a transcritical bifurcation at $R_0 = 1$.

Proof of Theorem 3.3: To prove the global stability of $D_0$, we follow the approach given by Castillo-Chavez et al. ^[52]. We rewrite the system (2.2) as follows

where $X = (S, R)\in \mathbb{R}^2$ signifies the number of uninfected individuals and $V = (E, I) \in \mathbb{R}^2$ signifies the number of infected individuals. Disease-free equilibrium $(D_0)$ is globally stable if the following two conditions are fulfilled:

(H1) For $\frac{dX}{dt} = F(X, V),$ $X^*$ is globally asymptotically stable,

(H2) $G(X, V) = MV - \hat{G}(X, V), \; \hat{G}(X, V) > 0$ for $(X, V)\in \mho$,

where $M = D_V G(X^*, 0)$ is an $M$-matrix. For the system (2.2), we have

It is obvious that the equilibrium $X^* = \left(\frac{A}{\mu}, 0 \right)$ is globally asymptotically stable of system (.2). Further, for system (2.2), we obtain

It is clear that $\hat{G}(X, V)\geq 0$. Hence, $D_0$ is globally stable.

Proof of Theorem 3.4: Define the following Lyapunov function:

Clearly, $L(D^*) = 0$ and $L(D)\neq 0$ when $D \neq D^*$. The upper right derivative of $L(S, E, I, R)$ can be calculated as:

In the above equation, there are 8 different types of situations depending on the size of $S$ and $S^*$, $E$ and $E^*$, $I$ and $I^*$, $R$ and $R^*$. It would be sufficient to show for $S > S^*, E > E^*, I > I^*, R > R^*$, similarly, one can also do for the other situations. Thus we obtain

Integrating the above inequality from $t_0$ to $t$ both sides, we obtain

Since $S, E, I, R$ are bounded and their derivatives are also bounded, therefore, $L$ is uniformly continuous on $[t_0, +\infty)$. From Barbalat's Lemma ^[53], $\lim\limits_{t \rightarrow +\infty}L(t) = 0$. Hence $D^{+}L < -\mu L < 0,$ which implies $D^*$ is globally stable.

Proof of Theorem 4.1: The characteristic equation of the Jacobian of system (4.1) at $D_0$ is given by

where

and $I$ is a $3 \times 3$ identity matrix. So the characteristic equation (.3) becomes in the following form

One eigenvalue is $-\mu$, which is negative and other eigenvalues are solution of the following equation

or

where

Let $\lambda = i \omega, (\omega > 0)$ be a root of Eq. (.5). Putting $\lambda = i \omega$ in Eq. (.5) and separating real and imaginary parts, we have

Squaring and adding above equations, we obtain

Let $\omega^2 = v$, then we have

By Descarte's rule of sign, if $V_{11}^{2}-V_{13}^{2}-2V_{12} > 0$ and $V_{12}^{2}-V_{14}^{2} > 0$, then Eq. (.6) has no positive roots. Hence all roots of Eq. (.6) have negative real parts.

If $\tau_{2} = 0$, then Eq. (.5) becomes in the following form

If $R_0 < 1$ then by Descarte's rule of sign, above equation has all its roots with negative real parts. Hence $D_0$ is locally asymptotically stable.

From $p(\lambda) = (\lambda + \mu +q) (\lambda +(\mu +\delta)+\frac{\gamma}{g(0)}e^{-\lambda \tau_{2}}) -\frac{\beta A q}{\mu f(0)}$, since $\lim\limits_{\lambda \rightarrow +\infty}p(\lambda) = +\infty$ and

Therefore, $p(\lambda) = 0$ has a positive real root. Hence, $D_0$ is unstable when $R_0 > 1$.

Proof of Lemma 4.3 and Theorem 4.4: If $\tau_{1} > 0, \tau_{2} = 0$, then the characteristic Eq. (4.3) becomes

Let $\lambda = i \omega_{1}, (\omega_{1} > 0)$ be a root of Eq. (.8). Putting $\lambda = i \omega_{1}$ in Eq. (.8) and separating imaginary and real parts, we have

Squaring and adding above equations, we obtain

where

Let $\omega_{1}^{2} = v_{1}$, then Eq. (.10) becomes

For the sign of the roots of Eq. (.10), by Descartes rule of sign, we have the following cases:

1. If $c_{11} > 0, c_{12} > 0, c_{13} > 0$, then Eq. (.10) has no positive roots. Hence all the roots of Eq. (.10) have negative real parts for $\tau_{1} \in (0, +\infty)$.

2. If $c_{13} < 0, c_{11} > 0, c_{12} > 0$ or $c_{13} < 0, c_{11} > 0, c_{12} < 0$ or $c_{13} < 0, c_{11} < 0, c_{12} < 0$, then Eq. (.10) has a unique positive root $\omega_{1}^{'} = \sqrt{v_{1}^{'}}$. Eliminating $\sin (\omega_{1} \tau_{1})$ from Eq. (.9), we obtain

3. If $c_{11} > 0$, $c_{12} < 0, c_{13} > 0$ or $c_{11} < 0$, $c_{12} > 0, c_{13} > 0$ or $c_{11} < 0$, $c_{12} < 0, c_{13} > 0$, then Eq. (.10) has two positive roots $\omega_{1}^{\pm} = \sqrt{v_{1}^{\pm}}$. Eliminating $\sin (\omega_{1} \tau_{1})$ from Eq. (.9), we obtain

4. If $c_{11} < 0$, $c_{12} > 0, c_{13} < 0$, then Eq. (.10) has three positive roots $\omega_{1}^{1, 2, 3} = \sqrt{v_{1}^{1, 2, 3}}$. Eliminating $\sin (\omega_{1} \tau_{1})$ from Eq. (.9), we obtain

Hence, Eq. (.8) has a pair of purely imaginary roots $\pm i\omega_{1}'$ with $\tau_{1} = \tau_{1_j}', j = 0, 1, 2, \dots$, a pair of purely imaginary roots $\pm i\omega_{1}^{\pm}$ with $\tau_{1} = \tau_{1_k}^{\pm}, k = 0, 1, 2, \dots$, and a pair of purely imaginary roots $\pm i\omega_{1}^{1, 2, 3}$ with $\tau_{1} = \tau_{1_l}^{1, 2, 3}, l = 0, 1, 2, \dots$. Let $\left\lbrace \tau_{1_j}\right\rbrace_{j = 0}^{+\infty} = \left\lbrace \tau_{1_j}'\right\rbrace_{j = 0}^{+\infty}$ such that $\tau_{1_0} < \tau_{1_1} < \tau_{1_2} < \dots < \tau_{1_j} < \dots$, $\left\lbrace \tau_{1_k}\right\rbrace_{k = 0}^{+\infty} = \left\lbrace \tau_{1_k}^{\pm}\right\rbrace_{k = 0}^{+\infty}$ such that $\tau_{1_0} < \tau_{1_1} < \tau_{1_2} < \dots < \tau_{1_k} < \dots$, where $\tau_{1_0} = \min\left\lbrace \tau_{1_0}^{+}, \tau_{1_0}^{-} \right\rbrace$, and $\left\lbrace \tau_{1_l}\right\rbrace_{l = 0}^{+\infty} = \left\lbrace \tau_{1_l}^{1, 2, 3}\right\rbrace_{l = 0}^{+\infty}$ such that $\tau_{1_0} < \tau_{1_1} < \tau_{1_2} < \dots < \tau_{1_l} < \dots$, where $\tau_{1_0} = \min\left\lbrace \tau_{1_0}^{1}, \tau_{1_0}^{2}, \tau_{1_0}^{3} \right\rbrace$.

Let $\lambda(\tau_{1}) = \alpha(\tau_{1})+i \omega_{1}(\tau_{1})$ be a root of Eq. (.8) near $\tau_{1} = \tau_{1_j}$ and $\alpha(\tau_{1_j}) = 0$, $\omega_{1}(\tau_{1_j}) = \omega_{1}'$, $j = 0, 1, 2, \dots$. Further, differentiating Eq. (.8) with respect to $\tau_{1}$ and substituting $\lambda = i\omega_{1}'$, $\tau_{1} = \tau_{1_j}$, $j = 0, 1, 2, \dots$, we obtain

where

Proof of Lemma 4.5 and Theorem 4.6: If $\tau_{1} = 0, \tau_{2} > 0$, then the characteristic Eq. (4.3) becomes

Let $\lambda = i \omega_{2}, (\omega_{2} > 0)$ be a root of Eq. (.16). Putting $\lambda = i \omega_{2}$ in Eq. (.16) and separating imaginary and real parts, we have

Squaring and adding above equations, we obtain

where

Let $\omega_{2}^{2} = v_{2},$ then Eq. (.18) becomes

1. If $d_{11} > 0, d_{12} > 0, d_{13} > 0$, then Eq. (.18) has no positive roots. Hence all roots of Eq. (.18) have negative real parts when $\tau_{2} \in (0, +\infty)$.

2. If $d_{13} < 0, d_{11} > 0, d_{12} > 0$ or $d_{13} < 0, d_{11} > 0, d_{12} < 0$ or $d_{13} < 0, d_{11} < 0, d_{12} < 0$, then Eq. (.18) has a unique positive root $\omega_{2}^{'} = \sqrt{v_{2}^{'}}$. Eliminating $\sin (\omega_{2} \tau_{2})$ from Eq. (.17), we obtain

3. If $d_{11} > 0$, $d_{12} < 0, d_{13} > 0$ or $d_{11} < 0$, $d_{12} > 0, d_{13} > 0$ or $d_{11} < 0$, $d_{12} < 0, d_{13} > 0$, then Eq. (.18) has two positive roots $\omega_{2}^{\pm} = \sqrt{v_{2}^{\pm}}$. Eliminating $\sin (\omega_{2} \tau_{2})$ from Eq. (.17), we obtain

4. If $d_{11} < 0$, $d_{12} > 0, d_{13} < 0$, then Eq. (.18) has two positive roots $\omega_{2}^{1, 2, 3} = \sqrt{v_{2}^{1, 2, 3}}$. Eliminating $\sin (\omega_{2} \tau_{2})$ from Eq. (.17), we obtain

Hence, Eq. (.18) has a pair of purely imaginary roots $\pm i\omega_{2}'$ with $\tau_{2} = \tau_{2_j}', j = 0, 1, 2, \dots$, a pair of purely imaginary roots $\pm i\omega_{2}^{\pm}$ with $\tau_{2} = \tau_{2_k}^{\pm}, k = 0, 1, 2, \dots$, and a pair of purely imaginary roots $\pm i\omega_{2}^{1, 2, 3}$ with $\tau_{2} = \tau_{2_l}^{1, 2, 3}, l = 0, 1, 2, \dots$. Let $\left\lbrace \tau_{2_j}\right\rbrace_{j = 0}^{+\infty} = \left\lbrace \tau_{2_j}'\right\rbrace_{j = 0}^{+\infty}$ such that $\tau_{2_0} < \tau_{2_1} < \tau_{2_2} < \dots < \tau_{2_j} < \dots$; $\left\lbrace \tau_{2_k}\right\rbrace_{k = 0}^{+\infty} = \left\lbrace \tau_{2_k}^{\pm}\right\rbrace_{k = 0}^{+\infty}$ such that $\tau_{2_0} < \tau_{2_1} < \tau_{2_2} < \dots < \tau_{2_k} < \dots$, where $\tau_{2_0} = \min\left\lbrace \tau_{2_0}^{+}, \tau_{2_0}^{-} \right\rbrace$; and $\left\lbrace \tau_{2_l}\right\rbrace_{l = 0}^{+\infty} = \left\lbrace \tau_{2_l}^{1, 2, 3}\right\rbrace_{l = 0}^{+\infty}$ such that $\tau_{2_0} < \tau_{2_1} < \tau_{2_2} < \dots < \tau_{2_l} < \dots$, where $\tau_{2_0} = \min\left\lbrace \tau_{2_0}^{1}, \tau_{2_0}^{2}, \tau_{2_0}^{3} \right\rbrace$.

Let $\lambda(\tau_{2}) = \alpha(\tau_{2})+i \omega_{2}(\tau_{2})$ be a root of Eq. (.18) near $\tau_{2} = \tau_{2_j}$ and $\alpha(\tau_{2_j}) = 0$, $\omega_{2}(\tau_{2_j}) = \omega_{2}'$, $j = 0, 1, 2, \dots$. Further, differentiating Eq. (.18) with respect to $\tau_{2}$ and substituting $\lambda = i\omega_{2}'$, $\tau_{2} = \tau_{2_j}$, $j = 0, 1, 2, \dots$, we obtain

where

Proof of Lemma 4.7 and Theorem 4.8: If $\tau_{1} = \tau_{2} = \tau$, then the characteristic Eq. (4.3) becomes

Let $\lambda = i \omega_{3}, (\omega_{3} > 0)$ be a root of Eq. (.24). Putting $\lambda = i \omega_{3}$ in Eq. (.24) and separating imaginary and real parts, we have

Squaring and adding above equations, we obtain

where

Let $\omega_{3}^{2} = v_{3},$ then Eq. (.26) becomes

1. If $e_{11} > 0, e_{12} > 0, e_{13} > 0$, then Eq. (.26) has no positive roots. Hence all roots of Eq. (.26) have negative real parts when $\tau \in (0, +\infty)$.

2. If $e_{13} < 0, e_{11} > 0, e_{12} > 0$ or $e_{13} < 0, e_{11} > 0, e_{12} < 0$ or $e_{13} < 0, e_{11} < 0, e_{12} < 0$, then Eq. (.26) has a unique positive root $\omega_{3}^{'} = \sqrt{v_{3}^{'}}$. Eliminating $\sin (\omega_{3} \tau)$ from Eq. (.25), we obtain

3. If $e_{11} > 0$, $e_{12} < 0, e_{13} > 0$ or $e_{11} < 0$, $e_{12} > 0, e_{13} > 0$ or $e_{11} < 0$, $e_{12} < 0, e_{13} > 0$, then Eq. (.26) has two positive roots $\omega_{3}^{\pm} = \sqrt{v_{3}^{\pm}}$. Eliminating $\sin (\omega_{3} \tau)$ from Eq. (.25), we obtain

4. f $e_{11} < 0$, $e_{12} > 0, e_{13} < 0$, then Eq. (.26) has two positive roots $\omega_{3}^{1, 2, 3} = \sqrt{v_{3}^{1, 2, 3}}$. Eliminating $\sin (\omega_{3} \tau)$ from Eq. (.25), we obtain

Hence, Eq. (.24) has a pair of purely imaginary roots $\pm i\omega_{3}'$ with $\tau = \tau_{j}', j = 0, 1, 2, \dots$, a pair of purely imaginary roots $\pm i\omega_{3}^{\pm}$ with $\tau = \tau_{k}^{\pm}, k = 0, 1, 2, \dots$, and a pair of purely imaginary roots $\pm i\omega_{3}^{1, 2, 3}$ with $\tau = \tau_{l}^{1, 2, 3}, l = 0, 1, 2, \dots$. Let $\left\lbrace \tau_{j}\right\rbrace_{j = 0}^{+\infty} = \left\lbrace \tau_{j}'\right\rbrace_{j = 0}^{+\infty}$ such that $\tau_{0} < \tau_{1} < \tau_{2} < \dots < \tau_{j} < \dots$, $\left\lbrace \tau_{k}\right\rbrace_{k = 0}^{+\infty} = \left\lbrace \tau_{k}^{\pm}\right\rbrace_{k = 0}^{+\infty}$ such that $\tau_{0} < \tau_{1} < \tau_{2} < \dots < \tau_{k} < \dots$, where $\tau_{0} = \min\left\lbrace \tau_{0}^{+}, \tau_{0}^{-} \right\rbrace$, and $\left\lbrace \tau_{l}\right\rbrace_{l = 0}^{+\infty} = \left\lbrace \tau_{l}^{1, 2, 3}\right\rbrace_{l = 0}^{+\infty}$ such that $\tau_{0} < \tau_{1} < \tau_{2} < \dots < \tau_{l} < \dots$, where $\tau_{0} = \min\left\lbrace \tau_{0}^{1}, \tau_{0}^{2}, \tau_{0}^{3} \right\rbrace$.

Let $\lambda(\tau) = \alpha(\tau)+i \omega_{3}(\tau)$ be a root of Eq. (.24) near $\tau = \tau_{j}$ and $\alpha(\tau_{j}) = 0$, $\omega_{3}(\tau_{j}) = \omega_{3}'$, $j = 0, 1, 2, \dots$. Further, differentiating Eq. (.24) with respect to $\tau$ and substituting $\lambda = i\omega_{3}'$, $\tau = \tau_{j}$, $j = 0, 1, 2, \dots$, we obtain

Proof of Lemma 4.9 and Theorem 4.10: If $\tau_{1} > 0, \tau_{2}\in (0, \tau_{2_j})$. In this case, we take $\tau_{1}$ as a bifurcation parameter and $\tau_{2}$ is in stable interval. Let $\lambda = \ i \omega_{4}, (\omega_{4} > 0)$ be a root of Eq. (4.3). Putting $\lambda = i \omega_{4}$ in Eq. (4.3) and separating imaginary and real parts, we have

From above equation, we obtain the following equation with respect to $\omega_{4}$

with

Assume that H1: Eq. (.33) has a positive root $\omega_{4}^{*}$. Eliminating $\sin(\omega_{4}\tau_{1})$ from Eq. (.32), we obtain

$j = 0, 1, 2, \dots.$ Now differentiating Eq. (4.3) with respect to $\tau_{1}$ and substituting $\lambda = i{\omega_{4}^{*}}$, $\tau_{1} = \tau_{1_j}^{*}$, $j = 0, 1, 2, \dots$, we obtain

where

Appendix B

Let $p_{1} = S-S^*, p_{2} = E-E^*, p_{3} = I-I^*$ and $\tau_{1} = \zeta + \tau_{1}^{*}$, where $\zeta \in \mathbb{R}$. Normalize the delay $\tau_{1}$ by the time-scaling $t\rightarrow t/\tau_{1}$. Then, the reduced system (4.1) becomes in the following functional differential form:

where $L_{\zeta}:C \rightarrow \mathbb{R}^{3}$ and $F:\mathbb{R} \times C \rightarrow \mathbb{R}^{3}$ are defined as

and

with

and

By the Riesz representation theorem ^[58], there exists a function $\vartheta (\theta, \zeta)$, whose components are of bounded variation for $\theta \in [-1, 0]$ such that

We can choose

with $\delta (\theta) = \left\lbrace \begin{array}{cc} 0, & \theta \neq 0, \\ 1, & \theta = 1. \end{array}\right.$

For $\varphi \in C([-1, 0], \mathbb{R}^{3})$, define

Then, the system (.36) is equivalent to the following system

For $\kappa \in C'([0, 1], (\mathbb{R}^{3})^{*})$, the adjoint operator $A^*$ of $A(0)$ is defined by

For $\varphi \in C([-1, 0], \mathbb{R}^{3})$ and $\kappa \in C'([0, 1], (\mathbb{R}^{3})^{*})$, we define the following bilinear inner form:

where $\vartheta(\theta) = \vartheta(\theta, 0)$.

From subsection 4.3, it can be seen that $\pm i \tau_{1}^{*}\omega_{4}^{*}$ are the eigenvalues of $A(0)$, therefore, $\pm i \tau_{1}^{*}\omega_{4}^{*}$ are also the eigenvalues of $A^{*}(0)$. Let $p(\theta) = (1, p_2, p_3)^{T}e^{i \tau_{1}^{*}\omega_{4}^{*}\theta}$ and $p^*(s) = D(1, p_2^*, p_3^*)^{T}e^{i \tau_{1}^{*}\omega_{4}^{*}s}$ be the eigenvectors of $A(0)$ and $A^{*}(0)$ corresponding to the eigenvalues $+ i \tau_{1}^{*}\omega_{4}^{*}$ and $- i \tau_{1}^{*}\omega_{4}^{*}$, respectively. Hence, we obtain

In order to show that $\langle p^*, \bar{p} \rangle = 0$ and $\langle p^*, p \rangle = 1$, we need to compute the value of $\bar{D}$. In view of Eq. (.38), we may choose

Using the same notation as in Hassard et al. ^[57], we compute the coordinates to describe the center manifold $C_0$ at $\zeta = 0$. Let $p_t$ be the solution of Eq. (.36), when $\zeta = 0$. Define

On the center manifold $C_0$, we have

where

and $z$ and $\bar{z}$ are local coordinates for center manifold $C_0$ in the direction of $p^*$ and $\bar{p}^*$. Note that $W$ is also real if $p_t$ is real, we consider only real solutions. For solutions $p_t \in C_0$ of Eq. (.36), since $\zeta = 0$,

where

then it follows from Eq. (.36), that

Hence, we have

After some algebraic calculation and comparing the coefficients with (.40), we obtain

with

$E_1$ and $E_2$ can be determined by the following equations

where

Hence, we can determine the following values

The above expressions regulate the behavior of bifurcating periodic solution at the critical value $\tau_{1}^{*}$. Specifically, $\mu_2$ regulates the direction of Hopf bifurcation: if $\mu_2 > 0\; (\mu_2 < 0)$ then the Hopf bifurcation is supercritical (subcritical) and bifurcating periodic solution exists for $\tau_{1} > \tau_{1}^{*} \; (\tau_{1} < \tau_{1}^{*})$. $\beta_ 2$ regulates the stability of the bifurcating periodic solution: the bifurcating periodic solution is stable (unstable) if $\beta_2 < 0 \; (\beta_2 > 0)$. $T_2$ regulates the period of bifurcating periodic solution: if $T_2 > 0$, the period increases, otherwise decreases.

Appendix C

Proof of Theorem 5.1: To prove this theorem, we use Theorem 4.3 of ^[60]. We note that disease-free equilibrium $(D_0)$ is unstable whenever $R_0 > 1$ for both systems (2.2) and (2.3). To show that (2.2) and (2.3) satisfy all the conditions of Theorem 4.3, when $R_0 > 1$, we choose $X = \mathbb{R}^{4}$ and $E = \mho.$ The maximal invariant set $M$ on the boundary $\partial \mho$ is the singleton ${D_0}$ and is isolated. Thus the hypothesis of Theorem 4.3 of ^[60] holds for both systems (2.2) and (2.3). We obtained that our systems have a unique endemic equilibrium whenever $R_0 > 1.$ Hence the necessary and sufficient condition for uniform persistence in Theorem 4.3 of ^[60] is equivalent to $D_0$ being unstable.