Stability analysis and backward bifurcation on an SEIQR epidemic model with nonlinear innate immunity

1.   Introduction

Infectious diseases can spread and turn into epidemics, taking thousands of lives within a matter of just a few days. Many scholars have studied the epidemiological mechanism of infectious diseases through various methods and given preventive strategies. Mathematical models can help us to gain insights into the dynamics of diseases and their control strategies ^[1,2,3]. Epidemic models have been developed by many scholars since the first epidemiological model created by Daniel Bernoulli ^[4]. In ^[5], considering both individual behavioral responses and governmental actions, Lin et al. proposed a conceptual model for the COVID-19 outbreak in Wuhan, China. In order to forecast the evolution of the COVID-19 outbreak in Mexico, Avila-Ponce de Le$\acute{o}$n et al. proposed an SEIARD mathematical model which incorporated the asymptomatic infected individuals ^[6]. In ^[7], a bacterial meningitis transmission dynamics was considered by Asamoah et al. The existence of backward bifurcation was discussed and optimal control problem was solved. In ^[8], an ordinary differential model of malaria was established. In this paper, stability of disease-free and endemic equilibria, bifurcation phenomena were investigated. In ^[9], a delay Ebola epidemic model was studied by Al-Darabsah et al. In ^[10], He et al. built an SEIR model for COVID-19 incorporating some general control strategies. In ^[11], based on the COVID-19 data in Ghana and Egypt, Asamoah et al. formulated a COVID-19 infection model to present the sensitivity assessment and optimal economic evaluation. In ^[12], Zhao et al. studied a stochastic switched SIRS epidemic model. The stationary distribution and extinction of the disease were discussed. In ^[13], Omame et al. considered a co-infection model for SARS-CoV-2 and ZIKV, which exhibited backward bifurcation. In ^[14], considering the infectious force in latent period and infected period, Zhao et al. lucubrated an SEIR epidemic model with discontinuous treatment strategy.

For the sake of the effective strategies to disease control and prevention, quarantine is the most effective way to reduce the transmission of the infected to the susceptible. For example, many countries have taken quarantine measures in the fight against COVID-19. Many scholars introduced the quarantine class into the epidemic models. In ^[15], Herbert et al. showed six epidemic models with quarantine and different forms of the incidence. In ^[16], a compartmental model incorporating asymptomatic class, quarantine and isolation was presented by Ali et al. The strategies for effective control of the epidemic were proposed by analyzing the model. In ^[17], Tulu et al. built a fractional-order model for Ebola with the strategies to vaccination and quarantine. They gained that vaccination and quarantine are effective control measures for Ebola.

Innate immunity is the body's natural immune defense function formed during germline development and evolution ^[18]. It is a series of defense mechanisms formed by organisms in the process of long-term evolution. In ^[19], Kabir et al. built an epidemic model for natural and artificial immunity with durability and imperfectness of protection. In ^[20,21], the authors considered the mathematical models with the nonlinear innate immunity rate. However, there currently no models that comprehensively consider the effect of the nonlinear innate immunity and quarantine.

This paper is organized as follows. In Section 2, we built an SEIQR epidemic model and describe it. Preliminaries, such as boundedness, positivity of system (2.2) are discussed in Section 3. In Section 4, we present the expressions of equilibria and the basic reproduction number. The global stability result of the disease-free equilibrium and the existence result of backward bifurcation for system (2.2) are derived in Section 5. In Section 6, the global stability of endemic equilibrium is shown via using geometric approach. We present some numerical examples to verify the theoretical results obtained in previous sections. Finally, we conclude this article with a brief discussion.

2.   Model formulation

In this work, therefore, we will investigate the effects of quarantine and nonlinear innate immunity. And we shall consider an SEIR model with quarantine and nonlinear innate immunity. Suppose that $N(t)$ denotes the number of total population. The total population is subdivided into five different classes, namely, the susceptible ($S(t)$) individuals, the exposed individuals ($E(t)$), the infected individuals $(I(t))$, the quarantined individuals $(Q(t))$ and the recovered individuals $(R(t))$. Hence, $N(t) = S(t) + E(t) + I(t) + Q(t) + R(t)$. In the present paper, we assume that the recovered individuals confer permanent immunity and they do not revert to the susceptible individuals, like varicella, measles, rubella ^[22,23]. The flowchart of the epidemiological SEIQR model is illustrated in Figure 1.

From the flowchart, we will establish the following model:

In model (2.1), the parameters $A$, $\mu$, $d$, $\sigma$, $\xi$, $\rho$, $\gamma$, $\phi$, $\beta$, $a$ and $k$ are nonnegative. $A$ is the replenishment rate of susceptibles; $\beta$ is the infection rate; $\mu$ is the natural death rate of the population; $d$ denotes the disease-related mortality rate; $\sigma$ denotes the transfer rate between the exposed and the infectious; $\xi$ is quarantine rate of the infected; $\rho$ is the recovery rate of the infectious individuals; $\phi$ is the recovery rate of the quarantine individuals. We assume that the exposed transformed into the susceptible with the nonlinear innate immunity rate $\frac{aE(t)}{1+kE(t)}$ ^[20,21]. For biological significance, we postulate that the initial conditions of system (2.1) satisfy: $S(0)\geq 0, \, \, E(0)\geq 0, \, \, I(0)\geq 0, \, \, Q(0)\geq 0\, \, \mbox{and} \, \, R(0) \geq 0$.

Since $R(t)$, the recovered population, is independent of the first four equations of system (2.1), the rest of the paper will consider only the following four-dimensional system:

The initial conditions associated with the system (2.2) are as follows

3.   Basic properties

In this section, we shall present some basic properties such as boundedness, positivity of system (2.2).

Theorem 3.1. For all $t > 0$, the solutions $(S(t), E(t), I(t), Q(t))$ of system (2.2) with initial condition (2.3) are positive.

Proof. Set $t_1 = \sup\{t > 0:S(0)\geq 0, E(0)\geq 0, I(0)\geq 0, Q(0)\geq 0\}.$ The following inequality is given by the first equation of system (2.2)

The above inequality can be rewritten as:

Thus

so that

Similar to the above method, we can obtain $E(t)\geq 0$, $I(t)\geq0$, $Q(t)\geq0$ for all time $t > 0$. Hence, for all $t > 0$, the solutions of system (2.2) satisfying the initial value condition (2.3) are positive. This completes the proof of Theorem 3.1.

Define

Theorem 3.2. The region $\mathcal{D}$ is invariant, which indicates that all solutions of system (2.2) with initial condition (2.3) in $\mathcal{D}$ remain in $\mathcal{D}$ for all $t > 0$.

Proof. Adding the two sides of system (2.1) respectively, we have

Hence

From above we obtain that

Therefore, if $N(0) < \frac{A}{\mu}$, then

Hence, for all $t > 0$,

That is, all orbits of system (2.2) with initial conditions $S(0)\geq 0, E(0)\geq 0, I(0)\geq 0, Q(0)\geq 0$ in $\mathcal{D}$ remain in $\mathcal{D}$ for all $t > 0$. Thus, the region $\mathcal{D}$ is positively-invariant. Furthermore, if $N(0) > \frac{A}{\mu}$, then either $N(t)$ approaches $\frac{A}{\mu}$ as $t\rightarrow \infty$ or the solution enters $\mathcal{D}$ in finite time. Hence, the region $\mathcal{D}$ attracts all solutions in $\mathbb{R}_+^4$.

Throughout this paper, we shall consider the dynamical behaviors of system (2.2) on the region $\mathcal{D}$.

4.   Equilibria

Setting the right-hand sides of system (2.2) to 0, we can get that system (2.2) has only one disease-free equilibrium, denoted by $P_0(\frac{A}{\mu}, 0, 0, 0)$.

Next, we shall calculate the basic reproduction number, denoted by $\mathcal{R}_{0}$, of system (2.2) by applying the next generation matrix method in ^[24] offered by van den Driessche et al.

Let $x = (I, E)^\top$. We can re-express system (2.2) as follows

where

We can obtain

Then the next generation matrix for system (2.2) is

The spectral radius of matrix $FV^{-1}$, denoted by $\rho(FV^{-1})$, equals to $\frac{A\sigma\beta}{\mu(\xi+\gamma+\mu+d)(\mu+\sigma+\rho+a)}$. According to Theorem 2 in literature ^[24], the basic reproductive rate of system (2.2) is spectral radius $\rho(FV^{-1})$, i.e.,

From ^[24], the local stability of the disease free equilibrium $P_0(\frac{A}{\mu}, 0, 0, 0)$ for system (2.2) can directly established. The result is listed as following.

Theorem 4.1. The disease-free equilibrium $P_0(\frac{A}{\mu}, 0, 0, 0)$ of system (2.2) is locally asymptotically stable when $\mathcal{R}_0 < 1$; and $P_0(\frac{A}{\mu}, 0, 0, 0)$ is unstable when $\mathcal{R}_0 > 1$.

In the following, we shall discuss the existence of the endemic equilibrium. Let $P^*(S^*, E^*, I^*, Q^*)$ be an arbitrary endemic equilibrium of system (2.2). Setting the right-hand sides of (2.2) to 0, we can get $E^* = \frac{\xi+\gamma+\mu+d}{\sigma}I^*$, $S^* = \frac{A(1+kE^*)+aE^*}{\beta I^*+\mu}$, $Q^* = \frac{\xi}{\phi+\mu}I^*$. Here $I^*$ satisfies the following quadratic equation

where

Denote

Solving the equation $\Delta = 0$, we can obtain that $\mathcal{R}_0 = \mathcal{R}^*$, where

The following equivalent relations are true:

Hence, for the existence of equilibria of system (2.2), the following conclusions are correct.

Theorem 4.2. System (2.2) always exists a disease free equilibrium $P_0$ and

(i) if $\mathcal{R}^* < \mathcal{R}_0$ or $\mathcal{R}^* = \mathcal{R}_0$ or $\mathcal{R}^* < \mathcal{R}_0 < 1$ and $A_2 > 0$, system (2.2) has no endemic equilibrium;

(ii) if $\mathcal{R}^* < \mathcal{R}_0 = 1$ and $A_2 < 0$ or $\mathcal{R}_0 > 1$ or $\mathcal{R}^* = \mathcal{R}_0 < 1$, system (2.2) has only one endemic equilibrium $P^*(S^*, E^*, I^*, Q^*)$;

(iii) if $\mathcal{R}^* < \mathcal{R}_0 < 1$ and $A_2 < 0$, system (2.2) has two unequal endemic equilibrium points denoted by $P_*(S_*, E_*, I_*, Q_*)$ and $P^*(S^*, E^*, I^*, Q^*)$, where

5.   Stability of disease-free equilibrium $P_0$

In this section we shall investigate the locally and globally asymptotical stability of disease-free equilibrium $P_0(\frac{A}{\mu}, 0, 0, 0)$ and endemic equilibrium $P^*$.

Theorem 5.1. The disease free equilibrium $P_0(\frac{A}{\mu}, 0, 0, 0)$ is globally asymptotically stable when $\mathcal{R}_1 < 1$, where $\mathcal{R}_1 = \frac{A\sigma\beta}{\mu(\xi+\gamma+\mu+d)(\sigma+\rho+\mu+\frac{\mu a}{\mu+kA})}$.

Proof. Consider Lyapunov function $V(t) = a_1E(t)+a_2I(t)+a_3Q(t)$, where $a_1$, $a_2$, $a_3$ are undetermined non-negative real numbers. Then the derivative of $V(t)$ along the solution curves of (2.2) has the following form,

Now we select the coefficients $a_1$, $a_2$ and $a_3$, with the zero coefficients of $E$ and $I$. Hence we obtain

Substituting the values of $a_1$, $a_2$ and $a_3$ to $V(t)$, the derivative of $V(t)$ can be expressed as

Clearly, $\frac{dV(t)}{dt}\leq 0$ when $\mathcal{R}_1 = \frac{A\sigma\beta}{\mu(\xi+r+\mu+d)(\sigma+\rho+\mu+\frac{\mu a}{\mu+kA})} < 1$. Furthermore, $\frac{dV(t)}{dt} = 0$ if and only if $E = I = Q = 0$.

Next, we shall examine the possibility of backward bifurcation for system (2.2). To do this, we firstly introduce the approach presented by Castillo-Chaves et al. (see ^[25]), which is based on the use of the general center manifold theory ^[26]. Considering a general system of ODEs with a parameter $\varsigma$:

Assume that $\chi = 0$ is an equilibrium for system (5.1) for all values of the parameter $\varsigma$, that is $\digamma(0, \varsigma)\equiv 0, \, \, \mbox{for}\, \, \mbox{all}\, \, \varsigma.$ Let $\mathcal{Q} = D_\chi \digamma(0, 0) = (\frac{ \partial \digamma_i}{ \partial \chi_i}(0, 0))$ be the Jacobian matrix of $\digamma(\chi, \varsigma)$ at point $(0, 0)$.

Lemma 5.1. (Castillo-Chavez and Song ^[25]) Assume:

(A1) 0 is a simple eigenvalue of $\mathcal{Q}$ and all other eigenvalues of $\mathcal{Q}$ have negative real parts;

(A2) Matrix $\mathcal{Q}$ has a (non-negative) right eigenvector ${\bf w} = (w_1, w_2, \cdots, w_n)^\top$ and a left eigenvector ${\bf v} = (v_1, v_2, \cdots, v_n)$ corresponding to the zero eigenvalue.

Let $\digamma_k$ denotes the $k^{th}$ component of $\digamma$ and,

Then the local dynamics of system (5.1) around $\chi = 0$ are totally determined by ${\bf a}$ and ${\bf b}$.

1) ${\bf a} > 0$, ${\bf b} > 0$. When $\varsigma < 0$, with $|\varsigma| \ll 1$, $\chi = 0$ is locally asymptotically stable and there exists a positive unstable equilibrium; when $0 < \varsigma\ll 1$, $\chi = 0$ is unstable and there exists a negative and locally asymptotically stable equilibrium;

2) ${\bf a} < 0$, ${\bf b} < 0$. When $\varsigma < 0$, with $|\varsigma| \ll 1$, $\chi = 0$ is unstable; when $0 < \varsigma\ll 1$, $\chi = 0$ is locally asymptotically stable and there exists a positive unstable equilibrium;

3) ${\bf a} > 0$, ${\bf b} < 0$. When $\varsigma < 0$, with $|\varsigma| \ll1$, $\chi = 0$ is unstable and there exists a locally asymptotically stable negative equilibrium; when $0 < \varsigma\ll1$, $\chi = 0$ is stable and a positive unstable equilibrium appears;

4) ${\bf a} < 0$, ${\bf b} > 0$. When $\varsigma$ changes from negative to positive, $\chi = 0$ changes its stability from stable to unstable. Correspondently, a negative unstable equilibrium becomes positive and locally asymptotically stable.

Remark 5.1. The requirement that ${\bf w}$ is non-negative is unnecessary (see ^[25]).

Introducing $S = \chi_1, \; E = \chi_2, \; I = \chi_3$, $Q = \chi_3$, we rewrite system (2.2) as

We will apply Lemma 5.1 to show that system (5.2) may exhibit a backward bifurcation when $\mathcal{R}_0 = 1$. We consider the parameter $\beta$ as bifurcation parameter. Corresponding to $\mathcal{R}_0 = 1$, we can get $\beta = \beta^* = \frac{\mu(\xi+r+\mu+d)(\mu+\sigma+\rho+a)}{A\sigma }.$

The Jacobi matrix of system (5.2) is

And the eigenvalues of $J(P_0, \beta^*)$ are given by $\lambda_1 = -\mu$, $\lambda_2 = 0$, $\lambda_3 = -(\mu+d+\gamma+\xi)$ and $\lambda_4 = -(\phi+\mu)$.

Obviously, the matrix $J(P_0, \beta^*)$ has a simple zero eigenvalue $\lambda_2 = 0$. And the other eigenvalues of $J(P_0, \beta^*)$ are negative real numbers. Therefore, we can use the center manifold theory to discuss the dynamics of system (2.2) when $\mathcal{R}_0 = 1$. Hence, the disease free equilibrium $P_0$ is a nonhyperbolic equilibrium.when $\beta = \beta^*$ (or equivalently when $\mathcal{R}_0 = 1$). Therefore, the assumption (A1) of Lemma 5.1 is then verified.

Now we will calculate a right eigenvector of the matrix $J(P_0, \beta^*)$ associated with the zero eigenvalue $\lambda_2 = 0$, denoted by ${\bf w} = (w_1, w_2, w_3, w_4)^\top$. It is found by

Thus, we can get

This implies $w_1 = \frac{\phi+\mu}{\mu}[\frac{a(\mu+d+\gamma+\xi)}{\sigma}-\frac{\beta A}{\mu}],$ $w_2 = \frac{(\phi+\mu)(\mu+d+\gamma+\xi)}{\sigma}, \; w_3 = \phi+\mu, \; w_4 = \xi.$ Therefore, the right eigenvector is

Furthermore, the left eigenvector ${\bf v} = (v_1, v_2, v_3, v_4)$ of the matrix $J(P_0, \beta^*)$ which satisfies ${\bf v}\cdot {\bf w} = 1$ is given by

Then, the left eigenvector v turns out to be

Calculating all the partial derivatives of $\digamma_i \; (i = 1, 2, 3, 4)$ with respect to $\chi_i \; (i = 1, 2, 3, 4)$ and $\beta$ at the disease-free equilibrium $P_0(\frac{A}{\mu}, 0, 0, 0)$, we get

and all the other second-order partial derivatives are equal to 0.

Thus, we can calculate the coefficients ${\bf a}$ and ${\bf b}$ defined in Lemma 5.1, i.e.,

Considering system (5.2) and taking into account of ${\bf a}$ and ${\bf b}$ only the nonzero derivatives of the terms $\frac{\partial^2\digamma_k}{\partial \chi_i\partial \chi_j}(P_0, \beta^*)$ and $\frac{\partial^2\digamma_k}{\partial \chi_i\partial \beta}(P_0, \beta^*)$, it follows that

and

From (5.3) and (5.4), we obtain

and

Obviously, the coefficient ${\bf b}$ is always positive. According to Lemma 5.1 we can conclude that the sign of ${\bf a}$ determines the local dynamics of the disease-free equilibrium $P_0$ when $\beta = \beta^*$.

Define $\mathcal{R}^*_1 = \frac{a(\mu+d+\gamma+\xi)}{\beta A \sigma}(\frac{1}{\mu}+\frac{\mu^2k}{\beta \sigma})$. Note that ${\bf a} < 0$ if $\mathcal{R}^*_1 < 1$, and ${\bf a} > 0$ if $\mathcal{R}^*_1 > 1$. Hence, from Lemma 5.1, we have the following results.

Theorem 5.2. (1) Assume that the basic reproductive rate $\mathcal{R}_0$ equals to 1. System (2.2) exhibits a backward bifurcation if $\mathcal{R}^*_1 < 1$. Otherwise, system (2.2) exhibits a forward bifurcation if $\mathcal{R}^*_1 > 1$.

(2) The endemic equilibrium is locally asymptotically stable when the basic reproductive rate $\mathcal{R}_0 > 1$ and close to one.

6.   Global stability of endemic equilibria

In the following, we will utilize the geometric approach to observe the global stability of the endemic equilibrium for system (2.2). Firstly, we present some preliminary results on the geometric approach to global dynamics, one can find them in ^[28,29,30]. Let $\mathcal{B}$ be the Euclidean unit ball in $\mathbb{R}^2$, and let $\partial \mathcal{B}$ and $\bar{\mathcal{B}}$ be its boundary and closure, respectively. Denote the set of Lipschitzian functions from $X$ to $Y$ by Lip$(X\rightarrow Y)$. We consider a function $\varphi\in$ Lip$(\bar{\mathcal{B}}\rightarrow \mathcal{D})$ as a simply connected rectifiable surface in $\mathcal{D}\subset \mathbb{R}^n$. A function $\psi\in$ Lip$(\partial \mathcal{B}\rightarrow \mathcal{D})$ is a closed rectifiable curve in $\mathcal{D}$ and is called simple if it is one-to-one. Define $\sum (\psi, \mathcal{D}) = \{\varphi\in \mbox{Lip}(\bar{\mathcal{B}}\rightarrow \mathcal{D})|\varphi_{\partial \mathcal{B}} = \psi\}$. If $\mathcal{D}$ is an open, simply connected set, then $\sum (\psi, \mathcal{D})$ is nonempty for each simple closed rectifiable curve $\psi$ in $\mathcal{D}$. Let $||\cdot||$ be a norm on $\mathbb{R}^ \left(\begin{array}{l} n \\ 2 \end{array}\right)$. A function $\mathcal{S}$ on surface in $\mathcal{D}$ is defined as follows:

where $u = (u_1, u_2), \; u \longmapsto \varphi(u)$ is Lipschitzian on $\bar{\mathcal{B}}$, $P$ is an $\left(\begin{array}{l} n \\ 2 \end{array}\right) \times \left(\begin{array}{l} n \\ 2 \end{array}\right)$ matrix such that $||P^{-1}||$ is bounded on $\varphi(\bar{\mathcal{B}})$, and the wedge product $\frac{\partial \varphi}{\partial u_1} \wedge \frac{\partial \varphi}{\partial u_2}$ is a vector in $\mathbb{R}^ \left(\begin{array}{l} n \\ 2 \end{array}\right)$. The following lemma develops on the results in ^[28] and ^[30].

Lemma 6.1. Suppose that $\psi$ is an arbitrary simple, closed and rectitiable curve in $\mathbb{R}^n$. Then there exists $\delta > 0$ such that $\mathcal{S}\psi \ge \delta$ for all $\varphi\in \sum (\psi, \mathbb{R}^n)$.

Let $x\longmapsto f(x)\in \mathbb{R}^n$ be a $C^1$ function for $x$ in a set $\mathcal{D} \subset \mathbb{R}^n$. We consider the autonomous equation in $\mathbb{R}^n$

For any surface $\varphi$, the new surface $\varphi_t$ is defined by $\varphi_t(u) = x(t, \varphi(u))$. If $\varphi_t(u)$ is reviewed as a function of $u$, $\varphi_t(u)$ is a time $t$ map determined by system (6.2). If $\varphi_t(u)$ is viewed as a function of $t$, $\varphi_t(u)$ is the solution of (6.2) passing through the initial point $(0, \varphi(u))$. The right-hand derivative of $\mathcal{S}\varphi_t$, denoted $D_+\mathcal{S}\varphi_t$ is defined by

where the matrix $Q = P_fP^{-1}+P \frac{\partial f}{\partial x}^{[2]}P^{-1}$. Here $P_f$ is the directional derivative of $P$ in the direction of the vector field $f$, $\frac{\partial f}{\partial x}^{[2]}$ is the second additive compound matrix (we can find its definition in ^[31]) of $\frac{\partial f}{\partial x}$, and $\varrho = P\cdot(\frac{\partial \varphi}{\partial u_1} \wedge \frac{\partial \varphi}{\partial u_2})$ is a solution to the differential equation

Then, the right-hand derivative $D_+\mathcal{S}\varphi_t$ is expressed as

At a general point $P(S, E, I, Q)$, the Jacobian matrix is given by

The second additive compound matrix of $\frac{\partial f}{\partial x}$ is the $6\times 6$ matrix given by

where

Let

Then

The derivative of $P$ in the direction of the vector field $f$ is

Then

and

where

We consider the following norm introduced in ^[34,35] for $\varrho = (\varrho_1, \varrho_2, \varrho_3, \varrho_4, \varrho_5, \varrho_6)^\top\in \mathbb{R}_6$, $\|\varrho\| = \max\{\mathcal{U}_1, \mathcal{U}_2\}$, where $\mathcal{U}_1(\varrho_1, \varrho_2, \varrho_3)$ has the following form

and $\mathcal{U}_2(\varrho_4, \varrho_5, \varrho_6)$ has the following form

Furthermore we use the following relations

and

Lemma 6.2. There exists a positive constant $\eta$, such that $D_+\|\varrho\|\leq -\eta \|\varrho\|$ for all $\varrho\in \mathbb{R}^4$ and all $S, E, I, Q > 0$, provided that the following inequation

holds for some $\eta$. Here $D_+\|\varrho\|$ the right-hand derivative of $\|\varrho\|$, $\varrho$ is the solution of $\frac{d\varrho}{dt} = Q\varrho$.

Proof. We show the existence of some $\eta > 0$ such that

for all $\varrho\in \mathbb{R}^4,$ where $\varrho$ is a solution of Eq (6.4). By linearity, if the above inequality is true for some $\varrho$, then it also holds for $-\varrho$. Based on the different orthants and the definition of the norm $\|\cdot\|$ in (6.6) and (6.7) within each orthant, the full calculation to demonstrate the proof involves 16 separate cases and subcases.

Case 1. $\mathcal{U}_1 > \mathcal{U}_2$, $\varrho_1, \varrho_2, \varrho_3 > 0$, and $|\varrho_1| > |\varrho_2|+|\varrho_3|.$ Then

so that

Case 2. $\mathcal{U}_1 > \mathcal{U}_2$, $\varrho_1, \varrho_2, \varrho_3 > 0$, and $|\varrho_1| < |\varrho_2|+|\varrho_3|.$ Then

so that

Case 3. $\mathcal{U}_1 > \mathcal{U}_2$, $\varrho_1 < 0$, $\varrho_2, \varrho_3 > 0$, and $|\varrho_1| > |\varrho_2|.$ Then

so that

Case 4. $\mathcal{U}_1 > \mathcal{U}_2$, $\varrho_1 < 0$, $\varrho_2, \varrho_3 > 0$, and $|\varrho_1| < |\varrho_2|.$ Then

so that

Case 5. $\mathcal{U}_1 > \mathcal{U}_2$, $\varrho_1$, $\varrho_2 > 0$, $\varrho_3 < 0$, and $|\varrho_2| > |\varrho_1|+|\varrho_3|.$ Then

so that

Case 6. $\mathcal{U}_1 > \mathcal{U}_2$, $\varrho_1$, $\varrho_2 > 0$, $\varrho_3 < 0$, and $|\varrho_2| < |\varrho_1|+|\varrho_3|.$ Then

so that

Case 7. $\mathcal{U}_1 > \mathcal{U}_2$, $\varrho_2$, $\varrho_3 > 0$, $\varrho_2 < 0$, and $|\varrho_1| = \max\{|\varrho_2|, |\varrho_3|\}.$ Then

so that

Case 8. $\mathcal{U}_1 > \mathcal{U}_2$, $\varrho_2$, $\varrho_3 > 0$, $\varrho_2 < 0$, and $|\varrho_2| = \max\{|\varrho_1|, |\varrho_3|\}.$ Then

so that

Case 9. $\mathcal{U}_1 > \mathcal{U}_2$, $\varrho_1$, $\varrho_3 > 0$, $\varrho_2 < 0$, and $|\varrho_3| > \max\{|\varrho_1|, |\varrho_2|\}.$ Then

so that

Case 10. $\mathcal{U}_1 < \mathcal{U}_2$, $\varrho_4$, $\varrho_5$, $\varrho_6 > 0$. Then

so that

Using $|\varrho_2+\varrho_3| < \mathcal{U}_1 < |\varrho_4|+|\varrho_5|+|\varrho_6|$, it follows

Case 11. $\mathcal{U}_1 < \mathcal{U}_2$, $\varrho_4$, $\varrho_5 > 0$, $\varrho_6 < 0$ and $|\varrho_5| > |\varrho_6|$. Then

so that

Using $|\varrho_2+\varrho_3| < \mathcal{U}_1 < |\varrho_4|+|\varrho_5|$, it follows

Case 12. $\mathcal{U}_1 < \mathcal{U}_2$, $\varrho_4$, $\varrho_5 > 0$, $\varrho_6 < 0$ and $|\varrho_5| < |\varrho_6|$. Then

so that

Using $|\varrho_2| < \mathcal{U}_1 < |\varrho_4|+|\varrho_5|$, it follows

Case 13. $\mathcal{U}_1 < \mathcal{U}_2$, $\varrho_4$, $\varrho_6 > 0$, $\varrho_5 < 0$ and $|\varrho_5| > |\varrho_4|+|\varrho_6|$. Then

so that

Using $|\varrho_3| < \mathcal{U}_1 < |\varrho_5|$, it follows

Case 14. $\mathcal{U}_1 < \mathcal{U}_2$, $\varrho_4$, $\varrho_6 > 0$, $\varrho_5 < 0$ and $|\varrho_5| < |\varrho_4|+|\varrho_6|$. Then

so that

Using $|\varrho_2| < \mathcal{U}_1 < |\varrho_4|+|\varrho_6|$, it follows

Case 15. $\mathcal{U}_1 < \mathcal{U}_2$, $\varrho_5$, $\varrho_6 > 0$, $\varrho_4 < 0$ and $|\varrho_5| < |\varrho_4|$. Then

so that

Using $|\varrho_2| < \mathcal{U}_1 < |\varrho_4|+|\varrho_6|$, it follows

Case 16. $\mathcal{U}_1 < \mathcal{U}_2$, $\varrho_5$, $\varrho_6 > 0$, $\varrho_4 < 0$ and $|\varrho_5| > |\varrho_4|$. Then

so that

Using $|\varrho_3| < \mathcal{U}_1 < |\varrho_5|+|\varrho_6|$, it follows

From Case 1 to Case 16 and the condition of Lemma 6.2, we can obtain

for all $\varrho\in \mathbb{R}^4.$ From ^[32,33] we know that the geometric approach can be applied to prove the globally asymptotic stability when the epidemic model has a unique endemic equilibrium. In this situation, there exists a compact absorbing set $\mathcal{D}$, and a surface remains in $\mathcal{D}$ for all time. From Section 4, we find that system (2.2) will possibly exhibit bistability. For this case, system (2.2) does not exist an absorbing set. Hence, we shall consider the following sequence of surface $\{\varphi^k\}$ in the following lemma.

Lemma 6.3. Let $\psi$ be a simple closed curve in $\mathcal{D}$. There exist a positive $\epsilon$ and a sequence of surface $\{\psi^k\}$ that minimize $\mathcal{S}$ given by (6.1) relative to $\sum (\psi, \mathcal{D})$ such that $\{\psi^k_t\subset \mathcal{D}\}$ for all $k = 1, 2, 3, \cdots$ and all $t\in [0, \epsilon]$.

Proof. Let $\zeta = \frac{1}{2}\min\{E, I:(S, E, I, Q)\in \psi\}$. It is easy to see $\zeta > 0$. From the second and third equations of (2.2), we can get the inequality $\frac{dE(t)}{dt}\geq -\sigma E(t)-\rho E(t)-\mu E(t)-\frac{aE(t)}{1+kE(t)}$, $\frac{dI(t)}{dt}\geq -(\xi+ \gamma+\mu+d) I$ in $\mathcal{D}$. Hence we can conclude that there is an $\epsilon > 0$ such that, if a solution satisfied $E(0)\geq \zeta, \; I(0)\geq \zeta,$ then it remains in $\mathcal{D}$ for $t\in [0, \epsilon]$. Thus we only require to prove there exists a sequence $\{\psi^k\}$, such that it minimizes $\mathcal{S}$ relative to $\Sigma (\psi, \tilde{\mathcal{D}})$, where $\tilde{\mathcal{D}} = \{(S, E, I, Q)\in\mathcal{\mathcal{D}} :E\geq \zeta, I\geq \zeta\}$.

Let $\varphi(u) = (S(u), E(u), I(u), Q(u))\in \Sigma (\psi, \mathcal{D})$. Define a new surface $\tilde{\varphi} (u) = (\tilde{S}(u), \tilde{E}(u), \; \tilde{I}(u), \tilde{Q}(u))$ by

From the above definition of $\tilde{\varphi}(u)$, it is not difficult to know $\tilde{\varphi}(u)\in\Sigma (\psi, \tilde{\mathcal{D}})$. Also $\mathcal{S} \tilde{\varphi} = \int_{\bar{\mathcal{B}}}|\frac{1}{\tilde{E}} (\frac{\partial \tilde{\varphi}}{\partial u_1} \wedge \frac{\partial \tilde{\varphi}}{\partial u_2})|du$ and $\mathcal{S} \varphi = \int_{\bar{\mathcal{B}}}|\frac{1}{E} (\frac{\partial \varphi}{\partial u_1} \wedge \frac{\partial \varphi}{\partial u_2})|du$. In the following, we will prove $\mathcal{S} \tilde{\varphi}\leq \mathcal{S} \varphi.$

According to the definition of wedge product, we obtain that

is a vector in $\mathbb{R}^6$ for each $u\in \mathcal{B}$. Denote $\frac{\partial \tilde{\varphi}}{\partial u_1} \wedge \frac{\partial \tilde{\varphi}}{\partial u_2} = (\tilde{x}_1, \tilde{x}_2, \tilde{x}_3, \tilde{x}_4, \tilde{x}_5, \tilde{x}_6)^\top$ and $\frac{\partial \varphi}{\partial u_1} \wedge \frac{\partial \varphi}{\partial u_2} = (x_1, x_2, x_3, x_4, x_5, x_6)^\top$. We will prove $|\tilde{x}_i|\leq |x_i| \; (i = 1, 2, 3, 4, 5, 6)$ in the following seven cases.

Case 1. If $E(u)\geq \zeta$, $I(u)\geq \zeta$, then $\tilde{\varphi} = \varphi$ and therefore, $|\tilde{x}_i|\leq |x_i| \; (i = 1, 2, 3, 4, 5, 6)$.

Case 2. If $E(u) < \zeta$, $I(u)\geq \zeta$ and $S(u)+E(u)+\zeta+Q(u)\leq \frac{A}{\mu},$ then $\tilde{\varphi}(u) = (S(u)+\zeta+I(u)+Q(u))$. Hence, we can obtain

Therefore, it follows $\tilde{x}_i = x_i \; (i = 1, 3, 5)$ and $\tilde{x}_i = 0 \; (i = 2, 4, 6)$. Hence $|\tilde{x}_i|\leq |x_i| \; (i = 1, 2, 3, 4, 5, 6)$.

Case 3. If $E(u) < \zeta$, $I(u)\geq \zeta$ and $S(u)+E(u)+\zeta+Q(u) > \frac{A}{\mu},$ then $\tilde{\varphi}(u) = (\frac{S}{S+Q}(\frac{A}{\mu}-2\zeta), \zeta, \zeta, \frac{S}{S+Q}(\frac{A}{\mu}-2\zeta))$. Therefore, $\frac{\partial \tilde{\varphi}}{\partial u_1} \wedge \frac{\partial \tilde{\varphi}}{\partial u_2} = 0$. Thus,

for $j = 1, 2$. Therefore, $\frac{\partial \tilde{\varphi}}{\partial u_1}$ and $\frac{\partial \tilde{\varphi}}{\partial u_2}$ have linear dependent. Hence, $\frac{\partial \tilde{\varphi}}{\partial u_1} \wedge \frac{\partial \tilde{\varphi}}{\partial u_2} = 0.$ Thus $\tilde{x}_i = 0 \; (i = 1, 2, 3, 4, 5, 6)$. Therefore $|\tilde{x}_i|\leq |x_i| \; (i = 1, 2, 3, 4, 5, 6)$.

Case 4. If $E(u)\geq \zeta$, $I(u) < \zeta$ and $S(u)+E(u)+\zeta+Q(u)\leq\frac{A}{\mu},$ then $\varphi(u) = (S(u), I(u), \zeta, Q(u))$. Hence,

almost everywhere. Therefore, $\tilde{x}_i = x_i \; (i = 1, 2, 4)$, $\tilde{x}_i = 0 \; (i = 3, 5, 6)$. Therefore $|\tilde{x}_i|\leq |x_i| \; (i = 1, 2, 3, 4, 5, 6)$.

Case 5. If $E(u)\geq \zeta$, $I(u) > \zeta$ and $S(u)+E(u)+\xi+Q(u) > \frac{A}{\mu},$ then $\tilde{\varphi}(u) = (\frac{S}{S+Q}(\frac{A}{\mu}-2\zeta), \zeta, \zeta, \frac{S}{S+Q}(\frac{A}{\mu}-2\zeta))$. Thus,

for $j = 1, 2$. Therefore, $\frac{\partial \tilde{\varphi}}{\partial u_1}$ and $\frac{\partial \tilde{\varphi}}{\partial u_2}$ have linear dependent. Hence, $\frac{\partial \tilde{\varphi}}{\partial u_1} \wedge \frac{\partial \tilde{\varphi}}{\partial u_2} = 0.$ Thus $\tilde{x}_i = 0 \; (i = 1, 2, 3, 4, 5, 6)$. Therefore $|\tilde{x}_i|\leq |x_i| \; (i = 1, 2, 3, 4, 5, 6)$.

Case 6. If $E(u) < \zeta$, $I(u) < \zeta$ and $S(u)+2\zeta+Q(u)\leq \frac{A}{\mu},$ then $\tilde{\varphi}(u) = (S(u), \zeta, \zeta, Q(u))$, thus,

Therefore, $\tilde{x}_i = x_i \; (i = 1)$ and $\tilde{x}_i = 0 \; (i = 2, 3, 4, 5, 6)$. Hence $|\tilde{x}_i|\leq |x_i| \; (i = 1, 2, 3, 4, 5, 6)$.

Case 7. If $E(u) < \zeta$, $I(u) < \zeta$ and $S(u)+E(u)+2\zeta > \frac{A}{\mu},$ then $\tilde{\varphi}(u) = (\frac{S}{S+Q}(\frac{A}{\mu}-2\zeta), \zeta, \zeta, \frac{S}{S+Q}(\frac{A}{\mu}-2\zeta))$. Thus,

for $j = 1, 2$. Therefore, $\frac{\partial \tilde{\varphi}}{\partial u_1}$ and $\frac{\partial \tilde{\varphi}}{\partial u_2}$ have linear dependent. Hence, $\frac{\partial \tilde{\varphi}}{\partial u_1} \wedge \frac{\partial \tilde{\varphi}}{\partial u_2} = 0.$ Thus $\tilde{x}_i = 0 \; (i = 1, 2, 3, 4, 5, 6)$. Therefore $|\tilde{x}_i|\leq |x_i| \; (i = 1, 2, 3, 4, 5, 6)$.

We also note that $\tilde{E}(u) = \max\{E(u), \zeta\}$ and $\tilde{I}(u) = \max\{I(u), \zeta\}$. Thus $\frac{1}{\tilde{E}(u)}\leq \frac{1}{E(u)}$, $\frac{1}{\tilde{I}(u)}\leq \frac{1}{I(u)}$. Let

by comparing the vector $P\cdot (\frac{\partial \varphi}{\partial u_1} \wedge \frac{\partial \varphi}{\partial u_2})$ and the vector $\tilde{P}\cdot(\frac{\partial \tilde{\varphi}}{\partial u_1} \wedge \frac{\partial \tilde{\varphi}}{\partial u_2})$, we have $|\frac{1}{\tilde{I}}\tilde{x}_i|\leq |\frac{1}{I}x_i| \; (i = 1, 2, 4)$ and $|\frac{1}{\tilde{E}}\tilde{x}_i|\leq |\frac{1}{E}x_i| \; (i = 3, 5, 6)$ for their corresponding component.

From above seven cases, we can obtain $\|\frac{\partial \tilde{\varphi}}{\partial u_1} \wedge \frac{\partial \tilde{\varphi}}{\partial u_2} \|\leq \|\frac{\partial \varphi}{\partial u_1} \wedge \frac{\partial \varphi}{\partial u_2}\|$ and $\mathcal{S} \tilde{\varphi}\leq \mathcal{S} \varphi$.

From Lemma 6.2, we can choose $\varpi = \inf\{\mathcal{S} \varphi:\varphi\in \Sigma(\psi, \mathcal{D}) \}$. Assume that $\{\varphi^k\}$ is a sequence of surfaces that minimizes $\mathcal{S}$ relative to $\Sigma(\psi, \mathcal{D})$, then $\lim\limits_{t\rightarrow \infty}\mathcal{S} \varphi^k = \varpi$. Let $\{\tilde{\varphi}^k\}$ be a sequence of surfaces that minimizes $\mathcal{S}$ relative to $\Sigma(\psi, \tilde{\mathcal{D}})$ defined by the above construction. Hence, for each $k$, $\mathcal{S}\tilde{\varphi}^k \leq \mathcal{S} \varphi^k$. On the other hand, since $\{\mathcal{S}\tilde{\varphi}^k\}$ is bounded, we assume that $\{\mathcal{S}\tilde{\varphi}^k\}$ is convergent without generality. From $\mathcal{S}\tilde{\varphi}^k \leq \mathcal{S} \varphi^k$ (for each $k$), we can obtain that $\lim\limits_{k\rightarrow \infty}\mathcal{S}\tilde{\varphi}^k \leq \varpi.$ From $\tilde{\varphi}^k\in \Sigma(\psi, \mathcal{D})$, we get $\tilde{\varphi}^k\in \Sigma(\psi, \tilde{\mathcal{D}})$. Then, for each $k$, $\lim\limits_{k\rightarrow \infty}\mathcal{S}\tilde{\varphi}^k \geq \varpi.$ Therefore, $\lim\limits_{k\rightarrow \infty}\mathcal{S}\tilde{\varphi}^k = \varpi.$ From $\mathcal{S}\tilde{\varphi}^k \leq \mathcal{S} \varphi^k$ (for each $k$), we can obtain

From $\tilde{\varphi}\in \Sigma(\psi, \mathcal{D})$, we have $\inf\{\mathcal{S} \tilde{\varphi}:\tilde{\varphi}\in\Sigma(\psi, \tilde{\mathcal{D}}) \}\geq \varpi$, which implies $\inf\{\mathcal{S} \tilde{\varphi}:\tilde{\varphi}\in\Sigma(\psi, \tilde{\mathcal{D}}) \} = \varpi$. At last, we can show that $\lim\limits_{k\rightarrow \infty}\mathcal{S}\tilde{\varphi}^k = \varpi = \inf\{\mathcal{S} \tilde{\varphi}:\tilde{\varphi}\in\Sigma(\psi, \tilde{\mathcal{D}}) \}$. It follows that $\tilde{\varphi}^k$ minimizes $\mathcal{S}$ relative to $\Sigma(\psi, \tilde{\mathcal{D}})$.

From Lemmas 6.2 and 6.3, we can obtain the following theorem.

Theorem 6.1. If the inequality (6.8) holds true, then any $\omega-$limit point of system (2.2) in the interior of $\mathcal{D}$ is an equilibrium, and so each positive semitrajectory of system (2.2) in $\bar{\mathcal{D}}$ limits to a single equilibrium.

The proof of Theorem 6.1 is similar to the proof of Corollary 5.4 in ^[27], we omit it.

From Theorem 6.1, we have the following theorem.

Theorem 6.2. Suppose that inequality (6.8) is satisfies, then

(1) when system (2.2) has only one disease-free equilibrium, then all solutions of system (2.2) limit to $P_0$;

(2) when $\mathcal{R}_0 > 1$, then all solutions of system (2.2) converge to the the endemic equilibrium $P^*$;

(3) when system (2.2) has two endemic equilibrium point, i.e., $\mathcal{R}^*_0 < \mathcal{R}_0 < 1$ and $A_2 < 0$, then solutions of system (2.2) either go to the disease-free equilibrium $P_0$, or tend to the upper equilibrium $P^*$.

7.   Numerical simulation

In this section, we present some numerical simulations by using the method of Runge-Kutta through Matlab 2018a software to corroborate the theoretical results and find the complex dynamics of system (2.2).

Case i.. Let $A = 2.5, \beta = 0.003, \mu = 0.006, \; a = 10.98, k = 15, \sigma = 0.84, \xi = 0.3,$$\gamma = 0.15, d = 0.009, \rho = 1.16, \phi = 0.32.$ We can get $A_2 = 0.00506 > 0$, $\mathcal{R}^* = 0.9973$, $\mathcal{R}_0 = 0.1739$ and $\mathcal{R}_0 < \mathcal{R}^* < 1$. It is easy to see that the conditions of Case (i) of Theorem 4.2 are satisfied. Hence, system (2.2) has only one disease free equilibrium $P_0(666.6667, 0, 0, 0)$. Furthermore, the condition of Theorem 4.1 is satisfied, $P_0(666.6667, 0, 0, 0)$ is locally asymptotically stable (See Figure 3).

Case ii. Set $A = 13, \beta = 0.003, \mu = 0.006, \; a = 10.98, k = 15, \sigma = 0.84, \xi = 0.3,$$\gamma = 0.15, d = 0.009, \rho = 1.16, \phi = 0.32.$ We can get $\mathcal{R}^* = -13.7347$, $\mathcal{R}_0 = 0.9042$ and $\mathcal{R}^* < \mathcal{R}_0 < 1$. The conditions of Case (iii) of Theorem 4.2 are satisfied. System (2.2) has a disease free equilibrium $P_0(2666.6667, 0, 0, 0)$, two endemic equilibria $P_*(2163.778095, 0.0086, \; 0.0156, 0.0144)$ and $P^*(395.3322, 5.2981,$$9.5708, 8.8075)$. $P_*$ is unstable and $P^*$ is locally asymptotically stable (See Figure 4).

Case iii. Suppose $A = 16, \beta = 0.003, \mu = 0.006, \; a = 10.98, k = 15, \sigma = 0.84, \xi = 0.3,$$\gamma = 0.15, d = 0.009, \rho = 1.16, \phi = 0.32.$ We can get $\mathcal{R}_0 = 1.1129 > 1$. The conditions of Case (ii) of Theorem 4.2 are satisfied. System (2.2) has a disease free equilibrium $P_0(2666.6667, 0, 0, 0)$ and an endemic equilibrium $P^*(389.7962, 6.8102, 12.3023, 11.3211)$. By the case (2) of Theorem 6.2, $P^*$ is locally asymptotically stable (See Figure 5).

Case iv. Set $A = 3.2, \beta = 0.0072, \mu = 0.006, k = 15, \sigma = 0.84, \xi = 0.3, \gamma = 0.15, d = 0.009, \rho = 1.16, \phi = 0.32.$ From Figure 6(a), we can find that Hopf bifurcation may occur around the endemic equilibrium when parameter $a$ changes. Set $A = 3.2, \beta = 0.0072, \mu = 0.006, a = 1.2, \sigma = 0.84, \xi = 0.3, \gamma = 0.15, d = 0.009, \rho = 1.16, \phi = 0.32.$ From Figure 6(b), we can find that system (2.2) exists Hopf bifurcation around the endemic equilibrium when parameter $k$ changes, which shows that nonlinear innate immunity rate can cause system to produce periodic solutions.

8.   Discussion

In this paper we considered an SEIQR epidemic model with nonlinear innate immunity. We found that system (2.2) might exist backward bifurcation under some conditions. The global stability of the disease free and endemic equilibria for system (2.2) was obtained. The global stability of the endemic equilibrium point for a four-dimensional nonlinear ordinary differential equation model is studied by using the method in ^[27,34]. At present, this approach is rarely used in studying the global stability of the endemic equilibrium for the epidemic models. This is a good method to study the global stability of the endemic equilibrium when the models exist backward bifurcation and the Lyapunov function is not well constructed. We find that Hopf bifurcation might occur by numerical simulation. In theory, we can further study the direction of bifurcation and the stability of the bifurcating periodic solutions. We'll leave the work for the future.

When the innate immunity is linear, i.e., $k = 0$, system (2.2) reduces to the following:

By calculating, we get that the basic reproduction number of (8.1) is

Hence $\mathcal{R}_{0} = \mathcal{R}_{0}^l$. That is to say the parameter $k$ does not effect the basic reproduction number. System (8.1) always exists a disease free equilibrium $P_0^l(\frac{A}{\mu}, 0, 0, 0)$. When $\mathcal{R}_{0} > 1$, system (8.1) has a unique endemic equilibrium $P^*_l(S^*_l, E^*_l, I^*_l, Q^*_l)$, where

By constructing the suitable Lyapunov functions we can conclude that (i) when $\mathcal{R}_{0} < 1$, the disease free equilibrium $P_0^l(\frac{A}{\mu}, 0, 0, 0)$ of system (8.1) is globally asymptotically stable; (ii) when $\mathcal{R}_{0} > 1$, the endemic equilibrium $P^*_l(S^*_l, E^*_l, I^*_l, Q^*_l)$ of system (8.1) is globally asymptotically stable. When the innate immunity rate is linear, backward bifurcation and Hopf bifurcation cannot exist. The dynamic properties of system (2.2) are more complex than those of system (8.1). Hence, the nonlinear innate immunity makes more difficult to eliminate the disease. In order to take appropriate prevention and control measures, we should pay attention to the impact of nonlinear innate immunity in disease control.

Some fascinating questions are well worth further investigation. For example, the phenomenon of stochastic disturbances is common in nature. The transmission coefficient is often randomly perturbed as the disease spreads. If we assume that the incidence rate in (2.1) is perturbed by white noise so that $\beta\rightarrow \beta+\nu \dot{B}(t)$, where $B(t)$ is a standard Brownian motion with intensity $\nu$. System (2.1) becomes the following stochastic model

System (8.2) is more reasonable than model (2.1). And the dynamical behaviors of system (8.2) are more complicated than system (2.1). We leave it in the future.

Acknowledgments

This work is sponsored by Natural Science Foundation of Henan (No. 222300420521), Research Project on Teacher Education Curriculum Reform in Henan Province in 2022 (No. 2022-JSJYYB-035) and Nanhu Scholars Program for Young Scholars of XYNU.

Conflict of interest

The authors declare there is no conflicts of interest.