Stability switches, periodic oscillations and global stability in an infectious disease model with multiple time delays

1.   Introduction

Mathematical models have been found a successful tool to understand the underlying dynamics of infectious diseases. This eventually helps in controlling epidemic threats and related consequences so that effective and applicable control measures can be developed ^[1,2,3]. A large number of mathematical models, after the fundamental work of Kermack and McKendrick ^[4], have been formulated and analyzed to comprehend the dynamics of disease transmission and progression ^{[1,2,3,5,6,7,8,9]}. It is important to mention here that in the last few years, various infectious diseases such as SARS, influenza, Ebola, Zika and most recent COVID-19 have posed serious challenges and socio-economic consequences worldwide. Therefore, researchers are attracted to study various control strategies so that limited resources can be best utilized. For that various control interventions: pharmaceutical and non-pharmaceutical and their suitable combination are an important area of study ^{[5,8,10,11,12]}. Education, awareness, information, isolation, social distancing etc. are considered as the non-pharmaceutical controls ^{[10,11,13,14,15,16]}.

Nowadays, during and after the outbreaks of diseases, it has been observed that individuals change their behaviour due to the effect of information generated via global connectivity of social media, educational campaigns etc. which eventually alters the progression of infectious diseases ^{[17,18,19,20,21,22]}. Therefore, researchers have intended to study the effect of behaviour influencing factors induced by information on the spread of infectious diseases. Researchers have studied the effect of information on model dynamics, using it either on the force of infection or by making a subclass of aware population ^{[14,20,21,23,24,25,26,27,28,29]}.

The time lags/delays are inherent in natural, biological as well as man made systems and they play a crucial role in the disease dynamics as well. Moreover, the effect of time delays, appeared in various interactions during the disease progression, has been well studied and explored in literature. Delay differential equation models have been used to explain different kinds of time lags in biological system for example delay in infection (disease transmission), maturation, waning immunity etc. ^{[30,31,32,33,34,35,36,37]}. These delay models show rich and complex dynamics leading to different kind of bifurcations, instability of equilibria, periodic solutions etc. ^{[33,35,36,38,39,40,41,42]}. However, study of multiple delay models is challenging as it increases the complexity of the system and analysis.

In 1999, an epidemic model which accounts for the effect of maturation delay in the growth of population was proposed and analyzed by Cooke et al. ^[30]. For $R_0 > 1$, authors found that the disease will remain endemic in the population for all the time, either at the equilibrium value or show oscillations around it. Further, in 2005, Greenhalgh et al. proposed an SIRS model with vaccination effect on the susceptible population and the effect of delay on waning the vaccine-induced immunity ^[31]. Analytically, the existence of Hopf bifurcation was established by the authors. In 2005, Kyrychko et al. proposed and analyzed a delay differential equation model in which they used a nonlinear force of infection along with delay effect in loss of vaccine immunity ^[34]. When the time lag in loss of vaccine immunity crosses a threshold quantity, the global stability of the infected equilibrium was established by authors. In 2008, a delayed SIRS model along with temporary immunity was formulated and explored by Wen et al. and the global stability of the infected equilibrium was shown by constructing a Lyapunov function ^[37]. Further, in 2010, various epidemiological models such as SIR, SIS, SEIR and SEI along with time delay effect and general force of infection term were proposed by Huang et al. and they studied their global stability properties under time delay effect ^[32].

There are some recent developments in this direction that are related to our study and some of them are discussed in the following ^{[43,44,45,46]}. In 2022, Mezouaghi and their coworkers proposed a delayed epidemic model in which they included a separate class of protected individuals ^[45]. They have established the global stability of disease free and endemic equilibrium points by constructing Lyapunov functions. Further, Yang and others formulated an epidemic model which quantifies a convex type of force of infection and a delay effect in recovery rate ^[46]. They proved the global exponential stability of endemic equilibrium point along with the existence of Hopf bifurcation and their stability and direction. Recently, Lv et al. developed an SVIR model which accounts for the effect of two delays related to COVID-19 booster vaccination and failure of antibody (waning immunity) ^[43]. The authors have explored the dynamic properties of the model consisting of local and Hopf bifurcation analysis. They estimated the model parameters from the real data set and numerically validated the obtained analytical results.

In 2017, Kumar et al. ^[11] have considered the growth of information as function of infective at instantaneous time. However, in practice, reporting the disease may not be possible instantaneously and hence, there is always a time lag in the report by media and health agencies after the outbreak. Thus, the dissemination/propagation of information about the disease prevalence may be affected by this time lag. So this delay in reporting these infected individuals may have an impact on generated information. In order to quantify this factor, Kumar et al. ^[47] have considered this delay effect in their model and observed important insights due to delay in reporting. As we are aware that vaccination and recovery provide a certain immunity to fight with infection but this immunity may not have long-term effects and hence, waning takes place. Therefore, the waning immunity is one of the crucial factors which leads to the high reoccurrence of infectious diseases such as measles, chickenpox etc. because recovered or immunized individuals will again join the susceptible class ^[48,49]. Moreover, it is very difficult to determine the duration of protection level through any kind of immunity either vaccine-induced or natural and it may vary from the case to case ^[50]. Therefore, time lag in waning the immunity becomes very important to quantify in the modeling process. In some of the studies, it is found that this kind of delay gives complex dynamics ^[48,49]. Keeping, these crucial factors in mind, in this study, we consider two discrete time delays in the model. First, we consider the delay effect in the growth of information but assume the growth function as a linear function of infective which is meaningful if we consider relatively low infection persistence. Second, we quantify another delay in the waning immunity of individuals in the recovered class related to loss of protective measures (such as vaccination, self-protection, responsive behaviour etc.) with the aim to explore its impact mathematically and biologically. Thus, the main objective of this study is to explore the impact of above mentioned time delays on the disease dynamics and we would like to mention here that to the best of the author's knowledge, this type of particular setup has not been explored much in the literature though delay models have been analyzed but not the impact of information and immunity loss together in delay context. Consideration of two delays in a model itself increases huge mathematical complexity which will be more challenging. One can refer to articles ^[51,52] for more details of mathematical difficulties that arise in two delay model systems.

The delay models concerning information have been studied in the literature in different modeling sets up. For example see references [19,21,42,53,54,55,56]. Misra et al. in 2011, proposed a model for the dynamics of infectious diseases in which susceptible individuals become aware due to awareness driven by media with delay effect ^[21]. Authors have shown the occurrence of oscillations via Hopf bifurcation due to the delay effect. Zhao et al. in 2014, considered an SIRS model which accounts for the media impact on the force of infection along with delay effect ^[42]. They found that the disease free equilibrium is locally stable and also observed that the periodic solutions bifurcate around the infected equilibrium via Hopf bifurcation at a critical value of the delay. Further global Hopf bifurcation is also established for the model system. Liu et al. in 2015, proposed and analyzed an SIS epidemic model which accounts for the effect of delayed behavioural response on the force of infection ^[35]. They found that different forms of the delayed behavioural response exhibit various kinds of characteristics such as oscillations via Hopf and epidemic bubble.

Greenhalgh et al. in 2015, formulated an epidemic model in which the effect of awareness programs on the dynamics of diseases has been incorporated ^[19]. Further, the authors have extended their model to the corresponding delay model and observed that the delay induces oscillations via Hopf bifurcation. Uniform persistence of the disease with delay effect has also been established by authors. In 2017, Zhao et al. formulated an SIR epidemic model by quantifying the delayed effect of awareness driven by media in force of infection term ^[56]. They performed the model analysis and found that at the critical value of the delay, an endemic equilibrium loses its stability and bifurcates into periodic oscillations via Hopf bifurcation. But these models have incorporated the effect of information differently. However, the results obtained in our study match with them in principle, i.e., we observe oscillations in the populations due to the presence of delays. Li and others, in 2020, proposed an epidemic model which accounts for the effect of media publicity and found that shortening the delay of media reports may gradually decrease the infection from the population ^[53]. Further, in 2021, an SIRS epidemic model along with delay effect in information-induced force of infection term was formulated by Yadav and Srivastava in which treatment is in saturation type ^[54]. They found that the model system exhibits the existence of oscillatory behaviour due to the delay effect whereas, the saturation in treatment leads to the existence of multiple endemic equilibrium points. Recently, in 2022, Zhang and others formulated a compartmental model in which the susceptible population was divided into three subclasses due to the effect of awareness along with delay in immunity loss ^[55]. Authors explored the local dynamical properties such as local stability and Hopf bifurcation.

The rest of the paper is presented in the following order: in the next section, a delay mathematical model is proposed which contains two discrete delays. Further, the model analysis is performed for different cases of time delays in consecutive sections followed by numerical validation. The local bifurcations analysis along with the global stability of the equilibrium point is established in each case.

2.   The proposed delay model

A mathematical model is formulated, in this section, which accounts for the effect of two different time delays: one in the growth of information and another in the waning of immunity. The model given by Kumar et al. ^[11] is modified to study the dynamics along with the modified rate equation of the information's dynamics. As in practice, there is always a time needed to lose the immunity related to protection. Therefore, in this model, the effect of delay is accounted for in the waning immunity of the recovered population. For this, a time delay $\tau_1 > 0$ is considered when recovered individuals will move to the susceptible class after losing the immunity. Also, there is always a significant time required for health agencies to measure the approximate cumulative density of the infective population and thus naturally a time delay comes in the dissemination of information generated by them. Therefore, during the outbreak, the effect of time lags or delays must be quantified in the growth factor of information to make the model dynamics more realistic. Therefore, a time delay $\tau_2 > 0$ is quantified in the growth of information as the information at time $t$ will be measured by the infective population present at time $t-\tau_2$. Keeping the above facts, the proposed delay differential equation model is given by

with the initial population size $S(\theta) = S_0\geq0, I(\theta) = I_0\geq0, R(\theta) = R_0\geq0$ and $Z(\theta) = Z_0\geq0$ where $\theta \in[- \tau, 0]$, where $\tau = max\{\tau_1, \tau_2\}$. Also, the initial condition lies in the Banach space of continuous functions $C([- \tau, 0], \mathbb{R}^4_+)$.

In this study, three sub-populations ($S(t)$-susceptible, $I(t)$-infective, $R(t)$-removed) of the entire population ($N(t)$) are considered to study at any given time $t$. Also, the density of the information is measured by the variable $Z(t)$ which depends on the cumulative count of the infective as well as other social media/activities ^[11,39]. Thus, the information variable $Z(t)$ dimensionally follows the number of infected individuals or prevalence of the disease ^[39]. Here, the growth of the susceptible population is given by parameter $\Lambda$. Natural mortality and disease related deaths are measured by the parameters $\mu$ and $\delta$ respectively whereas $\gamma$ represents the recovery of the infected population. A homogeneously mixed population and a mass action type interaction between susceptible and infective are assumed for the disease transmission with the rate $\beta$. The loss of immunity is represented by the parameter $\delta_0$ which consists of the loss of natural immunity as well as immunity of protective measures. The factor $u_1 dZ(t)S(t)$ is quantified as the information-induced behavioural response of susceptible individuals with response rate $u_1 d$. Where the parameter $u_1$ is the response intensity with $0\leq u_1\leq1$ and the $d$ denotes the information interaction rate that influences an individual's behaviour. The growth of information ($Z$) and their natural degradation are measured by the parameters $a$ and $a_0$ respectively.

3.   Qualitative analysis of the model

This section focuses to study the qualitative insight of the mathematical model by investigating the stability of equilibrium points, and the existence of bifurcations such as Hopf bifurcation and Hopf-Hopf bifurcation due to considered delays. For this purpose, first, the boundedness and positivity of the solutions are established.

3.1.   Positivity and boundedness

From the model system (2.1) for $t\in [0, \tau]$, the solution are given as

Clearly, from the second equation of (3.1), we have $I(t)\geq0$. It is also clear from the third equation, $R(t)$ may or may not be positive (refer to Claim 1 of ^[57]). Following similar argument from Lv et al. ^[43], during the outbreak, the count of $I(t)$ and $R(t)$ may be very large than the $R(t-\tau_1)$, thus it is reasonable to assume $\gamma I +u_1 dZ S-\delta_0 R(t-\tau_1)$ positive which leads the $R(t) > 0$ and hence, $R(t-\tau_1) > 0$. Using, the positivity of $R(t)$ and $I(t)$, we have $S(t) > 0$ and $Z(t) > 0$. Further, in order to establish the bound of the population, we consider from above that all the sub-populations will remain positive for all time when the initial populations are positive. Further, notice that the entire population $N(t)$, in the model system (2.1), follows the following differential equation

Therefore, we have $\limsup\limits_{t\rightarrow \infty} N(t)\leq\frac{\Lambda}{\mu}$. Clearly, all the sub-populations ($S(t), I(t)$ and $R(t)$) are bounded by $\frac{\Lambda}{\mu}$. Moreover, using the bound of $I$, the last equation of the model system (2.1) implies $\limsup\limits_{t\rightarrow \infty} Z(t)\leq \frac{a \Lambda}{a_0 \mu}$. Using the above facts, the following positive invariant set describes the biologically feasible region for the model system (2.1):

3.2.   Equilibrium points of the model

The basic reproduction number of the model system (2.1) (as computed in Kumar et al. ^[11]) is given below:

The model (2.1) have following equilibrium points:

(i) a disease free equilibrium $E_1 = \left(\frac{\Lambda}{\mu}, 0, 0, 0\right)$ which exists always and unconditionally, and

(ii) a unique endemic equilibrium $E_2 = (S_{*}, I_{*}, R_{*}, Z_{*})$, which exists if and only if $\mathcal{R}_0 > 1$. Here, $S_{*} = \frac{(\mu+\delta+\gamma)}{\beta}$, $R_* = \frac{I_*}{\mu+\delta_0}\left(\gamma +\frac{du_1 a(\mu+\delta+\gamma)}{a_0 \beta}\right)$, $Z_* = \frac{aI_*}{a_0}$ and $I_* = -\frac{C}{B}$, where $B = \frac{\mu(\mu+\delta+\gamma)+\delta_0(\mu+\delta)}{(\mu+\delta_0)}+\frac{\mu du_1 a(\mu+\delta+\gamma)}{\beta a_0(\mu+\delta_0)}$ and $C = \Lambda\left(\frac{1}{\mathcal{R}_0}-1\right)$.

The stability results obtained in Kumar et al. ^[11] are restated below for the case of no delay ($\tau_1 = \tau_2 = 0$).

Theorem 1. ^[11] For $\tau_1 = \tau_2 = 0$,

(i) the disease free equilibrium $E_1$ of the system (2.1) is locally asymptotically stable if $\mathcal{R}_0 < 1$ and is unstable if $\mathcal{R}_0 > 1$,

(ii) if $\mathcal{R}_0 > 1$, then the unique endemic equilibrium $E_2$ is locally asymptotically stable provided following conditions are satisfied:

Here, $P_1 = a_0+2 \mu+\delta_0+\beta I_*+du_1 Z_*,$ $P_2 = a_0(\mu+\delta_0)+(\mu+a_0)(\mu+\beta I_*+du_1 Z_*)+\delta_0(\mu+\beta I_*)+\beta^2 S_*I_*,$ $P_3 = \beta I_*((a_0+\mu)(\mu+\delta+\gamma)+ \delta_0 (\mu+\delta))+\mu a_0(\mu+\beta I_*+du_1 Z_*)+a_0\delta_0(\mu+\beta I_*)+a du_1 \beta S_* I_*$ and $P_4 = \beta a_0 I_*(\mu(\mu+\delta+\gamma)+\delta_0(\mu+\delta))+a\mu \beta du_1 S_* I_*$.

3.3.   Stability results of the delay model

Stability of the equilibrium points of the delay model (2.1) are established in the subsequent parts when delays are present.

3.3.1.   Stability properties of the disease free equilibrium $E_1$

Theorem 2. (i) The disease free equilibrium $E_1$ is locally asymptotically stable for any $\tau_1, \tau_2\geq 0$, if $\mathcal{R}_0 < 1$ and $\delta_0\leq\mu$.

(ii) The disease free equilibrium $E_1$ is locally asymptotically stable for $0\leq\tau_1 < \widehat{\tau_1}$ and $\tau_2\geq 0$, if $\mathcal{R}_0 < 1$ and $\delta_0 > \mu$. Here

(iii) If $\mathcal{R}_0 > 1$, the disease free equilibrium $E_1$ is unstable for any $\tau_1, \tau_2\geq 0$.

Proof. The characteristic equation at $E_1$ is given by,

where, $J_{E_1} = \left(\begin{array}{cccc} -\mu & -\beta \frac{\Lambda}{\mu} & 0 & -du_1\frac{\Lambda}{\mu} \\ 0 & \beta \frac{\Lambda}{\mu} -(\mu+\delta+\gamma) & 0 & 0 \\ 0 & \gamma & -\mu & du_1\frac{\Lambda}{\mu} \\ 0 & 0 & 0 & -a_0 \\ \end{array} \right)$, $J_2 = \left(\begin{array}{cccc} 0 & 0 & \delta_0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & -\delta_0 & 0\\ 0 & 0 & 0 & 0 \\ \end{array} \right)$, $J_3 = \left(\begin{array}{cccc} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0\\ 0 & a & 0 & 0 \\ \end{array} \right)$ and $I_4$ is the identity matrix of order four. Now

Hence, the eigenvalues are $\lambda = -a_0, -\mu, {(\mu+\delta+\gamma)(\mathcal{R}_0-1)}$ which are negative if $\mathcal{R}_0 < 1$. The remaining eigenvalues are solution to

When $\tau_1 = 0$, the solution is $\lambda = -\mu-\delta_0 < 0$. For $\tau_1 > 0$, let us suppose that $\lambda = i \omega (\omega > 0)$ is a pure imaginary root of the Equation 3.2. Further, the real and imaginary parts are separated as follows

which provides $\omega^2 = \delta_0^2-\mu^2 > 0$ if $\delta_0 > \mu$.

Corresponding to $\omega = \sqrt{\delta_0^2-\mu^2}$ for $\delta_0 > \mu$, the Eq (3.3) gives

When $\delta_0\leq\mu$, the Eq (3.3) has no positive solution. Hence, the theorem is proven.

Example 1. Here, we shall numerically support the aforesaid analytical results for the stability properties of the disease free equilibrium with the help of model parameters which are selected as: $\Lambda = 5, \beta = 0.003, \mu = 0.04,$ $d = 0.17, \delta = 0.5, \delta_0 = 0.02, a = 0.1, a_0 = 0.1, \gamma = 0.1, u_1 = 0.9$. For this set of model parameters, the delay model induces only disease free equilibrium $E_1 = (125, 0, 0, 0)$ as in this case $\mathcal{R}_0 = 0.585 < 1$. Clearly, $\mu > \delta_0$ which ensures the local stability of the disease free equilibrium $E_1$ for all time lags $\tau_1, \tau_2\geq 0$ (Theorem 2(i)) and the corresponding numerical result is given in Figure 1.

Further, we choose $\delta_0 = 0.4 \; (>\mu = 0.04)$ except with the same parameters as above. In this case, the model has the same disease free equilibrium ($E_1$) and $\mathcal{R}_0$ along with the threshold value of the time delay $\tau_1$ as $\widehat{\tau_1} = 4.198$. Clearly from Theorem 2(ii), when the delay $\tau_1 < \widehat{\tau_1}$, the disease free equilibrium $E_1$ is locally stable for all $\tau_2\geq0$ and the corresponding obtained results are shown in Figure 2. Whereas, when $\tau_1 > \widehat{\tau_1}$, the disease free equilibrium $E_1$ loses its stability and Figure 3 shows the corresponding result.

Remark 1. As we have assumed previously that the population is positive, so Figure 3(a) is not biologically feasible due to oscillations in the negative zone though it is mathematically sound to give the instability nature or property of the disease free equilibrium $E_1$.

3.3.2.   Stability properties of the endemic equilibrium $E_2$

Here, we shall establish the local stability of the endemic equilibrium $E_2$ when both the time lags are present in the system. For this, around endemic equilibrium $E_2$, we linearise the delay model (2.1) as.

Here, $J_4 = \left(\begin{array}{cccc} -(\mu+\beta I_*+du_1 Z_*) & -\beta S_* & 0 & -du_1 S_* \\ \beta I_* & 0 & 0 & 0 \\ du_1 Z_* & \gamma & -\mu & du_1 S_* \\ 0 & 0 & 0 & -a_0 \\ \end{array} \right)$,

$Y(t) = (S(t), I(t), R(t), Z(t))^T$. Here, $J_2$ and $J_3$ are defined as in the previous subsection.

For the linearized system (3.4), the characteristic equation is given by,

with four order identity matrix $I_4$. Further, it can be rewritten as

where

Notice that, infinitely many complex roots can be obtained from the characteristic Eq (3.5) due to its transcendental nature in $\lambda$. We consider a similar argument for the stability of $E_2$ as discussed in ^[58]. Thus, the sign of real parts of the roots of the characteristic Eq (3.5) will determine the stability nature of endemic equilibrium $E_2$. The negative real parts will lead to local stability whereas the instability will rise in the case when a root crosses the imaginary axis (existence of a purely imaginary root). Also, the transcendental nature of Eq (3.5) infers that it is very tedious to determine the sign of roots of the Eq (3.5). Thus, the sign will change when the roots of the Eq (3.5) cross the imaginary axis as per Rouche's Theorem and continuity in time delays ($\tau_1$ and $\tau_2$).

The subsequent sections will determine the stability properties of the endemic equilibrium point $E_2$ in different scenarios.

4.   Case-I: $\tau_1 > 0$ and $\tau_2=0$

In this section, the stability properties of $E_2$ under $\tau_1 > 0$ and $\tau_2 = 0$ will be examined and the corresponding characteristic Eq (3.5) is given by,

It is noted that when $\tau_1 = 0$, the Eq (4.1) possess roots with negative real parts (Theorem 1) under the parametric conditions $P_1 P_2 > P_3 \; \; \text{and} \; \; P_1(P_2P_3-P_1P_4) > P_3^2$ where $P_1 = A_1+B_1, P_2 = A_2+B_2, P_3 = A_3+B_3+C_1$ and $P_4 = A_4+B_4+C_2$. Further, in order to get the purely imaginary roots when $\tau_1 > 0$, we replace $\lambda = i\omega$ in the Eq (4.1) and the corresponding real and imaginary parts are listed below:

Further, we square and add both the Eqs (4.2) and (4.3), we have

Here, $A_{11} = A_1^2-2A_2-B_1^2, A_{12} = A_2^2+2(A_4+C_2)-2A_1(A_3+C_1)-B_2^2+2B_1B_3$, $A_{13} = (A_3+C_1)^2-2A_2(A_4+C_2)+2B_2B_4-B_3^2$ and $A_{14} = (A_4+C_2)^2-B_4^2$.

By substituting $m = \omega^2$ in Eq (4.4), we have

It is very clear that the Eq (4.5) will have all roots with negative real part if the Routh-Hurwitz criterion holds true for (4.5) which gives the following result.

Theorem 3. The unique endemic equilibrium $E_2$ of the delay system (2.1) will be locally asymptotically stable for all $\tau_1 > 0$ and $\tau_2 = 0$ provided following conditions hold

4.1.   Existence of Hopf bifurcation

In this part, the occurrence of Hopf bifurcation (periodic solutions) is established when an endemic equilibrium loses its stability. For this, the delay parameter $\tau_1$ is considered as a bifurcation parameter. There exists a Hopf bifurcation at a threshold value of the delay ${\tau_1}_0$ if

$(H_1)$ $\lambda_{1, 2}({\tau_1}_0) = \pm i{\omega_1}_0 ({\omega_1}_0 > 0)$ and all other eigenvalues are with negative real parts at $\tau = {\tau_1}_0$.

$(H_2)$ $\left.\left[Re\left(\frac{d\lambda_{1, 2}}{d\tau_1}\right)^{-1}\right]\right|_{\lambda = i{\omega_1}_0 }\neq0$.

For the $(H_1)$ condition, we require that there exists at least one positive root of the Eq (4.5). Descartes' rule of signs ^[59] will determine the conditions for at least one positive root of the Eq (4.5) in following result.

Lemma 1. The Eq (4.5) has

(i) at least one positive root (either one or three) if

(a) $A_{11} > 0, A_{12} < 0, A_{13} > 0, A_{14} < 0$,

(b) $A_{11} < 0, A_{12} < 0, A_{13} > 0, A_{14} < 0$,

(c) $A_{11} < 0, A_{12} > 0, A_{13} > 0, A_{14} < 0$,

(d) $A_{11} < 0, A_{12} > 0, A_{13} < 0, A_{14} < 0$.

(ii) exactly one positive root if

(a) $A_{11} < 0, A_{12} < 0, A_{13} < 0, A_{14} < 0$,

(b) $A_{11} > 0, A_{12} < 0, A_{13} < 0, A_{14} < 0$,

(c) $A_{11} > 0, A_{12} > 0, A_{13} < 0, A_{14} < 0$,

(d) $A_{11} > 0, A_{12} > 0, A_{13} > 0, A_{14} < 0$.

(iii) at most two positive roots if

(a) $A_{11} > 0, A_{12} < 0, A_{13} > 0, A_{14} > 0$,

(b) $A_{11} > 0, A_{12} > 0, A_{13} < 0, A_{14} > 0$,

(c) $A_{11} < 0, A_{12} > 0, A_{13} > 0, A_{14} > 0$.

Assume that the Eq (4.5) holds one of the conditions of Lemma 1, then, the Eq (4.5) has at least one positive root say ${m_1}_0 = {\omega_1}_0^2$. Then, for a threshold value of the delay $\tau_1$, there will be a pair of purely imaginary roots ($\pm i {\omega_1}_0$) for the Eq (4.1). We obtain the following threshold value for the delay $\tau_1$ with the help of Eqs (4.2) and (4.3) as:

where

Further, we differentiate the Eq (4.1) with respect to $\tau_1$ to get the transversality condition $(H_2)$ as:

Lemma 2. Let $i{\omega_1}_0$ be a purely imaginary root with $m_{10} = {\omega_1}_0 ^2$ such that $\psi({\omega_1}_0) = 0$ and $\psi^{'}({\omega_1}_0)\neq0$, then $\left.\left[Re\left(\frac{d\lambda}{d\tau_1}\right)^{-1}\right]\right|_{\lambda = i{\omega_1}_0 }\neq0$ and its sign is the same as $\psi^{'}({\omega_1}_0)\neq0$.

Above results ensure that the transversality condition holds true. The following result summarises the above discussion.

Theorem 4. The unique endemic equilibrium $E_2$ is locally asymptotically stable for $\tau_1 < {\tau_1}_0$ and is unstable for $\tau_1 > {\tau_1}_0$. At $\tau_1 = {\tau_1}_0$, a Hopf bifurcation occurs, i.e., a family of periodic solutions bifurcates from the endemic equilibrium $E_2$ as delay parameter $\tau_1$ crosses the threshold value ${\tau_1}_0$ ^[60,61].

Numerical validation

For numerical validation of the Hopf bifurcation result, a set of parameters are chosen as in Example 1 except $\beta = 0.0125, \mu = 0.02$ and $\delta_0 = 0.05$ so that the basic reproduction number is $\mathcal{R}_0 = 5.04 > 1$. In this case, the model has the unique endemic equilibrium $E_2 = (49.59, 1.47,162.04, 1.47)$ along with the unstable disease free equilibrium $E_1 = (250, 0, 0, 0)$. Further, we notice that the coefficients of the Eq (4.5) are $A_{11} > 0, A_{12} < 0, A_{13} > 0$ and $A_{14} > 0$, and hence, the condition $iii(a)$ of Lemma 1 holds. This gives that the Eqs (4.5) and (3.5) have a positive root ($0.04509$) and a pair of purely imaginary root ($\pm 0.212i$ with ${\omega_1}_0 = 0.212$) respectively. Further, by using Eq (4.6), the threshold value of the delay $\tau_1$ is calculated as ${\tau_1}_0 = 7.97$ along with transversality condition $\left.\left[Re\left(\frac{d\lambda}{d\tau_1}\right)^{-1}\right]\right|_{\lambda = i{\omega_1}_0 } = 2.2\times10^2 > 0$. Thus, Theorem 4 ensures that the delay system (2.1) will be stable for the delay range $\tau_1\in[0, {\tau_1}_0)$ and unstable for $\tau_1 > {\tau_1}_0$. Whereas periodic oscillations bifurcate at $\tau_1 = {\tau_1}_0 = 7.97$ around $E_2$ when $\tau_1$ crosses ${\tau_1}_0$.

Further, we numerically exhibit the stability and instability of the unique endemic equilibrium $E_2$ by solving the delay system (2.1) with the help of DDE23 in MATLAB. The delay model is solved for the delay $\tau_1 = 7 < {\tau_1}_0$ with initial population size $S(\theta) = 70, I(\theta) = 5, R(\theta) = 10$ and $Z(\theta) = 5$ for $\theta\in [-\tau_1, 0]$ and model parameters taken as above. Figure 4 describes the corresponding result which infers that the unique endemic equilibrium $E_2$ is asymptotically stable.

Furthermore, Theorem 4 infers that the unique endemic equilibrium $E_2$ loses its stability leading to the existence of a family of periodic solutions when delay $\tau_1$ crosses the threshold value ${\tau_1}_0 = 7.97$. To show this numerically, the delay system (2.1) is further solved for the delay parameter $\tau = 9 > {\tau_1}_0$. The corresponding periodic solutions are plotted in Figure 5 which mimics the obtained analytical result. Therefore, our study infers that the disease will survive in oscillatory nature due to the delay effect considered in the waning of the immunity. Moreover, if this delay takes more than eight days, then it will be very challenging to estimate the actual size of the epidemic due to oscillations which leads the difficulty in control implementation. In addition, in this case, disease elimination may be very critical and challenging due to oscillations and hence, the actual estimation of an epidemic cannot be measured. A similar kind of observation is found in the study of Barbarossa and others ^[48] though in this case, authors have taken waning immunity in a different context other than in our case. In addition, recently in 2022, Zhang and others also established the existence of such oscillations (Hopf bifurcation) when the delay in waning crosses a threshold value when the susceptible population is divided based on the constant effect of awareness. A similar kind of observation has also been encountered by Yadav and Srivastava when a constant effect of information with delay is quantified ^[54]. Whereas, in this proposed study, a dynamic effect of information is considered along with a delay in waning immunity. Thus, this proposed study gives more generalised biological and mathematical insights. Our obtained Hopf bifurcation result and study also have good agreement in results claimed by Lv et al. ^[43]. In 2023, Lv et al. investigated the delay effects related to the failure of antibody and COVID-19 booster vaccination dose, and they found the existence of periodic oscillations for a real data set.

In addition, the bifurcation diagrams are plotted to show the periodic orbits when the unique endemic equilibrium loses its stability at the threshold value of $\tau_1$. The delay parameter is varied in $ \tau_1\in^[7,11] $ to plot the diagrams and are shown in Figures 6–8. One can easily see from Figures 6 and 7 that when the delay parameter $\tau_1$ lies in the range $[7, 7.97)$, the unique endemic equilibrium $E_2$ is stable and as $\tau_1$ crosses the threshold value ${\tau_1}_0 = 7.97$ bifurcation takes place. In addition, we also plot the occurrence of periodic orbits as the delay parameter $\tau_1$ varies, and the corresponding bifurcation diagram is shown in Figure 8.

4.1.1.   Global stability of unique endemic equilibrium $E_2$

This particular part explores the global properties of endemic equilibrium $E_2$ by proposing a suitable Lyapunov function.

Theorem 5. If $\mathcal{R}_0 > 1$, $\tau_1 > 0$ and $\tau_2 = 0$, the unique endemic equilibrium $E_2$ of the model system (2.1) is globally asymptotically stable under following parametric conditions $a^2 < 2 a_0 (\mu+\delta),$ $\delta_0^2 < \frac{4\gamma \mu^2}{3(2\mu+\delta)},$ $(du_1)^2 < \frac{2a_0\gamma \mu^3}{3\Lambda^2 (2\mu+\delta)}$ and $((2\mu+\delta)du_1Z_* -2\mu\gamma)^2 < \frac{4}{3}\gamma \mu^2 (2\mu+\delta)$.

Proof. In order to establish the global stability, we construct the following positive definite function in $\sum$,

The positive constants $K_1$ and $K_2$ will be selected suitably later. Now, we shall differentiate the function $U_1$ with respect to $t$ along the solution trajectories of the model (2.1) and using parametric relationship of equilibrium points.

We have,

Further, selecting and replacing the values of positive constants as: $K_1 = \frac{2\mu+\delta}{\gamma}$ and $K_2 = \frac{2\mu+\delta}{\beta}$, we have

where $F_1((S-S_*)(R-R_*)) = \left[\mu(S-S_*)^2-\left(\frac{2\mu+\delta}{\gamma} du_1 Z_*-2\mu\right)(S-S_*)(R-R_*)+\frac{2\mu+\delta}{3\gamma} \mu (R-R_*)^2\right]$, $F_2((R-R_*)(R(t-\tau_1)-R_*)) = \left[\frac{2\mu+\delta}{3\gamma} \mu (R-R_*)^2+\frac{2\mu+\delta}{\gamma} \delta_0 (R-R_*)(R(t-\tau_1)-R_*)+\mu(R(t-\tau_1)-R_*)^2\right]$, $F_3((Z-Z_*)(R-R_*)) = \left[\frac{a_0}{2} (Z-Z_*)^2-\frac{2\mu+\delta}{\gamma}du_1S(Z-Z_*)(R-R_*)+\frac{2\mu+\delta}{3\gamma} \mu (R-R_*)^2\right]$ and $F_4((I-I_*)(Z-Z_*)) = \left[(\mu+\delta)(I-I_*)^2-a(I-I_*)(Z-Z_*)+\frac{a_0}{2} (Z-Z_*)^2\right]$. Note that within the feasible region $\sum$, $\dot{U_1}\leq0$ if following parametric conditions $a^2 < 2 a_0 (\mu+\delta),$ $\delta_0^2 < \frac{4\gamma \mu^2}{3(2\mu+\delta)},$ $(du_1)^2 < \frac{2a_0\gamma \mu^3}{3\Lambda^2 (2\mu+\delta)}$ and $((2\mu+\delta)du_1Z_* -2\mu\gamma)^2 < \frac{4}{3}\gamma \mu^2 (2\mu+\delta)$ hold true. It is observed that $\dot{U_1} = 0$ if and only if at endemic equilibrium $E_2$ and $\dot{U_1} < 0$ otherwise within the feasible region $\sum$, therefore the singleton set $\{E_2\}$ is the largest positively invariant set contained in $\{(S, I, R, Z)\in\sum:\dot{U_1} = 0\}$. Hence, in the interior of $\sum$, Lyapunov LaSalle's theorem ^[62] ensures that the endemic equilibrium $E_2$ is globally asymptotically stable.

Example 2. In this part, we shall numerically support the aforesaid global stability property of the unique endemic equilibrium with the help of model parameters which are selected as: $\Lambda = 5, \beta = 0.0125, \mu = 0.02, d = 0.0017, \delta = 0.5, \delta_0 = 0.005, a = 0.01, a_0 = 0.1, \gamma = 0.1, u_1 = 0.009$. For this set of model parameters, the delay model has the unique endemic equilibrium $E_2 = (49.6, 6.67, 26.73, 0.66)$ along with the unstable disease free equilibrium $E_1 = (250, 0, 0, 0)$ and in this case, $\mathcal{R}_0 = 5.0403 > 1$. Clearly, the parametric conditions for global stability $a^2-2 a_0 (\mu+\delta) = -0.1039 < 0,$ $\delta_0^2-\frac{4\gamma \mu^2}{3(2\mu+\delta)} = -7.3765\times 10^{-5} < 0,$ $(du_1)^2-\frac{2a_0\gamma \mu^3}{3\Lambda^2 (2\mu+\delta)} = -3.7165\times 10^{-9}$ and $((2\mu+\delta)du_1Z_* -2\mu\gamma)^2-\frac{4}{3}\gamma \mu^2 (2\mu+\delta) = -1.2844\times 10^{-5} < 0$ hold true. This confirms the global stability of the endemic equilibrium $E_2$ for all time lags $\tau_1, \tau_2 > 0$ (Theorem 5).

From biological point of view, our global stability result infers that irrespective of time lag in immunity loss (even if for larger delay), the disease will persist in stable sense and it can be eradicated by bringing the basic reproduction number below one. Thus, oscillations may appear locally and will disappear in larger context and this is one of the important biological as well as mathematical insight.

5.   Case-II: $\tau_1=0$ and $\tau_2 > 0$

The characteristic equation corresponding to this case is defined as follows:

Consider $P_1 = A_1+B_1, Q_1 = A_2+B_2, W_1 = A_3+B_3$ and $D_1 = A_4+B_4$. Using the similar argument as in Case-I, in the Eq (5.1), we further replace $\lambda = i\omega$ and the corresponding real and imaginary parts are given by.

Now, we square and add both the above Eqs (5.2) and (5.3), we have

Here, $B_{11} = P_1^2-2Q_1, B_{12} = Q_1^2+2D_1-2P_1W_1, B_{13} = W_1^2-2D_1Q_1-C_1^2$ and $B_{14} = D_1^2-C_2^2$. Further, we substitute $m = \omega^2$ in Eq (5.4), we get

Clearly, the Eq (5.5) will have all roots with negative real part if the Routh-Hurwitz criterion holds true for (5.5) which gives the following result.

Theorem 6. The unique endemic equilibrium $E_2$ of the delay system (2.1) will be locally asymptotically stable for all $\tau_2 > 0$ provided following conditions hold

5.1.   Existence of Hopf bifurcation

Following similarly as in Case-I, the delay parameter $\tau_2$ is taken as a bifurcation parameter and at the threshold value of delay $\tau = {\tau_2}_0$, the assumptions $(H_1)$ and $(H_2)$ must hold for Hopf bifurcation.

At least one positive root of the Eq (5.5) is needed to get assumption $(H_1)$ which can be obtained using Descartes' rule of signs and following the same result as discussed in Lemma 1 for $B_{ij}$ for $i, j = 1-4$ in place of $A_{ij}$ for $i, j = 1-4$. Assume that the Eq (5.5) satisfies one of the conditions given in Lemma 1 to have at least one positive root say ${m_2}_0 = {\omega_2}_0 ^2$. Hence, the Eq (5.1) will have a pair of purely imaginary root ($\pm i{\omega_2}_0$) at the threshold of the delay ${\tau_2}_0$. And the corresponding ${\tau_2}_0$ is calculated using Eqs (5.2) and (5.3) and is given by.

where,

Further, the Eq (5.1) is differentiated with respect to $\tau_2$ to get the transversality condition $(H_2)$ as.

Lemma 3. Let $i{\omega_2}_0$ be a purely imaginary root with $m_{20} = {\omega_2}_0 ^2$ such that $\psi({\omega_2}_0) = 0$ and $\Psi^{'}({\omega_2}_0)\neq0$, then $\left.\left[Re\left(\frac{d\lambda}{d\tau_2}\right)^{-1}\right]\right|_{\lambda = i{\omega_2}_0 }\neq0$ and its sign is the same as $\Psi^{'}({\omega_2}_0)\neq0$.

Above discussion ensures that the transversality condition holds true which leads the following result.

Theorem 7. The unique endemic equilibrium $E_2$ is locally asymptotically stable for $\tau_2 < {\tau_2}_0$ and is unstable for $\tau_2 > {\tau_2}_0$. At $\tau_2 = {\tau_2}_0$, a Hopf bifurcation occurs, i.e., a family of periodic solutions bifurcates from the endemic equilibrium $E_2$ as delay parameter $\tau_2$ crosses the threshold value ${\tau_2}_0$ ^[60,61].

Numerical validation

For the numerical experimentation of the Hopf bifurcation, a set of parameters as given in Example 1 are considered expect $\beta = 0.03$ so that the basic reproduction number becomes $\mathcal{R}_0 = 5.859 > 1$. For this set of parameters, the model has a unique endemic equilibrium $E_2 = (21.33, 1.49, 83.54, 1.49)$ along with the unstable disease free equilibrium $E_1 = (125, 0, 0, 0)$. For the numerical simulation, we consider the initial population size $S(\theta) = 30, I(\theta) = 5, R(\theta) = 85$ and $Z(\theta) = 1$ for $\theta \in [-\tau_2, 0]$. Notice that the coefficients $B_{11}, B_{12}$ of the Eq (5.5) are positive and $B_{13}, B_{14}$ are negative, and hence, the condition Lemma 1$ii$(c) holds true. Therefore, the Eqs (5.5) and (5.1) have exactly one positive root ($0.04005$) and a pair of purely imaginary root ($\pm 0.20012i$ with ${\omega_2}_0 = 0.20012$) respectively. Also, in this case, ${\tau_2}_0 = 1.537$ and $\left.\left[Re\left(\frac{d\lambda}{d\tau_2}\right)^{-1}\right]\right|_{\lambda = i{\omega_2}_0 } = 2.11 > 0$. Thus, Theorem 7 infers that the delay system (2.1) will remain stable for the delay range $\tau_2\in[0, {\tau_2}_0)$ and unstable for $\tau_2 > {\tau_2}_0$. Whereas, when $\tau_2$ crosses ${\tau_2}_0$, the $E_2$ loses stability and leads the occurrence of periodic oscillations at $\tau_2 = {\tau_2}_0 = 1.537$.

Further, the delay model is examined numerically for delay parameter $\tau_2 = 1 < {\tau_2}_0$ to check the stability of the unique endemic equilibrium $E_2$ and the results are plotted in Figure 10 which infers that $E_2$ is locally stable. In this case, the disease can be eradicated if the time lag in propagation of information ($\tau_2$) is less than couple of days after outbreak of the disease by bringing down the basic reproduction number $\mathcal{R}_0$ less than one (Theorem 2). From Theorem 7, as $\tau_2$ crosses ${\tau_2}_0 = 1.537$, the periodic solutions will appear around $E_2$. To show this, the delay model (2.1) is further solved for $\tau_2 = 3$ and the oscillatory trajectories are plotted in Figure 11. Clearly, due to the delay effect considered in information growth, the disease will survive oscillatory in behaviour within the population and hence the implementation of suitable control is complicated.

It is important to notice here that, if the propagation of information takes more than couple of days after the outbreaks to make aware the population, then the population will experience the oscillatory behaviour of the disease. Therefore, in this case, it will be very difficult to account the actual size of the epidemic and hence, the control will be very challenging due to high and low peaks of the disease. Therefore, in order to control the spread of the disease, health agencies and governments have to propagate the information very quickly after the disease outbreak.

Now to further get insight into the model numerically, we wish to plot the bifurcation diagram to depict the oscillations and instability around $E_2$. For this purpose, we vary $\tau_2\in[0.5, 3]$ and the corresponding bifurcation diagrams are shown in Figure 15. It is very clear from Figure 15 that the unique endemic equilibrium $E_2$ is stable for the range $\tau_2\in[0.5, 1.537)$ whereas it loses its stability at ${\tau_2}_0 = 1.537$ leads the occurrence of oscillations in the population for $\tau_2 > 1.537$. Some authors such as Misra et al. and Greenhalgh et al. ^[19,21] have investigated similar kind of results for Hopf bifurcation using different modeling techniques.

5.2.   Existence of Hopf-Hopf bifurcation

This particular section focuses to study a different kind of bifurcation named as a Hopf-Hopf bifurcation which appears at two different values of the delay. It can be easily seen that the Eq (5.5) may have more than one positive root depending on the signs of the coefficients $B_{11}, B_{12}, B_{13}$ and $B_{14}$. Let us consider that the Eqs (5.5) and (3.5) have two positive roots (${m_2}_1 = {\omega_2}_1^2$ and ${m_2}_2 = {\omega_2}_2^2$ and two pair of purely imaginary roots ($\pm i{\omega_2}_1$ and $\pm i{\omega_2}_2$) respectively. Also, the critical values of the delay $\tau_2$ are obtained as.

where $\Phi({\omega_2}_j) = \frac{{\omega_2}_j ^4(P_1C_1-C_2)+{\omega_2}_j ^2(C_2Q_1-C_1W_1)-C_2D_1} {C_2^2+C_1^2{\omega_2}_j ^2},$ $j = 1, 2.$ The corresponding transversality conditions are given by,

We state the following result using the results given in ^[33,63,64].

Lemma 4. Let $i{\omega_2}_j$ be the purely imaginary roots with ${m_2}_j = {\omega_2}_j^2$ such that $\psi({\omega_2}_j) = 0$ and $\Psi^{'}({\omega_2}_j)\neq0$, \; j = 1, 2, then $\left.\left[Re\left(\frac{d\lambda}{d\tau_2}\right)^{-1}\right]\right|_{\lambda = i{\omega_2}_1} > 0$ if $\Psi^{'}({\omega_2}_1) > 0$ and $\left.\left[Re\left(\frac{d\lambda}{d\tau_2}\right)^{-1}\right]\right|_{\lambda = i{\omega_2}_2} < 0$ if $\Psi^{'}({\omega_2}_2) < 0$.

Using the Lemma 4, the existence of Hopf-Hopf (double) bifurcation is ensured in the following.

Theorem 8. The delay model undergoes Hopf-Hopf (double) bifurcation at ${\tau_2}_1$ and ${\tau_2}_2$ respectively. The unique endemic equilibrium $E_2$ is locally asymptotically stable for $\tau_2 < {\tau_2}_1$. At $\tau_2 = {\tau_2}_1$ the system undergoes Hopf bifurcation and the endemic equilibrium $E_2$ loses its stability. Further as delay parameter $\tau_2$ increases, the endemic equilibrium $E_2$ regains its stability at $\tau_2 = {\tau_2}_2$ and again system undergoes Hopf bifurcation. The endemic equilibrium $E_2$ remains stable for $\tau_2 > {\tau_2}_2$.

Numerical validation

The Hopf-Hopf (double) bifurcation result is validated numerically in this part, for this, a set of parameters as given in Example 1 except $\beta = 0.03$ and $\delta_0 = 0.5$. In this case, we note $\mathcal{R}_0 = 5.859 > 1$ and the unique endemic equilibrium $E_2 = (21.33, 5.25, 32.73, 5.25)$ is obtained for the delay model along with unstable disease free equilibrium $E_1 = (125, 0, 0, 0)$. Notice that the Eq (5.5) consists two positive roots ${m_2}_1 = 0.0452, {m_2}_2 = 0.0097$ as $B_{11}$ and $B_{14}$ are positive, and $B_{12}$ and $B_{13}$ are negative. Hence, two pair of purely imaginary roots $\pm 0.212 i$ and $\pm 0.0985 i$ with ${\omega_2}_1 = 0.212, {\omega_2}_2 = 0.0985$ are found for the characteristic Eq (5.1). The corresponding critical values of the delay parameter $\tau_2$ are calculated by using (9.1) and given as: ${\tau_2}_1 = 6.18$ and ${\tau_2}_2 = 29.79$. Also The transversality conditions are determined as: $\left.\left[Re\left(\frac{d\lambda}{d\tau_2}\right)^{-1}\right]\right|_{\lambda = i{\omega_2}_1} = 1.75 > 0$ and $\left.\left[Re\left(\frac{d\lambda}{d\tau_2}\right)^{-1}\right]\right|_{\lambda = i{\omega_2}_2} = -0.792 < 0$. Hence Theorem 8 follows.

To show the stability and bifurcation, we simulate the delay system for the initial population size: $S(\theta) = 30, I(\theta) = 7, R(\theta) = 30$ and $Z(\theta) = 4$ for $\theta\in [-\tau_2, 0]$ along with above parameters. Figures 16 and 17 describe the corresponding Hopf-Hopf bifurcation result numerically. Notice from Figures 16 and 17 that the unique endemic equilibrium $E_2$ is locally stable for $\tau_2\in[0, {\tau_2}_1 = 6.18)$ and at $\tau_2 = {\tau_2}_1 = 6.18$ periodic oscillations appear through first Hopf bifurcation as $E_2$ loses its stability. Further, as the bifurcation parameter increases, at $\tau_2 = {\tau_2}_2 = 29.79$, periodic oscillations die out through second Hopf bifurcation and for $\tau_2 > {\tau_2}_2 = 29.79$, the endemic equilibrium $E_2$ regains its stability and remains stable. Thus, the delay in information propagation leads to the complex behaviour of the diseases.

Moreover, our finding suggests that Hopf-Hopf bifurcation gives a region for a delay that not only exhibits stability but also oscillations with different periods in population as shown in Figure 16. It is very clear from Figure 18 that $E_2$ is stable for $\tau_2 = 6 < {\tau_2}_1 = 6.18$ and $\tau_2 = 31 > {\tau_2}_2 = 29.79$ and also oscillations appear around $E_2$ for ${\tau_2}_1 < \tau_2 = 20 < {\tau_2}_2$. Therefore, multiple stability switches are observed due to the delay effect taken into account in the model system. We further infer that the system undergoes the existence of periodic orbits for $\tau_2\in[6.18, 29.79]$ as given in Figure 17 and hence infection will persist oscillatory within the population. In this case, it is clear from Figure 19 that the solution trajectories approach to periodic orbit for a fix $\tau_2 = 15$.

Thus, it is important to notice here that the control of disease spread will be very challenging when the delay in the dissemination of information lies $\tau_2\in[6.18, 29.79]$ due to oscillations (high and low peaks of an epidemic). Hence, if information propagates in the early phase of the epidemic, a reduction in disease transmission can be obtained. Therefore, biologically this Hopf-Hopf bifurcation result is very important which provides a range of specific delays where the infection will persist oscillatory and as discussed above this is a crucial range for delay parameter. In order to control the oscillations within the population, a reduction in delay value (propagation of information) is needed that is early dissemination of information is required. This kind of similar stability switches results were discussed in Barbarossa and others ^[48] but in a different context whereas a similar kind of observation is found in Kumar et al. ^[47]. In Liu et al. ^[35], authors obtained a closed loop shaped curve during Hopf-Hopf bifurcation which they named as 'endemic bubble'. In this section, we have also obtained a similar 'endemic bubble' by varying the bifurcation parameter. All the above studies include only a single delay effect in different contexts whereas, in this proposed model, we have considered two delays which gives the more generalised version of existing studies as two delays lead the high mathematical complexity.

5.2.1.   Global stability of unique endemic equilibrium $E_2$

Here, by formulating a suitable Lyapunov function, we shall explore the global properties of unique endemic equilibrium $E_2$.

Theorem 9. If $\mathcal{R}_0 > 1$, $\tau_1 = 0$ and $\tau_2 > 0$, the unique endemic equilibrium $E_2$ of the model system (2.1) is globally asymptotically stable under following parametric conditions $a < min \{\delta, 2a_0\},$ $(du_1)^2 < \frac{a_0\gamma \mu^2 (\mu+\delta)}{\delta \Lambda^2}$ and $(\delta du_1Z_* -2\mu\gamma)^2 < 2\gamma \mu \delta (\mu+\delta_0)$.

Proof. For this purpose, we consider the positive definite function in $\sum$ which is given by,

The positive constants $K_1$ and $K_2$ will be selected suitably later. Now, we shall differentiate the function $U_1$ with respect to $t$ along the solution trajectories of the model (2.1) and using the parametric relationship of equilibrium points. We have,

Now, choosing and replacing the values of positive constants as: $K_1 = \frac{\delta}{\gamma}$ and $K_2 = \frac{2\mu+\delta}{\beta}$, we get

where $G_1((S-S_*)(R-R_*)) = \left[\mu(S-S_*)^2-\left(\frac{\delta du_1 Z_*}{\gamma}-2\mu\right) (S-S_*)(R-R_*)+\frac{\delta(\mu+\delta_0)}{2\gamma} (R-R_*)^2\right]$, $G_2((Z-Z_*)(R-R_*)) = \left[\frac{a_0}{2} (Z-Z_*)^2-\frac{\delta du_1S}{\gamma}(Z-Z_*)(R-R_*) +\frac{\delta(\mu+\delta_0)}{2\gamma}(R-R_*)^2\right]$ and $G_3((I(t-\tau_2)-I_*)(Z-Z_*)) = \left[a(I(t-\tau_2)-I_*)^2-a(I(t-\tau_2)-I_*)(Z-Z_*)+\frac{a_0}{2} (Z-Z_*)^2\right]$. Clearly, within the feasible region $\sum$, if following parametric conditions $a < min \{\delta, 2a_0\},$ $(du_1)^2 < \frac{a_0\gamma \mu^2 (\mu+\delta)}{\delta \Lambda^2}$ and $(\delta du_1Z_* -2\mu\gamma)^2 < 2\gamma \mu \delta (\mu+\delta_0)$ hold true, then $\dot{U_2}\leq0$. We found that $\dot{U_2} = 0$ if and only if at endemic equilibrium $E_2$ and $\dot{U_2} < 0$ otherwise within the feasible region $\sum$, therefore the singleton set $\{E_2\}$ is the largest positively invariant set contained in $\{(S, I, R, Z)\in\sum:\dot{U_2} = 0\}$. Hence, in the interior of $\sum$, Lyapunov LaSalle's theorem ^[62] ensures that the endemic equilibrium $E_2$ is globally asymptotically stable.

Example 3. In this example, the aforesaid global stability result of the unique endemic equilibrium is validated numerically for a set of the same parameters as considered in the Example 2. Similarly, the delay model have the unique endemic equilibrium $E_2 = (49.6, 6.67, 26.73, 0.66)$ and a disease free equilibrium $E_1 = (250, 0, 0, 0)$ (unstable) along with the basic reproduction number $\mathcal{R}_0 = 5.0403 > 1$. The parametric conditions for global stability are defined as: $a-min \{\delta, 2a_0\} = -0.19 < 0,$ $(du_1)^2-\frac{a_0\gamma \mu^2 (\mu+\delta)}{\delta \Lambda^2} = -3.1977\times 10^{-7} < 0$ and $(\delta d u_1Z_* -2\mu\gamma)^2-2\gamma \mu \delta (\mu+\delta_0) = -3.4041\times 10^{-5} < 0$. Notice that all the conditions are satisfied for this set of parameters which insures that the endemic equilibrium $E_2$ is globally stable for all time lags $\tau_1, \tau_2 > 0$ (Theorem 9).

6.   Case-III: $\tau_1 > 0$ and $\tau_2 > 0$

Here, we study the stability of endemic equilibrium $E_2$ in the presence of both delays $\tau_1$ and $\tau_2$. In the following, we first state the result that infers about the sign of the real part of roots of the characteristic Eq (3.5).

Lemma 5. ^[65] If all roots of Eq (5.1) have negative real parts for $\tau_2 > 0$, then there exists a ${\tau_1}_0 (\tau_2) > 0$ such that all roots of Eq (3.5) have negative real parts when $\tau_1 < {\tau_1}_0 (\tau_2)$.

Proof. Proof follows similarly as of Lemma 6 of ^[65].

Theorem 10. Let $\mathcal{R}_0 > 1$ and ${\tau_2}_0 = \frac{1}{{\omega_2}_0 }[\arccos(\Phi({\omega_2}_0))]$, where $\Phi({\omega_2}_0)$ is given in Eq (5.7). Then for any $\tau_2 < {\tau_2}_0$ there exists a ${\tau_1}_0 (\tau_2) > 0$ such that the unique endemic equilibrium $E_2$ is locally asymptotically stable for $\tau_1 < {\tau_1}_0 (\tau_2)$ and $\tau_2 < {\tau_2}_0$.

Proof. It follows from Theorem 7 that the Eq (5.1) has all negative real parts for $\tau_2 < {\tau_2}_0$. Now using Lemma 5, the result follows.

Similarly, we state the following result for the delay $\tau_1$.

Theorem 11. Let $\mathcal{R}_0 > 1$ and ${\tau_1}_0 = \frac{1}{{\omega_1}_0 }[\arccos(\Upsilon({\omega_1}_0))]$, where $\Upsilon({\omega_1}_0)$ is given in Eq (4.7). Then for any $\tau_1 < {\tau_1}_0$ there exists a ${\tau_2}_0 (\tau_1) > 0$ such that the unique endemic equilibrium $E_2$ is locally asymptotically stable for $\tau_2 < {\tau_2}_0 (\tau_1)$.

Remark 2. From Theorem 10, for the stability, an explicit form for threshold value of the ${\tau_1}_0$ can not obtain rather it only informs the existence of such value. Using numerical experimentation, we shall find the region of stability in the presence of both delays.

In order to find the stability region when both delays are present, we simulate the delay model system for parametric values as considered in the pervious section (Hopf-Hopf bifurcation).

Further, we shall vary both delays as $0.25\leq\tau_1\leq 15$ and $0.25\leq\tau_2\leq 50$ to obtain the stability region numerically and the corresponding result is shown in Figure 20. The stable endemic equilibrium is plotted in '$\ast$' in $\tau_2\tau_1$-plane and $\diamond$ is showing otherwise.

This result is more precise for the above-considered set of parameters than the one established in Theorem 10. In addition, we would like to discuss here that, in general delay destabilises the system but from Theorems 10 and 11, it is clear that the system is still stable when both delays are present under some restrictions. Also, from the stability region (Figure 20), if we consider $\tau_1 = 2$ (small) and $\tau_2 = 32$ (large), one can find the system is stable whereas if we change it, the stability may change. Thus, stability or stability switches are important insights from the biological point of view as eradication of disease may be possible in the case of stability whereas, in the case of oscillations, it will be challenging.

6.1.   Global stability of unique endemic equilibrium $E_2$

In the previous sections, we have performed the local analysis of the delayed model system but due to mathematical complexity, it is very difficult to obtain the explicit parametric relation for local stability. Therefore, in this section, we shall establish the global stability of the unique endemic equilibrium $E_2$ by constructing a suitable Lyapunov function.

Theorem 12. If $\mathcal{R}_0 > 1$, $\tau_1 > 0$ and $\tau_2 > 0$, the unique endemic equilibrium $E_2$ of the model system (2.1) is globally asymptotically stable under following parametric conditions $a < min \{\delta, 2a_0\},$ $\delta_0^2 < \frac{4\gamma \mu^2}{3(2\mu+\delta)},$ $(du_1)^2 < \frac{2a_0\gamma \mu^2}{3\Lambda^2 (2\mu+\delta)}$ and $((2\mu+\delta)du_1Z_* -2\mu\gamma)^2 < \frac{4}{3}\gamma \mu^2 (2\mu+\delta)$.

Proof. For global stability, we consider the following positive definite function in $\sum$,

The positive constants $K_1$ and $K_2$ will be selected suitably later. Now, we shall differentiate the function $U_1$ with respect to $t$ along the solution trajectories of the model (2.1) and using the parametric relationship of equilibrium points. We have,

Now, selecting and replacing the values of positive constants as: $K_1 = \frac{2\mu+\delta}{\gamma}$ and $K_2 = \frac{2\mu+\delta}{\beta}$, we have

where $H_1((S-S_*)(R-R_*)) = \left[\mu(S-S_*)^2-\left(\frac{2\mu+\delta}{\gamma} du_1 Z_*-2\mu\right)(S-S_*)(R-R_*)+\frac{2\mu+\delta}{3\gamma} \mu (R-R_*)^2\right]$, $H_2((R-R_*)(R(t-\tau_1)-R_*)) = \left[\frac{2\mu+\delta}{3\gamma} \mu (R-R_*)^2+\frac{2\mu+\delta}{\gamma} \delta_0 (R-R_*)(R(t-\tau_1)-R_*)+\mu(R(t-\tau_1)-R_*)^2\right]$, $H_3((Z-Z_*)(R-R_*)) = \left[\frac{a_0}{2} (Z-Z_*)^2-\frac{2\mu+\delta}{\gamma}du_1S(Z-Z_*)(R-R_*)+\frac{2\mu+\delta}{3\gamma} \mu (R-R_*)^2\right]$ and $H_4((I(t-\tau_2)-I_*)(Z-Z_*)) = \left[a(I(t-\tau_2)-I_*)^2-a(I(t-\tau_2)-I_*)(Z-Z_*)+\frac{a_0}{2} (Z-Z_*)^2\right]$. Note that within the feasible $\sum$, $\dot{U_3}\leq0$ if following parametric conditions $a < min \{\delta, 2a_0\},$ $\delta_0^2 < \frac{4\gamma \mu^2}{3(2\mu+\delta)},$ $(du_1)^2 < \frac{2a_0\gamma \mu^2}{3\Lambda^2 (2\mu+\delta)}$ and $((2\mu+\delta)du_1Z_* -2\mu\gamma)^2 < \frac{4}{3}\gamma \mu^2 (2\mu+\delta)$ hold true. It is observed that $\dot{U_3} = 0$ if and only if at endemic equilibrium $E_2$ and $\dot{U_3} < 0$ otherwise within the feasible region $\sum$, therefore the singleton set $\{E_2\}$ is the largest positively invariant set contained in $\{(S, I, R, Z)\in\sum:\dot{U_3} = 0\}$. Hence, in the interior of $\sum$, Lyapunov LaSalle's theorem ^[62] ensures that the endemic equilibrium $E_2$ is globally asymptotically stable.

Example 4. Here, we perform the numerical validation of the aforesaid global stability result stated in Theorem 12 when the model has a unique endemic equilibrium when $\mathcal{R}_0 > 1$. For this purpose, we consider a set of the same parameters as mentioned in Example 2. For this set of parameters, the delay model has the unique endemic equilibrium $E_2 = (49.6, 6.67, 26.73, 0.66)$ and a disease free equilibrium $E_1 = (250, 0, 0, 0)$ (unstable) along with the basic reproduction number $\mathcal{R}_0 = 5.0403 > 1$. Now, we calculate the parametric conditions numerically and are given as: $a-min \{\delta, 2a_0\} = -0.19 < 0,$ $\delta_0^2-\frac{4\gamma \mu^2}{3(2\mu+\delta)} = -7.3765\times 10^{-5} < 0,$ $(du_1)^2-\frac{2a_0\gamma \mu^2}{3\Lambda^2 (2\mu+\delta)} = -1.9730\times 10^{-7} < 0$ and $((2\mu+\delta)du_1Z_* -2\mu\gamma)^2- \frac{4}{3}\gamma \mu^2 (2\mu+\delta) = -1.2844\times 10^{-5} < 0$. Clearly, for this set of parameters, all the parametric relations hold true which insures that the endemic equilibrium $E_2$ is globally stable for all time lags $\tau_1, \tau_2 > 0$ (Theorem 12).

Finally, a rich and complex dynamics are obtained for the model in the presence of time delays. Whereas, in the absence of one time delay either $\tau_1$ or $\tau_2$, the delay model system shows the occurrence of oscillations through Hopf bifurcation. Further, in the presence of $\tau_2$ only ($\tau_1 = 0$), the existence of Hopf-Hopf bifurcation is established at two different delays for the delay model. In addition, global analysis is also established for various cases by constructing suitable Lyapunov functions which is itself very mathematically challenging due to high mathematical complexity.

7.   Conclusions

The effects of delayed information ($\tau_2$) and the delay in immunity loss ($\tau_1$) are studied for infectious diseases by formulating a delay mathematical model. To describe the dynamics of the delayed information, a rate equation is considered in which the density of information depends on the infective population as well as social media. Model analysis is performed and stability of equilibrium points is carried out for different scenarios of delay parameters. We observed that when the rate of immunity loss is small, the disease free equilibrium will be always stable otherwise it becomes unstable at a large rate. This happens due to the delay effect considered in the waning of immunity, if the delay parameter is smaller than a threshold quantity for $\mathcal{R}_0 < 1$, the disease free equilibrium is stable whereas it loses its stability as the delay crosses the threshold quantity. Generally, such type of phenomenon does not appear in non-delay models. In addition, when the effect of both delays is small enough, the unique endemic equilibrium is found locally stable for $\mathcal{R}_0 > 1$ and under some parametric conditions. Moreover, when both the delays ($\tau_1$ and $\tau_2$) cross critical values (either in presence of a single delay only or both delays), the existence of periodic solutions (via Hopf bifurcation) is investigated in this study. This leads to the existence of oscillatory persistence of disease within the population and hence, execution of control measures becomes challenging due to high and low epidemic peaks. Such Hopf bifurcation results follow the observation posed by a few researchers such as Yadav and Srivastava ^[54], and Zhang et al. ^[55] (single delay cases with the constant effect of information) but in the more general context as our model considered the dynamic effect of information with two discrete time delays. Our findings also follow similar facts as discussed by Lv et al. ^[43] for the delay considered in COVID-19 booster dose and antibody failure. Furthermore, the existence of multiple stability switches, termed as Hopf-Hopf (double) bifurcation at two different values of delay, is also examined for the case when only $\tau_2$ is present ($\tau_1 = 0$). In this case, at the first critical value of the $\tau_2$, the system is destabilised and undergoes the first Hopf bifurcation (the endemic equilibrium loses stability) whereas as $\tau_2$ increases, the system again undergoes the second Hopf bifurcation (the endemic equilibrium regains the stability) at another critical value of the delay. Therefore, both delays (reporting of infective and delay in waning the immunity) not only destabilise the system but also retain its stability which is very important as the contraction of infection is possible. Here, we infer that our obtained stability switch result falls in line with the result discussed by Kumar and others ^[47] and by Barbarossa and others ^[48] but with only single delay cases and different biological context. Therefore, our study also establishes similar kinds of important results under a simple mathematical form of delay factors considered in the model system to the other studies. Moreover, the global stability of the unique endemic equilibrium point is established by constructing Lyapunov functions in various cases for time delays which is itself difficult and challenging under high mathematical complexity raised by two crucial delays. The global stability result infers that the disease will persist in an endemic sense irrespective of delay impacts and can be suppressed if the basic reproduction number is below one. Mostly, delay models with a single delay, in this line of work, explore the local bifurcation analysis, not global properties. Finally, our study infers that the model system exhibits rich and complex dynamics locally as well as globally, and provides important insight. As, during the outbreak, optimal usage and allocation of suitable control measures with minimum economic damage are very difficult in any delay environment, hence, the corresponding optimal control problems may be possible pathways to solve such problems. In coming future, we shall try to address such problems.

Acknowledgments

The financial support of Anuj Kumar is from TIET, Patiala project no. TU/DORSP/57/7284. Yasuhiro Takeuchi acknowledges the financial support from JSPS "Grand-in-Aid 20K03755".

Conflict of interest

The authors declare there is no conflict of interest.