Global Hopf bifurcation of a cholera model with media coverage

1.   Introduction

In modern society, media coverage has become an important strategy for controlling and preventing disease transmission. It can alter an individual's behavior, such as wearing protective masks, vaccination, self-isolating or avoiding gathering activities, and hence reduce the possibility of contracting the infection ^[1]. Additionally, the media coverage also influences the implementation of a public health policy intervention and control polices ^[2]. Therefore, how to quantify this impact through mathematical modeling is an important issue for epidemics control.

Many mathematical models have been developed to investigate the impact of media coverage on disease transmission and control. The most general approach to modeling the media impact is to change the transmission coefficient as a nonlinear decreasing function with respect to the number of infectious individuals $(I)$. For example, Cui et al. ^[3] used a contact transmission rate $\mu e^{-m I}$ to describe the impact of media coverage, where $\mu > 0$ is the probability of baseline transmission and $m > 0$ reflects the effect of the media. Li and Cui ^[4] employed $\frac{\beta_2 I}{m+I}$ to reflect the reduced amount of contact rate due to media coverage, producing the transmission rate of the form $\beta_1-\frac{\beta_2I}{m+I}\; (\beta_1 > \beta_2 > 0, m > 0)$. By considering the correlation between media impact and the number of infected individuals at different disease stages, Liu et al. ^[5] devised the media impact function $\beta e^{-a_1E-a_2I-a_3H}$ with $E$ denoting exposed individuals, $I$ infectives, $H$ hospitalized individuals and nonnegative constants $a_1, a_2, a_3$. However, it is noted that individuals can also change their behavior in response to the rate of change in case numbers. Xiao et al. ^[6] used an exponentially decreasing function $e^{-M(I, \frac{{\mathrm{d}} I}{{\mathrm{d}} t})}$ with $M(I, \frac{{\mathrm{d}} I}{{\mathrm{d}} t}) = \max\{0, p_1I+p_2 \frac{{\mathrm{d}} I}{{\mathrm{d}} t}\}$ to model such media impacts. For more studies on media impacts, we refer readers to ^{[1,2,7,8,9,10]} and references therein.

However, as mentioned in ^[10], impact of media coverage on the population transmission dynamics has lag, describing both the lag time of individuals' response to the media and the lag time of media reports about an infectious disease outbreak. Some studies have been conducted to investigate the impact of delayed mass media, see ^[1,8,10,11]. In this paper, we will consider the time delay of the media impact in the mathematical modeling of cholera, an acute intestinal infectious disease caused by the bacterium Vibrio cholerae (or, V. cholerae). Several mathematical models have been proposed to study cholera dynamics, including, but not limited to, ^{[12,13,14,15,16,17,18]}. However, to our knowledge, few of these have specifically taken into account the response delay in media impacts.

Our model is inspired by the disease transmission models in ^[10,13]. We divide the population into three subclasses: Susceptible $S(t)$, infectives $I(t)$ and recovered $R(t)$. The concentration of bacteria V. cholerae in contaminated water is indicated by $B(t)$. To study the lag effect of the impact of the media on cholera transmission, we incorporate a decreasing factor $e^{-m I(t-\tau)}$ into the direct and indirect incidence rate. That is,

Here, $\beta_1$ and $\beta_2$ are, respectively, the baseline direct and indirect transmission rate, $\tau$ represents the report delay and response time of individuals to the current infection and $m$ is the weight of media impact sensitive to the case number. The basic SIRB model under consideration is then

where $\Lambda$ stands for the influx rate of susceptible humans, $\mu$ is the natural death rate, $\sigma$ is the immunity waning rate, $\gamma$ is the recovery rate, $\xi$ denotes the rate of contribution to V. cholerae in the aquatic environment and $\delta = \delta_1-\delta_2 > 0$ is the net death rate of V. cholerae in the aquatic environment, where $\delta_1, \delta_2$ are the natural washout rate and the natural growth rate of V. cholerae in the aquatic environment, respectively. All the parameters of model (1.1) are assumed to be positive.

Note that the total population size $N(t) = S(t)+I(t)+R(t)$ satisfies

and has a unique equilibrium $N^* = \frac{\Lambda}{\mu}$ that is globally stable in $\mathbb{R}_+$. We then consider the following limiting system:

The remainder of this work is organized as follows. In Section 2, we first give the well-posedness of system (1.2) and then establish its threshold dynamics concerning the basic reproduction number. In Section 3, we are devoted to the stability and Hopf bifurcation analysis of the positive equilibrium. In Section 4, we investigate the global continuation of a local branch. In Section 5, we perform numerical simulations to illustrate our analytical results. The last section is a brief discussion.

2.   Well-posedness of the model and threshold dynamics

Let $\mathcal {C}: = C([-\tau, 0], \mathbb{R}^3)$ and $\mathcal {C}^{+}: = C([-\tau, 0], \mathbb{R}_+^3)$. Then $(\mathcal {C}, \mathcal {C}^{+})$ is an ordered Banach space equipped with the supremum norm $\|\phi\| = \sup_{-\tau\leq \theta \leq 0} |\phi(\theta)|$ for $\phi\in \mathcal {C}$. For any given continuous function $u = (u_1, u_2, u_3): [-\tau, \rho)\rightarrow \mathbb{R}^3$ with $\rho > 0$, we define $u_t\in \mathcal {C}$ by

for any $t\in [0, \rho)$. Set

We then have the following result.

Lemma 2.1. For any $\varphi\in \Gamma$, system (1.2) admits a unique nonnegative bounded solution $u(t, \varphi)$ on $[0, \infty)$ with $u_0 = \varphi$, and $u_{t}(\varphi) \in \Gamma, \; \forall t\geq0$.

Proof. For any $\varphi = (\varphi_1, \varphi_2, \varphi_3)\in \Gamma$, we define

Obviously, $f(\varphi)$ is continuous in $\varphi\in\Gamma$ and Lipschitz in $\varphi$ on each compact subset of $\Gamma$. It follows from [19, Theorems 2.2.3] that system (1.2) has a unique solution $u(t, \varphi)$ on its maximal existence interval $[0, t_\varphi)$ with $u_0 = \varphi$.

For any given $\varphi\in \Gamma$, one sees that if $\varphi_i(0) = 0$ for some $i\in\{1, 2, 3\},$ then $f_{i}(\varphi) \geq 0$. By [20, Theorem 5.2.1], we obtain $u(t, \varphi)\geq 0$ for all $t\in[0, t_\varphi)$. Let $H(t) = I(t)+R(t)$. Then one has

which implies that $I_t+R_t\leq N^*$ for $t\in [0, t_\varphi)$, and hence, $u_t(\varphi)\in \Gamma$ for all $t\in [0, t_\varphi)$. Moreover, by the third equation of (1.2), one easily sees that $u_3(t, \varphi)$ is bounded on $[0, t_\varphi)$. As a result, [19, Theorem 2.3.1] means that $t_{\varphi} = \infty$.

Note that the disease-free equilibrium of (1.2) is $E_0 = (0, 0, 0)$. Since the $R$ class do not participate in the transmission of cholera, we only consider the infected equations of the linearization of (1.2) at $E_0$:

Based on the next generation matrix method ^[21], we define

Hence, the basic reproduction number of system (1.2) is given by

where $r(FV^{-1})$ is the spectral radius of $FV^{-1}$, and $\mathcal {R}_1^{d}$ (resp. $\mathcal {R}_2^{i}$ is the basic reproduction number for the direct (resp. indirect) transmission in the absence of indirect (resp. direct) transmission.

Now we begin to study the local stability of the disease-free equilibrium. The Jacobian matrix of system (1.2) at $E_0$ has the form

and the corresponding characteristic equation is

where

When $\mathcal {R}_0 < 1$, one easily sees

Using the Routh-Hurwitz criterion, we know that all roots of (2.2) have negative real parts. Hence, $E_0$ is locally stable if $\mathcal {R}_0 < 1$. In the case where $\mathcal {R}_0 > 1$, the Eq (2.2) has at least one positive root since $d < 0$, and hence, $E_0$ is unstable. Now we proceed to the global stability of $E_0$.

Theorem 2.1. If $\mathcal {R}_0 < 1$, then $E_0$ is globally asymptotically stable in $\Gamma$.

Proof. From the second and third equations of (1.2), we see

which is equivalent to

Since the matrix $V^{-1}F$ is nonnegative, irreducible and $\mathcal {R}_0 = r(FV^{-1}) = r(V^{-1}F)$, there exists a positive left eigenvector $v$ such that

Define the functional $L:\; \Gamma\rightarrow \mathbb{R}$ as follows

For any solution $u(t, \varphi)$ of (1.2) with $u_0 = \varphi\in \Gamma$, with the help of (2.3), direct calculation yields

Thus, $L$ is a Lyapunov functional on $\Gamma$ relative to system (1.2).

Let $M$ be the largest compact invariant subset in the set $\{\phi\in \Gamma: \; \dot{L}_{(1.2)}(\phi) = 0\},$ where $\dot{L}_{(1.2)}(\phi)$ denotes the derivative of $L$ along the solution of (1.2). It is easy to see from (2.4) that $\dot{L}_{(1.2)}(\phi) = 0$ implies that $I(t) = B(t) = 0$ for any $t\geq 0$, and hence, $M = \{\phi\in \Gamma:\; \phi_1 = \phi_3 = 0\}$. It then follows from the LaSalle invariance principle (see, e.g., [22, Theorem 1]) that $\lim_{t\rightarrow \infty} (I(t), B(t)) = (0, 0)$ and further, $\lim_{t\rightarrow \infty} R(t) = 0$. Therefore, $E_0$ is globally attractive in $\Gamma$. This, together with the local stability of $E_0$, gives the desired result.

Next we are ready to show that the disease is uniformly persistent when $\mathcal{R}_0 > 1$. To this end, we need the following lemma.

Lemma 2.2. For any given $\varphi\in\Gamma$, if there exists some $t_0\geq 0$ such that $I(t_0, \varphi) > 0$ or $B(t_0, \varphi) > 0$, then $I(t, \varphi) > 0$ and $B(t, \varphi) > 0$ for all $t > t_0$.

Proof. If $I(t_0, \varphi) > 0$ for some $t_0$, then

By the third equation of (1.2), we have

Similarly, if $B(t_0, \varphi) > 0$ for some $t_0$, we obtain that $B(t, \varphi) > 0$ for all $t\geq t_0$. We will show that $H(t, \varphi): = I(t, \varphi)+R(t, \varphi) < N^*$ for all $t > t_0$. If the assertion was false, then there exists $t_1\in (t_0, \infty)$ such that

Clearly, $\frac{{\mathrm{d}} H(t, \varphi)}{{\mathrm{d}} t}\Big|_{t = t_1}\geq 0$ must hold. However, in view of the first two equations of (1.2), we find

This contradiction means that no such $t_1$ can exist. Accordingly, by the $I$-equation in (1.2), we have

for all $t > t_0$. Summarizing these two cases, we establish the desired result.

Theorem 2.2. If $\mathcal{R}_0 > 1$, then system (1.2) admits a unique componentwise positive equilibrium. Moreover, there exists $\eta > 0$ such that for any $\varphi\in \Gamma$ with $\varphi_1(0)\neq 0$ or $\varphi_3(0)\neq 0$, the solution $u(t, \varphi) = (I(t, \varphi), R(t, \varphi), B(t, \varphi))$ satisfies

Proof. Define the following two sets:

By the form of (1.2), it can be verified that both $\Gamma$ and $\Gamma_0$ are positively invariant. Clearly, $\partial \Gamma_0$ is relatively closed in $\Gamma$. Let $\Phi(t)$ be the solution maps of system (1.2), namely

Note that for each $t > \tau$, $\Phi(t)$ is continuous and compact (see [19, Theorem 3.6.1]). The ultimate boundedness of solutions implies that $\Phi(t)$ is point dissipative. It then follows from [23, Theorem 3.4.8] that $\Phi(t)$ has a global attractor $K$.

Let $\omega(\varphi)$ be the omega limit set of the forward orbit through $\varphi$ for $\Phi(t)$ and define

Then the following claim holds true.

Claim 1. $E_0$ is globally stable for $\Phi(t)$ in $M_\partial$.

For any given $\varphi \in M_{\partial}$, we have $\Phi(t) \varphi \in \partial \Gamma_{0}, \forall t\geq 0$. Hence, for each $t\geq0$, either $I(t, \varphi) = 0$ or $B(t, \varphi) = 0$. We further show

Suppose to the contrary that $\varphi\not\in M_0$. Then either $\varphi_1(0) = I(0, \varphi) > 0$ or $\varphi_3(0) = B(0, \varphi) > 0$. With the help of Lemma 2.2, we have that $I(t, \varphi) > 0$ and $B(t, \varphi) > 0$ for all $t > 0$, which contradicts the fact that $\varphi\in M_\partial$. Thus, (2.6) holds. It then follows that $I(t, \varphi) = 0$ and $B(t, \varphi) = 0$ for all $\varphi\in M_\partial$ and $t\geq 0$, and further, $\lim_{t\rightarrow \infty} R(t, \varphi) = 0$. Therefore, $\omega(\varphi) = E_0$ for any $\varphi\in M_\partial$. In other words, $E_0$ is globally attractive for $\Phi(t)$ in $M_\partial$. Note that system (1.2) restricted to $M_\partial$ becomes $\frac{{\mathrm{d}} R(t)}{{\mathrm{d}} t} = -(\mu+\sigma)R(t)$. In view of [24, Lemma 2.2.1], $E_0$ is locally Lyapunov stable for $\Phi(t)$ in $M_{\partial}$. So, the above claim is proved.

Since $\mathcal {R}_0 > 1$, we can choose $\epsilon > 0$ small enough such that

Claim 2. $\limsup _{t \rightarrow \infty}\|\Phi(t) \varphi-E_0\| \geq \epsilon, \; \; \forall \varphi\in \Gamma_0.$ Suppose the claim is not true. Then there exists some $\psi\in\Gamma_0$ such that

Thus, there is $\bar{t} > 0$ such that

It follows that for all $t\geq \bar{t}+\tau$, we have

Let $(x(t), y(t))$ be the solution of the following system:

and its Jacobian matrix is

Recall that the stability modulus of $M_\epsilon$ is defined by

Since $M_\epsilon$ is quasi-positive and irreducible, it follows from [20, Corollary 3.2] that $s(M_\epsilon)$ is a simple eigenvalue of $M_\epsilon$ with a strongly positive eigenvector $v_\epsilon$. In particular, $e^{s(M_\epsilon)t}v_\epsilon$ is a solution of (2.7). By using the proof of [21, Theorem 2], we know that $s(M_\epsilon) > 0$ if and only if $\mathcal {R}_0^\epsilon > 1$.

Since $(I(t, \psi), B(t, \psi))\gg(0, 0)$ for all $t > 0$, there exists $\kappa > 0$ such that

Hence by comparison,

In view of $s(M_\epsilon) > 0$, we see that $I(t, \psi)$ and $B(t, \psi)$ tend to infinity as $t\rightarrow \infty$. This contradicts the boundedness of solutions, and thus the claim holds.

Claim 1 shows that $E_0$ cannot form a cycle in $\partial\Gamma_0$, and claim 2 implies that $E_0$ is an isolated invariant set in $\Gamma$ and $W^s(E_0)\cap\Gamma_0 = \emptyset$, where $W^s(E_0)$ is the stable set of $E_0$ for $\Phi(t)$. According to the acyclicity theorem of the uniform persistence of maps (see, e.g., [24, Theorem 1.3.1,Remark 1.3.1 and Remark 1.3.2]), we conclude that $\Phi(t): \Gamma\rightarrow \Gamma$ is uniformly persistent with respect to $(\Gamma_0, \partial \Gamma_0)$. Moreover, by [25, Theorem 2.4], there exists a global attractor $A_0$ for $\Phi(t)$ in $\Gamma_0$ and system (1.2) has an stationary coexistence state $\bar{\varphi} = (\bar{\varphi}_1, \bar{\varphi}_2, \bar{\varphi}_3)\in \Gamma_0$, and $\Phi(t)\bar{\varphi} = \bar{\varphi}$ for all $t\geq 0$. Noting that $\bar{\varphi}$ is a constant function, we let $\bar{I} = \bar{\varphi}_1(0)$, $\bar{R} = \bar{\varphi}_2(0)$ and $\bar{B} = \bar{\varphi}_3(0)$. Then $\bar{R}\in \mathbb{R}_+$ and $(\bar{I}, \bar{B})\in {\rm int}(\mathbb{R}_+^2)$. We further claim that $\bar{R}\in \mathbb{R}_+\setminus \{0\}$. Suppose that $\bar{R} = 0$. By the $R$-equation in (1.2), we then get $0 = \gamma \bar{I}$, and hence, $\bar{I} = 0$, a contradiction. Therefore, $E_1 = (\bar{I}, \bar{R}, \bar{B})$ is a componentwise positive equilibrium of (1.2).

Next, we show the uniqueness of the positive equilibrium. Suppose that $(\hat{I}, \hat{R}, \hat{B})$ is the equilibrium of (1.2). Then $(\hat{I}, \hat{R}, \hat{B})$ satisfies

Set $\hat{S} = N^*-\hat{I}-\hat{R}$. A direct verification yields

This implies that $(\hat{S}, \hat{I}, \hat{R}, \hat{B})$ is the equilibrium of (1.1). In view of (2.8), we see that

Therefore, $\mathcal {F}(\hat{I})$ is strictly decreasing in $\hat{I}\in \mathbb{R}_+$, and $\mathcal {G}(\hat{I})$ is strictly increasing in $\hat{I}\in \mathbb{R}_+$. If $\mathcal {R}_0 > 1$, then $\mathcal {F}(0) > \mathcal {G}(0)$, which means that there is a unique intersection in $\mathbb{R}_+^2$ between $\mathcal {F}(\hat{I})$ and $\mathcal {G}(\hat{I})$, and thus $(\hat{I}, \hat{R}, \hat{B})$ is the unique positive equilibrium of (1.2). Moreover, the uniqueness of positive equilibrium also implies that of $E_1$.

Finally, to derive the practical persistence, we define a continuous function $p:\Gamma\rightarrow \mathbb{R}_+$ by

Clearly, $\Gamma_0 = p^{-1}(0, \infty)$ and $\partial\Gamma_0 = p^{-1}(0)$. Since $A_0$ is a compact subset of $\Gamma_0$, it follows that $\inf_{\varphi\in A_0} = \min_{\varphi\in A_0} p(\varphi) > 0$. Consequently, there exists $\eta > 0$ such that

Moreover, for any given $\varphi\in \Gamma$ with $\varphi_1(0)\neq 0$ or $\varphi_3(0)\neq 0$, Lemma 2.2 suggests that $I(t, \varphi) > 0$ and $B(t, \varphi) > 0$ for all $t > 0$. Fix a $t_0 > \tau$. By using the fact $\Phi(t)\varphi = \Phi(t-t_0)(\Phi(t_0)\varphi), \forall t > t_0$, we see that statement (2.5) holds true.

Remark 2.1. The uniqueness of the positive equilibrium $E_1$ can be obtained directly. In fact, from (1.2) we get

Solving the above equation with respect to $\bar{I}$ gives rise to

where the definition of Lambert W function is seen in ^[26]. It is positive (i.e. $\bar{I} > 0$) provided

which is equivalent to $\mathcal {R}_0 > 1$. In addition, this explicit expression makes the numerical approximation of $(\bar{I}, \bar{R}, \bar{B})$ more convenient.

3.   Stability analysis and local Hopf bifurcations

In this section, we will consider the global stability of the endemic equilibrium $E_1$ and the existence of local Hopf bifurcations of system (1.2).

3.1.   Stability of $E_1$ for system (1.2) without delay

The linearized form of system (1.2) at $E_1 = (\bar{I}, \bar{R}, \bar{B})$ is given by

where $A = e^{-m \bar{I}}(N^*-\bar{I}-\bar{R})$ and $C = \beta_1 e^{-m \bar{I}}\bar{I}+\beta_2 e^{-m \bar{I}}\bar{B}$. The characteristic equation of system (3.1) is

where

When $\tau = 0$, the equation (3.2) becomes

In view of $a_0+b_0 > 0$ and

where $X = \beta_2A\frac{\xi}{\delta}+\delta+C+m(\mu+\gamma) \bar{I}$, we deduce from the Routh-Hurwitz criterion that all roots of (3.3) have negative real parts. Hence, the following result holds true.

Theorem 3.1. Suppose that $\mathcal {R}_0 > 1$. Then the endemic equilibrium $E_1$ of system (1.2) is locally asymptotically stable when $\tau = 0$.

Now we present the geometric approach developed in ^[27] to study the global stability of $E_1$. Consider the following autonomous system

where $D\subset \mathbb{R}^n$ is a simply connected open set. Let $x(t, x_0)$ be the solution of (3.4) such that $x(0, x_0) = x_0$. The second compound equation with respect to $x(t, x_0)\in D$ is

where $\frac{\partial f}{\partial x}^{[2]}$ is the additive compound matrix of the Jacobian matrix $\frac{\partial f}{\partial x}$ (see ^[28]).

Let $P(x)$ be a $k\times k$ matrix-valued $C^1$ function, where $k = \frac{1}{2}n(n-1)$, and suppose $P^{-1}(x)$ exists for $x\in D$. Set

where $P_f$ is the matrix obtained by replacing each entry $p_{ij}$ in $P$ with its directional derivative in the direction of $f$. The Lozinskiǐ measure of $Q$ with respect to the matrix norm $|\cdot|$ in $\mathbb{R}^{n\times n}$ is defined as

Define a quantity $\bar{q}_2$ as

where $K$ denotes a compact absorbing set in $D$. The following result from [27, Theorem 3.5] is used to prove the global stability of $E_1$.

Lemma 3.1. Assume that system (3.4) has a unique equilibrium $x^*$ in $D$. Then $x^*$ is globally stable in $D$ if $\bar{q}_2 < 0$.

Theorem 3.2. If $\mathcal {R}_0 > 1$, $\tau = 0$ and $\mu+\gamma > \sigma$, then $E_1$ is globally asymptotically stable in $\Omega = \{(I, R, B)\in \mathbb{R}_+^3:\; I+R\leq N^*, \; I > 0, \; B > 0\}$.

Proof. Note that system (1.2) is equivalent to

Thus, under the condition $\mathcal {R}_0 > 1$, the global stability of $E_1$ is equivalent to that of the unique positive equilibrium $\hat{E}_1$ of system (3.5), where $\hat{E}_1 = (N^*-\bar{I}-\bar{R}, \bar{I}, \bar{B})$. Define

Clearly, $\tilde{\Omega}$ is a simply connected region in $\mathbb{R}^3$. Similar to the arguments in Theorem 2.2, one can prove that system (3.5) is uniformly persistent in $\tilde{\Omega}$. This, together with the boundedness of solutions to (3.5), implies the existence of a compact absorbing set $K\subset \tilde{\Omega}$. Therefore, by Lemma 3.1, it suffices to show that there exists a matrix-valued function $P(x)$ such that $\bar{q}_2 < 0$.

Define the diagonal matrix $P$ as

Then $P$ is $C^1$ and nonsigular in $\tilde{\Omega}$. Let $f$ denote the vector field of (3.5). Then

where $'$ is the derivative with respect to time $t$. The Jacobian matrix $J$ associated with a general solution $(S(t), I(t), B(t))$ to (3.5) is

where $g(I, B) = \beta_1e^{-mI}I+\beta_2e^{-mI}B$. The second additive compound matrix of $J = (j_{ik})_{3\times 3}$, denoted by $J^{[2]}$, is expressed by

As a result, the matrix $Q = P_{f} P^{-1}+PJ^{[2]} P^{-1}$ can be written in block form as follows:

where

The vector norm $|\cdot|$ in $\mathbb{R}^3$ is taken as

As described in ^[27], we have the following estimate

with

Here, $\left|Q_{12}\right|$ and $\left|Q_{21}\right|$ are matrix norms induced by the $l_1$ norm, and $\mu_1$ denotes the Lozinskiǐ measure concerning the $l_1$ norm. Moreover, using the method in ^[29], we calculate $\mu_1(Q_{22})$ as

where $H(S, I, B) = \beta_1e^{-mI}S-mSg(I, B)+\sigma.$ In view of

we have

Define the function

Along each solution $(S(t), I(t), B(t))$ to (3.5) such that $(S(0), I(0), B(0))\in K$, we find that

Here we have used the boundedness of solutions to (3.5). Observe that

then $\bar{q}_2$ has two possibilities:

for all $(S(0), I(0), B(0))\in K$. Therefore, the desired result holds.

3.2.   Local Hopf bifurcations

From Theorem 3.1, we know that if $\mathcal {R}_0 > 1$ and $\tau = 0$, $E_1$ is locally asymptotically stable. Further, 0 cannot be an eigenvalue of (3.2) due to $p(0, \tau) = a_2+b_2 > 0$ for any $\tau\geq 0$. Therefore, the stability of $E_1$ changes only when at least a pair of eigenvalues of (3.2) cross the imaginary axis to the right. We thus suppose that $\lambda = {\rm i}\omega$ ($\omega > 0$) is a purely imaginary solution of (3.2) for some $\tau > 0$, that is,

Separating the real and imaginary parts, we obtain

Equivalently,

Squaring and adding both equations of (3.6) yields

with

Let $x = \omega^2$. Then equation (3.7) becomes

Therefore, if ${\rm i} \omega$ is a purely imaginary root of (3.2), then the equation

has a positive root $x = \omega^2$. Define the set $W$ as

By using the properties for the general cubic equations given in [30, Lemma A.2], we have $(p, q, r)\in W$ if and only if one of the following conditions holds:

(C1) $\Delta = (pq-9r)^2-4(p^2-3q)(q^2-3pr) > 0, r < 0$;

(C2) $r = 0, q = 0, p < 0$;

(C3) $\Delta < 0, q > 0, p > 0, r < 0$.

In the case where $(p, q, r)\in W$, let $\omega_0 = \sqrt{x^*}$. Solving (3.6) for $\tau$ and $\omega = \omega_0$, we obtain

Then we have the following result.

Theorem 3.3. Let $\mathcal {R}_0 > 1$. For $(p, q, r)\in W$, $E_1$ is locally asymptotically stable for $\tau\in[0, \tau_0)$ and unstable for $\tau > \tau_0$. Besides, system (1.2) undergoes Hopf bifurcation at $E_1$ when $\tau = \tau_n$, $n = 0, 1, 2, \cdots$.

Proof. Differentiating both sides of equation (3.2) with respect to $\tau$ yields

Substituting $\tau = \tau_n$ into the above equality, we obtain

Since $\omega^2[(b_2-b_0\omega^2)^2+b_1^2\omega^2] > 0$, we have

Therefore, the transversality condition holds, and so the desired result follows.

4.   Global Hopf bifurcation

Theorem 3.3 states that periodic solutions can bifurcate from $E_1$ when $\tau$ is near the local Hopf bifurcation values $\tau_n$, $n = 0, 1, 2, \cdots$. In this section, we use the global Hopf bifurcation theorem (see [31, Therorem 3.2]) to study the global continuation of these locally bifurcating periodic solutions. In the remainder of this section, we always assume that $\mathcal {R}_0 > 1$, $\mu+\gamma > \sigma$ and $(p, q, r)\in W$.

Our arguments are similar to those in ^{[1,10,30,32,33,34,35]}. Let $z(t) = (I(\tau t), R(\tau t), B(\tau t))^{\mathrm{T}}.$ Rewrite system (1.2) as the following functional differential equation

where $z_t\in Y: = C([-1, 0], \mathbb{R}_+^3)$ with $z_t(\theta) = z(t+\theta)$ for $\theta\in[-1, 0],$ parameter $T$ is the period of the non-constant periodic solution of (4.1), and

with $z_t = (z_{1t}, z_{2t}, z_{3t})\in Y$. Restricting $F$ to the subspace of $Y$, we get a restricted function

Clearly, $\tilde{F}$ is twice continuously differentiable, that is, the assumption (A1) in ^[31] holds.

By Theorems 2.1 and 2.2, the set of stationary solutions of system (4.1) is given by

For any stationary solution $(\tilde{z}, \tau, T)\in N(F)$, the characteristic matrix is

where Id is the 3 $\times$ 3 identity matrix, $\tilde{A} = e^{-m \tilde{z}_1}(N^*-\tilde{z}_1-\tilde{z}_2)$ and $\tilde{C} = \beta_1 e^{-m \tilde{z}_1}\tilde{z}_1+\beta_2 e^{-m \tilde{z}_1}\tilde{z}_3$. Thus, for any stationary solution $(\tilde{z}, \tau, T)$, the characteristic equation reads

where

Under the condition $\mathcal {R}_0 > 1$, 0 cannot be an eigenvalue of any stationary solution of (4.1). Therefore, the condition (A2) in ^[31] holds. From Eq (4.2), it can be easily verified that the smoothness condition (A3) in ^[31] is also valid.

As defined in ^[36], the stationary solution $(\tilde{z}, \tilde{\tau}, \tilde{T})$ of (4.1) is called a center if ${\rm det} \Delta_{(\tilde{z}, \tilde{\tau}, \tilde{T})}({\rm i} k\frac{2\pi}{\tilde{T}}) = 0$ for some positive integer $k$. A center is isolated if no other center exists in some neighborhood of $(\tilde{z}, \tilde{\tau}, \tilde{T})$ and there are only finite pure imaginary eigenvalues of the form ${\rm i}k\frac{2\pi}{\tilde{T}}$. Let $J(\tilde{z}, \tilde{\tau}, \tilde{T})$ represent the set of all such positive integers $k$. Theorem 3.3 implies that if $(p, q, r)\in W$, for any integer $n\geq0$, $(E_1, \tau_n, \frac{2\pi}{\omega_0 \tau_n})$ is an isolated center of (4.1), and there is only one pure imaginary root of the form ${\rm i}k\frac{2\pi}{\tilde{T}}$ with $k = 1$ and $\tilde{T} = \frac{2\pi}{\omega_0 \tau_n}$. Therefore,

Moreover, it follows from Theorem 3.3 that the crossing number at each of these center is

Thus the condition (A4) in ^[31] holds.

Next, we define a closed subset $\Sigma(F)$ of $Y\times (0, \infty) \times \mathbb{R}_+$ by

For each integer $n\geq 0$, let $C\left(E_1, \tau_n, \frac{2\pi}{\omega_0 \tau_n}\right)$ denote the connected branch of $C\left(E_1, \tau_n, \frac{2\pi}{\omega_0 \tau_n}\right)$ in $\Sigma(F)$. Theorem 3.3 guarantees that $C\left(E_1, \tau_n, \frac{2\pi}{\omega_0 \tau_n}\right)$ is a nonempty subset of $\Sigma(F)$. The global bifurcation theorem [31, Theorem 3.4] means that one of the following holds:

(ⅰ) $C\left(E_1, \tau_n, \frac{2\pi}{\omega_0 \tau_n}\right)$ is unbounded in $Y\times (0, \infty) \times \mathbb{R}_+$;

(ⅱ) $C\left(E_1, \tau_n, \frac{2\pi}{\omega_0 \tau_n}\right)$ is bounded, $C\left(E_1, \tau_n, \frac{2\pi}{\omega_0 \tau_n}\right) \cap N(F)$ is finite and

for all $k = 1, 2, 3, \cdots$, where $\gamma_{k}(\tilde{z}, \tau, T)$ is the $k$-th crossing number of $(\tilde{z}, \tau, T)$ if $k \in J(\tilde{z}, \tau, T)$, otherwise, $\gamma_{k}(\tilde{z}, \tau, T)$ is zero.

For each $n = 1, 2, \cdots,$ $(\tilde{z}, \tau, T)\in C\left(E_1, \tau_n, \frac{2\pi}{\omega_0 \tau_n}\right),$ based on (4.3) and (4.4), we see

Hence, (ⅱ) fails and (ⅰ) holds. The following two lemmas help us confirm the boundedness of the projections of $C\left(E_1, \tau_n, \frac{2\pi}{\omega_0 \tau_n}\right)$ on $z$-space and $T$-space.

Lemma 4.1. For initial value $\phi = \left(\phi_{1}, \phi_{2}, \phi_{3}\right) \in Y$ with $\phi_1(\theta)+\phi_2(\theta)\leq\frac{\Lambda}{\mu}, \; \forall \theta\in [-1, 0]$, all periodic solutions of system (4.1) are uniformly bounded.

The proof is a direct result of Lemma 2.1, and hence is omitted.

Lemma 4.2. If $\mathcal {R}_0 > 1$ and $\mu+\gamma > \sigma$, then system (4.1) has no periodic solutions of period $1$.

Proof. Suppose that $z(t) = (z_1(t), z_2(t), z_3(t))$ is the periodic solution of system (4.1) with period $1$, then $z(t)$ is a periodic solution of the following ordinary differential equations

By Theorem 3, system (4.5) has a unique positive equilibrium that is globally stable, and thus no periodic solutions appear. This leads to a contradiction.

Lemma 4.1 shows that for any integer $n\geq0$, the projection of $C\left(E_1, \tau_n, \frac{2\pi}{\omega_0 \tau_n}\right)$ onto $z$-space is bounded. Lemma 4.2 implies that system (4.1) also has no periodic solutions of period $\frac{1}{n+1}$ for any $n\geq0$, Moreover, with the help of (3.8), we obtain

Therefore, the projection of $C\left(E_1, \tau_n, \frac{2\pi}{\omega_0 \tau_n}\right)$ onto $T$-space is bounded. Accordingly, the projection of $C\left(E_1, \tau_n, \frac{2\pi}{\omega_0 \tau_n}\right)$ onto $\tau$-space is unbounded.

Summarizing the above discussion, we arrive at the following result.

Theorem 4.1. Assume that $\mathcal {R}_0 > 1$, $\mu+\gamma > \sigma$ and $(p, q, r) \in W$, then for any $\tau > \tau_1$ system (1.2) has at least one nontrivial periodic solution.

5.   Numerical simulation

In this section, we carry out numerical simulations to demonstrate the theoretical results. Particularly, the global Hopf branches are computed by a Matlab package DDE-BIFTOOL developed by Engelborghs et al. ^[37,38].

For illustrative purpose, we choose the parameters of system (1.1) as follows:

Direct calculation gives $\mathcal {R}_0 = 12.1964 > 1$ and $(p, q, r)\in W.$ By using (3.8), we further obtain $\tau_0 = 5.5131$, $\tau_1 = 18.3517$, $\tau_2 = 31.1903, \cdots$. Figure 1 shows that the endemic equilibrium $E_1 = (17.8880, 1.6115, 44.7201)$ is asymptotically stable when $\tau = 4 < \tau_0$, while Figure 2 displays that $E_1$ loses stability and a Hopf bifurcation occurs when $\tau = 8 > \tau_0$. These numerical results agree with the result of Theorem 3.3.

Moreover, we depict the global Hopf branches of periodic solutions emanating from the Hopf bifurcation points $\tau_0$, $\tau_1$ and $\tau_2$. As seen in Figure 3, when $\tau_0 < \tau < \tau_1$, system (1.2) has only one periodic solution originating from $\tau_0$. When $\tau$ lies between $\tau_1$ and $\tau_2$, the periodic solutions from $\tau_0$ and $\tau_1$ coexist. With the further increase of $\tau$ and $\tau > \tau_2$, three periodic solutions originating from $\tau_0$, $\tau_1$ and $\tau_2$ coexist.

Finally, we use the delay as the bifurcation parameter to plot the bifurcation diagram. Figure 4 demonstrates the onset and global continuation of Hopf bifurcations as $\tau$ varies.

6.   Discussion

Models related to the impact of media coverage on disease spread have shown great popularity in recent years. To study the effect of media coverage on cholera transmission, we considered a cholera model with delayed media impact. We showed that the basic reproduction number $\mathcal {R}_0$ is an epidemic threshold parameter that determines the extinction and uniform persistence of the disease. However, this threshold phenomenon is not influenced by the delayed media impact, since $\mathcal {R}_0$ is independent of $m$ and $\tau$. This observation motivates us to further explore the impact of media coverage. We proved that the positive equilibrium $E_1$ is locally asymptotically stable when $\mathcal {R}_0 > 1$ and $\tau\in [0, \tau_0)$, and is unstable when $\tau > \tau_0$ (see Theorem 3.3). Furthermore, system (1.2) undergoes a Hopf bifurcation at $E_1$ along the sequence $\tau_n,$ $n = 0, 1, 2, \cdots$. To examine the onset and termination of periodic solutions bifurcated from $E_1$, we use delay as the bifurcation parameter and establish the existence of global bifurcation (see Theorem 4.1).

The local stability and the Hopf bifurcation analysis of the positive equilibrium $E_1$ are critically dependent on the existence and distribution of the roots of the cubic equation $h(x) = 0$ ^[39]. In this paper, we considered only case where $h(x) = 0$ has only one simple positive root. Therefore, it is exciting to study properties of the global Hopf branches that accompany the stability switch when $h(x) = 0$ has exactly two or three simple positive roots. Another possible project is to consider a spatial version of model (1.1). For such a model with multiple compartments, the investigation of Hopf bifurcation is challenging. We leave these topics for future research.

Use of AI tools declaration

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

Acknowledgments

This research was supported by the National Natural Science Foundation of China (No. 11971369) and the Fundamental Research Funds for the Central Universities (No. QTZX23004).

Conflict of interest

The authors declare that there is no conflict of interest.