Stability and Hopf bifurcation of an SIR epidemic model with density-dependent transmission and Allee effect

1.   Introduction

For the past few years, with the acceleration of global warming and the processes of industrialization, urbanization and globalization, the natural environment has become unbalanced to some extent, and new epidemics such as COVID-19 continue to emerge. Increasingly more people believe that infectious diseases threaten the survival of human beings (see ^{[1,2,3,4,5,6,7,8]} and the references cited therein). Understanding the transmission process of disease is of great significance for studying the dynamics of disease ^[9]. By analyzing mathematical models, the dynamics of epidemic diseases transmission can be obtained and the future development trend of a disease can be further grasped (the disease disappears or becomes endemic) in order to formulate effective control measures. While mathematical modeling is a significant element in understanding the epidemiological characteristics of diseases, many such models are based on very simple demographic processes. It has been reported that new dynamic behaviors emerge from synergies between epidemiological and demographic processes that do not emerge in epidemiological models that assume constant population magnitude ^[10]. Therefore, studies of the transmission dynamics of epidemic diseases should not only consider the effects of simple demographic processes. Over a period of time, a non-negligible gap exists between the lifespan of animal populations infected with fatal diseases and their life expectancy ^[11]. For infectious diseases with a more extended incubation period, the birth and mortality rates of the host population must be considered ^[5], which may be influenced by strong Allee effects.

The Allee effect takes into account the social characteristics, mutual aggregation and cooperation of members in an animal population ^[12]. When the population density falls below a critical level, the population reproduction rate will decrease, a phenomenon known as the Allee effect ^[13]. The strong Allee effect typically occurs in endangered and small-scale populations, possibly due to genetic inbreeding, difficulty finding mates at small-scale populations, population randomness or reduced cooperative interactions ^{[14,15,16,17]}. The principal result of the Allee effect is that when the density of a population or social group falls below a threshold, the population has the possibility of extinction ^[15]. Hilker et al. ^[11] demonstrated that rich and fascinating dynamical behaviors can arise from strong Allee effects, such as persistent oscillation and multiple steady states, even in simple epidemiological models.

In recent years, increasing numbers of researchers have focused on the dynamic behavior of predator-prey interactions with the Allee effect ^{[18,19,20,21,22,23,24,25,26,27,28]}. At the same time, some researchers have paid attention to modeling and mathematical analysis of the interaction between epidemics and the Allee effect (see ^{[5,11,13,29,30,31]} and references therein). The study by Hilker et al. ^[11] suggested that infectious diseases could have essential effects on populations with strong Allee effects, as any slight reduction in population density could bring population sizes below critical levels and lead to population extinction. Studies have shown that some species with the Allee effect are also plagued by infectious diseases ^{[32,33,34,35,36,37,38]}. Based on these conditions, we have sufficient reasons to reveal and analyze the interaction between epidemics and the Allee effect. Kang and Castillo-Chavez ^[39] argued that biological conservation theories could assess the vulnerability of populations with the Allee effect to experiencing infectious disease outbreaks. Populations experiencing the Allee effect, being subjected to extinction or both, must be managed as efficiently as possible ^[11,40]. Deredec and Courchamp ^[29] and Hilker et al. ^[6] indicated that the integration of the Allee effect and parasitism raises the possibility of extinction. Thieme et al. ^[31] and Hilker et al. ^[11] investigated the impact of the density-dependent spread of infectious diseases on host populations and showed that host populations might be at risk of extinction. Friedman and Yakubu ^[30] demonstrated that understanding the interaction between disease epidemiological models and the Allee effect may have significant influences on epidemiology, ecology and conservation biology. Lv et al. ^[41] performed a theoretical analysis of an SI infectious disease model with additive Allee effect and delay, where the time delay is represented due to the latency period of the disease. The elementary mathematical characteristic of the model, as well as the equilibrium point stability and the bifurcation behavior with and without time delay, were investigated. The results show that the SI model with additive Allee effect and delay can produce rich dynamic behavior such as saddle-node bifurcation, supercritical bifurcation, Hopf bifurcation and limit cycle, which is an exciting study. Yin et al. ^[42] considered a Leslie-Gower predator-prey system with Allee effect and prey refuge, investigated the bifurcation behavior by taking the prey refuge constant as a parameter, including saddle-node bifurcation, Hopf bifurcation and Bogdanov-Takens bifurcation. More notably, it is found that the model has two limit cycles. In addition, the Allee effect has been recently investigated in several references, yielding many rich and exciting results to varying degrees (see ^[43,44,45]).

For the modeling of epidemiological models, in addition to considering the impact of demographic processes on the population, it is also essential to consider the incidence function of infectious diseases. Infectious disease transmission involves both frequency and density dependence. Density dependence presumes that the number of contacts is proportional to population size, while frequency dependence means there is no correlation between the two. Kang and Castillo-Chavez ^[46] reported that a density-dependent pathogen infecting a single host had a low probability of causing host extinction, but it did lead to the genocide of the pathogen. As mentioned by Hilker ^[47], density-dependent transmission is a necessary condition for achieving cyclic dynamics. They also found that if the spread of infectious diseases were density dependent, it would likely bring about more elaborate dynamics ^[11]. Hence, we consider density-dependent transmission as the form of disease transmission and establish an SIR model.

In Section 2, we describe the fundamental assumptions of disease transmission and build a mathematical model. In Sections 3 and 4, we obtain the general properties of the mathematical model, e.g., boundedness of solutions and the emergence of equilibrium points, and then the basic reproduction number is defined. In Sections 5 and 6, we perform theoretical analyses of the stability of equilibrium points as well as prove the emergence of the Hopf bifurcation. We simulate the stability of equilibrium points, and the emergence and stability of the limit cycle, by choosing different parameter combinations in Section 7. In the last section, we give a concise biological explanation and conclusions.

2.   Model description

We let $S\left( t \right)$ means the susceptible, $I(t)$ indicates the infective and $R(t)$ refers to the recovered, with $N\left( t \right)$ is the sum of the three part, denoting the population size. First of all, we assume that there is horizontal transmission, namely, the infected do not transmit the virus to their descendant. Meanwhile, we presume that the infectious disease transmission is density-dependent and transferred as bilinear incidence rate $\beta SI$, with $\beta > 0$. The path of disease transmission in the population is as follows:

According to reference ^[47], the net per capita growth rate is quadratic in the Allee effect:

The parameter $r$ represents the intrinsic rate of growth. $M$ indicates the maximal capacity of the population. $m$ suggests the Allee threshold. In most biologically meaningful cases, we have $0 < m < M$ and $m\ll M$.

Hence, the mathematical model can be expressed as

parameter $\delta$ and $\mu$ are disease-related mortality rate and the natural death rate, respectively. $\alpha$ indicates the recovery rate.

3.   The positive and bounded of the solution

In the first place, we need to verify the solutions of model (2) with positive initial value will keep positive for $t > 0$. We will prove Theorem 1 to obtain the conclusion.

Theorem 1. Let $S(0) > 0$, $I(0) \geqslant 0$, $R(0) \geqslant 0$, for $t > 0$, and the solution $S(t)$, $I(t)$, $R(t)$ of the model (2) is positive.

Proof. Since the equation on the right-hand side of the model (2) is a continuous smooth function on ${\mathbb{R}}^{3} = \left\{S\left(t\right), I\left(t\right), R\left(t\right):S\left(t\right) > 0, I\left(t\right)\ge 0, R\left(t\right)\ge 0\right\}$, we have the following result for the first equation of model

then we can get

Similarly, from the other equations of model (2), we also can get

for all $t \geqslant 0$.

Therefore, we can see that the solution $S(t) > 0$, $I(t) \geqslant 0$ and $R(t) \geqslant 0$ for all $t \geqslant 0$. Then we also can get the bounded of the solution from the follow theorem.

Theorem 2. The model (2) of the initial condition in $R_ + ^3$ is positively invariant in

Proof. Add all equations in model (2) to get:

we have

And

Then according to comparison principle we can easily deduce the solution of the equation

when $t \to \infty$, we have

Let $\Omega = \left\{ {\left. {(S, I, R) \in R_ + ^3} \right|0 \leqslant S(t) + I(t) + R(t) \leqslant \frac{{mMr}}{{4\mu \left( {m + M} \right)}}} \right\}$, therefore, we can examine the dynamics behavior of model (2) in $\Omega$.

4.   The emergence of equilibrium and the basic reproduction number

In this section, we will consider whether the model (2) has equilibrium points and determine the basic reproduction number. Let the right part of the equation of the model (2) equal to zero, for the disease-free states, we can make $I = 0$, then we can obtain $R = 0$. Hence, we just need to consider the first equation of the model (2):

By calculation, we can get that the model has an extinction state $O(0, 0, 0)$. And when ${D_1} > 0$, the model has two disease-free equilibrium points ${A_1}({S_{01}}, 0, 0)$ and ${A_2}({S_{02}}, 0, 0)$; when ${D_1} = 0$, the model has a disease-free equilibrium point ${A_3}({S_{00}}, 0, 0)$, where

In order to obtain the conditions for disease prevalence, the basic reproduction number is defined as

If ${D_1} > 0$, compare the relationship between ${R_0}$ and 1, we can deduce the following results:

(1) ${R_0} > 1 \Leftrightarrow {S_{01}} > \frac{{\mu + \delta + \alpha }}{\beta } > {S_{02}};$

(2) ${R_0} < 1 \Leftrightarrow {S_{02}} > \frac{{\mu + \delta + \alpha }}{\beta }$ or ${S_{01}} < \frac{{\mu + \delta + \alpha }}{\beta }$;

(3) ${R_0} = 1 \Leftrightarrow \frac{{\mu + \delta + \alpha }}{\beta } = {S_{01}}$ or $\frac{{\mu + \delta + \alpha }}{\beta } = {S_{02}}$.

If ${D_1} = 0$, compare the relationship between ${R_0}$ and 1, we can deduce the following results:

(1) ${R_0} < 1 \Leftrightarrow \frac{{\mu + \delta + \alpha }}{\beta } < \frac{{M + m}}{2}$ or $\frac{{\mu + \delta + \alpha }}{\beta } > \frac{{M + m}}{2}$;

(2) ${R_0} = 1 \Leftrightarrow \frac{{\mu + \delta + \alpha }}{\beta } = \frac{{M + m}}{2}.$

If ${R_0} > 1$, we make the right part of the equation of the model (2) equal to zero, then, there is a unique endemic equilibrium $A({S^*}, {I^*}, {R^*})$, with

5.   Local stability of the equilibrium

After determining the equilibrium points and ${R_0}$ of the model (2), it is necessary to grasp the equilibrium points' behavior to understand the model's dynamic behavior. When the disease-free equilibrium point is stable, the number of infections gradually tends to zero and the disease disappears, whereas when the endemic equilibrium point is stable, the disease becomes endemic and the number of infections does not tend to zero. Obviously, we prefer to keep the disease-free equilibrium point of the model stable. Therefore, we will discuss the stability of the extinction state $O$, the disease-free equilibrium ${A_1}$, ${A_2}$, ${A_3}$ and the endemic equilibrium $A$ by linear stability analysis in this section, i.e., using the sign of the characteristic roots of the characteristic equation of the Jacobi matrix at the equilibrium point to determine the stability of the equilibrium point.

Theorem 3. (a) If ${R_0} < 1$ and ${S_{01}} > {S_2}$, the disease-free equilibrium ${A_1}({S_{01}}, 0, 0)$ is locally asymptotically stable. If ${R_0} < 1$ and ${S_1} < {S_{01}} < {S_2}$, ${A_1}({S_{01}}, 0, 0)$ is unstable;

(b) If ${R_0} < 1$ and ${S_{02}} > {S_2}$, the disease-free equilibrium ${A_2}({S_{02}}, 0, 0)$ is locally asymptotically stable. If ${R_0} < 1$ and ${S_1} < {S_{02}} < {S_2}$, ${A_2}({S_{02}}, 0, 0)$ is unstable, where

(c) the disease-free equilibrium ${A_3}({S_{00}}, 0, 0)$ is a saddle node, ${A_3}$ is a repelling saddle-node when $\beta {S_{00}} - \left( {\mu + \delta + \alpha } \right) > 0$ and an attracting saddle-node when $\beta {S_{00}} - \left( {\mu + \delta + \alpha } \right) < 0$;

(d) the extinction state $O\left( {0, 0, 0} \right)$ is always stable.

Proof. (a) The Jacobian at ${A_1}({S_{01}}, 0, 0)$ is given

The characteristic equation is

the eigenvalues are ${\lambda _1} = - \frac{{3r}}{{mM}}{S_{01}}^2 + \frac{{2r\left( {m + M} \right)}}{{mM}}{S_{01}} - \left( {r + \mu } \right)$, ${\lambda _2} = \beta {S_{01}} - \left( {\mu + \delta + \alpha } \right)$, ${\lambda _3} = - \mu$.

Next, we discuss the eigenvalues of ${A_1}({S_{01}}, 0, 0)$ are greater than zero or not.

Firstly, for ${\lambda _1} = - \frac{{3r}}{{mM}}{S_{01}}^2 + \frac{{2r\left( {m + M} \right)}}{{mM}}{S_{01}} - \left( {r + \mu } \right)$, let

thus

It is clear to see that

Then we have

Next, we discuss that there are always ${S_1} < {S_{01}}$.

There is no doubt that

so ${S_1} < {S_{01}}$. Therefore, when ${S_1} < {S_{01}} < {S_2}$, we can get ${\lambda _1} > 0$, then in this case the disease-free equilibrium ${A_1}({S_{01}}, 0, 0)$ is unstable. When ${S_{01}} > {S_2}$, we have ${\lambda _1} < 0$, then we need to further verify whether ${\lambda _2} = \beta {S_{01}} - (\mu + \delta + \alpha )$ is greater than zero or less than zero.

Finally, for ${\lambda _2} = \beta {S_{01}} - (\mu + \delta + \alpha )$, according to the equation

we can get ${\lambda _2} < 0$ if ${R_0} < 1$ and ${\lambda _2} > 0$ if ${R_0} > 1$.

So, if ${R_0} < 1$ and ${S_{01}} > {S_2}$, the roots of the characteristic equation are ${\lambda _1} = - \frac{{3r}}{{mM}}{S_{01}}^2 + \frac{{2r\left( {m + M} \right)}}{{mM}}{S_{01}} - \left( {r + \mu } \right) < 0$, ${\lambda _2} = \beta {S_{01}} - \left( {\mu + \delta + \alpha } \right) < 0$, ${\lambda _3} = - \mu < 0$, hence, there is a locally asymptotically stable equilibrium ${A_1}({S_{01}}, 0, 0)$.

(b) The Jacobian at ${A_2}({S_{02}}, 0, 0)$ is given

The characteristic equation is

the eigenvalues are ${\lambda _1} = - \frac{{3r}}{{mM}}{S_{02}}^2 + \frac{{2r\left( {m + M} \right)}}{{mM}}{S_{02}} - \left( {r + \mu } \right)$, ${\lambda _2} = \beta {S_{02}} - \left( {\mu + \delta + \alpha } \right)$, ${\lambda _3} = - \mu$.

Similar to (a), next, we discuss that there are always ${S_1} < {S_{02}}$.

we can also get if ${S_1} < {S_{02}} < {S_2}$, then ${\lambda _1} > 0$, and if ${S_{02}} > {S_2}$, then ${\lambda _1} < 0$.

For ${\lambda _2} = \beta \left[ {{S_{02}} - \frac{{\mu + \delta + \alpha }}{\beta }} \right]$, we can know if ${R_0} < 1$, then ${\lambda _2} < 0$, and if ${R_0} > 1$, then ${\lambda _2} > 0$. So, if ${R_0} < 1$ and ${S_{02}} > {S_2}$, the roots of the characteristic equation are ${\lambda _1} = - \frac{{3r}}{{mM}}{S_{02}}^2 + \frac{{2r\left( {m + M} \right)}}{{mM}}{S_{02}} - \left( {r + \mu } \right) < 0$, ${\lambda _2} = \beta {S_{02}} - \left( {\mu + \delta + \alpha } \right) < 0$, ${\lambda _3} = - \mu < 0$. Therefore, there is a locally asymptotically stable equilibrium ${A_2}({S_{02}}, 0, 0)$.

(c) Similarly, the Jacobian at ${A_3}({S_{00}}, 0, 0)$ is given

The characteristic equation is

the eigenvalues are ${\lambda _1} = 0$, ${\lambda _2} = \beta {S_{00}} - \left( {\mu + \delta + \alpha } \right)$, ${\lambda _3} = - \mu$.

By calculation, it can be concluded that ${A_3}({S_{00}}, 0, 0)$ is a saddle node, and the detailed proof is shown in the Appendix B.

(d) The Jacobian at $O\left( {0, 0, 0} \right)$ is given by

The characteristic equation is

the eigenvalues are ${\lambda _1} = - \left( {r + \mu } \right)$, ${\lambda _2} = - \left( {\mu + \delta + \alpha } \right)$, ${\lambda _3} = - \mu$. Hence, the extinction state $O\left( {0, 0, 0} \right)$ of the model (2) is always stable.

Theorem 4. If $\frac{{\mu + \delta + \alpha }}{\beta } > \frac{{m + M}}{2}$ and ${R_0} > 1$, the endemic equilibrium $A\left( {{S^*}, {I^*}, {R^*}} \right)$ is locally asymptotically stable. If $\frac{{m + M}}{2} > \frac{{\mu + \delta + \alpha }}{\beta }$ and ${R_0} > 1$, there is an unstable endemic equilibrium $A\left( {{S^*}, {I^*}, {R^*}} \right)$.

Proof. If ${R_0} > 1$ (that is ${S_{01}} > \frac{{\mu + \delta + \alpha }}{\beta } > {S_{02}}$), the Jacobian at the endemic equilibrium $A\left( {{S^*}, {I^*}, {R^*}} \right)$ is given by

Let

The characteristic equation of ${J_1}$ is

where

Hence, we get if $\frac{{\mu + \delta + \alpha }}{\beta } > \frac{{m + M}}{2}$, then $\sigma < 0$; if $\frac{{\mu + \delta + \alpha }}{\beta } < \frac{{m + M}}{2}$, then $\sigma > 0$.

Therefore, the characteristic equation of the Jacobian ${J_{A\left( {{S^*}, {I^*}, {R^*}} \right)}}$ is $\left( {{\lambda ^2} - \sigma \lambda + \Delta } \right)\left( {\lambda + \mu } \right) = 0$.

Clearly, if $\frac{{\mu + \delta + \alpha }}{\beta } > \frac{{m + M}}{2}$ and ${R_0} > 1$ (that is, ${S_{01}} > \frac{{\mu + \delta + \alpha }}{\beta } > \frac{{m + M}}{2}$), then $\sigma < 0$ and $\Delta > 0$. Thus the real part of all the characteristic roots are less than zero, that is, there is a locally asymptotically stable equilibrium $A\left( {{S^*}, {I^*}, {R^*}} \right)$. If $\frac{{\mu + \delta + \alpha }}{\beta } < \frac{{m + M}}{2}$ and ${R_0} > 1$ (that is, $\frac{{m + M}}{2} > \frac{{\mu + \delta + \alpha }}{\beta } > {S_{02}}$), then $\sigma > 0$ and $\Delta > 0$. Hence, the real part of all the characteristic roots is more than zero and equilibrium $A\left( {{S^*}, {I^*}, {R^*}} \right)$ is unstable.

6.   Bifurcation analysis

Because the first two equations of the model (2) are independent of the third equation, we consider to analyze the following model:

Clearly, the model (41) have a positively invariant set ${\mathrm{\Omega }}_{1}$, where

We perform an analysis of the model (41) in ${\mathrm{\Omega }}_{1}$.The model (41) always has a extinction state $P\left( {0, 0} \right)$, the disease-free equilibrium points ${B_1}\left( {{S_{01}}, 0} \right)$, ${B_2}\left( {{S_{02}}, 0} \right)$, ${B_3}\left( {{S_{00}}, 0} \right)$ If ${R_0} > 1$, the model has a unique endemic equilibrium $E\left( {{S^*}, {I^*}} \right)$.

Through the above qualitative analysis, it is easy to obtain that the Jacobian at $E\left( {{S^*}, {I^*}} \right)$ is given by

Then we can find that the characteristic equation is ${\lambda ^2} - \sigma \lambda + \Delta = 0$, where $\sigma$ and $\Delta$ are the same as mentioned above. And we have if $\frac{{\mu + \delta + \alpha }}{\beta } > \frac{{m + M}}{2}$, then $\sigma < 0$; if $\frac{{\mu + \delta + \alpha }}{\beta } < \frac{{m + M}}{2}$, then $\sigma > 0$; if $\frac{{\mu + \delta + \alpha }}{\beta } = \frac{{m + M}}{2}$, that is, $\beta = \frac{{2\left( {\mu + \delta + \alpha } \right)}}{{m + M}}$, then $\sigma = 0$.

Choosing $\beta$ as the bifurcation parameter. Next, we will analyze that model (41) goes through a Hopf bifurcation. Prior to this, we introduce a lemma.

Lemma 1 ^[48]: Any general one-parameter, two-dimensional system $\dot x = f\left( {x, \alpha } \right)$, $x\in {\mathbb{R}}^{2}$, $\alpha \in {\mathbb{R}}^{1}$, with smooth $f$, having for all sufficiently small $\left| \alpha \right|$ the equilibrium $x = 0$ with eigenvalues ${\lambda _{1, 2}}\left( \alpha \right) = \varphi \left( \alpha \right) \pm i\omega \left( \alpha \right),$ $\varphi \left( 0 \right) = 0,$ $\omega \left( 0 \right) = {\omega _0} > 0$ and at $\alpha = 0$ the equilibrium $x = 0$ with eigenvalues ${\lambda _{1, 2}}\left( 0 \right) = \pm i{\omega _0}$, is locally topologically corresponding around the origin to one of the following normal forms:

The genericity conditions assumed in Lemma 1 are as follows:

(B.1) ${l_1} = \frac{1}{{2{\omega _0}^2}}{Re} \left( {i{g_{20}}{g_{11}} + {\omega _0}{g_{21}}} \right) \ne 0$, where ${l_1}$ is the first Lyapunov coefficient, and

(B.2) ${\left. {\frac{{d\varphi }}{{d\alpha }}} \right|_{\alpha = 0}} \ne 0.$

Theorem 5. Model (41) has a stable limit cycle for $\beta = {\beta _0}$, where ${\beta _0} = \frac{{2\left( {\mu + \delta + \alpha } \right)}}{{m + M}}$.

Proof. If ${R_0} > 1$, model (41) has an endemic equilibrium point $E\left( {{S^*}, {I^*}} \right)$, where

It is not hard to find that the Jacobian at $E\left( {{S^*}, {I^*}} \right)$ is given

the characteristic equation is

where $\sigma = \sigma \left( \beta \right) = trA\left( \beta \right)$, $\Delta = \Delta \left( \beta \right) = \det A\left( \beta \right)$, so, ${\lambda _{1, 2}}\left( \beta \right) = \frac{1}{2}\left( {\sigma \left( \beta \right) \pm \sqrt {{\sigma ^2}\left( \beta \right) - 4\Delta \left( \beta \right)} } \right)$.

The Hopf bifurcation condition implies that the characteristic equation has a pair of pure imaginary roots, i.e., it satisfies the condition that $\sigma \left( 0 \right) = 0$, $\Delta \left( 0 \right) = {\omega _0}^2 > 0$.

For small $\left| \beta \right|$ we can introduce $\varphi \left( \beta \right) = \frac{1}{2}\sigma \left( \beta \right)$, $\omega \left( \beta \right) = \frac{1}{2}\sqrt {4\Delta \left( \beta \right) - {\sigma ^2}\left( \beta \right)}$, and therefore obtain the following eigenvalues: ${\lambda _1}\left( \beta \right) = \lambda \left( \beta \right)$, ${\lambda _2}\left( \beta \right) = \overline {\lambda \left( \beta \right)}$, where

and thus

We have $\varphi \left( {{\beta _0}} \right) = 0$ for ${\beta _0} = \frac{{2\left( {\mu + \delta + \alpha } \right)}}{{M + m}}$.

Moreover,

Further on,

Therefore, at $\beta = {\beta _0}$, the characteristic equation has eigenvalues ${\lambda _1} = \lambda \left( {{\beta _0}} \right) = i\omega \left( {{\beta _0}} \right)$, ${\lambda _2} = \overline {\lambda \left( {{\beta _0}} \right)} = - i\omega \left( {{\beta _0}} \right)$, a Hopf bifurcation emerges. Meanwhile, equilibrium $E\left( {{S^*}, {I^*}} \right)$ is unstable for $\beta > {\beta _0}$ and stable for $\beta < {\beta _0}$.

To discuss the bifurcation behavior around $E\left( {{S^*}, {I^*}} \right)$, we need to verify whether all the qualifications of Lemma 1 are sufficient. It is easy to check the condition (B.2):

Let the parameter $\beta$ be fixed at ${\beta _0}$, we have ${I^{\left( 0 \right)}} = \frac{{m + M}}{{2\left( {\mu + \delta + \alpha } \right)}}\left[ {\frac{{r{{\left( {m + M} \right)}^2}}}{{4mM}} - \left( {r + \mu } \right)} \right]$, and ${S^{\left( 0 \right)}} = \frac{{m + M}}{2}$.

Transform the origin of the coordinates to $\left( {{S^{\left( 0 \right)}}, {I^{\left( 0 \right)}}} \right)$ by the conversion of variables

This transforms model (41) into

This model can be represented as

where $A{\text{ }} = {\text{ }}A\left( {{\beta _0}} \right)$, and the functions $B$ and $C$ carry out the planar vectors $\xi = {\left( {{\xi _1}, {\text{ }}{\xi _2}} \right)^T}$, $\eta = {\left( {{\eta _1}, {\text{ }}{\eta _2}} \right)^T}$ and $\zeta {\text{ }} = {\text{ }}{\left( {{\zeta _1}, {\text{ }}{\zeta _2}} \right)^T}$ the values

and

We can write $A\left( {{\beta _0}} \right)$ in the form

where ${\omega _0}^2$ is given by Eq (50). Obviously, it is convenient to obtain the following complex vectors

Let $q\left(\alpha \right)\in {\mathbb{C}}^{2}$ be an eigenvector of $A\left( \alpha \right)$ belonging to $\lambda \left( \alpha \right)$ :

let $p\left(\alpha \right)\in {\mathbb{C}}^{2}$ be an eigenvector of ${A^T}\left( \alpha \right)$ belonging to $\overline {\lambda \left( \alpha \right)}$ :

Get $p$, $q$ that satisfies the vital normalization $\left\langle {p, q} \right\rangle = \overline {{p_1}} {q_1} + \overline {{p_2}} {q_2} = 1$. To attain the vital normalization $\left\langle {p, q} \right\rangle = \overline {{p_1}} {q_1} + \overline {{p_2}} {q_2} = 1$, we can take

Now we can simply calculate

the first Lyapunov coefficient can be obtained by calculation is

Hence, (B.1) of Lemma 1 is satisfied. Therefore, there is a unique and stable limit cycle for $\beta = {\beta _0}$.

The proof of existence and stability of the limit cycle for the model (41) is obtained by first translating the model to the origin to obtain a new model with a linear representation of the $A$, $B$, $C$ matrix, and then taking the eigenvectors $p$ and $q$ corresponding to the eigenvalues of the matrix $A$ when the condition $\left\langle {p, q} \right\rangle = \overline {{p_1}} {q_1} + \overline {{p_2}} {q_2} = 1$ is satisfied. Further, letting $\varphi \left( \beta \right) = \frac{{tr\left( A \right)}}{2}$, we can obtain ${\left. {\frac{{d\varphi }}{{d\beta }}} \right|_{\beta = {\beta _0}}} \ne 0$. Eventually, ${l_1} = \frac{1}{{2{\omega _0}^2}}{Re} \left( {i{g_{20}}{g_{11}} + {\omega _0}{g_{21}}} \right)$ is obtained by computing ${g_{20}}$, ${g_{11}}$ and ${g_{21}}$. In case ${l_1} < 0$, the model has a stable limit cycle; in case ${l_1} > 0$, the model has an unstable limit cycle; in case ${l_1} = 0$, it is necessary to further prove whether the model has more than one limit cycle. It is obvious that ${l_1} < 0$ in the proposed model, so that a stable limit cycle exists when $\beta = {\beta _0}$.

Remark: We find that the model without the third equation of the model (2) has a stable limit cycle, because the first two equations of the model (2) are independent of the third equation, and we only discuss the existence of the limit cycle of the simplified model, i.e., there is a unique stable limit cycle when $\beta {\text{ }} = {\text{ }}{\beta _0}$. In fact, there is a stable limit cycle of the model (2) when the same conditions are satisfied.

In a biological sense, when a stable limit cycle is generated at the endemic equilibrium point, a recurrence of disease is implied, a situation that has greater challenges for disease control.

7.   Numerical simulation

In this section, to confirm the above conclusions, we mainly use Matlab draw numerical results. We select $r$, $m$, $M$, $\mu$, $\delta$ and $\alpha$ are fixed, to give detailed numerical simulations of the stability of the equilibrium points, as well as the emergence and stability of limit cycle.

7.1.   Stability of the extinction state $O\left( {0, 0, 0} \right)$

We choose $r = 0.04386$, $M = 100$, $m = 2$, $\beta = 0.05$, $\mu = 0.01246$, $\delta = 0.915$ and $\alpha = 0.6$, then the extinction state $O\left( {0, 0, 0} \right)$ is locally asymptotically stable in Figure 2.

7.2.   Stability of the disease-free equilibrium ${A_1}\left( {{S_{01}}, 0, 0} \right)$

We choose $r = 0.04386$, $M = 100$, $m = 2$, $\beta = 0.008$, $\mu = 0.01246$, $\delta = 0.915$ and $\alpha = 0.6$, then ${R_0} = 0.5213$, ${S_{01}} = 99.4168$, ${S_2} = 66.7169$, we have ${R_0} < 1$ and ${S_{01}} > {S_2}$, hence, Figure 3 shows a locally asymptotically stable equilibrium ${A_1}\left( {{S_{01}}, 0, 0} \right)$.

7.3.   Stability of the endemic equilibrium $A\left( {{S^*}, {I^*}, {R^*}} \right)$

We choose $r = 0.04386$, $M = 100$, $m = 2$, $\beta = 0.0256$, $\mu = 0.01246$, $\delta = 0.915$ and $\alpha = 0.6$, then ${R_0} = 1.3633$, ${S_{01}} = 99.4168$, $\frac{{\mu + \delta + \alpha }}{\beta } = 59.6664$, $\frac{{m + M}}{2} = 51$, we have ${S_{01}} > \frac{{\mu + \delta + \alpha }}{\beta } > \frac{{m + M}}{2}$, thus $A\left( {{S^*}, {I^*}, {R^*}} \right)$ is locally asymptotically stable in Figure 4.

7.4.   Stability of the limit cycle

We choose $r = 0.04386$, $M = 100$, $m = 2$, $\beta = 0.029950197$, $\mu = 0.01246$, $\delta = 0.915$ and $\alpha = 0.6$, the initial values as ${y_0} = \left( {30, 5, 1} \right)$, ${y_0} = \left( {51, 12, 3} \right)$ and ${y_0} = \left( {52, 17.5, 1500} \right)$. We have ${R_0} = 1.4842$, then the initial value ${y_0} = \left( {30, 5, 1} \right)$ and ${y_0} = \left( {51, 12, 3} \right)$ from under the limit cycle (red line) is moving towards the limit cycle; in addition, so as to better understand the emergence and stability of the limit cycle from a mathematical point of view, we find an initial value ${y_0} = \left( {52, 17.5, 1500} \right)$ above the limit cycle, which will move towards the red line from above it. Therefore, there is a stable limit cycle in Figure 5.

8.   Discussion

In the study of infectious diseases, we should pay attention not only to the dynamic behavior of the disease, but also to the demographic processes. The Allee effect, an ecological phenomenon, can influence population dynamics ^[49]. Considering the significance of demographic processes for grasping disease transmission, we constructed an SIR epidemic model with density-dependent transmission and a strong Allee effect. We demonstrated that fascinating dynamics could arise in a simple epidemiological model. The local asymptotical stability of each equilibrium of the model (2) is discussed and the emergence and uniqueness of the limit cycle were analyzed.

We find that the model's dynamics depend not only on the basic reproduction number. The extinction state $O(0, 0, 0)$ is always stable. The disease-free equilibrium ${A_3}({S_{00}}, 0, 0)$ is an attracting saddle-node if $\beta {S_{00}} < \mu + \delta + \alpha$, i.e., ${S_{00}} < \frac{{\mu + \delta + \alpha }}{\beta }$. The disease-free equilibrium ${A_2}({S_{02}}, 0, 0)$ exists and is locally asymptotically stable if ${R_0} < 1$ and ${S_{02}} > {S_2}$, i.e., ${S_{02}} > \frac{{\mu + \delta + \alpha }}{\beta }$ and ${S_{02}} > {S_2}$. The disease-free equilibrium ${A_1}({S_{01}}, 0, 0)$ emerges and is locally asymptotically stable if ${S_{01}} > {S_2}$ and ${R_0} < 1$, i.e., $\frac{{\mu + \delta + \alpha }}{\beta } > {S_{01}} > {S_2}$. If $\beta$ increases or $\mu + \delta + \alpha$ decreases so that ${S_{01}} > \frac{{\mu + \delta + \alpha }}{\beta } > \frac{{m + M}}{2}$ and thus ${R_0} > 1$, then the endemic equilibrium $A({S^*}, {I^*}, {R^*})$ emerges and is locally asymptotically stable. Therefore, the development of a disease as an endemic disease can be controlled by reducing contact between infected and susceptible individuals, or by increasing the mortality rate of infected individuals, which is also in line with the current policy related to quarantine measures to control further spread of a disease. If $\frac{{m + M}}{2} > \frac{{\mu + \delta + \alpha }}{\beta } > {S_{02}}$ so that ${R_0} > 1$ holds, then a locally asymptotically stable limit cycle emerges near the unstable equilibrium $A\left( {{S^*}, {I^*}, {R^*}} \right)$. The Hopf bifurcation caused by density-dependent transmission is verified in detail, uniqueness of the Hopf bifurcation is obtained. Our numerical simulations confirmed the result as well.

As mentioned earlier, the synergy between the Allee effect and epidemic disease could induce interesting disease transmission dynamics, such as bistability and sustained oscillations. Similar models with logistic growth do not emerge in these phenomena ^[10]. Hence, these complex dynamics are primarily caused by the Allee effect.

Furthermore, the propagation is supposed to be density dependent and is transferred as the bilinear incidence rate $\beta SI$, where $\beta$ the quantity of contacts, increases linearly with population density. As Hilker et al. ^[47] mentioned, the essential condition of achieving cyclic dynamics is density-dependent transmission. For the infectious disease transmission model with the Allee effect, density-dependent transmission displays the appearance and disappearance of stable oscillations, while frequency-dependent transmission (also called normal incidence) manifests no such phenomenon. We obtained the emergence and stability of the limit cycle, consistent with the conjecture of Hilker et al. ^[11].

In terms of mathematical and biological applications, the emergence of a limit cycle and bistability in this epidemic model is both exciting and meaningful. On one hand, from the perspective of biological conservation, the existence of a limit cycle can explain the periodic existence of some populations, and in order to maintain the diversity of biological populations, we hope that limit cycles can exist; from the perspective of infectious diseases, the existence of a limit cycle explains the periodic recurrence and disappearance of some diseases, and if a certain endemic equilibrium point is perturbed by parameters to produce a stable limit cycle, then the disease might periodically recur, and this result is not desired and must be avoided, similar results have been presented in references ^[11,13]. On the other hand, when the model reaches bistability, small perturbations in the parameters may cause the model to reach extinction rapidly because the extinction state in the model is always stable. Thus, bistability may imply more cautious measures for infectious disease control because some control measures may lead to perturbations of key parameters and thus affect the stability of the model equilibrium point. Therefore, the above conclusions may be informative for biological conservation and infectious disease control.

In summary, adding an Allee effect to epidemiological models with density-dependent transmission allows a rich spectrum of dynamical behavior, including infectious diseases showing periodic oscillations. It will be fascinating to consider a delay differential equation model in future mathematical analyses of disease.

Acknowledgments

This work was supported by the Gansu Science and Technology Fund (20JR5RA512), the Research Fund for Humanities and Social Sciences of the Ministry of Education (20XJAZH006) and the Fundamental Research Funds for the Central Universities (31920220066).

Conflict of interest

The authors declare that there are no personal or competing interest relationships to the best of their knowledge that could have influenced the study results.

Appendix A

In this section, we discuss the method of defining the basic reproduction number. For the unique endemic equilibrium point $A\left( {{S^*}, {I^*}, {R^*}} \right)$, we obtain the characteristic equation of the Jacobi matrix at this point is ${\lambda ^2} - \sigma \lambda + \Delta = 0$, where $\sigma$ and $\Delta$ denote the same meaning as in Eq (38), i.e.,

and

By simplifying $\Delta$, we obtain

Whether or not the endemic equilibrium point is stable, it must be ensured that $\Delta > 0$, that is,

Then

furthermore,

Therefore, we let the basic reproduction number

If ${R_0} > 1$ and $\sigma < 0$, there is a locally asymptotically stable endemic equilibrium point. If ${R_0} > 1$ and $\sigma > 0$, the endemic equilibrium point is unstable.

Appendix B

In this section we discuss the stability of the disease-free equilibrium ${A_3}\left( {{S_{00}}, 0, 0} \right)$. Since the first two equations of the model (2) are independent of the last one, the first two equations constitute model (41), which is equivalent to model (2), so here we discuss the stability of the disease-free equilibrium point ${B_3}\left( {{S_{00}}, 0} \right)$ of the model (41).

In order to discuss the stability of ${B_3}$, we first do a translational transformation of ${B_3}$ to the origin. Let $s = S - {S_{00}}, i = I$, the model becomes

The Jacobi matrix of this model at the origin is also the Jacobi matrix of model (41) at ${B_3}$. The calculation yields

Then ${J_{\left( {{B_3}} \right)}}$ has two eigenvalues, which are ${\lambda _1} = 0$ and ${\lambda _2} = \beta {S_{00}} - \left( {\mu + \delta + \alpha } \right)$.

Further, the construction of the transformation matrix

After passing the transformation

the model becomes

where

Introduce a new time variable $\tau = \left[ {\beta {S_{00}} - \left( {\mu + \delta + \alpha } \right)} \right]t$, we obtain the model

where

It can be seen that the coefficient of ${S^2}$ is ${c_2} = \frac{{r\left( {m + M} \right)}}{{2mM\left[ {\left( {\mu + \delta + \alpha } \right) - \beta {S_{00}}} \right]}}$. From the Theorem 7.1 in chapter 2 in ^[50], if ${\lambda _2} = \beta {S_{00}} - \left( {\mu + \delta + \alpha } \right) > 0$, then ${c_2} < 0, \tau > 0$, and ${B_3}$ is a repelling saddle-node. If ${\lambda _2} = \beta {S_{00}} - \left( {\mu + \delta + \alpha } \right) < 0$, then ${c_2} > 0, \tau < 0$, and ${B_3}$ is an attracting saddle-node.

Remark: Model (2) has three characteristic roots of the characteristic equation of the Jacobi at the disease-free equilibrium ${A_3}$, which are ${\lambda _1} = 0$, ${\lambda _2} = \beta {S_{00}} - \left( {\mu + \delta + \alpha } \right)$, ${\lambda _3} = - \mu$. It is important to emphasize that we have $\beta {S_{00}} - \left( {\mu + \delta + \alpha } \right) > 0$ or $\beta {S_{00}} - \left( {\mu + \delta + \alpha } \right) < 0$ when ${R_0} < 1$. The disease-free equilibrium ${A_3}$ of the model (2) is a saddle-node. Further, ${A_3}$ is a repelling saddle-node when ${\lambda _2} > 0$ and an attracting saddle-node when ${\lambda _2} < 0$.