Codimension-two bifurcation analysis at an endemic equilibrium state of a discrete epidemic model

1.   Introduction

Humans have suffered from a number of contagious diseases over the ages, including cholera, influenza, and plague ^[1,2,3]. Contagious diseases have long been ranked alongside conflicts and a food shortage as key threats to human advancement and existence. The transmission of contagious diseases in populations, as well as how to prevent and eradicate them, are crucial and essential topics. Mathematical analysis and modelling is an important part of infectious diseases epidemiology. Applications of mathematical models to disease surveillance data can be used to address both scientific hypotheses and disease control policy questions. The mathematical description of disease epidemics immediately leads to several useful results, including the expected size of an epidemic and the critical level that is needed for an interaction to achieve effective disease control. There are several mathematical models proposed by eminent mathematicians to investigate what transpires when a population becomes infected, and under what conditions depending on the circumstances, the disease will be eradicated. Medical professionals have created vaccinations against a variety of viruses and suggested numerous epidemic preventive strategies as part of humanity's fight against contagious diseases. Furthermore, mathematicians have significantly aided in the effort to stop the spread of disease. From the standpoint of mathematical models, it is possible to roughly determine the duration between the contagious disease's outbreak and containment, the number of individuals infected at a given point in the disease's development who require quarantine, and the number of individuals in contact with the disease at any given time. After examining the number of plague cases and patient survival days, Kermack and Mckendrick ^[4] developed the SIR model, a ground-breaking mathematical epidemiology model. The great majority of research that have examined contagious diseases mathematically up to this point have operated under the shadow of this model. Many academics used the contagious disease dynamics model as a foundational research tool in the study and forecasting of the COVID-19 epidemic. They also suggested a number of enhancement techniques to enhance the classical contagious disease dynamics model, allowing it to more accurately describe the real context of the epidemic's transmission and produce more plausible forecasts regarding the epidemic's development trend. Numerous practical recommendations for governance, control, and prevention have been made ^[5,6,7]. Dynamical system's research has gained significant interest in the last several years. Dynamical systems is an interdisciplinary area that has several applications, including predator-prey models and tumor models, in addition to the study of epidemic diseases ^[8,9]. We typically use mathematical models to characterize the rich dynamic behavior of epidemic diseases while some mathematical models are tailored to specific diseases, where the majority are appropriate for broad investigations into the principles underlying different epidemic diseases ^[10,11]. Standard differential equations are used in the majority of epidemic disease models. Nonetheless, discrete-time dynamic models are far simpler, and more computationally efficient than continuous models. In addition to being straightforward analogs of continuous epidemic models, the majority of discrete-time epidemic models also exhibit intricate dynamic characteristics that the corresponding continuous models are unable to display. Even in a one-dimensional instance, the discrete-time model can produce incredibly complex dynamics ^[12], and for dynamical properties of higher-dimensional epidemic models we refer the reader to work of eminent researchers ^{[13,14,15,16,17]}. Therefore, from many years, mathematical infectious models have been a popular and interesting topic ^{[18,19,20,21,22]}. From an epidemiological point of view, epidemic dynamics are a extremely transited topic of investigation. The majority of researchers investigated bifurcation phenomena when a single systemic parameter changes. Indeed, many practical models include a number of systemic parameters when more than single systemic parameter is altered simultaneously, and so it is probable that complicated bifurcation such as codimension-two bifurcations are likely to occur. Nevertheless, because of the impact of higher-order nonlinear components, codimension-two bifurcations remain extremely difficult to understand. Even elementary dynamical systems have complex dynamics that cannot be satisfactorily illustrated by theoretical analysis. A few scholars have investigated simulating the dynamics within the local parameter fluctuation using computers. As a result, the numerical techniques allow to better illustrate and comprehend the dynamics of the model in addition to validating our analytical results ^{[23,24,25,26]}. On the other hand, in recent years, many researchers have investigated the codimension-one and codimension-two bifurcations of discrete model by bifurcation theory. For instance, Ruan and Wang ^[27] have studied Bogdanov-Takens bifurcation for the model:

where $d, \gamma$ and $v$, respectively denote death rate, recovery rate and removed individuals rate whereas $\nu$ is a nonnegative constant. Eskandari and Alidousti ^[28] have examined codimension-two bifurcations of the following discrete model:

where $d, \lambda, A$ and $r$ denote the natural death rate, the bilinear incidence rate, the recruitment rate of the population, and the recovery rate of the infective individuals, respectively. Ruan et al. ^[29] have examined codimension-two bifurcations of the following discrete model:

where $I_{t}$ and $S_{t}$ denote infected and susceptible individuals, respectively. Abdelaziz et al. ^[30] have examined codimension-two bifurcations of the following discrete model:

where the parameter $h > 0$ is the time step size. Chen et al. ^[31] have examined codimension-two bifurcations of the following discrete model:

where $\gamma, \bigwedge, \ \mu, \ \beta$ and $\phi$, respectively denote disease related death rate, recruitment rate, natural death rate, disease transmission coefficient, and rate at with individuals $I$ return to class $S$, and $\delta > 0$ is a integral step size. Liu et al. ^[32] have examined codimension-two bifurcations of following discrete model:

where $A, \ d, \ r$, respectively denote recruitment rate, natural death rate, recovery rate, and finally, $\lambda$ denotes bilinear incidence rate. Yi et al. ^[33] have examined codimension-two bifurcations of the model:

where $\eta, \ d, N_{0}, \ \beta > 0$ denote effective contact rate, death rate, rate of recruitment and recovery rate, respectively. Ma and Duan ^[34] have explored codimension-two bifurcations of a two-dimensional discrete time Lotka-Volterra predator-prey model. Yousef et al. ^[35] have explored codimension-one and codimension-two bifurcations in a discrete Kolmogorov type predator-prey model. Eskandari et al. ^[36] have explored codimension-two bifurcations of a discrete game model. Guo et al. ^[37] have examined hopf bifurcations of a bioeconomic model. Inspired by the aforementioned research, in this paper, we aim to examine codimension-two bifurcations of the following discrete epidemic model with vital dynamics and vaccination ^[38]:

where $\mu_3$, $\mu_1$, $\mu_4$ and $\mu_6$, respectively denote natural death rate, contact rate, cure rate, and rate of immunity loss while $\mu_7$ and $\mu_5$ are rates of vaccination in individuals $S_t$ and newcomers. More precisely, our goal of this paper is to examine the existence of codimension-two bifurcation sets, and codimension-two bifurcation at endemic equilibrium state (EES) associated with $1:2$, $1:3$ and $1:4$ strong resonances of a discrete epidemic model (1.8). Furthermore, our theoretical results are confirmed by numerical simulation.

The organization of the paper is as follows: The existence of codimension-two bifurcation sets at EES of a discrete epidemic model (1.8) are identified in Section 2 whereas Section 3 is about the study of codimension-two bifurcation at EES. In order to confirm theoretical results, simulations are presented in Section 4 whereas conclusion is given in Section 5.

2.   Codimension-two bifurcation sets at endemic equilibrium state

In this section, we examine codimension-two bifurcation sets at EES for the discrete epidemic model (1.8). For this, the simple calculation shows that if $\mu_6 > \frac{\left(\mu_3+\mu_7\right)\left(\mu_3+\mu_4\right)-\left(1-\mu_5\right)\mu_3 \mu_1}{\mu_1 - \left(\mu_3 +\mu_4\right)}$ then model (1.8) has EES $\left(\frac{\left(\left(1-\mu_5\right)\mu_3 +\mu_6\right)\mu_1 \mu_2 -\left(\mu_3 +\mu_7 +\mu_6\right)\left(\mu_3 +\mu_4 \right)\mu_2}{\mu_1\left(\mu_3 +\mu_6\right)}, \frac{\left(\mu_3 +\mu_4\right)\mu_2}{\mu_1}\right)$ with basic reproduction number is $\mathfrak{R}_0: = \frac{\left(\left(1-\mu_5\right)\mu_3+\mu_6\right)\mu_1}{\left(\mu_3+\mu_6+\mu_7\right)\left(\mu_3+\mu_4\right)} > 1$. Now variation matrix $V|_{\mathrm{EES}}$ of the linearized system of model (1.8) at EES is

with characteristic equation is

where

Setting

where

Finally, the roots of (2.2) are

where

Hereafter, following three cases are to considered in order to get codimension-two bifurcations sets:

Case 1. If $H = 4 = -G$, then from (2.6) one gets $\lambda_{1, 2} = -1$ with $\mu_7 = \frac{4\left(\mu_3 +\mu_6\right)-\left(\mu_3 +\mu_6\right)^2 -4}{\left(\mu_3 +\mu_6\right)}$ and $\mu_5 = \frac{\mu_1\left(\mu_3 +\mu_6\right)^2 -4\left(\mu_3 +\mu_6\right)\left(\mu_3 +\mu_4\right)+4\left(\mu_4 -\mu_6\right)}{\mu_1 \mu_3 \left(\mu_3 +\mu_6\right)}$. Therefore, at EES, model (1.8) may undergoes codimension-two bifurcation with $1:2$ strong resonance where

Case 2. If $H = 3 = -G$ then from (2.6) one gets $\lambda_{1, 2} = \frac{-1\pm\sqrt{3} \iota}{2}$ with $\mu_7 = \frac{3\left(\mu_3 +\mu_6\right)-\left(\mu_3 +\mu_6\right)^2 -3}{\mu_3 +\mu_6}$ and $\mu_5 = \frac{\mu_1\left(\mu_3 +\mu_6\right)^2 -3\left(\mu_3 +\mu_6\right)\left(\mu_3 +\mu_4\right)+3\left(\mu_4 -\mu_6\right)}{\mu_1 \mu_3 \left(\mu_3 +\mu_6\right)}$. Therefore, at EES, model (1.8) may undergoes codimension-two bifurcation with $1:3$ strong resonance where

Case 3. If $H = 2 = -G$ then from (2.6) one gets $\lambda_{1, 2} = \pm \iota$ with $\mu_7 = \frac{2\left(\mu_3 +\mu_6\right)-\left(\mu_3 +\mu_6\right)^2 -2}{\mu_3 +\mu_6}$ and $\mu_5 = \frac{\mu_1\left(\mu_3 +\mu_6\right)^2 -2\left(\mu_3 +\mu_6\right)\left(\mu_3 +\mu_4\right)+2\left(\mu_4 -\mu_6\right)}{\mu_1 \mu_3 \left(\mu_3 +\mu_6\right)}$. Therefore, at EES, model (1.8) may undergoes codimension-two bifurcation with $1:4$ strong resonance where

Remark 1. If $H = 0 = G$ then $\lambda_{1, 2} = 1$. However, from (2.5), it is noted that $G$ and $H$ are never zero and so, at EES, model (1.8) does not undergo codimension-two bifurcation with $1:1$ resonance.

3.   Codimension-two bifurcation at EES

The codimension-two bifurcations at EES of model (1.8) will be examined in this section by bifurcation theory ^{[39,40,41,42,43]}.

3.1.   $1:2$ strong resonance at EES: $\lambda_{1, 2}=-1$

From (2.8), if $\mu_7 = \frac{4\left(\mu_3 +\mu_6\right)-\left(\mu_3 +\mu_6\right)^2 -4}{\mu_3 +\mu_6}$ and $\mu_5 = \frac{\mu_1\left(\mu_3 +\mu_6\right)^2 -4\left(\mu_3 +\mu_6\right)\left(\mu_3 +\mu_4\right)+4\left(\mu_4 -\mu_6\right)}{\mu_1 \mu_3 \left(\mu_3 +\mu_6\right)}$ then calculation shows that $\lambda_{1, 2}|_{\mu_7 = \frac{4\left(\mu_3 +\mu_6\right)-\left(\mu_3 +\mu_6\right)^2 -4}{\left(\mu_3 +\mu_6\right)}, \ \mu_5 = \frac{\mu_1\left(\mu_3 +\mu_6\right)^2 -4\left(\mu_3 +\mu_6\right)\left(\mu_3 +\mu_4\right)+4\left(\mu_4 -\mu_6\right)}{\mu_1 \mu_3 \left(\mu_3 +\mu_6\right)}} = -1$ which implies that if $\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7 \right)\in \mathcal{F}_{12}|_{\mathrm{EES}}$ then at $\mathrm{EES}$ model (1.8) may undergoes codimension-two bifurcation with $1:2$ strong resonance, by choosing $\mu_5$ and $\mu_7$ as bifurcation parameters. Now using following transformation, EES of model (1.8) transform to $\left(0, 0\right)$:

where $I^* = \frac{\left(\left(1-\mu_5\right)\mu_3 +\mu_6\right)\mu_1 \mu_2 -\left(\mu_3 +\mu_7 +\mu_6\right)\left(\mu_3 +\mu_4 \right)\mu_2}{\mu_1\left(\mu_3 +\mu_6\right)}$ and $S^* = \frac{\left(\mu_3 +\mu_4\right)\mu_2}{\mu_1}$. In view of (3.1), one write the model (1.8) as follows:

Now on expanding (3.2) at $\left(0, 0\right)$ up to order-$2^{nd}$, one gets:

Now at EES, (3.3) becomes

where

From (3.4), if one denotes

then

where $\varrho = \left(\mu_5, \mu_7 \right)$ with characteristic roots are $\lambda_{1, 2} = -1$. Furthermore, eigenvector and generalized eigenvector of $A_{0}$ corresponding to characteristic roots $-1$, respectively are $q_{0} = \left(\begin{array}{cc} 1\\ -\frac{\mu_3 +\mu_6}{2} \end{array} \right)$ and $q_{1} = \left(\begin{array}{cc} 1\\ -\frac{\mu_3 +\mu_6}{4} \end{array} \right)$. Additionally, eigenvector and generalized eigenvector of $A_{0}^T$ corresponding to characteristic roots $-1$, respectively are $p_{1} = \left(\begin{array}{cc} 2\\ \frac{4}{\mu_3 +\mu_6} \end{array} \right)$ and $p_{0} = \left(\begin{array}{cc} -1\\ -\frac{4}{\mu_3 +\mu_6} \end{array} \right)$, where $p_i, \ q_i\ (i = 0, 1)$ satisfying the following relations:

Now if

with $x = \left(\begin{array}{cc} u_{t}\\ v_{t} \end{array} \right)$ then straightforward calculation yields

Now in coordinates $\left(n_{t}, m_{t}\right)$, the model general representation of (3.4) is

where

and

From (3.9), (3.12) and (3.13), the calculation yields

and

Moreover, the calculation shows that $a\left(\varrho_{0}\right) = b\left(\varrho_{0}\right) = c\left(\varrho_{0}\right) = d\left(\varrho_{0}\right) = 0$. Now denote

with the following non-singular coordinate transformation:

model (3.17) can be writen as

where

and

Now if

then $\beta_{1}\left(\varrho_{0}\right) = \beta_{2}\left(\varrho_{0}\right) = 0$, So, (3.18) along with (3.19)–(3.21) becomes

where

and

$g_{jk} = h_{jk} = 0\ \forall\ j, k > 0$ and $j+k = 3$. Furthermore, by employing the transformation:

with

one gets the following $1:2$ resonance normal form:

where

Based on above analysis, one has the following result:

Theorem 3.1. If $\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7\right)\in \mathcal{F}_{12}|_{\mathrm{EES}}$ and the discriminatory quantities, which are depicted in (3.28), that is, $C\not = 0, \ D+3 C = \left(\frac{\mu_1}{\mu_2}\right)^2\left(4\times \left(\mu_3 +\mu_6\right)^2 -9\left(\mu_3 +\mu_6\right)+5\right)\neq0$, then at EES, model (1.8) undergoes $1:2$ strong resonance. Additionally, EES is elliptic (respectively, saddle) if $C > 0$ (respectively, $C < 0$), and near $1:2$ point $D+3 C\neq0$ defines the following bifurcation curves:

(i) Pitchfork bifurcation curve

and furthermore, nontrivial equilibrium state exists for $\beta_1 < 0$;

(ii) Heteroclinic bifurcation curve:

(iii) Non-degenerate N-S bifurcation curve:

(iv) Homologous bifurcation curve:

Remark 1. The occurrence of codimension-two bifurcation associated with $1:2$ strong resonance at EES of model (1.8) indicates the complex dynamical behavior if $\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7\right)\in \mathcal{F}_{12}|_{\mathrm{EES}}$. An important biological consequence of the non-degenerate N-S bifurcation is the existence of periodic or quasiperiodic oscillations between individuals $I$ and $S$ if $\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7 \right)\in \mathcal{F}_{12}|_{\mathrm{EES}}$. Furthermore, periodic oscillations or the homoclinic structure may cause long-period oscillations or even chaos in the $I$ and $S$ individuals.

3.2.   $1:3$ strong resonance at EES: $\bf \lambda_{1, 2}=\frac{{-1}\pm\sqrt{3} \iota}{2}$

From (2.9), if $\mu_7 = \frac{3\left(\mu_3 +\mu_6\right)-\left(\mu_3 +\mu_6\right)^2 -3}{\mu_3 +\mu_6}$ and $\mu_5 = \frac{\mu_1\left(\mu_3 +\mu_6\right)^2 -3\left(\mu_3 +\mu_6\right)\left(\mu_3 +\mu_4\right)+3\left(\mu_4 -\mu_6\right)}{\mu_1 \mu_3 \left(\mu_3 +\mu_6\right)}$ then from (2.1) one gets:

with $\lambda_{1, 2} = \frac{{-1}\pm\sqrt{3} \iota}{2}$ which implies that if $\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7 \right)\in \mathcal{F}_{13}|_{\mathrm{EES}}$ then at $\mathrm{EES}$, model (1.8) may undergoes codimension-two bifurcation with $1:3$ strong resonance. Furthermore, eigenvector and adjoint eigenvector of $A_{0}$ corresponding to eigenvalues $\frac{-1\pm\sqrt{3 }\iota}{2}$, respectively are $q = \left(\begin{array}{cc} \frac{3}{ \mu_3 +\mu_6}\\\frac{3-\sqrt{3} \iota}{2} \end{array} \right)$ and $p = \left(\begin{array}{cc} \frac{(\mu_3 +\mu_6)(1+\sqrt{3} \iota)}{6 }) \\ \frac{-\sqrt{3 }\iota}{3} \end{array} \right)$ satisfying

Now if $Z = \left(\begin{array}{cc} u_{t}\\ v_{t} \end{array} \right) \in \mathbb{R}^2$ where it can be represented by $Z_t = z_tq+\bar{z}_t\bar{q}$ then (3.4) becomes

where

with

and

Now following transformation

along with its inverse transformation is utilized in order to eliminate quadratic terms from (3.35), where it becomes

with

Now, it should be noted that quadratic terms of (3.40) should be vanished if

Furthermore, the following transformation yields to annihilate cubic terms:

Using (3.43) along with inverse transformation, from (3.40), one gets:

where

On setting

the ingredients from (3.45) should be zero except $r_{21}$. Therefore, the desired normal form of $1:3$ resonance is

where

Finally, let

then one has the following theorem for the codimension-two bifurcation with $1:3$ strong resonance:

Theorem 3.2. If the discriminatory quantities, which are depicted in (3.49), that is, $C_{1}\neq 0$ and $\Re\left(D_{1}\right)\neq 0$ then at EES, model (1.8) undergoes codimension-two bifurcation with $1:3$ strong resonance as the parameters $\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7\right)\in \mathcal{F}_{13}|_{\mathrm{EES}}$. Additionally, if $\Re\left(D_{1}\right)\neq 0$ examine the bifurcation behavior then at EES, model (1.8) has the following dynamical characteristics:

(i) If $\Re\left(D_{1}\right) > 0$ then invariant closed curve occurs at $1:3$ resonance point is unstable;

(ii) If $\Re\left(D_{1}\right) < 0$ then invariant closed curve occurs at $1:3$ resonance point is stable;

(iii) At trivial equilibrium state of (3.35) one has the non-degenerate N-S bifurcation.

Remark 2. If $\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7\right)\in \mathcal{F}_{13}|_{\mathrm{EES}}$ then there exists codimension-two bifurcation associated with $1:3$ strong resonance. From a biological perspective, periodic or quasi-periodic oscillation may occur in individuals $I$ and $S$ as a result of non-degenerate N-S bifurcation.

3.3.   $1:4$ strong resonance : $\bf \lambda_{1, 2}=\pm\iota$

From (2.10), if $\mu_7 = \frac{2\left(\mu_3 +\mu_6\right)-\left(\mu_3 +\mu_6\right)^2 -2}{\left(\mu_3 +\mu_6\right)}$ and $\mu_5 = \frac{\mu_1\left(\mu_3 +\mu_6\right)^2 -2\left(\mu_3 +\mu_6\right)\left(\mu_3 +\mu_4\right)+2\left(\mu_4 -\mu_6\right)}{\mu_1 \mu_3 \left(\mu_3 +\mu_6\right)}$ then from (2.1) one gets:

with $\lambda_{1, 2} = \pm\iota$ which implies that if $\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7 \right)\in \mathcal{F}_{14}|_{\mathrm{EES}}$ then at $\mathrm{EES}$, model (1.8) may undergoes codimension-two bifurcation with $1:4$ strong resonance. Furthermore, eigenvector and adjoint eigenvector of $A_{0}$ corresponding to eigenvalues $\pm\iota$, respectively are $q = \left(\begin{array}{cc}-\frac{2}{\mu_3 +\mu_6}\\ 1-\iota \end{array} \right)$ and $p = \left(\begin{array}{cc}\frac{(\mu_3 +\mu_6)(-1-\iota)}{4}\\ \frac{-\iota}{ 2} \end{array} \right)$ satisfying

Now if $Z = \left(\begin{array}{cc} u_{t}\\ v_{t} \end{array} \right) \in \mathbb{R}^2$ where it can be represented by $Z_t = z_tq+\bar{z}_t\bar{q}$ then (3.4) becomes

where

with

and

Now the transformation, which is depicted in (3.39) along its inverse transformation, is utilized to eliminate quadratic terms from (3.52), where it becomes

with same $\sigma$'s as in (3.41). Now, it should be noted that quadratic terms of (3.56) should be vanished if

Using (3.43) along with inverse transformation, from (3.56), one gets:

where

On setting

the ingredients from (3.59) should be zero except $r_{21}$ and $r_{03}$. Therefore, the desired normal form of $1:4$ resonance is

where

Let

Now if $D_{1}(\mu_5, \mu_7)\neq 0$ and

then one has the following theorem for the codimension-two bifurcation with $1:4$ strong resonance:

Theorem 3.3. If the discriminatory quantity, which is depicted in (3.64), such that, $\Re\left(B(\mu_5, \mu_7)\right)\neq0$ and $\Im\left(B(\mu_5, \mu_7)\right)\neq0$ then at EES, model (1.8) undergoes codimension-two bifurcation with $1:4$ strong resonance as $\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7\right)\in \mathcal{F}_{14}|_{\mathrm{EES}}$. Additionally, $B(\mu_5, \mu_7)$ determine the bifurcation behavior near EES of model (1.8) and so, there are two-parameter families of equilibrium state of order four bifurcation from EES near it. Depending on the choices made for $\mu_5$ and $\mu_7$, one of these families contains unstable, attracting, or repelling invariant circles. Furthermore, in a sufficiently small neighborhood of $(\mu_5, \mu_7)$, there exist numerous complex codimension-one bifurcation curves of (3.61).

(i) At trivial equilibrium state of (3.61) there is a N-S bifurcation. Moreover, there is an invariant circle if $\lambda = -\iota$ and invariant circle will disappear if $\lambda = \iota$;

(ii) If $\left|B(\mu_5, \mu_7)\right| > 1$ then there exists eight equilibrium states that disappear or appear in pairs via fold bifurcation at $\mu_5$ and $\mu_7$;

(iii) At eight equilibrium states, there exists N-S bifurcations. Furthermore, four small invariant circles bifurcate from equilibrium states, and vanish near the homoclinic loop bifurcation curve.

Remark 3. The presence of a non-degenerate N-S bifurcation is indicated by the occurrence of codeminsion-two bifurcation associated with $1:4$ strong resonance. In a specific parametric region, it is also feasible to generate an invariant cycle of period-$4$ orbit. In biology, the non-degenerate N-S bifurcation may give rise to periodic or quasi-periodic oscillations in individuals $I$ and $S$.

4.   Numerical simulations

Example 1. In this Example, it is proved numerically that if $\mu_1 = 3.9$, $\mu_2 = 0.5$, $\mu_3 = 0.28$, $\mu_4 = 0.22$, $\mu_6 = 0.9$ and varying $\mu_5\in[0.2, 2.9]$, $\mu_7\in[-0.6, 1.9]$ with $(I_0, S_0) = (0.03, 0.6)$ then at EES, model (1.8) undergoes codimension-two bifurcation with $1:2$ strong resonance. For this, in the following, first one need to prove the eigenvalues criterion for the existence of $1:2$ strong resonance holds. For this, if $\mu_1 = 3.9$, $\mu_2 = 0.5$, $\mu_3 = 0.28$, $\mu_4 = 0.22$, $\mu_6 = 0.9$ then from (2.8) one gets: $\mu_5 = 0.271900415968213, \mu_7 = -0.5698305084745758$. Therefore, if $\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7\right) = \left(3.9, 0.5, 0.28, 0.22, 0.271900415968213, 0.9, -0.5698305084745758\right)$ then model (1.8) has $\mathrm{EES} = \left(0.43459365493263796, 0.06410256410256411\right)$ and moreover, from (2.1) one gets:

with $\lambda_{1, 2} = -1$ and so $\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7\right) = \left(3.9, 0.5, 0.28, 0.22, 0.271900415968213, 0.9, -0.5\right. \\ \left. {}698305084745758\right)\in \mathcal{F}_{12}|_{\mathrm{EES} = (0.43459365493263796, 0.06410256410256411)}$, and finally, from (3.28) one gets $C = -19.207188000000006\not = 0$ and $D = 54.555228\not = 0$, $D+3 C = -3.0663359999999864\not = 0$ which imply that model (1.8) undergoes a codimension-two bifurcation with 1:2 strong resonance. Therefore, the simulation agree with the conclusion of Theorem 3.1. Hence, codimension-two bifurcation diagrams with $1:2$ strong resonance are drawn in Figure 1.

Example 2. In this Example, it is proved numerically that if $\mu_1 = 3.9$, $\mu_2 = 0.2$, $\mu_3 = 0.3985$, $\mu_4 = 0.75$, $\mu_6 = 0.95$ and varying $\mu_5\in[0.612332, 2.9]$, $\mu_7\in[-0.6, 0.855]$ with $(I_0, S_0) = (0.0669333, 0.0794103)$ then at EES, model (1.8) undergoes codimension-two bifurcation with $1:3$ strong resonance. For this, first one need to prove the eigenvalues criterion for the existence of $1:3$ strong resonance holds. So, if $\mu_1 = 3.9$, $\mu_2 = 0.2$, $\mu_3 = 0.3985$, $\mu_4 = 0.75$, $\mu_6 = 0.95$ then from (2.9) one gets: $\mu_5 = 0.8806815166411702$, $\mu_7 = -0.5731941045606225$. Therefore, if $\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7\right) = \left(3.9, 0.2, 0.3985, 0.75, 0.8806815166411702, 0.95, -0.5731941045606225\right)$ then model (1.8) has $\mathrm{EES} = \left(0.11408687715695504, 0.05889743589743591\right)$ and moreover, from (2.1) one gets:

with $\lambda_{1, 2} = \frac{{-1}\pm\sqrt{3} \iota}{2}$ and so $\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7\right) = \left(3.9, 0.2, 0.3985, 0.75, 0.8806815166411702, 0.95, -0.5731941045606225\right)\in \mathcal{F}_{13}|_{EES = (0.11408687715695504, 0.05889743589743591)}$, and finally, from (3.48) one gets $C = 25.470383759733036 + 27.228503386338268 \iota$ and $D = 527.560465001899 + 1209.8979374673806 \iota$. From (3.49), one gets $C_1 = 32.53615127919911 - 107.01675321970966 \iota \not = 0$ and $D_1 = -1818.3291845871258 - 3185.4892003729983\iota$ which shows that $\Re(D_{1}) = -1818.3291845871258 < 0$ which imply that model (1.8) undergoes a codimension-two bifurcation with 1:3 strong resonance. Therefore, the simulation agree with the conclusion of Theorem 3.2. Hence, codimension-two bifurcation diagrams with $1:3$ strong resonance are drawn in Figure 2.

Example 3. Finally, it is proved numerically that if $\mu_1 = 3.9$, $\mu_2 = 0.5$, $\mu_3 = 0.3985$, $\mu_4 = 0.45$, $\mu_6 = 0.95$ and varying $\mu_5\in[0.95607, 2.9]$, $\mu_7\in[-0.9, 0.855]$ with $(I_0, S_0) = (0.0446222, 0.0794103)$ then at EES, model (1.8) undergoes codimension-two bifurcation with $1:4$ strong resonance. For this, in the following, first one need to prove the eigenvalues criterion for the existence of $1:4$ strong resonance holds. So, if $\mu_1 = 3.9$, $\mu_2 = 0.5$, $\mu_3 = 0.3985$, $\mu_4 = 0.45$, $\mu_6 = 0.95$ then from (2.10) one gets: $\mu_5 = 1.8148732738022662, \mu_7 = -0.8316294030404152$. Therefore, if $\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7\right) = \left(3.9, 0.5, 0.3985, 0.45, 1.8148732738022662, 0.95, -0.8316294030404152\right)$ then model (1.8) has $\mathrm{EES} = \left(0.1901447952615917, 0.1087820512820513\right)$ and moreover, from (2.1) one gets:

with $\lambda_{1, 2} = \pm \iota$ and so $\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7\right) = \left(3.9, 0.5, 0.3985, 0.45, 1.8148732738022662, 0.95, -0.83162940304041\right)\in \mathcal{F}_{14}|_{\mathrm{EES} = (0.1901447952615917, 0.1087820512820513)}$. On the other hand, if $\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7\right)\in \mathcal{F}_{14}|_{EES = (0.1901447952615917, 0.1087820512820513)}$ then from (3.62) and (3.63), one gets $C = 0.7573645541146341 - 54.00102260452535 \iota$, $D = -51.459886377274955 + 9.38011362272504 \iota$, $C_1 = -216.00409041810144 - 3.0294582164584933\iota \not = 0$ and $D_1 = 37.52045449090016 + 205.8395455090998\iota \not = 0$. Finally, from (3.64), one gets $B = -1.032370271795126 - 0.014478997116506187 \iota$ as $\Re(B) = -1.032370271795126\not = 0$ and $\Im (B) = - 0.014478997116506187 \not = 0$ which imply that model (1.8) undergoes a codimension-two bifurcation with 1:4 strong resonance. Therefore, the simulation agree with the conclusion of Theorem 3.3. Hence, codimension-two bifurcation diagrams with $1:4$ strong resonance are drawn in Figure 3.

5.   Conclusions

The work is about the codimension-two bifurcation analysis of a discrete epidemic model (1.8) in the region $\mathbb{R}_+^2 = \left\{(I, S): I, \ S\ge 0\right\}$. It is proved that if $\mu_6 > \frac{\left(\mu_3+\mu_7\right)\left(\mu_3+\mu_4\right)-\left(1-\mu_5\right)\mu_3 \mu_1}{\mu_1 - \left(\mu_3 +\mu_4\right)}$ then model (1.8) has EES $\left(\frac{\left(\left(1-\mu_5\right)\mu_3 +\mu_6\right)\mu_1 \mu_2 -\left(\mu_3 +\mu_7 +\mu_6\right)\left(\mu_3 +\mu_4 \right)\mu_2}{\mu_1\left(\mu_3 +\mu_6\right)}, \frac{\left(\mu_3 +\mu_4\right)\mu_2}{\mu_1}\right)$. At EES of model (1.8), we first identified the codimension-two bifurcations sets associated with (i) $1:2$ strong resonance $\mathcal{F}_{12}|_{\mathrm{EES}}: = \left\{\left(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6, \mu_7\right):\mu_7 = \frac{4\left(\mu_3 +\mu_6\right)-\left(\mu_3+\mu_6\right)^2 -4}{\mu_3 +\mu_6}, \mu_5 = \frac{\mu_1\left(\mu_3 +\mu_6\right)^2 -4\left(\mu_3 +\mu_6\right)\left(\mu_3 +\mu_4\right)+4\left(\mu_4 -\mu_6\right)}{\mu_1 \mu_3 \left(\mu_3 +\mu_6\right)}\right\}$, (ii) $1:3$ strong resonance

(iii) $1:4$ strong resonance

and then we have studied detailed codimension-two bifurcations with $1:2$, $1:3$, and $1:4$ by bifurcation theory and series of affine transformations. Furthermore, we have also given biological interpretations of theoretical results. Finally, theoretical results are carried out numerically.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

The Authors Adil Jhangeer and Muhammad Bilal Riaz are highly thankful to Ministry of Education, Youth and Sports of the Czech Republic for their support through the e-INFRA CZ (ID: 90254).

Conflict of interest

The authors declare that they have no conflict of interest regarding the publication of this paper.