Transmission dynamics of cholera with hyperinfectious and hypoinfectious vibrios: mathematical modelling and control strategies

1.   Introduction

In the modeling of infectious diseases, incidence rate and recovery rate are two important factors that affect the dynamical behavior of the models. The bilinear incidence rate $\beta SI$ are frequently used to describe the epidemic transmission process, where $\beta$ is the transmission rate, $S$ and $I$ represent the numbers of susceptible and infective respectively. But, it is unrealistic when the number of infective individuals is large. To study the cholera epidemic spread in Bari in 1973, Capasso and Serio ^[1] proposed a saturated incidence rate:

where $\beta I$ describes the infection force of the disease, $\frac{1}{1+\alpha I}$ measures the inhibition effect from the behavioral change of the susceptible individuals as the number of infected gradually increases.

In classical epidemic models, the recovery rate is usually assumed to be proportional to the number of infected, which means that medical resources are plentiful for the infectious diseases ^[2]. In fact, when diseases break out, medical resources tend to be scarce ^[3,4,5]. In order to study the influence of hospital beds on infectious diseases, the following nonlinear recovery rate function ^[5] is proposed by Shan and Zhu:

where $\gamma_{0}$ and $\gamma_{1}$ are the minimum and maximum per-capita recovery rates respectively, $b$ denotes the number of hospital beds. An SIR model with (1.2) is studied, and it is found that backward bifurcation, saddle-node bifurcation, Hopf bifurcation and cusp type of Bogdanov-Takens bifurcation of codimension 3 will happen with different values of $b$ ^[5].

Cui et al. studied an SIR epidemic model (system (1.3)) with saturate incidence rate and nonlinear recovery rate ^[6]. It is proved that when the number of hospital beds is small enough, system (1.3) can take place backward bifurcation, saddle-node bifurcation and Hopf bifurcation.

where $S(t), I(t), R(t)$ are the numbers of susceptible, infected, recovered at time $t$ respectively, $\gamma(b, I)$ is defined in (1.2). Besides, $\Lambda$ is the birth rate; the interpretations of $\beta$ and $\alpha$ are the same as in ^[1]; $\mu$ denotes the natural death rate; $\epsilon$ represents the disease-induced death rate.

To describe the heterogeneity of contact, complex networks are therefore incorporated into the epidemic models. Li and Yousef studied a network-based SIR epidemic model with saturated treatment function ^[7]. Huang and Li studied the complex dynamical properties of a network-based SIS epidemic model with saturated treatment function ^[8]. A condition which can determine the direction of bifurcation at $\mathit{\boldsymbol{R_{0}}} = 1$ is derived in ^[7,8]. Wang and Yang ^[9] employed an edge-compartmental approach which can reduce the dimension of a mean-field vector-borne model to research the global dynamics of the system. Wang and Yang ^[10] proposed a degree-edge-mixed SIS model and used a novel geometric method to completely investigate the stability of feasible equilibrium. Based on the SIR model ^[7], Yang et al. adopted the edge-compartmental approach to prove the existence of backward bifurcation and deeply analyzed the local stability of each endemic equilibrium ^[11].

To our knowledge, few people combined saturated incidence rate and nonlinear recovery rate to discuss the dynamics of network-based epidemic models. So the research of this paper is still valuable, the highlights are summarized as follows:

$\bullet$ A new network-based SIR epidemic model is established and then reduced to a degree-edge-mixed model by an edge-compartmental approach;

$\bullet$ A necessary and sufficient condition which determines the existence of backward bifurcation at $\mathit{\boldsymbol{R_{0}}} = 1$ is derived;

$\bullet$ The stability of the feasible endemic equilibrium is proved by a geometric approach.

We know that time delay exists in the epidemic spreading, Jiang and Zhou ^[12] studied the influence of time delay on epidemic spreading under limited resources. They found that time delay can induce first-order, continuous and hybrid phases and a small resources amount can effectively control the spread of infectious diseases if the delay exceeds a threshold. The stability of nonlinear systems with time delay has been studied in the literature ^{[13,14,15,16,17]}. Yang et al. obtained the existence of a positive periodic solution for neutral-type integral differential equation arising in epidemic model with time-varying delay ^[18]. Of course, we do not consider the delay factor in our paper, which will be a worthy of attention and research work.

The structure of this paper is as follows. In Section 2, a new SIR epidemic model is proposed on complex networks. Besides, the dimension of the system is reduced by the edge-compartmental method. In Section 3, the existence and stability of the disease-free equilibrium and endemic equilibrium are studied. The results of numerical simulations are given in Section 4. Lastly, conclusions and discussions are presented in Section 5.

2.   Model description

In this section, we extend model (1.3) to the network. Furthermore, disease-induced death rate is not considered and the susceptible individuals will not be infected once vaccinated. The other parameters are exactly the same as in system (1.3) and we show the interpretation of the parameters in Table 1 for better visualization. All the nodes are classified into $n$ groups and the nodes in the same group have the same degree. Suppose that $S_{k}(t), I_{k}(t), R_{k}(t)$ be the densities of susceptible, infected, recovered with degree $k$ at time $t$ respectively. So, the specific mean-field equations of network-based SIR model are as follows:

with $\theta(t)$ is the probability that any given edge is connected to an infected node, assume the network is uncorrelated, then

where $p(k)$ is the probability that a node is connected to $k$ other nodes and $\langle k \rangle = \sum\limits_{k = 1}\limits^{n}k p(k)$ denotes the mean degree of the network. It is supposed that newborns are balanced by deaths, that is $\Lambda = \mu.$

Table 1.  Definition of parameters.

Parameter	Definition
$\mu$	The birth rate / The death rate
$\beta$	Transmission rate of infected individuals
$\alpha$	The saturated coefficient
$b$	The number of hospital beds
$\gamma_{0}$	The minimum per-capita recovery rate
$\gamma_{1}$	The maximum per-capita recovery rate
$\eta$	Proportion that has been vaccinated

---

 | Show Table

DownLoad: CSV

Inspired by the edge-compartmental method ^[9,10,11], we can rewrite system (2.1) as the following degree-edge-mixed model:

$\mathbf{Remark\ 1.}$ The dimension of degree-edge-mixed model (2.2) is much lower than system (2.1). So the process of analysis presented in this paper will be relatively concise.

$\mathbf{Lemma\ 1.}$ If $\theta(0) > 0, S_{k}(0) > 0,$ then for all $k\in N$ and $t\in R_{+},$ we have $S_{k}(t) > 0, \theta(t) > 0.$

Proof. Firstly, we prove that $\theta(t) > 0$ for all $t > 0.$ From the $(n+1)$th equation of system (2.2), we have $\theta(t) = \theta(0) e^{\int_{0}^{t} z(\tau) d \tau}, z(t) = \frac{\beta}{\langle k\rangle (1+\alpha \theta(t))}\sum\limits_{k = 1}^{n} k^{2} p(k) S_{k}(t)-\mu-\gamma_{0}-\frac{b(\gamma_{1}-\gamma_{0})}{b+\theta(t)}.$ Since $\theta(0) > 0$, we can deduce $\theta(t) > 0$ for all $t > 0.$

Note that $S_{k}(0) > 0,$ in view of the continuity of $S_{k}(t)$, we can find a small $\delta > 0,$ such that $S_{k}(t) > 0$ for $t\in (0, \delta)$. Now we demonstrate that $S_{k}(t) > 0$ for all $t > 0.$ If not, we may assume that there exists $t_{0}\geq \delta > 0,$ s.t. $S_{k}(t_{0}) = 0$ and $S_{k}(t) > 0$ for $t\in(0, t_{0}).$ Thus, from the first $n$ equation of system (2.2), we can obtain $S_{k}^{\prime}(t_{0}) = \mu > 0.$ This leads to $S_{k}(t) < 0$ for some $t\in (0, t_{0}),$ which is apparently a contradiction. Therefore, $S_{k}(t) > 0$ for all $t > 0.$

$\mathbf{Remark\ 2.}$ From the result of Lemma 1 and the fact $S_{k}(t)+I_{k}(t)+R_{k}(t) = 1,$ we can assert that $\Omega: = \left\{(S_{1}, S_{2}, \cdots, S_{n}, \theta)\mid 0 < S_{k}(t) \leq 1, 0 \leq \theta(t) < 1\right\}$ is a positive invariant set of system (2.2).

3.   Existence and stability of equilibrium

System (2.2) has a unique disease-free equilibrium given by $E_{0} = (\frac{\mu}{\mu+\eta}, \frac{\mu}{\mu+\eta}, \cdots, \frac{\mu}{\mu+\eta}, 0)\in R^{n+1}.$ Linearizing the $(n+1)$th equation of system (2.2) at $E_{0}$ yields

Solving Eq (3.1) by a constant variation method, we can build the following equation:

Using the approach introduced in ^[19], we can calculate the basic reproduction number $\mathit{\boldsymbol{R_{0}}}$ by taking Laplace transform:

where $\langle k^{2}\rangle = \sum\limits_{k = 1}\limits^{n}k^{2}p(k).$

3.1.   Stability of the disease-free equilibrium

The local and global stability of the disease-free equilibrium $E_{0}$ will be proved in this subsection.

$\mathbf{Theorem\ 1.}$ If $\mathit{\boldsymbol{R_{0}}} < 1$, the disease-free equilibrium $E_{0}$ of system (2.2) is locally asymptotically stable, whereas if $\mathit{\boldsymbol{R_{0}}} > 1, E_{0}$ is unstable.

Proof. Linearizing system (2.2) at $E_{0}$ yields to

Furthermore, the Jacobian matrix of Eq (3.2) is

Obviously, $\lambda = -(\mu+\eta)$ is $n$-multiple negative eigenvalue. So, the stability of $E_{0}$ is determined by the following eigenvalue:

Therefore, all eigenvalues of the Jacobian matrix $J(E_{0})$ are negative when $\mathit{\boldsymbol{R_{0}}} < 1$, and hence, $E_{0}$ is locally asymptotically stable. By contrast, if $\mathit{\boldsymbol{R_{0}}} > 1, E_{0}$ is unstable.

Next, we will prove the global stability of the disease-free equilibrium $E_{0}$.

$\mathbf{Theorem\ 2.}$ Denote $\hat{\mathit{\boldsymbol{R_{0}}}}$ $= \frac{\beta(b+1)\langle k^{2}\rangle}{(\mu+b\mu+\gamma_{0}+b\gamma_{1})\langle k\rangle}$, if $\hat{\mathit{\boldsymbol{R_{0}}}} \leq 1,$ then $E_{0}$ is globally asymptotically stable.

Proof. Let a Lyapunov function $V(t) = \theta(t),$ then

if $\hat{\mathit{\boldsymbol{R_{0}}}}$ $< 1$, $V^{\prime}(t)\leq 0;$ if $\hat{\mathit{\boldsymbol{R_{0}}}}$ $= 1$, then $V^{\prime}(t)\leq-\frac{b(\gamma_{1}-\gamma_{0})(1-\theta(t))}{(b+1)(b+\theta(t))}\theta(t)\leq0.$ The equality holds if and only if $\theta(t) = 0.$ Hence, we can conclude that the disease-free equilibrium $E_{0}$ is globally asymptotically stable if $\hat{\mathit{\boldsymbol{R_{0}}}}\leq 1$.

$\mathbf{Remark\ 3.}$ Because $\gamma_{0} < \gamma_{1},$ so $\mathit{\boldsymbol{R_{0}}} = \frac{\mu(\mu+b\mu+\gamma_{0}+b\gamma_{1})}{(\mu+\eta)(\mu+b\mu+\gamma_{1}+b\gamma_{1})}\hat {\mathit{\boldsymbol{R_{0}}}} < \hat{\mathit{\boldsymbol{R_{0}}}}.$

$\mathbf{Remark\ 4.}$ If $\hat{\mathit{\boldsymbol{R_{0}}}} < 1,$ that is $\mathit{\boldsymbol{R_{0}}} < \frac{\mu(\mu+b\mu+\gamma_{0}+b\gamma_{1})}{(\mu+\eta)(\mu+b\mu+\gamma_{1}+b\gamma_{1})} < 1$, then $E_{0}$ is globally asymptotically stable.

$\mathbf{Remark\ 5.}$ It seems that the condition $\mathit{\boldsymbol{R_{0}}} < 1$ is not sufficient to ensure the global asymptotic stability of $E_{0}.$ So system (2.2) may exist backward bifurcation at $\mathit{\boldsymbol{R_{0}}} = 1$ which will be proved later.

3.2.   Existence of endemic equilibrium

The existence of feasible positive equilibrium will be proved in this subsection.

$\mathbf{Theorem\ 3.}$ If $\mathit{\boldsymbol{R_{0}}} > 1,$ system (2.2) admits at least an endemic equilibrium.

Proof. Assume that $E^{\ast} = (S_{k}^{\ast}, \theta^{\ast}), k = 1, 2, \ldots, n,$ is an endemic equilibrium of system (2.2), then $E^{\ast}$ should satisfy

Then we can obtain $S_{k}^{\ast} = \frac{\mu(1+\alpha \theta^{\ast})}{\beta k \theta^{\ast}+(\mu+\eta)(1+\alpha \theta^{\ast})}, k = 1, 2, \ldots, n.$ Substituting $S_{k}^{\ast}$ into the $(n+1)$th equation of (3.3), then we can get a self-consistency equation:

Since $\mathit{\boldsymbol{R_{0}}} > 1, F(0) < 0, F(1) > 0,$ using the Immediate Value Theorem, we can conclude that Eq (3.4) has at least one positive root in $(0, 1).$ This implies that system (2.2) admits at least an endemic equilibrium.

$\mathbf{Remark\ 6.}$ If $\mathit{\boldsymbol{R_{0}}}\leq 1,$ then $F(0)\geq 0, F(1) > 0,$ we can not conclude whether system (2.2) exists endemic equilibrium by the method of the Immediate Value Theorem. However, we already know that $E_{0}$ is not globally asymptotically stable, so system (2.2) may also exist endemic equilibrium when $\mathit{\boldsymbol{R_{0}}}\leq 1.$ In other words, the system may take place backward bifurcation at $\mathit{\boldsymbol{R_{0}}} = 1$ which will be proved in Theorem 4.

Next, we will derive a necessary and sufficient condition which determines the existence of backward bifurcation at $\mathit{\boldsymbol{R_{0}}} = 1.$

$\mathbf{Theorem\ 4.}$ System (2.2) exists backward bifurcation at $\mathit{\boldsymbol{R_{0}}} = 1$ if and only if $b < \hat{b};$ the system undergoes forward bifurcation if and only if $b > \hat{b},$ where

Proof. The endemic equilibrium should satisfy Eq (3.4). Replacing $\beta$ by $\frac{(\mu+\eta)(\mu+\gamma_{1})\langle k \rangle \mathit{\boldsymbol{R_{0}}}}{\mu\langle k^{2}\rangle}$ and we can obtain the following equation:

where $Q_{0} = (\mu+\eta)(\mu+\gamma_{1}) \langle k\rangle k \theta^{\ast}\mathit{\boldsymbol{R_{0}}}+(\mu+\eta)(1+\alpha \theta^{\ast})\mu\langle k^{2}\rangle.$

If we keep in mind that $\theta^{\ast}$ is a function of $\mathit{\boldsymbol{R_{0}}}$, then the direction of bifurcation is depended on the sign of ${\frac{\partial\theta^{\ast}}{\partial \mathit{\boldsymbol{R_{0}}}}}\big|_{(\mathit{\boldsymbol{R_{0}}}, \theta^{\ast}) = (1, 0)}.$ More exactly, if ${\frac{\partial\theta^{\ast}}{\partial \mathit{\boldsymbol{R_{0}}}}}\big|_{(\mathit{\boldsymbol{R_{0}}}, \theta^{\ast}) = (1, 0)} < 0,$ then backward bifurcation occurs at $\mathit{\boldsymbol{R_{0}}} = 1.$ Conversely, the forward bifurcation happens.

Next, taking the derivative of Eq (3.5) associated with $\mathit{\boldsymbol{R_{0}}}$ by the implicit function theorem, then yields to the following equation:

where $Q_{1} = \mu(\mu+\eta)\big[(\mu+\gamma_{1}) \langle k\rangle k \theta^{\ast}+(\mu+\gamma_{1})\langle k\rangle k \frac{\partial\theta^{\ast}}{\partial \mathit{\boldsymbol{R_{0}}}}\mathit{\boldsymbol{R_{0}}}+\mu\alpha\langle k^{2}\rangle\frac{\partial\theta^{\ast}}{\partial \mathit{\boldsymbol{R_{0}}}}\big].$

Substituting $(\mathit{\boldsymbol{R_{0}}}, \theta^{\ast}) = (1, 0)$ into Eq (3.6), we can obtain the following equation by simple calculation and arrangement:

Thus if

The same holds true also for the other proposition.

$\mathbf{Remark\ 7.}$ It is a fact that the number of hospital beds does play a key role in determining whether there is a backward bifurcation at $\mathit{\boldsymbol{R_{0}}} = 1.$ More precisely, if $b$ is small enough to satisfy the result of Theorem 4, backward bifurcation will occur. In this case, there will also exist endemic equilibrium even if $\mathit{\boldsymbol{R_{0}}} < 1$ (see Figure 2(b)). That is to say, $\mathit{\boldsymbol{R_{0}}} < 1$ is not sufficient to eradicate the diseases from the population, which is unfavorable to the control of infectious diseases.

Figure 1.  The time evolutions of $\theta(t)$ with different initial conditions $\theta(0) = 0.001\times i, i = 1, 2, \ldots, 20.$ (a) $\hat{\mathit{\boldsymbol{R_{0}}}} = 0.8933 < 1.$ (b) $\mathit{\boldsymbol{R_{0}}} = 1.5008 > 1.$.

DownLoad: Full-Size Img PowerPoint

Figure 2.  (a) The time evolutions of $\theta(t)$ with different initial conditions $\theta(0) = 0.001\times i, i = 1, 2, \ldots, 20.$ (b) The backward bifurcation diagrams in the $\mathit{\boldsymbol{R_{0}}}\theta$-plane with different values of $b.$.

DownLoad: Full-Size Img PowerPoint

Otherwise, when $b$ is big enough, forward bifurcation will take place at $\mathit{\boldsymbol{R_{0}}} = 1.$ In this case, $E_{0}$ will be globally asymptotically stable when $\mathit{\boldsymbol{R_{0}}} < 1.$ This also shows the existence of the threshold for the number of hospital beds. When the relevant departments provide enough hospital beds such that $b > \hat{b},$ infectious diseases can be completely eliminated if $\mathit{\boldsymbol{R_{0}}} < 1.$

From the above discussions, we can summarize the following results.

$\mathbf{Corollary\ 1.}$ For system (2.2),

(1) if $\mathit{\boldsymbol{R_{0}}} < 1$ and

(a) $b < \hat{b},$ there exists two endemic equilibria;

(b) $b > \hat{b},$ there is no endemic equilibrium.

(2) when $\mathit{\boldsymbol{R_{0}}} = 1,$ if $b < \hat{b},$ then there exists a unique endemic equilibrium; otherwise, there is no endemic equilibrium.

(3) if $\mathit{\boldsymbol{R_{0}}} > 1,$ the endemic equilibrium always exist (see Theorem 3) and the number of endemic equilibrium is associated with the value of $b$ (see Figure 3).

Figure 3.  Bifurcation diagrams in the $\mathit{\boldsymbol{R_{0}}}\theta$-plane with different values of $b.$ Forward bifurcation takes place at $\mathit{\boldsymbol{R_{0}}} = 1$ for each subgraph.

DownLoad: Full-Size Img PowerPoint

3.3.   Stability of endemic equilibrium

In this subsection, we will use a geometric approach to prove the local stability of the endemic equilibrium. Linearizing system (2.2) at $E^{*}$ yields

And hence the characteristic equation can be expressed as follows:

where $(x_{k}, y)$ is the eigenvector that correspond to eigenvalue $\lambda.$ From the first equation of (3.8), we can obtain

Substituting $x_{k}$ into the second equation of (3.8) leads to

where $H(\lambda) = \frac{\beta}{\langle k\rangle} \sum\limits_{k = 1}^{n} \frac{k^{3} p(k)\beta S_{k}^{*}\theta^{*}}{\beta k \theta^{*}(1+\alpha\theta^{*})^{2}+(\lambda+\mu+\eta)(1+\alpha\theta^{*})^{3}}+\lambda+\mu+\gamma_{0}+\frac{b^{2}(\gamma_{1}-\gamma_{0})} {(b+\theta^{*})^{2}}-\frac{\beta}{\langle k\rangle}\sum\limits_{k = 1}^{n}\frac{k^{2} p(k) S_{k}^{*}}{(1+\alpha\theta^{*})^{2}}.$

If $y = 0,$ then $\lambda = -\mu-\eta-\frac{\beta k \theta^{*}}{1+\alpha\theta^{*}} < 0.$ Otherwise, when $y\neq 0,$ then $H(\lambda) = 0.$ Since

and hence

Besides, taking the derivative of $F$ associated with $\theta^{*},$ we can have:

Therefore, we can obtain that

$\mathbf{Theorem\ 5.}$ When $F'(\theta^{*}) < 0,$ then $E^{*}$ is unstable; when $F'(\theta^{*}) > 0,$ then $E^{*}$ is locally asymptotically stable.

Proof. When $F'(\theta^{*}) < 0,$ since $H(0) = \theta^{*}F'(\theta^{*}) < 0, \lim\limits_{\lambda\rightarrow +\infty}H(\lambda) = +\infty,$ by the Intermediate Value Theorem, we can conclude that the equation $H(\lambda) = 0$ has at least one positive real solution, so $E^{*}$ is unstable.

When $F'(\theta^{*}) > 0,$ we rewrite the equation $H(\lambda) = 0$ as the following equation:

where $L(\lambda, k, \theta^{*}) = \beta k \theta^{*}(1+\alpha\theta^{*})+(\lambda+\mu+\eta)(1+\alpha\theta^{*})^{2}.$ Suppose Eq (3.13) has a solution $\lambda^{*}$ with $Re{\lambda^{*}}\geq 0,$ then since

if $F'(\theta^{*}) > 0,$ we can conclude that $|\widehat{H}(\lambda^{*})| > 0,$ this contradicts with the fact $\widehat{H}(\lambda^{*}) = 0.$ So Eq (3.13) has no solution with nonnegative real parts. Therefore, when $F'(\theta^{*}) > 0, E^{*}$ is locally asymptotically stable.

$\mathbf{Corollary\ 2.}$ If $\mathit{\boldsymbol{R_{0}}} < 1$ and $b < \hat{b},$ system (2.2) exists two endemic equilibria, the one with smaller value is unstable; the other one with larger value is locally asymptotically stable.

Proof. If $\mathit{\boldsymbol{R_{0}}} < 1$ and $b < \hat{b},$ there exists two endemic equilibria (see Corollary 1). For the endemic equilibrium with smaller value, since $F'(\theta^{*}) < 0$ (see Figure 2(b)), then $E^{*}$ is unstable; and for the endemic equilibrium with larger value, since $F'(\theta^{*}) > 0,$ so $E^{*}$ is locally asymptotically stable.

4.   Simulations

Next, numerical simulations will be presented to validate the previous theoretical results. Our simulations are based on the scale-free network with $p(k) = \xi k^{-2.5}, k = 1, 2, \ldots, 100,$ and the constant $\xi$ is chosen so that the equation $\sum\limits_{k = 1}^{100}p(k) = 1$ holds. For better visualization, we present the values of the parameters in Table 2.

Table 2.  Parameter values of numerical simulations.

Parameter	Figure 1(a)	Figure 1(b)	Figure 2(a)	Figure 2(b)	Figure 3	Figure 4(a)	Figure 4(b)
$\mu$	0.02	0.02	0.01	0.01	0.01	0.02	0.02
$\beta$	0.006	0.03	0.009	—	—	0.03	0.03
$\alpha$	10	10	10	10	10	0/3/6/9/12	10
$b$	0.03	0.03	0.003	0.003/0.005/0.007/0.009/0.011	0.0165/0.017/0.0175/0.018	0.03	0/0.003/0.006/0.009/0.012
$\gamma_{0}$	0.03	0.03	0.005	0.005	0.005	0.03	0.03
$\gamma_{1}$	0.09	0.09	0.08	0.08	0.08	0.09	0.09
$\eta$	0.008	0.008	0.008	0.008	0.008	0.008	0.008

---

 | Show Table

DownLoad: CSV

Firstly, we verify the stability of the disease-free equilibrium $E_{0}.$ According to the selected parameter values, we can calculate $\mathit{\boldsymbol{R_{0}}} = 0.3002, \hat{\mathit{\boldsymbol{R_{0}}}} = 0.8933.$ It is known that $E_{0}$ is globally asymptotically stable from the result of Theorem 2. As seen in Figure 1(a), $\lim\limits_{t \rightarrow +\infty}\theta(t) = 0.$ Therefore, this does support the global stability of $E_{0}$.

In Figure 1(b), we can obtain $\mathit{\boldsymbol{R_{0}}} = 1.5008 > 1$ and hence model (2.2) exists endemic equilibrium. Observing the simulation results, we can see that all trajectories move towards to a positive constant and the unique endemic equilibrium could be stable. However, when $\mathit{\boldsymbol{R_{0}}} > 1,$ model (2.2) can exhibit more complex dynamics (see Figure 3).

Secondly, we will verify that system (2.2) may exist endemic equilibrium when $\mathit{\boldsymbol{R_{0}}} < 1.$ The parameter values as shown in Table 2, since $\mathit{\boldsymbol{R_{0}}} = 0.4280 < 1, \hat{\mathit{\boldsymbol{R_{0}}}} = 4.5545 > 1,$ therefore $E_{0}$ is not globally asymptotically stable from the result of Theorem 2. From Figure 2(a), it can be observed that the trajectories partially converge to zero or partially move towards to a positive level which corresponds to different initial conditions. Therefore, the initial infectious density must be controlled to a lower level for the control of the diseases. The result of the co-existence of two locally asymptotically stable equilibria is named as the bistable phenomenon.

Different values of $b$ are taken to show the bifurcation diagrams. In Figure 2(b), backward bifurcation occurs at $\mathit{\boldsymbol{R_{0}}} = 1$ when $b < \hat{b} = 0.016.$ Besides, the larger the value of $b,$ the smaller the depth of backward bifurcation. In Figure 3, the parameters are the same as those given in Figure 2, except for the value of $b.$ It can be seen that forward bifurcation occurs at $\mathit{\boldsymbol{R_{0}}} = 1$ for each subgraph. The reason is that although the values of $b$ are different, they are all bigger than $\hat{b}.$ Hence, the numerical simulation results are consistent with the theoretical results of Theorem 4.

At the same time, we can observe that the shapes of the bifurcation curves are different from Figure 3. For example, the endemic equilibrium curve has an "S shape" when $b = 0.017.$ That is to say, backward bifurcation takes place at a certain value of $\mathit{\boldsymbol{R_{0}}}.$ And the situation results in the existence of three endemic equilibria. So, the number of endemic equilibrium is uncertain when $\mathit{\boldsymbol{R_{0}}} > 1.$

Lastly, we will investigate the effect of $\alpha$ and $b$ on diseases transmission through numerical simulations. The parameters in Figure 4 are the same as those in Figure 1(b) and $\mathit{\boldsymbol{R_{0}}} = 1.5008 > 1.$ From Figure 4, we can observe that the larger $\alpha$ or $b,$ the value of $\theta(t)$ at steady state will be smaller, which means that the final density of infected individuals will also be smaller. This indicates that the behavioral change of individuals can indeed suppress the spread of the diseases. Besides, enriching the adequate hospital beds is also important for the control of the diseases.

Figure 4.  The time evolutions of $\theta(t)$ with initial condition $\theta(0) = 0.01.$ (a) $\alpha = 3\times(i-1),$ (b) $b = 0.003\times(i-1), i = 1, 2, \ldots, 5$.

DownLoad: Full-Size Img PowerPoint

5.   Conclusions and discussion

Many scholars proved that the number of hospital beds plays an important role in the occurrence of bifurcation ^{[5,6,20,21,22]}, but most of them studied the homogeneously mixed models. In ^[6], Cui et al. obtained a quadratic equation with respect to $I$ and discussed the existence and classification of equilibrium from the algebraic perspective. The condition for the existence of backward bifurcation is proved by reducing system (1.3) into the center manifold. In this paper, a network-based SIR epidemic model which is the extension of model (1.3) is proposed. Using the edge-compartmental approach, we reduced the dimension of the system. It is proved that when ${\frac{\partial\theta^{\ast}}{\partial \mathit{\boldsymbol{R_{0}}}}}\Big|_{(\mathit{\boldsymbol{R_{0}}}, \theta^{\ast}) = (1, 0)} < 0 \Leftrightarrow b < \hat{b},$ backward bifurcation will occur at $\mathit{\boldsymbol{R_{0}}} = 1.$ The method is different from the approach in ^[6]. Moreover, the condition which determines the direction of bifurcation at $\mathit{\boldsymbol{R_{0}}} = 1$ is related to the network structure.

It is also possible that system (2.2) exists multiple endemic equilibria when $\mathit{\boldsymbol{R_{0}}} > 1$ (see Figure 3). In conclusion, the condition $\mathit{\boldsymbol{R_{0}}} < 1$ is not enough to guarantee the global asymptotic stability of $E_{0},$ and $\mathit{\boldsymbol{R_{0}}} > 1$ does not confirm the uniqueness of endemic equilibrium. In addition, we established a relationship between the local stability of endemic equilibrium and the tangent slope of the epidemic curve. It is observed that the larger $\alpha$ or $b$ can reduce the endemic level (see Figure 4). Therefore, the effect from the individual behavioral change should not be ignored and the hospital beds should be adequately provided in the process of controlling the infectious diseases.

However, whether system (2.2) can exhibit more complex dynamical behavior such as Hopf, saddle-node, transcritical and Bogdanov-Takens bifurcation are not discussed. At present, the bifurcation theory on networks is still less. It will be interesting and important to solve the problem and we leave this for future work.

Use of AI tools declaration

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

Acknowledgments

Research Project Supported by National Nature Science Foundation of China (Grant Nos. 12071445, 12271519), the Fundamental Research Program of Shanxi Province (Grant Nos. 20210302123031, 202203021211091), the Science and Technology Innovation Program of higher education in Shanxi Province (Grant No. 2023L436).

Conflict of interest

The authors declare there is no conflict of interest in this paper.