Bogdanov–Takens Bifurcation of Kermack–McKendrick Model with Nonlinear Contact Rates Caused by Multiple Exposures

Abstract

In this paper, we consider the influence of a nonlinear contact rate caused by multiple contacts in classical SIR model. In this paper, we unversal unfolding a nilpotent cusp singularity in such systems through normal form theory, we reveal that the system undergoes a Bogdanov-Takens bifurcation with codimension 2. During the bifurcation process, numerous lower codimension bifurcations may emerge simultaneously, such as saddle-node and Hopf bifurcations with codimension 1. Finally, employing the Matcont and Phase Plane software, we construct bifurcation diagrams and topological phase portraits. Additionally, we emphasize the role of symmetry in our analysis. By considering the inherent symmetries in the system, we provide a more comprehensive understanding of the dynamical behavior. Our findings suggest that if this occurrence rate is applied to the SIR model, it would yield different dynamical phenomena compared to those obtained by reducing a 3-dimensional dynamical model to a planar system by neglecting the disease mortality rate, which results in a stable nilpotent cusp singularity with codimension 2. We found that in SIR models with the same occurrence rate, both stable and unstable Bogdanov-Takens bifurcations occur, meaning both stable and unstable limit cycles appear in this system.

1. Introduction

Kermack–McKendrick models of SIR type are always used to describe infectious diseases where the immunity of infected individuals can persist,

$\begin{array}{cc}\hfill \frac{dS}{dt}& =A-\beta SI-dS,\hfill \\ \hfill \frac{dI}{dt}& =\beta SI-(d+\mu +\delta )I,\hfill \\ \hfill \frac{dR}{dt}& =\mu I-dR,\hfill \end{array}$

such as measles, smallpox, and polio, et al. In the model, the population is divided into three categories: susceptible (S), infected (I), and removed (R). A and d are the birth rate and natural death rate of people in the area, $\mu$ is the removal rate, $\delta$ denotes the disease-related death rate, $\beta$ denotes the effective contact rate between susceptible and infected person, which is the average probability of infection per infected person per contact with a susceptible person, and $\beta S$ is the the number of susceptible people who are infected by an infected person. The bilinear incidence rate $\beta SI$ is considered a mass action incidence, and it is considered to be correct in the early stages of an outbreak of diseases or a huge number of the population. During the recent COVID-19 pandemic, this incidence rate was widely adopted [1]. However, this is insufficient for many epidemics.

For example, the effective contact rate sometimes depends on the proportion of infected persons. To overcome the deficiencies, many nonlinear occurrence rates are used in the epidemic model [2,3,4,5]. However, these nonlinear incidence rates often imply more complex dynamical behaviors, such as Hopf bifurcation and nilpotent singularity bifurcations, et al. [6,7,8,9]. Moreover, in multi-timescale scenarios, infectious disease models can exhibit new bifurcation phenomena, such as canard explosions and relaxation oscillations, as discussed in works like [10,11], among others. Multiple exposures will increase the possibility of infection, saturation exposure, etc. These factors will significantly change the contact rate. Through experiment and analysis, P. van den Driessche et al. [12] proposed a nonlinear incidence rate that has multiple exposure effects, as follows:

(1)$\beta SI(1+\nu I).$

This paper [13] generalizes this incidence rate as

(2)$\beta SI(1+f(I,v\left)I\right).$

This paper finds that SIV models with this incidence may have a bi-stable phenomenon. By calculating the first Lyapunov coefficient, it is proved that Hopf bifurcation can occur in the system. However, the more degenerate bifurcation of singularities is not covered. Subsequently, the incidence rate (2) is applied to an SIR model with both nonlinear incidence and nonlinear clearance in [14]. It is proved that the system can generate Bautin bifurcation and Bogdanov–Takens bifurcation.

Susceptible individuals become infected due to multiple exposures to certain infected individuals.This incidence rate is used to describe the infection rate of multiple exposures. Compared with experimental results, this incidence rate can better simulate the incidence data of mumps. This incidence rate is frequently used in infectious disease models. Recently, in many COVID-19 models, this multi-exposure transmission rate is often employed to describe the infection of symptomatic and asymptomatic infectious. Ref. [15] provides the basic reproduction number under this transmission rate in a COVID-19 model. Meanwhile, Ref. [16] demonstrates that models of this type can exhibit bistability and even periodic phenomena in an SIRS model. Additionally, Ref. [17] establishes an infectious disease model in a two-patch environment with this incidence rate, analyzing various complex bifurcations in the system; especially, cusps with codimension 2 and codimension 3 are also found in this model. Other information about nilpotent singularity in the infectious disease model can be found in [18,19] et al.

The dynamic properties of the SIR model are explicit in the study of classical epidemic models. By applying the regenerative matrix method [20], it is easy to obtain the basic reproductive number $R_0=\frac{A{\beta }_1}{d(d+\mu +\delta )}.$ When $R_0<1$, the system only has disease-free equilibrium, which is globally stable;when $R_0>1$, there is a globally stable positive equilibrium, whereas the disease-free equilibrium becomes unstable. However, if we adopt the incidence rate (1) in the SIR model, will the system generate other interesting dynamic phenomena? To address this, we proceed to investigate the following model in the subsequent sections of this paper:

(3)$\begin{array}{c}\frac{dS}{dt}=A-dS-({\beta }_1+{\beta }_{0}I)SI,\hfill \\ \frac{dI}{dt}=({\beta }_1+{\beta }_{0}I)SI-(d+\mu +\delta )I,\hfill \\ \frac{dR}{dt}=\mu I-dR.\hfill \end{array}$

Note that the first two equations of the above system are independent, so we consider the following derived system:

(4)$\begin{array}{c}\frac{dS}{dt}=A-dS-({\beta }_1+{\beta }_{0}I)SI,\hfill \\ \frac{dI}{dt}=({\beta }_1+{\beta }_{0}I)SI-(d+\mu +\delta )I,\hfill \end{array}$

where $\rho =d+\mu +\delta$. In this paper, we will analyze its dynamic properties completely.

This paper is organized as follows: in the next section, we obtain the existence of the equilibria of the system and analyze the type of equilibria. We have obtained the existence of disease-free and epidemic equilibrium. In the third part, we analyze the bifurcation of the system, including the saddle-node bifurcation, Hopf bifurcation, and Bogdanov–Takens bifurcations. In addition, the germs of the system are universally unfolded, and the global phase diagrams on each parameter interval are obtained. Finally, we conducted numerical simulations of this system and discussed related biological significance.

2. Existence and Types of Equilibria

It is easy to obtain that there is a positive invariant set of system (4): $D=\{(S,I)\in [0,+\infty )\times [0,+\infty )|S+I<\frac{A}{d}\left)\right\}$. In view of the biological significance of the model, we will discuss the dynamic properties of the model in the set in the following proofs. This means that if $S\left(0\right)>0, I\left(0\right)\ge 0$, then the solution $\left(S\right(t),I(t\left)\right)$ of model (4) is positive for all $t\ge 0$, and all non-negative solutions of model (4) are ultimately uniformly bounded in forward time.

Firstly, to find the equilibria, we let system (4) have zero rights.A disease-free equilibrium $E_0=(\frac{A}{d},0)$ always exists. Furthermore, the other equilibrium satisfies

$A-dS-\rho I=0,({\beta }_1+{\beta }_{0}I)S-\rho =0.$

By direct calculation, the coordinates of the non-negative equilibria should satisfy

$S^{*}=\frac{\rho }{{\beta }_{0}I+{\beta }_1},$

and the coordinates I are to be the positive roots of the quadratic equation

$\rho {\beta }_0{I}^2+(\rho {\beta }_1-A{\beta }_0)I+(d\rho -A{\beta }_1)=0.$

If the discriminant of the equation above is positive, the equation may have two roots

$I_1^{*}=\frac{A{\beta }_0-\rho {\beta }_1-\sqrt{\Delta }}{2\rho {\beta }_0},I_2^{*}=\frac{A{\beta }_0-\rho {\beta }_1+\sqrt{\Delta }}{2\rho {\beta }_0},$

where

$\Delta =4\rho {\beta }_0(A{\beta }_1-d\rho )+{(A{\beta }_0-\rho {\beta }_1)}^2,$

which imply that the system can have two different positive equilibria $E_1^{*}=(S_1^{*},I_1^{*})$ and $E_2^{*}=(S_2^{*},I_2^{*})$. Although $\Delta =0$, i.e., $A=\frac{\rho \left(2\sqrt{{\beta }_{0}d}-{\beta }_1\right)}{{\beta }_0}$ and ${\beta }_1<\sqrt{{\beta }_{0}d}$, these two equilibria collide into a unique positive equilibrium with multiplicity 2, $E^{*}=(\frac{\rho }{\sqrt{{\beta }_{0}d}},\frac{\sqrt{{\beta }_{0}d}-{\beta }_1}{{\beta }_0})$.

Theorem 1. 

For system (4),

• (1) 

  If $A<\frac{\rho \left(2\sqrt{{\beta }_{0}d}-{\beta }_1\right)}{{\beta }_0}$, only disease-free equilibrium $E_0$ exists.

• (2) 

  If $A=\frac{\rho \left(2\sqrt{{\beta }_{0}d}-{\beta }_1\right)}{{\beta }_0}$, and

• (i) 

  ${\beta }_1<\sqrt{{\beta }_{0}d}$, epidemic equilibrium $E^{*}$ and disease-free equilibrium $E_0$ exist.

• (ii) 

  ${\beta }_1\ge \sqrt{{\beta }_{0}d}$, only $E_0$ exists.

• (3) 

  If $A>\frac{\rho \left(2\sqrt{{\beta }_{0}d}-{\beta }_1\right)}{{\beta }_0}$, and

• (i) 

  ${\beta }_0\le \frac{A^{2}d}{{\rho }^2},{\beta }_1<\frac{d\rho }{A}$ or ${\beta }_1<\frac{A{\beta }_0}{\rho },{\beta }_0>\frac{A^{2}d}{{\rho }^2}$, system (4) has two epidemic equilibria $E_1^{*}$, $E_2^{*}$ and and disease-free equilibrium $E_0$.

• (ii) 

  ${\beta }_1>\frac{d\rho }{A},$ there is an epidemic equilibrium $E_2^{*}$ and a disease-free equilibrium $E_0$.

• (iii) 

  $\frac{A{\beta }_0}{\rho }<{\beta }_1<\frac{d\rho }{A},{\beta }_0>\frac{A^{2}d}{{\rho }^2},$ only epidemic equilibrium $E_0$ exists.

Types of Equilibria

In this section, and throughout the subsequent discussions, we will employ the Central Limit Theorem and the Reduction Principle, presented in lemma form. For further details, please refer to the citation [21].

Consider a continuous-time dynamical system defined by

(5)$\frac{dx}{dt}=f\left(x\right),x\in R^n,$

where $f\in C^{\infty }$ and $f\left(0\right)=0$. ${\lambda }_1,{\lambda }_2....{\lambda }_n$ are the eigenvalues of the Jacobian matrix A at the equilibrium point $x_0=0$. If the eigenvalues with zero real part, the equilibrium is not hyperbolic. Suppose there are $n_{+}$ eigenvalues with Re $\lambda >0$, $n_0$ eigenvalues with Re $\lambda =0$, and $n_{-}$ eigenvalues with Re $\lambda <0$. $T^c$ is the linear eigenspace of A corresponding to the collection of the $n_0$ eigenvalues on the imaginary axis, and ${\varphi }^t$ denotes the flow associated with Equation (5). Given the assumptions outlined above, the following lemma is valid.

Lemma 1. 

There is a locally defined smooth $n_0$ dimensional invariant manifold $W_{loc}^c\left(0\right)$ of (5) that is tangent to $T^c$ at $x=0$. Moreover, there is a neighborhood U of $x_0=0,$ such that if ${\varphi }^{t}x\in U$ for all $t\ge 0(t\le 0),$ then ${\varphi }^{t}x\to W_{loc}^c\left(0\right)$ for $t\to +\infty (t\to -\infty )$.

$W_{loc}^c\left(0\right)$ is called the center manifold. Applying the lemma above, we can derive the following conclusion.

Theorem 2. 

For system (4), the disease-free equilibrium $E_0(\frac{A}{d},0)$

• (1) 

  if ${\beta }_1>\frac{d\rho }{A}$ is a saddle;

• (2) 

  if ${\beta }_1<\frac{d\rho }{A}$ is a stable node;

• (3) 

  if ${\beta }_1=\frac{d\rho }{A}$ is a saddle-node. Moreover,

• (i) 

  the codimension is 1 if ${\beta }_0=\frac{d{\rho }^2}{A^2}$,

• (i) 

  the codimension is 2 if ${\beta }_0=\frac{d{\rho }^2}{A^2}$.

Proof. 

With the change of variables $(S,I)\to (x,y)$ by $x=S-A/d,y=I$. The eigenvalues of Jacobian at $E_0$ are $-d$, $A{\beta }_1/d-\rho$. Hence, these are the first two conclusions of Theorem 2. If ${\beta }_1=\frac{d\rho }{A}$, the system (4) can be rewritten as

(6)$\begin{array}{c}\frac{dx}{dt}=-dx+\frac{\rho (\rho -d)}{A}xy+\frac{(d-\rho )(d{\rho }^2-{\beta }_0{A}^2)}{d^{2}A}{y}^2+(\frac{\rho }{d}-1){\beta }_{0}x{y}^2\hfill \\ +\frac{(d-\rho )\rho {\beta }_0}{d^2}{y}^3:=X(x,y),\hfill \\ \frac{dy}{dt}=\frac{d\rho }{A}xy+(\frac{A{\beta }_0}{d}-\frac{{\rho }^2}{A})y^2+{\beta }_{0}x{y}^2-\frac{\rho {\beta }_0}{d}{y}^3:=Y(x,y).\hfill \end{array}$

By the Lemma 1, there is a unique function $x=x\left(y\right)$, such that $x\left(0\right)=0$ and $X(x,x(y\left)\right)=0$. In fact, the solution can be given as

$x\left(y\right)=\frac{(d-\rho )\left(A^2{\beta }_0-d{\rho }^2\right)}{A{d}^2}{y}^2+\frac{\rho (d-\rho )\left(A^2{\beta }_0(\rho -2d)+d{\rho }^2(d-\rho )\right)}{A^2{d}^3}{y}^3+{\mathcal{O}\left(\right|y|}^4).$

Substituting this expression into Equation (6), we can obtain

$\frac{dy}{dt}=(\frac{A{\beta }_0}{d}-\frac{{\rho }^2}{A})y^2+\frac{\rho \left(\frac{d{\rho }^2(d-\rho )}{A^2}+{\beta }_0(\rho -2d)\right)}{d^2}{y}^3+{\mathcal{O}\left(\right|y|}^4).$

From [22], $E_0$ is a saddle-node of codimension 1 if ${\beta }_0\ne \frac{d{\rho }^2}{A^2}$. If ${\beta }_0=\frac{d{\rho }^2}{A^2}$, the system has been changed into

$\frac{dy}{dt}=-\frac{{\rho }^3}{A^2}{y}^3+{\mathcal{O}\left(\right|y|}^4),$

hence, the codimension is 2. □

Remark 1. 

The existence of a saddle-node $E^{*}$ reflects the critical conditions for disease occurrence. When the infection rate ${\beta }_1$ is relatively small, the disease-free equilibrium is stable, indicating that the disease will not become prevalent. However, if ${\beta }_1$ increases, under certain parameter conditions, the disease will spiral out of control, leading to its sustained existence.

Theorem 3. 

For system (4), if $A=\frac{\rho \left(2\sqrt{{\beta }_{0}d}-{\beta }_1\right)}{{\beta }_0}$, ${\beta }_1<\sqrt{{\beta }_{0}d}$, the unique epidemic equilibrium $E^{*}$ is a saddle-node.

Proof. 

At $E^{*}$, by $X=S-S^{*},Y=I-I^{*}$, we translate this equilibrium to zero and linearize the system (4), and obtain the following system:

$\begin{array}{cc}\hfill \frac{dx}{dt}& =d\left(\frac{{\beta }_1}{\sqrt{{\beta }_{0}d}}-2\right)x+\rho \left(\frac{{\beta }_1}{\sqrt{{\beta }_{0}d}}-2\right)y+({\beta }_1-2\sqrt{{\beta }_{0}d})xy\hfill \\ \hfill & -\frac{\rho \sqrt{{\beta }_{0}d}}{d}{y}^2+{O\left(\right|x,y|}^3),\hfill \\ \hfill \frac{dy}{dt}& =\left(d-\frac{{\beta }_{1}d}{\sqrt{{\beta }_{0}d}}\right)x+\left(\rho -\frac{{\beta }_1\rho }{\sqrt{{\beta }_{0}d}}\right)y+\left(2\sqrt{{\beta }_{0}d}-{\beta }_1\right)xy\hfill \\ \hfill & +\frac{\rho \sqrt{{\beta }_{0}d}}{d}{y}^2+{O\left(\right|x,y|}^3).\hfill \end{array}$

Next, we make the transformation $x_1=d\left(\frac{{\beta }_1}{\sqrt{{\beta }_{0}d}}-2\right)+\rho \left(\frac{{\beta }_1}{\sqrt{{\beta }_{0}d}}-2\right)y,y_1=y;$ then, the system can be changed into

(7)$\begin{array}{cc}\hfill \frac{d{x}_1}{dt}& =a_{1,0}{x}_1+a_{1,1}{x}_1{y}_1+a_{0,2}{y}_1^2+O\left(\right|x_1,y_1{|}^3),\hfill \\ \hfill \frac{d{y}_1}{dt}& =b_{1,1}{x}_1{y}_1+b_{0,2}{y}_1^2+O\left(\right|x_1,y_1{|}^3),\hfill \end{array}$

where

$\begin{array}{cc}\hfill a_{1,0}& =\frac{{\beta }_1(d-\rho )}{\sqrt{{\beta }_{0}d}}-2d+\rho ,a_{1,1}=\frac{\left({\beta }_1-2\sqrt{{\beta }_{0}d}\right)(d-\rho )}{d},\hfill \\ \hfill a_{0,2}& =-\frac{\rho \left(-{\beta }_1^3+4{\left({\beta }_{0}d\right)}^{3/2}+5{\beta }_1^2\sqrt{{\beta }_{0}d}-8{\beta }_0{\beta }_{1}d\right)(d-\rho )}{d\left(2{\beta }_{0}d-{\beta }_1\sqrt{{\beta }_{0}d}\right)},\hfill \\ \hfill b_{1,1}& =\sqrt{{\beta }_{0}d},b_{0,2}=\rho \left(\sqrt{{\beta }_{0}d}-{\beta }_1\right).\hfill \end{array}$

If ${\beta }_1\ne \frac{\sqrt{{\beta }_{0}d}(2d-\rho )}{d-\rho },$ then $a_{1,0}\ne 0.$ Applying the center manifold theorem to system Lemma 1, the system has a center manifold

$x_1=-\frac{a_{0,2}}{a_{1,0}}{y}_1^2+\frac{a_{0,2}{a}_{1,1}-a_{0,3}{a}_{1,0}}{a_{1,0}^2}{y}_1^3+O\left(y_1^4\right).$

We can obtain the equation on the center manifold,

$\frac{dy}{dt}=\rho \left(\sqrt{{\beta }_{0}d}-{\beta }_1\right)y^2+{\mathcal{O}\left(\right|y|}^3).$

Hence, if ${\beta }_1<\sqrt{{\beta }_{0}d}$, the equilibrium is saddle-node with codimension 1. □

As for $E_1^{*}$ and $E_2^{*}$, we have the following conclusion:

Theorem 4. 

For system (4), if simple endemic equilibrium $E_1^{*}$ exists, it is a hyperbolic saddle. If $E_2^{*}$ exists, it is an anti-saddle.

Proof. 

Denoting any endemic equilibrium as $\tilde{E}(\tilde{S},\tilde{I})$, and noting $\tilde{S}({\beta }_0\tilde{I}+{\beta }_1)=\rho$, then the Jacobian matrix of (4) is given by

$\left(\begin{array}{cc}-d-\tilde{I}({\beta }_1+{\beta }_0\tilde{I})& -\tilde{S}(2\tilde{I}{\beta }_0+{\beta }_1)\\ \tilde{I}({\beta }_1+{\beta }_0\tilde{I})& {\beta }_0\tilde{S}\tilde{I}\end{array}\right),$

the eigenvalues of this matrix are roots of ${\lambda }^2+C_1\lambda +C_0=0$, where

$C_1=d+\tilde{I}({\beta }_1+{\beta }_0\tilde{I}+{\beta }_0\tilde{S}),C_0=\tilde{S}\tilde{I}\left({\tilde{I}}^2{\beta }_0^2+{\beta }_1^2-{\beta }_0(d-2\tilde{I}{\beta }_1)\right).$

Notice $\tilde{S}=\frac{\rho }{{\beta }_1+{\beta }_0\tilde{I}}$, hence,

$C_0\left(\tilde{I}\right)=\rho \tilde{I}\left({\beta }_1+{\beta }_0\tilde{I}-\frac{d{\beta }_0}{{\beta }_1+{\beta }_0\tilde{I}}\right).$

On the other hand, discriminant $\Delta >0$ is equivalent to ${\beta }_1>2\sqrt{d{\beta }_0}-\frac{A{\beta }_0}{\rho }$. Through this inequality, we can easily prove

$\begin{array}{c}\hfill {\beta }_1+{\beta }_0{\tilde{I}}_1^{*}-\frac{d{\beta }_0}{{\beta }_1+{\beta }_0{\tilde{I}}_1^{*}}<0,\\ \hfill {\beta }_1+{\beta }_0{\tilde{I}}_2^{*}-\frac{d{\beta }_0}{{\beta }_1+{\beta }_0{\tilde{I}}_2^{*}}>0,\end{array}$

i.e., $C_0\left({\tilde{I}}_1^{*}\right)<0$, $C_0\left({\tilde{I}}_2^{*}\right)>0$. That is to say, $E_1^{*}$ is a hyperbolic saddle and $E_2^{*}$ is an anti-saddle. Especially, if $C_1=0$, $E_2^{*}$ is a weak focus or center; if $C_1>0$, $E_2^{*}$ it is attracting a node or focus. If $C_1<0$, $E_2^{*}$, it is repelling a node or focus. □

If we denote $S_0^{*}=\frac{\rho }{\sqrt{{\beta }_{0}d}},I_0^{*}=\frac{d^2}{\sqrt{{\beta }_{0}d}(\rho -d)}$, then the system has a nilpotent singularity. For further information on this topic, readers are referred to the literature [23]. For system (4), we have the following:

Theorem 5. 

Suppose $A=\frac{d{\rho }^2}{\sqrt{{\beta }_{0}d}(\rho -d)}$, ${\beta }_1=\frac{\sqrt{{\beta }_{0}d}(\rho -2d)}{\rho -d}$, $\rho >2d$ and $\rho \ne 3d$ system (4) has a unique positive equilibrium $E_0^{*}(S_0^{*},I_0^{*})$, which is a codimension 2 cusp-type nilpotent singularity, and the system at $E_0^{*}$ is topologically equivalent to system

(8) $\begin{array}{cc}\hfill \frac{dx}{d\tau }& =y,\hfill \\ \hfill \frac{dy}{d\tau }& =x^2\pm xy+{O\left(\right|x,y|}^3),\hfill \end{array}$

where if $2d<\rho <3d$, then the $“\pm ”$ above takes a positive sign, but if $\rho >3d$, then the $“\pm ”$ above takes a negative sign. Hence, $E_0^{*}$ is a codimension 2 cusp-type nilpotent singularity.

Proof. 

If $A=\frac{d{\rho }^2}{\sqrt{{\beta }_{0}d}(\rho -d)}$, ${\beta }_1=\frac{\sqrt{{\beta }_{0}d}(2d-\rho )}{d-\rho }$, $\rho >2d,$ the existence of $E_0^{*}$ can be proved easily, so we omit here, and only prove this equilibrium is a degenerate nilpotent singularity, and the system at $E_0^{*}$ is topologically equivalent to system (8). After translating the equilibrium to origin and performing Taylor expansion, we obtain

(9)$\begin{array}{cc}\hfill \frac{d{x}_1}{dt}& =\sum _{i+j=1}^2{\tilde{a}}_{i,j}{x}_1^i{y}_1^j+\mathcal{O}\left(\right|x_1,y_1{|}^3),\hfill \\ \hfill \frac{d{y}_1}{dt}& =\sum _{i+j=1}^2{\tilde{b}}_{i,j}{x}_1^i{y}_1^j+\mathcal{O}\left(\right|x_1,y_1{|}^3),\hfill \end{array}$

where

$\begin{array}{c}{\tilde{a}}_{1,0}=\frac{d\rho }{d-\rho },{\tilde{a}}_{0,1}=\frac{{\rho }^2}{d-\rho },{\tilde{a}}_{1,1}=\frac{\rho \sqrt{{\beta }_{0}d}}{d-\rho },{\tilde{a}}_{0,2}=-\frac{\rho \sqrt{{\beta }_{0}d}}{d},\hfill \\ {\tilde{b}}_{1,0}=-\frac{d^2}{d-\rho },{\tilde{b}}_{0,1}=\frac{d\rho }{\rho -d},{\tilde{b}}_{1,1}=-\frac{\rho \sqrt{{\beta }_{0}d}}{d-\rho },{\tilde{b}}_{0,2}=\frac{\rho \sqrt{{\beta }_{0}d}}{d}.\hfill \end{array}$

Let $x_2=x_1$ and $y_2=\frac{d\rho }{d-\rho }{x}_1++\frac{{\rho }^2}{d-\rho }{y}_1,$ then system (9) can be changed into

(10)$\begin{array}{cc}\hfill \frac{d{x}_2}{dt}& =y_2+\sum _{i+j=2}{\tilde{c}}_{i,j}{x}_2^i{y}_2^j+\mathcal{O}\left(\right|x_2,y_2{|}^3),\hfill \\ \hfill \frac{d{y}_2}{dt}& =\sum _{i+j=2}{\tilde{d}}_{i,j}{x}_2^i{y}_2^j+\mathcal{O}\left(\right|x_2,y_2{|}^3),\hfill \end{array}$

where

$\begin{array}{c}{\tilde{c}}_{2,0}=-\frac{d^2\sqrt{{\beta }_{0}d}}{d\rho -{\rho }^2},{\tilde{a}}_{1,1}=\frac{\sqrt{{\beta }_{0}d}(2d-\rho )}{{\rho }^2},{\tilde{c}}_{0,2}=-\frac{\sqrt{{\beta }_{0}d}{(d-\rho )}^2}{d{\rho }^3},\hfill \\ {\tilde{d}}_{2,0}=-\frac{d^2\sqrt{{\beta }_{0}d}}{d-\rho },{\tilde{d}}_{1,1}=\frac{\sqrt{{\beta }_{0}d}(2d-\rho )}{\rho },{\tilde{d}}_{0,2}=-\frac{\sqrt{{\beta }_{0}d}{(d-\rho )}^2}{d{\rho }^2}.\hfill \end{array}$

Make near-identity transformation

$x_3=x_2,y_3=y_2+\sum _{i+j=2}{\tilde{c}}_{i,j}{x}_2^i{y}_2^j+\mathcal{O}\left(\right|x_2,y_2{|}^3),$

in system (10) to eliminate the non-resonant terms. Then, we obtain

(11)$\begin{array}{cc}\hfill \frac{d{x}_3}{dt}& =y_3,\hfill \\ \hfill \frac{d{y}_3}{dt}& =\sum _{i+j=1}^2{\tilde{e}}_{i,j}{x}_3^i{y}_3^j+O\left(\right|x_3,y_3{|}^3),\hfill \end{array}$

where

$\begin{array}{c}{\tilde{e}}_{2,0}=\frac{d^2\sqrt{{\beta }_{0}d}}{\rho -d}>0,\hfill \\ {\tilde{e}}_{1,1}=\frac{\sqrt{{\beta }_{0}d}(\rho -3d)}{d-\rho }\ne 0,\hfill \\ {\tilde{e}}_{0,2}=\frac{\sqrt{{\beta }_{0}d}\left(d^2+d\rho -{\rho }^2\right)}{d{\rho }^2}.\hfill \end{array}$

At the last step, by scale transformation

$x_3=\frac{{\tilde{e}}_{2,0}}{{\tilde{e}}_{1,1}^2}{x}_4,y_3=\frac{{\tilde{e}}_{2,0}^2}{{\tilde{e}}_{1,1}^3}{y}_4,t=\frac{{\tilde{e}}_{1,1}}{{\tilde{e}}_{2,0}}\tau ,$

the system is $C^{\infty }$ equivalent to system (8). Thus, we prove that there exist smooth coordinate changes which take system (4) into the system above. From [22,24,25,26], the germ in the system will produce a nilpotent singularity bifurcation of codimension 2. This completes the proof of Theorem 5. □

Remark 2. 

Theorem 5 discusses the case when $\rho >2d$ and $\rho \ne 3d$, but if $\rho =3d$, then the situation is very complicated; the system may have a degenerate shoot of higher codimension, which will be discussed in the following works.

Remark 3. 

The nilpotent singularity serves as the organizational core of the model, and its existence reflects the complexity of the model.

3. Bifurcation

To investigate the dynamics implied in system (4), we first consider the bifurcation direction and codimension of the Hopf bifurcation.

3.1. Hopf Bifurcation

Recall $S_1^{*}=\frac{A-\rho I_1^{*}}{d}$ and $A\left(E_1^{*}\right)=-d-I_2^{*}({\beta }_1+{\beta }_0{I}_1^{*}+{\beta }_0{S}_1^{*})=-\frac{A}{S_1^{*}}+\rho -{\beta }_1{S}_1^{*},$ and we have the following:

Theorem 6. 

A generic Hopf bifurcation could occur if $A\left(E_1^{*}\left({\beta }_0\right)\right)=0$.

Proof. 

We only need to verify the transversal condition. The real part of the eigenvalue of system (4) is $-A/2$. In the following, we will prove $\partial A/\partial {\beta }_0\ne 0.$

$\begin{array}{cc}\hfill \frac{\partial A\left({\beta }_0\right)}{\partial {\beta }_0}& =-\frac{\partial S_1^{*}}{\partial {\beta }_0}({\beta }_1+\frac{A}{S_1^{*2}})\hfill \\ \hfill & =\frac{\rho ({\beta }_1+\frac{A}{S_1^{*2}})}{2d{\beta }_0^2\sqrt{A^2{\beta }_0^2+{\rho }^2{\beta }_1^2+2\rho {\beta }_0(A{\beta }_1-2d\rho )}}f\left({\beta }_0\right),\hfill \end{array}$

here $f\left({\beta }_0\right)={\beta }_0(A{\beta }_1-2d\rho )+{\beta }_1\left(\rho {\beta }_1-\sqrt{A^2{\beta }_0^2+{\rho }^2{\beta }_1^2+2\rho {\beta }_0(A{\beta }_1-2d\rho )}\right)$, and $f\left({\beta }_0\right)=0$ holds if and only if ${\beta }_0=0$, contradicting the range of ${\beta }_0$. Thus, Theorem 6 holds. □

3.2. Cusp-Type Nilpotent Bifurcation with Codimension 2

If $A^{*}=\frac{d{\rho }^2}{\sqrt{{\beta }_{0}d}(\rho -d)}, {\beta }_1^{*}=\frac{\sqrt{{\beta }_{0}d}(2d-\rho )}{d-\rho },$ $\rho >2d$ and $\rho \ne 3d,$ $E_0^{*}$ is a cuspital point from Theorem 5. In this section, we will unfold the singularity. First, we chose ${\beta }_1$ and d as bifurcation parameters and let $(A,{\beta }_1)=(A^{*},{\beta }_1^{*})$, satisfying the conditions in Theorem 5. To prove that these parameters can unfold the singularity, we perturb the parameters and eliminate ${\beta }_1$ and d. Let $A=A^{*}+{ϵ}_1,{\beta }_1={\beta }_1^{*}+{ϵ}_2$. Next, we will study the bifurcation of the perturbed system

(12)$\begin{array}{cc}\hfill \frac{dS}{dt}& =A+{ϵ}_1-dS-({\beta }_1^{*}+{ϵ}_2+{\beta }_{0}I)SI,\hfill \\ \hfill \frac{dI}{dt}& =({\beta }_1^{*}+{ϵ}_2+{\beta }_{0}I)SI-\rho I.\hfill \end{array}$

By

$x_1=S-S_0^{*}, y_1=I-I_0^{*},$

we have

(13)$\begin{array}{cc}\hfill \frac{d{x}_1}{dt}& =\sum _{i+j=0}^2{\widehat{a}}_{i,j}{x}_1^i{y}_1^j+\mathcal{O}\left(\right|x_1,y_1{|}^3),\hfill \\ \hfill \frac{d{y}_1}{dt}& =\sum _{i+j=0}^2{\widehat{b}}_{i,j}{x}_1^i{y}_1^j+\mathcal{O}\left(\right|x_1,y_1{|}^3),\hfill \end{array}$

where

$\begin{array}{c}{\widehat{a}}_{0,0}={ϵ}_1+\frac{d\rho {ϵ}_2}{{\beta }_0(d-\rho )}+O(|{ϵ}_1,{ϵ}_2{|}^2),{\widehat{a}}_{1,0}=\frac{d\rho }{d-\rho }+\frac{d^2{ϵ}_2}{\sqrt{{\beta }_{0}d}(d-\rho )}+O\left(\right|{ϵ}_1,{ϵ}_2{|}^2),\hfill \\ {\widehat{a}}_{0,1}=\frac{{\rho }^2}{d-\rho }-\frac{\rho {ϵ}_2}{\sqrt{{\beta }_{0}d}}+O\left(\right|{ϵ}_1,{ϵ}_2{|}^2),{\widehat{a}}_{2,0}=0,\hfill \\ {\widehat{a}}_{1,1}=\frac{\rho \sqrt{{\beta }_{0}d}}{d-\rho }-{ϵ}_2+O\left(\right|{ϵ}_1,{ϵ}_2{|}^2),{\widehat{a}}_{0,2}=-\frac{\rho \sqrt{{\beta }_{0}d}}{d},\hfill \\ {\widehat{b}}_{0,0}=\frac{d\rho {ϵ}_2}{{\beta }_0(\rho -d)}+O\left(\right|{ϵ}_1,{ϵ}_2{|}^2),\hfill \\ b_{1,0}=\frac{d^2}{\rho -d}+\frac{d^2{ϵ}_2}{\sqrt{{\beta }_{0}d}(\rho -d)}+O\left(\right|{ϵ}_1,{ϵ}_2{|}^2),\hfill \\ {\widehat{b}}_{0,1}=-\frac{d\rho }{d-\rho }+\frac{\rho {ϵ}_2}{\sqrt{{\beta }_{0}d}}+O\left(\right|{ϵ}_1,{ϵ}_2{|}^2),\hfill \\ {\widehat{b}}_{1,1}=-\frac{\rho \sqrt{{\beta }_{0}d}}{d-\rho }+{ϵ}_2+O\left(\right|{ϵ}_1,{ϵ}_2{|}^2),{\widehat{b}}_{0,2}=\frac{\rho \sqrt{{\beta }_{0}d}}{d}.\hfill \end{array}$

By the transformation

$x_1=x_2-\frac{{\widehat{a}}_{0,0}}{{\widehat{a}}_{1,0}},y_1=y_2.$

System (13) is $C^{\infty }$ equivalent to

(14)$\begin{array}{cc}\hfill \frac{d{x}_2}{dt}& =\sum _{i+j=1}^2{\widehat{c}}_{i,j}{x}_2^i{y}_2^j+O\left(\right|x_2,y_2{|}^3),\hfill \\ \hfill \frac{d{y}_2}{dt}& =\sum _{i+j=0}^2{\widehat{d}}_{i,j}{x}_2^i{y}_2^j+O\left(\right|x_2,y_2{|}^3),\hfill \end{array}$

where

$\begin{array}{cc}\hfill {\widehat{c}}_{1,0}& ={\widehat{a}}_{1,0},{\widehat{c}}_{0,1}={\widehat{a}}_{0,1}-\frac{{\widehat{a}}_{0,0}{\widehat{a}}_{1,1}}{{\widehat{a}}_{1,0}},{\widehat{c}}_{1,1}={\widehat{a}}_{1,1},{\widehat{c}}_{0,2}={\widehat{a}}_{0,2},\hfill \\ \hfill {\widehat{d}}_{0,0}& ={\widehat{b}}_{0,0}-\frac{{\widehat{a}}_{0,0}{\widehat{b}}_{1,0}}{{\widehat{a}}_{1,0}},{\widehat{d}}_{1,0}={\widehat{b}}_{1,0},{\widehat{d}}_{0,1}={\widehat{b}}_{0,1}-\frac{{\widehat{a}}_{0,0}{\widehat{b}}_{1,1}}{{\widehat{a}}_{1,0}},{\widehat{d}}_{1,1}={\widehat{b}}_{1,1},{\widehat{d}}_{0,2}={\widehat{b}}_{0,2}.\hfill \end{array}$

Next, we perform the following transmission,

$x_3=x_2,y_3={\widehat{c}}_{1,0}{x}_2+{\widehat{c}}_{0,1}{y}_2,$

then the system can be changed into

(15)$\begin{array}{cc}\hfill \frac{d{x}_3}{dt}& =y_3+\sum _{i+j=2}{\widehat{e}}_{i,j}{x}_3^i{y}_3^j+O\left(\right|x_3,y_3{|}^3),\hfill \\ \hfill \frac{d{x}_3}{dt}& =\sum _{i+j=0}^2{\widehat{f}}_{i,j}{x}_3^i{y}_3^j+O\left(\right|x_3,y_3{|}^3),\hfill \end{array}$

where

$\begin{array}{cc}\hfill {\widehat{e}}_{2,0}& =\frac{{\widehat{c}}_{0,2}{\widehat{c}}_{1,0}^2}{{\widehat{c}}_{0,1}^2}-\frac{{\widehat{c}}_{1,0}{\widehat{c}}_{1,1}}{{\widehat{c}}_{0,1}},{\widehat{e}}_{1,1}=\frac{{\widehat{c}}_{1,1}}{{\widehat{c}}_{0,1}}-\frac{2{\widehat{c}}_{0,2}{\widehat{c}}_{1,0}}{{\widehat{c}}_{0,1}^2},{\widehat{e}}_{0,2}=\frac{{\widehat{c}}_{0,2}}{{\widehat{c}}_{0,1}^2},\hfill \\ \hfill {\widehat{f}}_{0,0}& ={\widehat{c}}_{0,1}{\widehat{d}}_{0,0},{\widehat{f}}_{1,0}={\widehat{c}}_{0,1}{\widehat{d}}_{1,0}-{\widehat{c}}_{1,0}{\widehat{d}}_{0,1},{\widehat{f}}_{0,1}={\widehat{c}}_{1,0}+{\widehat{d}}_{0,1},\hfill \\ \hfill {\widehat{f}}_{2,0}& =\frac{{\widehat{c}}_{0,2}{\widehat{c}}_{1,0}^3-{\widehat{c}}_{0,1}{\widehat{c}}_{1,0}\left({\widehat{c}}_{1,0}\left({\widehat{c}}_{1,1}-{\widehat{d}}_{0,2}\right)+{\widehat{c}}_{0,1}{\widehat{d}}_{1,1}\right)}{{\widehat{c}}_{0,1}^2},\hfill \\ \hfill {\widehat{f}}_{1,1}& =\frac{{\widehat{c}}_{0,1}\left({\widehat{c}}_{1,0}\left({\widehat{c}}_{1,1}-2{\widehat{d}}_{0,2}\right)+{\widehat{c}}_{0,1}{\widehat{d}}_{1,1}\right)-2{\widehat{c}}_{0,2}{\widehat{c}}_{1,0}^2}{{\widehat{c}}_{0,1}^2},{\widehat{f}}_{0,2}=\frac{{\widehat{c}}_{0,1}{\widehat{d}}_{0,2}+{\widehat{c}}_{0,2}{\widehat{c}}_{1,0}}{{\widehat{c}}_{0,1}^2}.\hfill \end{array}$

Next, letting

$x_4=x_3,y_4=y_3+\sum _{i+j=2}{\widehat{e}}_{i,j}{x}_3^i{y}_3^j+O\left(\right|x_3,y_3{|}^3),$

then system (15) becomes

(16)$\begin{array}{cc}\hfill \frac{d{x}_4}{dt}& =y_4,\hfill \\ \hfill \frac{d{y}_4}{dt}& =\sum _{i+j=0}^2{\widehat{g}}_{i,j}{x}_4^i{y}_4^j+O\left(\right|x_4,y_4{|}^3),\hfill \end{array}$

where

$\begin{array}{cc}\hfill {\widehat{g}}_{0,0}& ={\widehat{f}}_{0,0},{\widehat{g}}_{1,0}={\widehat{e}}_{1,1}{\widehat{f}}_{0,0}+{\widehat{f}}_{1,0},{\widehat{g}}_{0,1}=2{\widehat{e}}_{0,2}{\widehat{f}}_{0,0}+{\widehat{f}}_{0,1},\hfill \\ \hfill {\widehat{g}}_{2,0}& =-2{\widehat{e}}_{0,2}{\widehat{e}}_{2,0}{\widehat{f}}_{0,0}-{\widehat{e}}_{2,0}{\widehat{f}}_{0,1}+{\widehat{e}}_{1,1}{\widehat{f}}_{1,0}+{\widehat{f}}_{2,0},\hfill \\ \hfill {\widehat{g}}_{1,1}& ={\widehat{e}}_{0,2}\left(2{\widehat{f}}_{1,0}-2{\widehat{e}}_{1,1}{\widehat{f}}_{0,0}\right)+2{\widehat{e}}_{2,0}+{\widehat{f}}_{1,1},\hfill \\ \hfill {\widehat{g}}_{0,2}& =-2{\widehat{e}}_{0,2}^2{\widehat{f}}_{0,0}+{\widehat{e}}_{0,2}{\widehat{f}}_{0,1}+{\widehat{e}}_{1,1}+{\widehat{f}}_{0,2}.\hfill \end{array}$

Make the following transformation

$x_5=x_4,y_5=\frac{y_4}{\sqrt{{\widehat{g}}_{2,0}}},\tau ={\int }_0^t\sqrt{{\widehat{g}}_{2,0}}ds,$

then system (16) can be rewritten as

(17)$\begin{array}{cc}\hfill \frac{d{x}_5}{dt}& =y_5,\hfill \\ \hfill \frac{d{y}_5}{dt}& =h_{0,0}+h_{1,0}{x}_5+x_5^2+y_5\left(h_{0,1}+h_{1,1}{x}_5+\mathcal{O}\left(x_5^2\right)\right)+y_5^{2}O\left(\right|x_5,y_5\left|\right),\hfill \end{array}$

where

$\begin{array}{c}\hfill h_{0,0}=\frac{{\widehat{g}}_{0,0}}{{\widehat{g}}_{2,0}},h_{1,0}=\frac{{\widehat{g}}_{1,0}}{{\widehat{g}}_{2,0}},h_{0,1}=\frac{{\widehat{g}}_{0,1}}{{\widehat{g}}_{2,0}},h_{1,1}=\frac{{\widehat{g}}_{1,1}}{{\widehat{g}}_{2,0}}.\end{array}$

At the last step, let $u=x_5+\frac{h_{1,0}}{2},v=h_{1,1}{y}_5,t=h_{1,1}\tau ,$ and system (17) can be changed into

(18)$\begin{array}{cc}\hfill \frac{du}{d\tau }& =v,\hfill \\ \hfill \frac{dv}{d\tau }& ={\nu }_1({ϵ}_1,{ϵ}_2)+{\nu }_2({ϵ}_1,{ϵ}_2)v+x^2+xy+O\left({|u,v|}^3,{ϵ}_1,{ϵ}_2\right),\hfill \end{array}$

where

$\begin{array}{cc}\hfill {\nu }_1& =h_{0,0}-\frac{h_{1,0}^2}{4}\hfill \\ \hfill & =-\frac{\rho }{d\sqrt{{\beta }_{0}d}}{ϵ}_1-\frac{{\rho }^2\sqrt{{\beta }_{0}d}}{{\beta }_0^2{d}^2}{ϵ}_2+O\left(\right|{ϵ}_1,{ϵ}_2{|}^2),\hfill \\ \hfill {\nu }_2& =\frac{h_{0,1}}{h_{1,1}}-\frac{h_{1,0}}{2}\hfill \\ \hfill & =\frac{4d^4-2d^3\rho -d{\rho }^3+{\rho }^4}{2d^2{\rho }^2(3d-\rho )}{ϵ}_1+\frac{4d^4-4d^3\rho +4d^2{\rho }^2-3d{\rho }^3+{\rho }^4}{2{\beta }_0{d}^2\rho (3d-\rho )}{ϵ}_2+O\left(\right|{ϵ}_1,{ϵ}_2{|}^2).\hfill \end{array}$

and if $2d<\rho <3d,$

(19)$\frac{D({\nu }_1,{\nu }_2)}{D({ϵ}_1,{ϵ}_2)}{|}_{({ϵ}_1,{ϵ}_2)=(0,0)}=\frac{{\beta }_0\rho {(d-\rho )}^2}{{\left({\beta }_{0}d\right)}^{5/2}(3d-\rho )}\ne 0.$

Hence, we have the following conclusion.

Theorem 7. 

For $ϵ=({ϵ}_1,{ϵ}_2)$, which is sufficiently small, system (12) is a universal unfolding of the cusp singularity of codimension 2.

Remark 4. 

The existence of unstable periodic solutions provides a threshold condition. If the initial values lie within the limit cycle, the number of infected individuals stabilizes at a constant after a period of oscillation. However, if the initial values lie outside the limit cycle, the number of infected individuals tends to die out after oscillation.

3.3. Bifurcation Diagram and the Phase Diagrams

In this section, we will conduct numerical simulations based on the previous analysis. Firstly, we employ the MatCont7P4 program to generate bifurcation diagrams. Taking $A=0.864$, $d=0.2$, and $\mu =0.2$, $\delta =0.15$, then $\rho =d+\mu +\delta =0.55<3d.$ We keep these parameters constant and only vary the values of ${\beta }_0$ and ${\beta }_1$, thus obtaining the bifurcation diagram Figure 1a. In this figure, the blue curve represents the saddle-node bifurcation curve, whereas the red curve represents the Hopf bifurcation curve. BT denotes the Bogdanov–Takens point with codimension 2. These three curves divide the parameter plane ${\beta }_1-{\beta }_0$ into three parts. In Figure 2, we use the Phase Plane software to provide the global phase portraits in these three parts. Combining Figure 1a and Figure 2, we can clearly observe the variation process of the phase portraits with the changes in parameters ${\beta }_0$ and ${\beta }_1$.

Next, we consider the case where $\rho >3d$, taking $A=0.98$, $d=0.2$, $\mu =0.3$, and $\delta =0.2$. Then, $\rho =0.7>3d$. In this case, the Hopf bifurcation curve and the saddle-node bifurcation curve divide the entire region into five parts (see Figure 3). Unlike Figure 1a, the Hopf bifurcation curve exhibits self-intersections. This complicates the dynamical properties of the system. In fact, through the phase portraits, we discover that multiple limit cycles may emerge outside the positive singularity point $E_2^{*}$ in regions $III$ and V, and the stability of these two limit cycles may be opposite. To thoroughly understand the relative positions and stability of these two limit cycles, it may be necessary to compute higher-order focal quantities, which is inherently very challenging. We will discuss this in future work.

Remark 5. 

In Figure 3c,e, although there are two concentric limit cycles around $E_2^{*}$, the stability of the larger and smaller cycles may be opposite. We did not provide a formal proof for this in the paper, but we observed this phenomenon through numerical simulations.

4. Discussion

This paper presents an analysis of the dynamics of an epidemic model featuring a nonlinear incidence rate, originally proposed by Van den Driessche and Watmough [12], and subsequently studied by Jin [7]. Jin’s findings indicate that in the SIRS model, after neglecting the disease-induced mortality rate, the three-dimensional model can be reduced to a planar polynomial system. Utilizing qualitative theory, it was discovered that the system may harbor a stable cusp-type singularity, leading to the emergence of stable Bogdanov–Takens bifurcations accompanied by stable limit cycles. Here, we extend the same incidence rate to the SIR model. Consequently, the stable nodes in the system become unstable, resulting in the emergence of unstable limit cycles through bifurcation. Dynamically, if the limit cycle is stable, trajectories near it tend towards the cycle. Conversely, if the limit cycle is unstable, the system exhibits a peculiar bistable structure, with the limit cycle itself acting as a boundary. Trajectories starting within the limit cycle converge towards a positive equilibrium point, whereas those outside tend towards a boundary equilibrium point, indicating disease extinction.

In practical applications, our results underscore that the most effective disease control strategy involves regulating ${\beta }_1$, representing the probability of contact between susceptible and infected individuals. This highlights the paramount importance of controlling the effective contact rates between these two groups when managing infectious diseases. Additionally, we observe that elevated multiple contact rates lead to unstable periodic solutions in the system. Consequently, the system’s dynamics become less governed by the basic reproduction number $R_0$ and entail a more intricate mechanism. Importantly, the disease’s progression depends not only on system parameters but also on initial conditions. Trajectories initiated within the limit cycle stabilize at a constant level of infected individuals, whereas those outside lead to disease eradication. This underscores the necessity of controlling the populations of susceptible and infected individuals.

It is worth noting that whereas this paper identifies a Bogdanov–Takens bifurcation with codimension 2, further analysis of higher codimension bifurcations and the dynamic impact of immune factors and social distancing measures remain avenues for future research in our model.

Author Contributions

Formal analysis, M.M.; investigation, M.M.; writing–original draft preparation, J.L.; writing–review and editing, M.M.; funding acquisition, J.L. and M.M. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Bifurcation diagram of system (4). “Hopf” represents the Hopf bifurcation curve, “Saddle-Node” refers to the saddle-node bifurcation curve, “DH” represents the degenerate Hopf bifurcation point, and “BT” represents the location of the Bogdanov–Takens bifurcation.

Figure 2. Global phase diagram of system (4) in Figure 1a when [Forumla omitted. See PDF.]. (a) In region I, the system has only one disease-free equilibrium [Forumla omitted. See PDF.], which is globally stable; (b) in region [Forumla omitted. See PDF.], there are two positive equilibrium points [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] in the system, and near [Forumla omitted. See PDF.] there exists an unstable limit cycle; (c) in region [Forumla omitted. See PDF.], there are two positive equilibrium points [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] in the system. [Forumla omitted. See PDF.] is a stable focus.

Figure 2. Global phase diagram of system (4) in Figure 1a when [Forumla omitted. See PDF.]. (a) In region I, the system has only one disease-free equilibrium [Forumla omitted. See PDF.], which is globally stable; (b) in region [Forumla omitted. See PDF.], there are two positive equilibrium points [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] in the system, and near [Forumla omitted. See PDF.] there exists an unstable limit cycle; (c) in region [Forumla omitted. See PDF.], there are two positive equilibrium points [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] in the system. [Forumla omitted. See PDF.] is a stable focus.

Figure 3. Global phase diagram of system (4) in Figure 1b when [Forumla omitted. See PDF.]. (a) In region I, the system has only one disease-free equilibrium E0, which is globally stable; (b) in region [Forumla omitted. See PDF.], there are two positive equilibrium points [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] in the system, and [Forumla omitted. See PDF.] is a stable focus or node, whereas [Forumla omitted. See PDF.] is a saddle; (c) in region [Forumla omitted. See PDF.] the system may exhibit two limit cycles; (d) in region [Forumla omitted. See PDF.] [Forumla omitted. See PDF.] is a stable focus or node, and [Forumla omitted. See PDF.] is a saddle; (e) in region [Forumla omitted. See PDF.] the system may exhibit two limit cycles.

References

1. Msmali, A.; Zico, M.; Mechai, I.; Ahmadini, A. Modeling and simulation: A study on predicting the outbreak of COVID-19 in Saudi Arabia. Discret. Dyn. Nat. Soc.; 2021; 2021, 5522928. [DOI: https://dx.doi.org/10.1155/2021/5522928]

2. Liu, W.M.; Levin, S.A.; Iwasa, Y. Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological models. J. Math. Biol.; 1986; 23, pp. 187-204. [DOI: https://dx.doi.org/10.1007/BF00276956] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/3958634]

3. Ruan, S.; Wang, W. Dynamical behavior of an epidemic model with a nonlinear incidence rate. J. Differ. Equ.; 2003; 188, pp. 135-163. [DOI: https://dx.doi.org/10.1016/S0022-0396(02)00089-X]

4. Tang, Y.; Huang, D.; Ruan, S.; Zhang, W. Coexistence of limit cycles and homoclinic loops in a SIRS model with a nonlinear incidence rate. SIAM J. Appl. Math.; 2008; 69, pp. 621-639. [DOI: https://dx.doi.org/10.1137/070700966]

5. Li, G.; Wang, W. Bifurcation analysis of an epidemic model with nonlinear incidence. Appl. Math. Comput.; 2009; 214, pp. 411-423. [DOI: https://dx.doi.org/10.1016/j.amc.2009.04.012]

6. Lu, M.; Huang, J.; Ruan, S.; Yu, P. Bifurcation analysis of an SIRS epidemic model with a generalized nonmonotone and saturated incidence rate. J. Differ. Equ.; 2019; 267, pp. 1859-1898. [DOI: https://dx.doi.org/10.1016/j.jde.2019.03.005] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32226129]

7. Jin, Y.; Wang, W.; Xiao, S. An SIRS model with a nonlinear incidence rate. Chaos Solitons Fractals; 2007; 34, pp. 1482-1497. [DOI: https://dx.doi.org/10.1016/j.chaos.2006.04.022]

8. Lu, M.; Xiang, C.; Huang, J. Bogdanov-Takens bifurcation in a SIRS epidemic model with a generalized nonmonotone incidence rate. Discret. Contin. Dyn. Syst.-S; 2020; 13, pp. 3125-3138. [DOI: https://dx.doi.org/10.3934/dcdss.2020115]

9. Huang, C.; Jiang, Z.; Huang, X.; Zhou, X. Bifurcation analysis of an SIS epidemic model with a generalized non-monotonic and saturated incidence rate. Int. J. Biomath.; 2024; 17, 2350033. [DOI: https://dx.doi.org/10.1142/S179352452350033X]

10. Li, C.; Li, J.; Ma, Z.; Zhu, H. Canard phenomenon for an SIS epidemic model with nonlinear incidence. J. Math. Anal. Appl.; 2014; 420, pp. 987-1004. [DOI: https://dx.doi.org/10.1016/j.jmaa.2014.06.035]

11. Li, J.; Li, S.; Wang, X. Canard, homoclinic loop, and relaxation oscillations in a Lotka–Volterra system with Allee effect in predator population. Chaos Interdiscip. J. Nonlinear Sci.; 2023; 33, 073130. [DOI: https://dx.doi.org/10.1063/5.0152946] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37459220]

12. van den Driessche, P.; Watmough, J. A simple SIS epidemic model with a backward bifurcation. J. Math. Biol.; 2000; 40, pp. 525-540. [DOI: https://dx.doi.org/10.1007/s002850000032] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10945647]

13. Alexander, M.E.; Moghadas, S.M. Bifurcation analysis of an SIRS epidemic model with generalized incidence. SIAM J. Appl. Math.; 2005; 65, pp. 1794-1816. [DOI: https://dx.doi.org/10.1137/040604947]

14. Moghadas, S.M.; Alexander, M.E. Bifurcations of an epidemic model with non-linear incidence and infection-dependent removal rate. Math. Med. Biol. J. IMA; 2006; 23, pp. 231-254. [DOI: https://dx.doi.org/10.1093/imammb/dql011] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16648145]

15. Ali, Z.; Rabiei, F.; Rashidi, M.M.; Khodadadi, T. A fractional-order mathematical model for COVID-19 outbreak with the effect of symptomatic and asymptomatic transmissions. Eur. Phys. J. Plus; 2022; 137, 395. [DOI: https://dx.doi.org/10.1140/epjp/s13360-022-02603-z] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35368740]

16. Wang, D.; Zhao, Y.; Luo, J.; Leng, H. Simplicial SIRS epidemic models with nonlinear incidence rates. Chaos Interdiscip. J. Nonlinear Sci.; 2021; 31, 053112. [DOI: https://dx.doi.org/10.1063/5.0040518] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34240944]

17. Lu, M.; Gao, D.; Huang, J.; Wang, H. Relative prevalence-based dispersal in an epidemic patch model. J. Math. Biol.; 2023; 86, 52. [DOI: https://dx.doi.org/10.1007/s00285-023-01887-8] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36877332]

18. Shan, C.; Yi, Y.; Zhu, H. Nilpotent singularities and dynamics in an SIR type of compartmental model with hospital resources. J. Differ. Equ.; 2016; 260, pp. 4339-4365. [DOI: https://dx.doi.org/10.1016/j.jde.2015.11.009]

19. Cui, W.; Zhao, Y. Saddle-node bifurcation and Bogdanov-Takens bifurcation of a SIRS epidemic model with nonlinear incidence rate. J. Differ. Equ.; 2024; 384, pp. 252-278. [DOI: https://dx.doi.org/10.1016/j.jde.2023.11.030]

20. Van den Driessche, P.; Watmough, J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci.; 2002; 180, pp. 29-48. [DOI: https://dx.doi.org/10.1016/S0025-5564(02)00108-6]

21. Kuznetsov, Y.A.; Kuznetsov, I.A.; Kuznetsov, Y. Elements of Applied Bifurcation Theory; Springer: Berlin/Heidelberg, Germany, 1998; Volume 112.

22. Zhang, Z. Qualitative Theory of Differential Equations; American Mathematical Soc.: Providence, RI, USA, 1992; Volume 101.

23. Chow, S.N.; Li, C.; Wang, D. Normal Forms and Bifurcation of Planar Vector Fields; Cambridge University Press: Cambridge, UK, 1994.

24. Dumortier, F.; Fiddelaers, P.; Li, C. Generic unfolding of the nilpotent saddle of codimension four. Global Analysis of Dynamical Systems; CRC Press: Boca Raton, FL, USA, 2001; pp. 137-172.

25. Dumortier, F.; Roussarie, R.; Sotomayor, J. Generic 3-parameter families of vector fields on the plane, unfolding a singularity with nilpotent linear part. The cusp case of codimension 3. Ergod. Theory Dyn. Syst.; 1987; 7, pp. 375-413. [DOI: https://dx.doi.org/10.1017/S0143385700004119]

26. Dumortier, F.; Roussarie, R.; Sotomayor, J.; Zoladek, H. Bifurcations of Planar Vector Fields: Nilpotent Singularities and Abelian Integrals; Springer: Berlin/Heidelberg, Germany, 2006.