Stability, Hopf Bifurcation and Optimal Control of Multilingual Rumor-Spreading Model with Isolation Mechanism

Abstract

The propagation of rumors on online social networks (OSNs) brings an awful lot of trouble to people’s life and society. Aiming at combating rumors spreading on OSNs, two novel rumor-propagation models without and with time delays are proposed, which combine with the influence of the immune mechanism, isolation mechanism and network structure. Firstly, we analyze the existence of rumor equilibria and obtain some existence conditions of backward bifurcation. Secondly, the local stabilities of rumor-free and rumor equilibria are proved by using the Jacobian matrix method, and some critical conditions for the existence of Hopf bifurcation are acquired by selecting critical parameters and delays as bifurcation parameters. Furthermore, an optimal control method is proposed, which can prevent the spread of rumors within an expected time period and minimize the cost of control. Finally, some numerical simulations are provided to verify the effectiveness of the proposed theoretical results.

1. Introduction

The current information age has witnessed the emergence of various communication platforms and online social network software, which provide convenience for users to obtain all kinds of information [1]. As an attractive information statement, a rumor catches people’s attention naturally. It spreads rapidly and widely on OSNs since the content is close to people’s daily life [2,3,4]. For example, the rumors that ‘eating garlic can cure cancer’ and ‘smearing ginger can cure hair loss’ spread widely on OSNs. As a result, some users believe these rumors and miss the best treatment period. Another well-known example is about the treatment of COVID-19, such as the claims that drinking Banlangen, fumigation vinegar and Shuanghuanglian oral liquid can prevent COVID-19, and the thicker the mask, the better the antiviral effect. This rumor not only reduced the public’s vigilance to the virus but also intensified the anxiety and panic to a certain extent. Hence, the spread of rumor has a negative impact on people’s lives [5,6]. It is of great significance to explore spreading rules and propose some feasible control methods to prevent the spread of rumors on OSNs.

The research on the dynamics of rumor propagation is mainly inspired by the infectious disease model due to the similarity of transmission mechanism. The systematic study of rumor spreading model began in 1965. As a milestone of the rumor-spreading model, the DK [7,8] model and MT [9] model divide the population into three categories. Each person can only be in one of the following three states: ignorant (the user who does not know the rumor); spreader (the one who propagates the rumor); and removed (the one who knows the rumor but does not spread it). However, these models [7,8,9] are only suitable for describing the spread of rumor by word of mouth. With the rapid development of social science and technology, the spread of rumors on social media has gradually replaced oral communication and become a new mainstream way. Therefore, some new research methods have been proposed by scholars. In 2001, Zanette [10] first applied complex network theory to the study of rumor spreading. After that, an increasing number of scholars began to investigate the rumor propagation model on small world networks [11,12,13], uniform networks [14,15,16] and scale-free networks [17,18].

With the widespread use of new Internet media and the occurrence of frequent emergencies in the past few years, research about rumor spreading on OSNs has also been widely carried out [19,20,21]. Meanwhile, a growing number of attention has been paid to the influence of group psychological characteristics [22,23,24,25,26], government refutation [27,28,29], the immune mechanism [30] and the educational mechanism [31]. In addition, the authors studied a rumor spreading model with homogeneous network, and introduced expert intervention strategy into the model in [32]. In [33], an impulsive rumor-blocking strategy was proposed, in which the rumor users were isolated for a limited period and the optimal control method was applied to save cost. Different from the above articles, considering the difference in regulatory capabilities of network platforms, a kind of forced silence mechanism was introduced to suppress the spreading of rumors in [34,35]. In practical application, the forced silence mechanism can be realized by forbidding the rumor disseminator. Unfortunately, these useful models were all proposed in a single language environment, rarely involving multiple language environments.

Multilingual OSNs refer that there are multiple languages on social platforms. The emergence of a large number of social network platforms has greatly facilitated people’s communication. At the same time, online translation and other software also provide convenience for information exchange between multiple languages. Therefore, it is of great practical significance to study the spread of rumors on multilingual OSNs. There are some results [36,37,38,39] about rumor spreading in multilingual environment. In [36], an I2S2R rumor-spreading model was proposed in homogeneous complex networks, and the global stability of the rumor-free and rumor equilibria was proved. In [37], a multilingual SIR rumor-spreading model with a cross-transmitted mechanism was established. In [38,39], two rumor-propagation models with heterogeneous networks were proposed. In the process of rumor spreading, it is an effective way to isolate the rumor disseminator or forbid the users. Therefore, it is meaningful to study the rumor-transmission model with the isolation mechanism. However, as far as we know, there is no research on the modeling and control of rumor propagation with the isolation mechanism in a multilingual environment.

Based on the above-mentioned factors, in this paper, we propose two new rumor-propagation models with an isolation mechanism on multilingual OSNs. In the model, all users are divided into six categories, including two types of ignorants, two types of spreaders, quarantined users and recovered users. We abbreviate this type of model as the 2I2SQR model. Moreover, we study the dynamic behaviors of the model, and propose an optimal control by using the Pontriagin maximum principle to minimize the control cost. The main contributions of this paper are as follows. Firstly, based on some basic assumptions, a rumor-spreading model with quarantine control is designed on multilingual OSNs. The existence of a rumor equilibrium is analyzed, and some conditions for the existence of backward bifurcation are obtained. Secondly, we prove the local stability of the rumor-free and rumor equilibria without and with time delays. Concurrently, by selecting the critical parameters and the time delays, some conditions for the existence of Hopf bifurcation are given. In addition, using the optimal control theory, an optimal control is obtained, which can minimize the objective function and suppress the spread of rumor.

The rest of this paper is arranged as follows. In Section 2, a novel rumor-propagation model without time delay is proposed, and the dynamic analysis of the model is given. An optimal control is proposed to optimize the isolation mechanism in Section 3. In Section 4, a rumor-propagation model with time delays is given, and the dynamics of the model is studied. The validity of the theoretical results is verified by several numerical simulations in Section 5. Section 6 concludes this paper.

2. 2I2SQR Rumor Model without Time-Delay

In this section, a rumor spreading model without time-delay is proposed on OSNs with multilingual environment, and the dynamics of the model is carefully analyzed.

2.1. Model Formulation

Users may speak multiple languages on an OSN when they come from different regions or countries. According to the language ability, we classify users into two groups. The first group is the ones who can speak both the official language and unofficial language, but they are more likely to browse the information published in the unofficial language when they are skimming through the web, and they also prefer to edit it in the unofficial language. The second group is the ones who can only speak the official language. In our model, we divide users into six categories. Ignorants-1($I_1(t)$) and Ignorants-2($I_2(t)$) represent the ratios of users who are not aware of the rumor in the first and second groups. Spreaders-1($S_1(t)$) and Spreaders-2($S_2(t)$) denote the ratios of users who know the rumor and spread it in the first and second groups, respectively. Quarantined ($Q(t)$) denotes the ratio of the user whose communicator is temporarily banned. Recovered ($R(t)$) represents the ratio of the individual who knows the rumor and no longer spreads it. Moreover, the densities satisfy $I_1(t)+I_2(t)+S_1(t)+S_2(t)+Q(t)+R(t)=1$.

Based on the mention above, the following detailed assumptions are presented before we establish the rumor-spreading model on OSNs with a multilingual environment:

(1) Assume that the immigration rates are denoted by $B_i(i=1,2)$ for the first and second groups, respectively. The removed rate is represented by d.

(2) In the process of rumor spreading, $I_1(t)$(or $I_2(t)$) will turn into $S_1(t)$ (or $S_2(t)$) with conversion rate ${\alpha }_1$ (or ${\alpha }_2$). In addition, $S_1$ will change to $S_2$ with a certain probability $\rho$ to expand the spread range of the rumor.

(3) When $I_1(t)$ (or $I_2(t))$ glance over a rumor, they will choose to disbelieve the rumor and transform into $R(t)$ with probability ${\mu }_1$(or ${\mu }_2$) because of the level of their education (identifying the rumor) or not being interested in the rumor. After spreading the rumor for a period of time, $S_1(t)$ (or $S_2(t)$) will change into $R(t)$ with probability ${\beta }_1$ (or ${\beta }_2$) on account of losing interest or other factors. Similarly, the $Q(t)$ will turn into $R(t)$ with probability c after passing the compulsory prohibition.

(4) When a rumor spreads on OSNs, the network regulatory authorities will prohibit the propagation of the rumor according to the number of communicators with a certain proportion. We describe it by a quarantine function $h_i(S_i)=\frac{r_i{S}_i(t)}{{\sigma }_i+S_i(t)}$, $(i=1,2)$, where $r_i>0$ is the probability of partition, and ${\sigma }_i>0$ is a half-saturation constant that measures the isolation delay index. If ${\sigma }_i$ is small, then the efficiency is high.

On the basis of these facts, $B_1$, $B_2$, ${\alpha }_1$, ${\alpha }_2$, ${\mu }_1$, ${\mu }_2$, ${\beta }_1$, ${\beta }_2$, d are all positive constants. The process of rumor spreading is shown in Figure 1. The dynamics of the rumor-propagation model in a multilingual environment is established as follows:

$$\begin{array}{c}\hfill {\displaystyle \left\{\begin{array}{cc}\hfill & \frac{\mathrm{d}{I}_1(t)}{\mathrm{d}t}=B_1-\langle k\rangle {\alpha }_1{I}_1(t)S_1(t)-{\mu }_1{I}_1(t)-d{I}_1(t),\hfill \\ \hfill & \frac{\mathrm{d}{I}_2(t)}{\mathrm{d}t}=B_2-\langle k\rangle {\alpha }_2{I}_2(t)S_2(t)-{\mu }_2{I}_2(t)-d{I}_2(t),\hfill \\ \hfill & \frac{\mathrm{d}{S}_1(t)}{\mathrm{d}t}=\langle k\rangle {\alpha }_1{I}_1(t)S_1(t)-\frac{r_1{S}_1(t)}{{\sigma }_1+S_1(t)}-{\beta }_1{S}_1(t)-\rho S_1(t)-d{S}_1(t),\hfill \\ \hfill & \frac{\mathrm{d}{S}_2(t)}{\mathrm{d}t}=\langle k\rangle {\alpha }_2{I}_2(t)S_2(t)+\rho S_1(t)-\frac{r_2{S}_2(t)}{{\sigma }_2+S_2(t)}-{\beta }_2{S}_2(t)-d{S}_2(t),\hfill \\ \hfill & \frac{\mathrm{d}Q(t)}{\mathrm{d}t}=\frac{r_1{S}_1(t)}{{\sigma }_1+S_1(t)}+\frac{r_2{S}_2(t)}{{\sigma }_2+S_2(t)}-cQ(t)-dQ(t),\hfill \\ \hfill & \frac{\mathrm{d}R(t)}{\mathrm{d}t}={\beta }_1{S}_1(t)+{\beta }_2{S}_2(t)+{\mu }_1{I}_1(t)+cQ(t)+{\mu }_2{I}_2(t)-dR(t),\hfill \end{array}\right.}\end{array}$$

(1)

where $\langle k\rangle ={\sum }_{k=1}^{n}kP(k)$ denotes the average degree of the individual, $P(k)$ stands for the probability of an individual with degree k and satisfies ${\sum }_{k=1}^{n}P(k)=1$. The initial conditions satisfy

$$\begin{array}{c}\hfill I_1(0)\ge 0,I_2(0)\ge 0,S_1(0)>0,S_2(0)>0,Q(0)\ge 0,R(0)\ge 0.\end{array}$$

Thus, according to the above conditions, the positive invariant set of system (1) is

$$\begin{array}{c}\hfill \Omega =\{(I_1,I_2,S_1,S_2,Q,R)\in R_6^{+}:I_1+I_2+S_1+S_2+Q+R=1\}.\end{array}$$

Since $Q(t)$ and $R(t)$ do not affect other equations of system (1), we simplify system (1) as follows:

$$\begin{array}{c}\hfill {\displaystyle \left\{\begin{array}{cc}\hfill & \frac{\mathrm{d}{I}_1(t)}{\mathrm{d}t}=B_1-\langle k\rangle {\alpha }_1{I}_1(t)S_1(t)-{\mu }_1{I}_1(t)-d{I}_1(t),\hfill \\ \hfill & \frac{\mathrm{d}{I}_2(t)}{\mathrm{d}t}=B_2-\langle k\rangle {\alpha }_2{I}_2(t)S_2(t)-{\mu }_2{I}_2(t)-d{I}_2(t),\hfill \\ \hfill & \frac{\mathrm{d}{S}_1(t)}{\mathrm{d}t}=\langle k\rangle {\alpha }_1{I}_1(t)S_1(t)-\frac{r_1{S}_1(t)}{{\sigma }_1+S_1(t)}-{\beta }_1{S}_1(t)-d{S}_1(t)-\rho S_1(t),\hfill \\ \hfill & \frac{\mathrm{d}{S}_2(t)}{\mathrm{d}t}=\langle k\rangle {\alpha }_2{I}_2(t)S_2(t)+\rho S_1(t)-\frac{r_2{S}_2(t)}{{\sigma }_2+S_2(t)}-{\beta }_2{S}_2(t)-d{S}_2(t),\hfill \end{array}\right.}\end{array}$$

(2)

with the initial conditions $I_1(0)\ge 0$, $I_2(0)\ge 0$, $S_1(0)>0$, and $S_2(0)>0$.

Remark 1. 

Although some rumor-spreading models on OSNs have been proposed in [19,20,21,34,35], they are only used to describe the spread of rumors in a single language environment. In this paper, an new multilingual rumor propagation model is established. In particular, when $I_1(t)=S_1(t)=0$ in system (1), it will be simplified into a language rumor propagation model. Therefore, the model proposed in this paper has a wider application.

2.2. Existence of Equilibria

It is easily to verify that system (2) always has a rumor-free equilibrium, which is given by

$$\begin{array}{c}\hfill E_0=\left\{\frac{B_1}{d+{\mu }_1},\frac{B_2}{d+{\mu }_2},0,0\right\}.\end{array}$$

By using the next generation matrix method [37], we calculate the basic reproduction number ${\Re }_0$, which is defined as

$$\begin{array}{c}\hfill {\Re }_0=max\{{\Re }_{01},{\Re }_{02}\},\end{array}$$

(3)

where ${\displaystyle {\Re }_{01}=\frac{\langle k\rangle {\alpha }_1{B}_1{\sigma }_1}{(d+{\mu }_1)[r_1+({\beta }_1+\rho +d){\sigma }_1]}}$, ${\displaystyle {\Re }_{02}=\frac{\langle k\rangle {\alpha }_2{B}_2{\sigma }_2}{(d+{\mu }_2)[r_2+({\beta }_2+d){\sigma }_2]}}$.

Next, we will determine the existence of a rumor equilibrium. Supposing that $E^{*}=(I_1^{*},I_2^{*},S_1^{*},S_2^{*})$ is a rumor equilibrium, then it satisfies

$$\begin{array}{c}\hfill {\displaystyle \left\{\begin{array}{cc}\hfill & B_1-\langle k\rangle {\alpha }_1{I}_1^{*}{S}_1^{*}-{\mu }_1{I}_1^{*}-d{I}_1^{*}=0,\hfill \\ \hfill & B_2-\langle k\rangle {\alpha }_2{I}_2^{*}{S}_2^{*}-{\mu }_2{I}_2^{*}-d{I}_2^{*}=0,\hfill \\ \hfill & \langle k\rangle {\alpha }_1{I}_1^{*}{S}_1^{*}-\frac{r_1{S}_1^{*}}{{\sigma }_1+S_1^{*}}-{\beta }_1{S}_1^{*}-d{S}_1^{*}-\rho S_1^{*}=0,\hfill \\ \hfill & \langle k\rangle {\alpha }_2{I}_2^{*}{S}_2^{*}+\rho S_1^{*}-\frac{r_2{S}_2^{*}}{{\sigma }_2+S_2^{*}}-{\beta }_2{S}_2^{*}-d{S}_2^{*}=0.\hfill \end{array}\right.}\end{array}$$

Therefore, it has

$$\begin{array}{ccc}\hfill & I_1^{*}=\frac{B_1}{\langle k\rangle {\alpha }_1{S}_1^{*}+{\mu }_1+d},\hfill & \hfill I_2^{*}=\frac{B_2}{\langle k\rangle {\alpha }_2{S}_2^{*}+{\mu }_2+d}.\end{array}$$

Moreover, $S_1^{*}$ and $S_2^{*}$ are determined by the following equations:

$$\begin{array}{cc}\hfill & a_2{(S_1^{*})}^2+a_1{S}_1^{*}+a_0=0,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill & b_3{(S_2^{*})}^3+b_2{(S_2^{*})}^2+b_1{S}_2^{*}+b_0=0,\hfill \end{array}$$

(5)

where

$$\begin{array}{cc}\hfill a_0& =({\mu }_1+d)[r_1+{\sigma }_1({\beta }_1+d+\rho )](1-{\Re }_{01}),\hfill \\ \hfill a_1& =r_1{\alpha }_1\langle k\rangle +{\alpha }_1{\sigma }_1({\beta }_1+d+\rho )\langle k\rangle +({\mu }_1+d)({\beta }_1+d+\rho )-B_1{\alpha }_1\langle k\rangle ,\hfill \\ \hfill a_2& ={\alpha }_1({\beta }_1+d+\rho )\langle k\rangle ,\hfill \\ \hfill b_0& =\rho ({\mu }_2+d){\sigma }_2{S}_1^{*},\hfill \\ \hfill b_1& =[r_2+({\beta }_2+d){\sigma }_2]({\mu }_2+d)({\Re }_{02}-1)+[\rho {\alpha }_2{\sigma }_2\langle k\rangle +\rho (\mu +d)]S_1^{*},\hfill \\ \hfill b_2& ={\alpha }_2{B}_2\langle k\rangle +{\alpha }_2\rho S_1^{*}\langle k\rangle -{\alpha }_2\langle k\rangle [r_2+({\beta }_2+d){\sigma }_2]-({\beta }_2+d)({\mu }_2+d),\hfill \\ \hfill b_3& =-\langle k\rangle {\alpha }_2({\beta }_2+d).\hfill \end{array}$$

Firstly, we will solve the positive solutions of Equation (4) with respect to $S_1^{*}$. For convenience, we denote

$$\begin{array}{cc}\hfill & {\widehat{\Re }}_{01}=1-\frac{{[(r_1-B_1){\alpha }_1\langle k\rangle +({\alpha }_1{\sigma }_1\langle k\rangle +{\mu }_1+d)({\beta }_1+d+\rho )]}^2}{4{\alpha }_1\langle k\rangle ({\beta }_1+d+\rho )({\mu }_1+d)[r_1+{\sigma }_1({\beta }_1+d+\rho )]}.\hfill \end{array}$$

By simple calculation, the following results can be obtained.

(i).

If ${\Re }_{01}>1$, Equation (4) has a positive solution.

(ii).

If ${\Re }_{01}=1$, Equation (4) is transformed into

$$\begin{array}{c}\hfill a_2{(S_1^{*})}^2+a_1{S}_1^{*}=0.\end{array}$$

(6)

Obviously, $S_1^{*}=0$ or $S_1^{*}=-\frac{a_1}{a_2}>0$ if $a_1<0$.

(iii).

If ${\Re }_{01}<1$ and $a_1<0$, the following results can be easily verified.

(1). If ${\Re }_{01}>{\widehat{\Re }}_{01}$, Equation (4) has two positive roots.

(2). If ${\Re }_{01}={\widehat{\Re }}_{01}$, Equation (4) has two equal positive roots $-\frac{a_1}{2a_2}$.

(3). If ${\Re }_{01}<{\widehat{\Re }}_{01}$, Equation (4) has no positive root.

Next, we will consider the positive solutions of Equation (5) with respect to $S_2^{*}$.

Case (1). It is clear that $b_0=0$ when $S_1^{*}=0$, then Equation (5) is transformed into

$$\begin{array}{c}\hfill b_3{(S_2^{*})}^3+{\tilde{b}}_2{(S_2^{*})}^2+{\tilde{b}}_1{S}_2^{*}=0,\end{array}$$

(7)

where ${\tilde{b}}_1=[r_2+({\beta }_2+d){\sigma }_2]({\mu }_2+d)({\Re }_{02}-1)$, ${\tilde{b}}_2={\alpha }_2{B}_2\langle k\rangle -{\alpha }_2\langle k\rangle [r_2+({\beta }_2+d){\sigma }_2]-({\beta }_2+d)({\mu }_2+d)$, $b_3=-\langle k\rangle {\alpha }_2({\beta }_2+d)$.

Apparently, $S_2^{*}=0$ or $b_3{(S_2^{*})}^2+{\tilde{b}}_2(S_2^{*})+{\tilde{b}}_1=0$. Let us discuss the solutions of the following equation:

$$\begin{array}{c}\hfill b_3{(S_2^{*})}^2+{\tilde{b}}_2(S_2^{*})+{\tilde{b}}_1=0.\end{array}$$

(8)

By calculation, it has

$$\begin{array}{c}\hfill {\tilde{b}}_1>0\iff {\Re }_{02}>1,{\tilde{b}}_1=0\iff {\Re }_{02}=1,{\tilde{b}}_1<0\iff {\Re }_{02}<1.\end{array}$$

We denote

$$\begin{array}{c}\hfill {\widehat{\Re }}_{02}=1-\frac{{[{\alpha }_2{B}_2\langle k\rangle -{\alpha }_2\langle k\rangle [r_2+({\beta }_2+d){\sigma }_2]-({\beta }_2+d)({\mu }_2+d)]}^2}{4\langle k\rangle {\alpha }_2({\beta }_2+d)[[r_2+({\beta }_2+d){\sigma }_2]({\mu }_2+d)]}.\end{array}$$

Then, the following conclusions are obtained.

(i).

If ${\Re }_{02}>1$, Equation (8) has a positive solution.

(ii).

If ${\Re }_{02}=1$, Equation (8) has a solution $S_2^{*}=0$ or $S_2^{*}=-\frac{{\tilde{b}}_2}{b_3}>0$ if and only if ${\tilde{b}}_2>0$.

(iii).

If ${\Re }_{02}<1$ and ${\tilde{b}}_2>0$, one of the following three cases holds

(1). If ${\Re }_{02}>{\widehat{\Re }}_{02}$, Equation (8) has two positive roots.

(2). If ${\Re }_{02}={\widehat{\Re }}_{02}$, Equation (8) has two identical positive roots $S_2^{*}=-\frac{{\tilde{b}}_2}{2b_3}$.

(3). If ${\Re }_{02}<{\widehat{\Re }}_{02}$, Equation (8) does not have a positive root.

Case (2). When $S_1^{*}>0$, it has $b_0>0$. We denote

$$\begin{array}{c}\hfill G(S_2^{*})=b_3{(S_2^{*})}^3+b_2{(S_2^{*})}^2+b_1{S}_2^{*}+b_0=0.\end{array}$$

(9)

Then, the discriminant of the cubic polynomial $G(S_2^{*})$ is given by

$$\begin{array}{c}\hfill D=b_1^2{b}_2^2+18b_0{b}_1{b}_2{b}_3-4b_0{b}_2^3-4b_1^2{b}_3-27b_0^2{b}_3^2.\end{array}$$

By simple calculation, $G(S_2^{*})$ has one real root and two complex roots if $D<0$. $G(S_2^{*})$ has three real roots (at least two of which are equal) if $D=0$, and G has three distinct real roots if $D>0$. In order to analyze the roots of Equation (5), the following lemmas are given:

Lemma 1. 

Suppose that ${\alpha }_2{\sigma }_2\langle k\rangle >{\mu }_2+d$, then $b_2<0$ if $b_1<0$.

Proof. 

According to ${\alpha }_2{\sigma }_2\langle k\rangle >{\mu }_2+d$, we have

$$\begin{array}{cc}\hfill b_2=& {\alpha }_2{B}_2\langle k\rangle +{\alpha }_2\rho S_1^{*}\langle k\rangle -{\alpha }_2\langle k\rangle [r_2+({\beta }_2+d){\sigma }_2]-({\beta }_2+d)({\mu }_2+d)\hfill \\ \hfill =& {\alpha }_2{B}_2\langle k\rangle +{\alpha }_2\rho S_1^{*}\langle k\rangle -{\alpha }_2{r}_2\langle k\rangle -({\beta }_2+d)({\sigma }_2{\alpha }_2\langle k\rangle +{\mu }_2+d),\hfill \\ \hfill b_1=& {\alpha }_2{\sigma }_2{B}_2\langle k\rangle -[r_2+({\beta }_2+d){\sigma }_2]({\mu }_2+d)+[\rho {\alpha }_2{\sigma }_2\langle k\rangle +\rho (\mu +d)]S_1^{*}\hfill \\ \hfill =& {\alpha }_2{\sigma }_2{B}_2\langle k\rangle -r_2({\mu }_2+d)-{\sigma }_2({\beta }_2+d)({\mu }_2+d)+\rho {\sigma }_2{\alpha }_2\langle k\rangle S_1^{*}+\rho (\mu +d)S_1^{*}\hfill \\ \hfill \ge & {\sigma }_2[{\alpha }_2{B}_2\langle k\rangle -r_2{\alpha }_2\langle k\rangle -({\beta }_2+d)({\mu }_2+d)+\rho {\alpha }_2{S}_1^{*}\langle k\rangle ]+\rho (\mu +d)S_1^{*}\hfill \\ \hfill \ge & {\sigma }_2[{\alpha }_2{B}_2\langle k\rangle -r_2{\alpha }_2\langle k\rangle -({\beta }_2+d)({\alpha }_2{\sigma }_2\langle k\rangle +{\mu }_2+d)+\rho {\alpha }_2{S}_1^{*}\langle k\rangle ]+\rho (\mu +d)S_1^{*}\hfill \\ \hfill \ge & {\sigma }_2{b}_2.\hfill \end{array}$$

Since ${\sigma }_2>0$, it is easy to obtain that $b_2<0$ if $b_1<0$. □

Lemma 2. 

If ${\alpha }_2{\sigma }_2\langle k\rangle >{\mu }_2+d$ holds, Equation (9) has one positive root.

Proof. 

By deriving Equation (5), we can obtain

$$\begin{array}{c}\hfill G{(S_2^{*})}^{{}^{\prime }}=3b_3{(S_2^{*})}^2+2b_2{S}_2^{*}+b_1.\end{array}$$

Denote the roots of $G{(S_2^{*})}^{{}^{\prime }}=0$ by $x^{\pm }$, then one of the following three cases holds:

(i).

$D<0$. Using the fact that $G(0)=b_0>0$ and ${lim}_{S_2^{*}\to \infty }G(S_2^{*})=-\infty$, it follows that the real root is positive.

(ii).

$D=0$. Solving the roots $x^{+}$ and $x^{-}$, we find that Equation (5) has one positive root.

(iii).

$D>0$. The analysis method is similar to (ii), and Equation (5) has one positive root.

Note that there is only one positive root in the previous three cases. Therefore, Equation (5) only has one positive solution $S_2^{*}$ if ${\alpha }_2{\sigma }_2\langle k\rangle >{\mu }_2+d$ and $b_0>0$.

Based on the above analysis, the following theorem is given to ensure the existence of the rumor equilibria. □

Theorem 1. 

For system (2), ${\Re }_{01}$ and ${\Re }_{02}$ are defined as Equation (3). The following conclusions about the existence of the rumor equilibria are presented.

(1). For the case of $b_0>0$ and ${\alpha }_2{\sigma }_2\langle k\rangle >{\mu }_2+d$:

If ${\Re }_{01}>1$, there is a unique rumor equilibrium point $E_1^{*}$.

If ${\Re }_{01}=1$ and $a_1<0$, system (2) admits a unique rumor equilibrium point $E_1^{*}$.

If ${\widehat{\Re }}_{01}<{\Re }_{01}<1$ and $a_1<0$, there are two rumor equilibrium points $E_1^{*}$ and $E_2^{*}$ in system (2), and no rumor equilibrium point when $a_1\ge 0$.

If ${\widehat{\Re }}_{01}={\Re }_{01}<1$ and $a_1<0$, system (2) admits a unique rumor equilibrium point $E_1^{*}$, and no rumor equilibrium point when $a_1\ge 0$.

If ${\Re }_{01}<{\widehat{\Re }}_{01}<1$, system (2) has no rumor equilibrium point. (2). For the case of $S_1^{*}=0(b_0=0)$:

If ${\Re }_{02}>1$, there is a unique rumor equilibrium point $E_1^{*}$ in system (2).

If ${\Re }_{02}=1$ and ${\tilde{b}}_2>0$, system (2) admits a unique rumor equilibrium points $E_1^{*}$.

If ${\widehat{\Re }}_{02}<{\Re }_{02}<1$ and ${\tilde{b}}_2>0$, there are two rumor equilibria $E_1^{*}$ and $E_2^{*}$ in system (2), and no rumor equilibrium point when ${\tilde{b}}_2\le 0$.

If ${\widehat{\Re }}_{02}={\Re }_{02}<1$ and ${\tilde{b}}_2>0$, system (2) admits a unique rumor equilibrium point $E_1^{*}$, and no rumor equilibrium point when ${\tilde{b}}_2\le 0$.

If ${\Re }_{02}<{\widehat{\Re }}_{02}<1$, system (2) has no rumor equilibrium.

We find that the conditions for the existence of a rumor equilibrium are very complicated. Under different conditions, system (2) may have different equilibria. Based on the above theorem, it can be seen that there exist rumor equilibria in system (2) when ${\Re }_{01}<1$ and ${\Re }_{02}<1$. For instance, when ${\widehat{\Re }}_{01}=0.89$ and ${\widehat{\Re }}_{02}=0.86$, the backward bifurcation graph of the system (2) is shown in Figure 2.

2.3. Stability and Hopf Bifurcation of the Equilibria

In order to study the evolution tend of the solution of system (2), we need to analyze the stability of the equilibrium. Based on Theorem 1, system (2) may have multiple equilibria. It is very difficult to analyze the global stabilities of the rumor-free and rumor equilibria. Therefore, we will discuss the local stability and Hopf bifurcation of the equilibria.

Theorem 2. 

For system (2), the rumor-free equilibrium $E_0$ is locally asymptotically stable when ${\Re }_0<1$ and unstable when ${\Re }_0>1$.

Proof. 

The Jacobian matrix of system (2) at $E_0$ is given as

$$\begin{array}{c}\hfill \mathcal{J}(E_0)=\left(\begin{array}{cccc}-({\mu }_1+d)& 0& -\frac{{\alpha }_1{B}_1\langle k\rangle }{d+{\mu }_1}& 0\\ 0& -({\mu }_2+d)& 0& -\frac{{\alpha }_2{B}_2\langle k\rangle }{d+{\mu }_2}\\ 0& 0& \frac{{\alpha }_1{B}_1\langle k\rangle }{d+{\mu }_1}-\frac{[r_1+{\sigma }_1({\beta }_1+d+\rho )]}{{\sigma }_1}& 0\\ 0& 0& \rho & \frac{{\alpha }_2{B}_2\langle k\rangle }{d+{\mu }_2}-\frac{[r_2+{\sigma }_2({\beta }_2+d)]}{{\sigma }_2}\end{array}\right).\end{array}$$

Then, the characteristic equation is equivalent to

$$\begin{array}{cc}\hfill [\lambda & -\frac{[r_1+{\sigma }_1({\beta }_1+d+\rho )]}{{\sigma }_1}({\Re }_{01}-1)]\left[\lambda -\frac{[r_2+{\sigma }_2({\beta }_2+d)]}{{\sigma }_2}({\Re }_{02}-1)\right]\hfill \\ \hfill & \times (\lambda +{\mu }_1+d)(\lambda +{\mu }_2+d)=0.\hfill \end{array}$$

Thus, ${\lambda }_1=-({\mu }_1+d)$, ${\lambda }_2=-({\mu }_2+d)$, ${\lambda }_3=\frac{[r_1+{\sigma }_1({\beta }_1+d+\rho )]}{{\sigma }_1}({\Re }_{01}-1)$, ${\lambda }_4=\frac{[r_2+{\sigma }_2({\beta }_2+d)]}{{\sigma }_2}$ $\times ({\Re }_{02}-1)$. Since ${\Re }_0<1$, it has ${\lambda }_3<0$ and ${\lambda }_4<0$. Based on the stability theory, we can conclude that $E_0$ is locally asymptotically stable. Otherwise, $E_0$ is unstable if ${\Re }_0>1$. □

Lemma 3

([40]). For any $a>0$ and $b>0$, if $\frac{\mathrm{d}x(t)}{\mathrm{d}t}\ge b-ax(t)$ for $t\ge 0$ and $x(0)>0$, it has ${lim}_{t\to +\infty }infx(t)\ge \frac{b}{a}$; if $\frac{\mathrm{d}x(t)}{\mathrm{d}t}\le b-ax(t)$ for $t\ge 0$ and $x(0)>0$, it has ${lim}_{t\to +\infty }supx(t)\le \frac{b}{a}$.

Theorem 3. 

For system (2), the rumor-free equilibrium $E_0$ is globally asymptotically stable if ${\Re }_{01}+\frac{r_1}{(1+{\sigma }_1)[r_1+{\sigma }_1({\beta }_1+\rho +d)]}<1$ and ${\Re }_{02}+\frac{r_2}{(1+{\sigma }_2)[r_2+{\sigma }_2({\beta }_2+d)]}<1$.

Proof. 

Based on system (2), it has

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{I}_1(t)}{\mathrm{d}t}=& B_1-\langle k\rangle {\alpha }_1{I}_1(t)S_1(t)-{\mu }_1{I}_1(t)-d{I}_1(t)\hfill \\ \hfill \le & B_1-({\mu }_1+d)I_1(t).\hfill \end{array}$$

By Lemma 3, one can obtain that

$$\begin{array}{c}\hfill \underset{t\to +\infty }{lim}sup{I}_1(t)\le \frac{B_1}{{\mu }_1+d}.\end{array}$$

(10)

Hence, for any ${\epsilon }_1>0$, there exists $t_1>0$ such that $I_1(t)\le \frac{B_1}{{\mu }_1+d}+{\epsilon }_1$ for $t>t_1$. Then, when $t>t_1$, one has

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{S}_1(t)}{\mathrm{d}t}\le & \left(\langle k\rangle {\alpha }_1(\frac{B_1}{{\mu }_1+d}+{\epsilon }_1)-(\frac{r_1}{{\sigma }_1+1}+{\beta }_1+\rho +d)\right)S_1(t)\hfill \\ \hfill =& \frac{r_1+({\beta }_1+d+\rho ){\sigma }_1}{{\sigma }_1}[{\Re }_{01}-1+\frac{r_1}{(1+{\sigma }_1)[r_1+{\sigma }_1({\beta }_1+\rho +d)]}\hfill \\ \hfill & +\frac{{\sigma }_1\langle k\rangle {\alpha }_1{\epsilon }_1}{r_1+({\beta }_1+d+\rho ){\sigma }_1}]S_1(t).\hfill \end{array}$$

Denote $\omega ({\epsilon }_1)={\Re }_{01}-1+\frac{r_1}{(1+{\sigma }_1)[r_1+{\sigma }_1({\beta }_1+\rho +d)]}+\frac{{\sigma }_1\langle k\rangle {\alpha }_1{\epsilon }_1}{r_1+({\beta }_1+d+\rho ){\sigma }_1}$. According to comparison theorem, one can obtain that

$$\begin{array}{c}\hfill S_1(t)\le S_1(0)exp\left\{\frac{r_1+({\beta }_1+d+\rho ){\sigma }_1}{{\sigma }_1}\omega ({\epsilon }_1)t\right\}.\end{array}$$

Because ${\Re }_{01}+\frac{r_1}{(1+{\sigma }_1)[r_1+{\sigma }_1({\beta }_1+\rho +d)]}<1$ and ${\epsilon }_1>0$ is an arbitrarily small real number, we can chose ${\epsilon }_1>0$ such that $\omega ({\epsilon }_1)<0$. Therefore, ${lim}_{t\to +\infty }{S}_1(t)=0$.

Next, we prove that ${lim}_{t\to +\infty }{S}_2(t)=0$. Similarly, for any ${\epsilon }_2>0$, there exists $t_2>0$ such that $I_2(t)\le \frac{B_2}{{\mu }_2+d}+{\epsilon }_2$ and $S_2(t)\le {\epsilon }_2$ for $t>t_2$. Then, for $t>t_2$, it follows that

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{S}_2(t)}{\mathrm{d}t}\le & \left(\langle k\rangle {\alpha }_2(\frac{B_2}{{\mu }_2+d}+{\epsilon }_2)+\rho {\epsilon }_2-(\frac{r_2}{{\sigma }_2+1}+{\beta }_2+d)\right)S_2(t)\hfill \\ \hfill =& \frac{r_2+({\beta }_2+d){\sigma }_2}{{\sigma }_2}[{\Re }_{02}-1+\frac{r_2}{(1+{\sigma }_2)[r_2+{\sigma }_2({\beta }_2+d)]}\hfill \\ \hfill & +\frac{{\sigma }_2(\langle k\rangle {\alpha }_2+\rho ){\epsilon }_2}{r_2+({\beta }_2+d){\sigma }_2}]S_2(t).\hfill \end{array}$$

Denote $\varpi ({\epsilon }_2)={\Re }_{02}-1+\frac{r_2}{(1+{\sigma }_2)[r_2+{\sigma }_2({\beta }_2+d)]}+\frac{{\sigma }_2(\langle k\rangle {\alpha }_2+\rho ){\epsilon }_2}{r_2+({\beta }_2+d){\sigma }_2}$. By using comparison theorem, one can obtain that

$$\begin{array}{c}\hfill S_2(t)\le S_2(0)exp\left\{\frac{r_2+({\beta }_2+d){\sigma }_2}{{\sigma }_2}\varpi ({\epsilon }_2)t\right\}.\end{array}$$

Since ${\Re }_{02}+\frac{r_2}{(1+{\sigma }_2)[r_2+{\sigma }_2({\beta }_2+d)]}<1$ and ${\epsilon }_2>0$ is an arbitrarily small real number, then we can chose ${\epsilon }_2>0$ such that $\varpi ({\epsilon }_2)<0$. Therefore, ${lim}_{t\to +\infty }{S}_2(t)=0$.

Due to ${lim}_{t\to +\infty }{S}_1(t)=0$, then for any ${\epsilon }_3>0$, there exists $t_3>0$ such that $S_1(t)\le {\epsilon }_3$ for $t>t_3$. Therefore, it has

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{I}_1(t)}{\mathrm{d}t}=& B_1-\langle k\rangle {\alpha }_1{I}_1(t)S_1(t)-{\mu }_1{I}_1(t)-d{I}_1(t)\hfill \\ \hfill \ge & B_1-(\langle k\rangle {\alpha }_1{\epsilon }_3+{\mu }_1+d)I_1(t).\hfill \end{array}$$

By Lemma 3, one can get that ${lim}_{t\to +\infty }{I}_1(t)\ge \frac{B_1}{\langle k\rangle {\alpha }_1{\epsilon }_3+{\mu }_1+d}$. Let ${\epsilon }_3\to 0$, it follows that

$$\begin{array}{c}\hfill \underset{t\to +\infty }{lim}inf{I}_1(t)\ge \frac{B_1}{{\mu }_1+d}.\end{array}$$

(11)

In combination with (10) and (11), it is clear that ${lim}_{t\to +\infty }{I}_1(t)=\frac{B_1}{{\mu }_1+d}$. Similar to the analysis of $I_1(t)$, we can obtain that ${lim}_{t\to +\infty }{I}_2(t)=\frac{B_2}{{\mu }_2+d}$. This proves that the equilibrium $E_0$ is globally asymptotically stable. □

In order to facilitate the subsequent analysis process, we denote

$$\begin{array}{cc}\hfill H_1(S_1^{*})=& {\alpha }_1{S}_1^{*}\langle k\rangle +\frac{r_1{\sigma }_1}{{({\sigma }_1+S_1^{*})}^2}+{\Theta }_1+{\Theta }_3-{\alpha }_1{I}_1^{*}\langle k\rangle ,\hfill \\ \hfill G_1(S_1^{*})=& (\frac{r_1{\sigma }_1}{{\sigma }_1+S_1^{*}}+{\Theta }_1)({\alpha }_1{S}_1^{*}\langle k\rangle +{\Theta }_3)-{\alpha }_1{I}_1^{*}{\Theta }_3\langle k\rangle ,\hfill \\ \hfill H_2(S_2^{*})=& {\alpha }_2{S}_2^{*}\langle k\rangle +\frac{r_2{\sigma }_2}{{({\sigma }_2+S_2^{*})}^2}+{\Theta }_2+{\Theta }_4-{\alpha }_2{I}_2^{*}\langle k\rangle ,\hfill \\ \hfill G_2(S_2^{*})=& (\frac{r_2{\sigma }_2}{{\sigma }_2+S_2^{*}}+{\Theta }_2)({\alpha }_2{S}_2^{*}\langle k\rangle +{\Theta }_4)-{\alpha }_2{I}_2^{*}{\Theta }_4\langle k\rangle ,\hfill \end{array}$$

where ${\Theta }_1={\beta }_1+d+\rho$, ${\Theta }_2={\beta }_2+d$, ${\Theta }_3={\mu }_1+d$, and ${\Theta }_4={\mu }_2+d$.

Theorem 4. 

Suppose that $E_1^{*}$ is the rumor equilibrium of system (2). If $H_i(S_i^{*})>0$ and $G_i(S_i^{*})>0$, for $i=1,2$, then $E_1^{*}$ is locally asymptotically stable. Otherwise, $E_1^{*}$ is unstable.

Proof. 

The Jacobian matrix of system (2) at $E_1^{*}$ is expressed as

$$\begin{array}{c}\hfill \mathcal{J}(E_1^{*})=\left(\begin{array}{cccc}-{\alpha }_1{S}_1^{*}\langle k\rangle -{\Theta }_3& 0& -{\alpha }_1{I}_1^{*}\langle k\rangle & 0\\ 0& -{\alpha }_2{S}_2^{*}\langle k\rangle -{\Theta }_4& 0& -{\alpha }_2{I}_2^{*}\langle k\rangle \\ {\alpha }_1{S}_1^{*}\langle k\rangle & 0& {\alpha }_1{I}_1^{*}\langle k\rangle -\frac{r_1{\sigma }_1}{{({\sigma }_1+S_1^{*})}^2}-{\Theta }_1& 0\\ 0& {\alpha }_2{S}_2^{*}\langle k\rangle & \rho & {\alpha }_2{I}_2^{*}\langle k\rangle -\frac{r_2{\sigma }_2}{{({\sigma }_2+S_2^{*})}^2}-{\Theta }_2\end{array}\right).\end{array}$$

The characteristic equation for matrix $\mathcal{J}(E_1^{*})$ is as follows:

$$\begin{array}{c}\hfill ({\lambda }^2+H_1(S_1^{*})\lambda +G_1(S_1^{*}))({\lambda }^2+H_2(S_2^{*})\lambda +G_2(S_2^{*}))=0.\end{array}$$

(12)

Denote

$$\begin{array}{c}\hfill {\lambda }^2+H_1(S_1^{*})\lambda +G_1(S_1^{*})=0,\end{array}$$

(13)

$$\begin{array}{c}\hfill {\lambda }^2+H_2(S_2^{*})\lambda +G_2(S_2^{*})=0.\end{array}$$

(14)

According to the Routh–Hurwitz criterion, the rumor equilibrium $E_1^{*}$ is locally asymptotically stable when $H_i(S_i^{*})>0$ and $G_i(S_i^{*})>0$, $i=1,2$. Otherwise, it is unstable. □

Remark 2. 

In Theorem 4, we give a criterion for the local stability of the rumor equilibrium $E_1^{*}$. From Theorem 1, it can be seen that there are many cases above the rumor equilibrium in system (2). In either case, the criterion of Theorem 4 is applicable since it is of a general form. When system (2) has two equilibria, each point may be locally stable under different conditions. Therefore, it is very meaningful to give the threshold from stable to unstable for the equilibrium of system (2).

Based on the definition of $H_i(S_i^{*})$ and $G_i(S_i^{*})$, $i=1,2$, we let

$$\begin{array}{cc}\hfill {\widehat{r}}_1=& \frac{({\alpha }_1{I}_1^{*}\langle k\rangle -{\alpha }_1{S}_1^{*}\langle k\rangle -{\Theta }_1-{\Theta }_3){({\sigma }_1+S_1^{*})}^2}{{\sigma }_1},\hfill \\ \hfill {\check{r}}_1=& (\frac{{\alpha }_1{I}_1^{*}\langle k\rangle {\Theta }_3}{{\alpha }_1{S}_1^{*}\langle k\rangle +{\Theta }_3}-{\Theta }_1)\frac{({\sigma }_1+S_1^{*})}{{\sigma }_1},\hfill \\ \hfill {\widehat{r}}_2=& \frac{({\alpha }_2{I}_2^{*}\langle k\rangle -{\alpha }_2{S}_2^{*}\langle k\rangle -{\Theta }_2-{\Theta }_4){({\sigma }_2+S_2^{*})}^2}{{\sigma }_2},\hfill \\ \hfill {\check{r}}_2=& (\frac{{\alpha }_2{I}_2^{*}\langle k\rangle {\Theta }_4}{{\alpha }_2{S}_2^{*}\langle k\rangle +{\Theta }_4}-{\Theta }_2)\frac{({\sigma }_2+S_2^{*})}{{\sigma }_2},\hfill \end{array}$$

and

$$\begin{array}{c}\hfill {\overline{r}}_1=max\{{\widehat{r}}_1,{\check{r}}_1\},{\overline{r}}_2=max\{{\widehat{r}}_2,{\check{r}}_2\}.\end{array}$$

Then, the following theorem is given.

Theorem 5. 

For system (2), the following conclusions are true.

(1). If $r_1>{\overline{r}}_1$ and $r_2>{\overline{r}}_2$, the rumor equilibrium $E_1^{*}$ is locally asymptotically stable, otherwise, $E_1^{*}$ is unstable.

(2). System (2) has a Hopf bifurcation at $r_1={\overline{r}}_1$ or $r_2={\overline{r}}_2$.

Proof. 

Let $H_i(S_i^{*})=0$ for $i=1,2$, then one can get $r_i={\widehat{r}}_i$. By simple calculation, we can obtain that $H_i(S_i^{*})>0$ if $r_i>{\widehat{r}}_i$. Similarly, let $G_i(S_i^{*})=0$ for $i=1,2$, then it follows $r_i={\check{r}}_i$, and $G_i(S_i^{*})>0$ if $r_i>{\check{r}}_i$.

For the case of ${\widehat{r}}_1>{\check{r}}_1$, consider $r_1$ in the neighborhood of ${\widehat{r}}_1$, then Equation (13) has two positive roots, which are given by

$$\begin{array}{c}\hfill {\lambda }_{1,2}=\frac{-H_1(S_1^{*})\pm {(H_1{(S_1^{*})}^2-4G_1(S_1^{*}))}^{\frac{1}{2}}}{2}.\end{array}$$

The derivative of $r_1$ can be obtained:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\lambda }_{1,2}}{\mathrm{d}{r}_1}=-\frac{{\sigma }_1}{2{({\sigma }_1+S_1^{*})}^2}\pm \frac{1}{4\sqrt{H_1{(S_1^{*})}^2-4G_1(S_1^{*})}}[2H_1(S_1^{*})\frac{{\sigma }_1}{{({\sigma }_1+S_1^{*})}^2}-4\frac{{\sigma }_1({\alpha }_1{S}_1^{*}\langle k\rangle +{\Theta }_3)}{{\sigma }_1+S_1^{*}}].\end{array}$$

When $r_1={\widehat{r}}_1$, one has

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\lambda }_{1,2}}{\mathrm{d}{r}_1}=-\frac{{\sigma }_1}{2{({\sigma }_1+S_1^{*})}^2}\pm \frac{1}{2i\sqrt{G_1(S_1^{*})}}\frac{{\sigma }_1({\alpha }_1{S}_1^{*}\langle k\rangle +{\Theta }_3)}{{\sigma }_1+S_1^{*}}.\end{array}$$

Obviously,

$$\begin{array}{c}\hfill \frac{\mathrm{d}Re({\lambda }_{1,2})}{\mathrm{d}{r}_1}{|}_{r_1={\widehat{r}}_1}=Re(\frac{\mathrm{d}{\lambda }_{1,2}}{\mathrm{d}{r}_1}{|}_{r_1={\widehat{r}}_1})=-\frac{{\sigma }_1}{2{({\sigma }_1+S_1^{*})}^2}\ne 0.\end{array}$$

Therefore, there is a Hopf bifurcation at $r_1={\widehat{r}}_1$.

Similarly, for the case of ${\widehat{r}}_2>{\check{r}}_2$, we can obtain

$$\begin{array}{c}\hfill \frac{\mathrm{d}Re({\lambda }_{3,4})}{\mathrm{d}{r}_2}{|}_{r_2={\widehat{r}}_2}=Re(\frac{\mathrm{d}{\lambda }_{3,4}}{\mathrm{d}{r}_2}{|}_{r_2={\widehat{r}}_2})=-\frac{{\sigma }_2}{2{({\sigma }_2+S_2^{*})}^2}\ne 0.\end{array}$$

Hence, there is a Hopf bifurcation at $r_2={\widehat{r}}_2$.

For the case of ${\check{r}}_1>{\widehat{r}}_1$, consider $r_1$ in the neighborhood of ${\check{r}}_1$, so Equation (13) has one positive root, which is presented as follows:

$$\begin{array}{c}\hfill {\lambda }_5=\frac{-H_1(S_1^{*})+{(H_1{(S_1^{*})}^2-4G_1(S_1^{*}))}^{\frac{1}{2}}}{2}.\end{array}$$

The derivative of $r_1$ can be obtained:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\lambda }_5}{\mathrm{d}{r}_1}=-\frac{{\sigma }_1}{2{({\sigma }_1+S_1^{*})}^2}+\frac{1}{4\sqrt{H_1{(S_1^{*})}^2-4G_1(S_1^{*})}}[2H_1(S_1^{*})\frac{{\sigma }_1}{{({\sigma }_1+S_1^{*})}^2}-4\frac{{\sigma }_1({\alpha }_1{S}_1^{*}\langle k\rangle +{\Theta }_3)}{{\sigma }_1+S_1^{*}}].\end{array}$$

When $r_1={\check{r}}_1$, one has $G_1(S_1^{*})=0$, $H_1(S_1^{*})>0$,

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\lambda }_5}{\mathrm{d}{r}_1}{|}_{r_1={\check{r}}_1}=-\frac{{\sigma }_1}{H_1(S_1^{*})({\sigma }_1+S_1^{*})}({\alpha }_1{S}_1^{*}\langle k\rangle +{\Theta }_3)\ne 0.\end{array}$$

Hence, there exists a Hopf bifurcation at $r_1={\check{r}}_1$.

Similarly, for the case of ${\check{r}}_2>{\widehat{r}}_2$, we can obtain

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\lambda }_6}{\mathrm{d}{r}_2}{|}_{r_2={\check{r}}_2}=-\frac{{\sigma }_2}{H_2(S_2^{*})({\sigma }_2+S_2^{*})}({\alpha }_2{S}_2^{*}\langle k\rangle +{\Theta }_4)\ne 0.\end{array}$$

Therefore, there exists a Hopf bifurcation at $r_2={\check{r}}_2$. □

Remark 3. 

In this paper, we introduce the isolation mechanism to the rumor-spreading model and study the local stability of the equilibrium. Our purpose is to give the condition of rumor extinction by choosing the control parameters. As we all know, the basic reproduction number ${\Re }_0$ is an important reference index. However, we find that system (2) may contain a rumor equilibrium and it is locally stable when ${\Re }_0<1$. Therefore, it is difficult to give the criterion of rumor extinction under any initial condition.

3. Optimal Control Model

In system (2), the isolation mechanism is applied to suppress the spread of rumor, in which the isolation rates $r_1$ and $r_2$ are constants. In this section, we will consider the time-varying isolation mechanism. In order to reduce the control cost and achieve the desired control effect, we use the optimal control method to optimize the control cost. The model with the isolation mechanism is as follows:

$$\begin{array}{c}\hfill {\displaystyle \left\{\begin{array}{cc}\hfill & \frac{\mathrm{d}{I}_1(t)}{\mathrm{d}t}=B_1-\langle k\rangle {\alpha }_1{I}_1(t)S_1(t)-{\mu }_1{I}_1(t)-d{I}_1(t),\hfill \\ \hfill & \frac{\mathrm{d}{I}_2(t)}{\mathrm{d}t}=B_2-\langle k\rangle {\alpha }_2{I}_2(t)S_2(t)-{\mu }_2{I}_2(t)-d{I}_2(t),\hfill \\ \hfill & \frac{\mathrm{d}{S}_1(t)}{\mathrm{d}t}=\langle k\rangle {\alpha }_1{I}_1(t)S_1(t)-\frac{r_1(t)S_1(t)}{{\sigma }_1+S_1(t)}-{\beta }_1{S}_1(t)-d{S}_1(t)-\rho S_1(t),\hfill \\ \hfill & \frac{\mathrm{d}{S}_2(t)}{\mathrm{d}t}=\langle k\rangle {\alpha }_2{I}_2(t)S_2(t)+\rho S_1(t)-\frac{r_2(t)S_2(t)}{{\sigma }_2+S_2(t)}-{\beta }_2{S}_2(t)-d{S}_2(t),\hfill \\ \hfill & \frac{\mathrm{d}Q(t)}{\mathrm{d}t}=\frac{r_1(t)S_1(t)}{{\sigma }_1+S_1(t)}+\frac{r_2(t)S_2(t)}{{\sigma }_2+S_2(t)}-cQ(t)-dQ(t),\hfill \\ \hfill & \frac{\mathrm{d}R(t)}{\mathrm{d}t}={\beta }_1{S}_1(t)+{\beta }_2{S}_2(t)+{\mu }_1{I}_1(t)+cQ(t)+{\mu }_2{I}_2(t)-dR(t),\hfill \end{array}\right.}\end{array}$$

(15)

where the control $r_1(t)$ and $r_2(t)$ are isolation rates of the social platform. Considering the number of rumor spreaders and control cost, the following objective function is established:

$$\begin{array}{c}\hfill J(r_1(t),r_2(t))={\int }_0^T[u_1{S}_1(t)+u_2{S}_2(t)+u_3{r}_1^2(t)+u_4{r}_2^2(t)]\mathrm{d}t.\end{array}$$

(16)

where ${\mu }_1,{\mu }_2,{\mu }_3$ and ${\mu }_4$ are positive numbers, which are trade-off parameters among these items. T is the expected control time.

The feasible region of $r_1(t)$ and $r_2(t)$ is $U=\{(r_1(t),r_2(t))|0\le r_1(t)\le 1,0\le r_2(t)\le 1,t\in (0,T]\}$. Optimal control $r_1^{*}$ and $r_2^{*}$ satisfy $J(r_1^{*},r_2^{*})=min\{J(r_1(t),r_2(t)):(r_1(t),r_2(t))\in U\}$.

In order to obtain the optimal control, we construct the Lagrangian function:

$$\begin{array}{c}\hfill L(S_1(t),S_2(t),r_1(t),r_2(t))=u_1{S}_1(t)+u_2{S}_2(t)+u_3{r}_1^2(t)+u_4{r}_2^2(t).\end{array}$$

The Hamiltonian function is defined as

$$\begin{array}{cc}\hfill H(& I_i(t),S_i(t),Q(t),R(t),r_i(t),{\lambda }_j(t))=L(S_i(t),r_i(t))+\lambda (t)P(I_i(t),S_i(t),Q(t),R(t),r_i(t)),\hfill \end{array}$$

where $i=1,2$, $j=1,2,\dots ,6$. $\lambda (t)={({\lambda }_1(t),{\lambda }_2(t),\dots ,{\lambda }_6(t))}^T$,$P(I_i(t),S_i(t),Q(t),R(t),$ $r_i(t))={(\frac{\mathrm{d}{I}_1(t)}{\mathrm{d}t},\frac{\mathrm{d}{I}_1(t)}{\mathrm{d}t},\frac{\mathrm{d}{I}_1(t)}{\mathrm{d}t},\frac{\mathrm{d}{I}_2(t)}{\mathrm{d}t},\frac{\mathrm{d}{S}_1(t)}{\mathrm{d}t},\frac{\mathrm{d}{S}_2(t)}{\mathrm{d}t})}^T.$

Using Pontryagin’s maximum principle, the following theorem is given.

Theorem 6. 

Let ${\widehat{I}}_1$, ${\widehat{I}}_2$, ${\widehat{S}}_1$, ${\widehat{S}}_2$, $\widehat{Q}$ and $\widehat{R}$ be the optimal states of system (15) under the optimal control $(r_1^{*},r_2^{*})$. Then, there exist adjoint variables ${\lambda }_j(t)$, $j=1,\cdots ,6$ satisfying

$$\begin{array}{c}\hfill {\displaystyle \left\{\begin{array}{cc}\hfill \frac{\mathrm{d}{\lambda }_1(t)}{\mathrm{d}t}=& {\lambda }_1(t)(\langle k\rangle {\alpha }_1{\widehat{S}}_1+{\mu }_1+d)-{\lambda }_3(t)\langle k\rangle {\alpha }_1{\widehat{S}}_1-{\lambda }_6(t){\mu }_1,\hfill \\ \hfill \frac{\mathrm{d}{\lambda }_2(t)}{\mathrm{d}t}=& {\lambda }_2(t)(\langle k\rangle {\alpha }_2{\widehat{S}}_2+{\mu }_2+d)-{\lambda }_4(t)\langle k\rangle {\alpha }_2{\widehat{S}}_2-{\lambda }_6(t){\mu }_2,\hfill \\ \hfill \frac{\mathrm{d}{\lambda }_3(t)}{\mathrm{d}t}=& -u_1+{\lambda }_1\langle k\rangle {\alpha }_1{\widehat{I}}_1-{\lambda }_3(t)[\langle k\rangle {\alpha }_1{\widehat{I}}_1-\frac{r_1(t){\sigma }_1}{{({\sigma }_1+{\widehat{S}}_1)}^2}-d-{\beta }_1-\rho ]\hfill \\ \hfill & -{\lambda }_4(t)\rho -{\lambda }_5(t)\frac{r_1(t){\sigma }_1}{{({\sigma }_1+{\widehat{S}}_1)}^2}-{\lambda }_6(t){\beta }_1,\hfill \\ \hfill \frac{\mathrm{d}{\lambda }_4(t)}{\mathrm{d}t}=& -u_2+{\lambda }_2\langle k\rangle {\alpha }_2{\widehat{I}}_2-{\lambda }_4(t)[\langle k\rangle {\alpha }_2{\widehat{I}}_2-\frac{r_2(t){\sigma }_2}{{({\sigma }_2+{\widehat{S}}_2)}^2}-{\beta }_2-d]\hfill \\ \hfill & -{\lambda }_5(t)\frac{r_2(t){\sigma }_2}{{({\sigma }_2+{\widehat{S}}_2)}^2}-{\lambda }_6(t){\beta }_2,\hfill \\ \hfill \frac{\mathrm{d}{\lambda }_5(t)}{\mathrm{d}t}=& {\lambda }_5(t)(c+d)-{\lambda }_6(t)c,\hfill \\ \hfill \frac{\mathrm{d}{\lambda }_6(t)}{\mathrm{d}t}=& {\lambda }_6(t)d,\hfill \end{array}\right.}\end{array}$$

with the transversality conditions ${\lambda }_j(T)=0$. Optimal control $r_1^{*}$ and $r_2^{*}$ are given by

$$\begin{array}{cc}\hfill & r_1^{*}=max\{min\{\frac{({\lambda }_3-{\lambda }_5){\widehat{S}}_1}{2u_3({\sigma }_1+{\widehat{S}}_1)},0\},r_1^{max}\},\hfill \\ \hfill & r_2^{*}=max\{min\{\frac{({\lambda }_4-{\lambda }_5){\widehat{S}}_2}{2u_4({\sigma }_2+{\widehat{S}}_2)},0\},r_2^{max}\}.\hfill \end{array}$$

Proof. 

Let $I_1(t)={\widehat{I}}_1$, $I_2(t)={\widehat{I}}_2$, $S_1(t)={\widehat{S}}_1$, $S_2(t)={\widehat{S}}_2$, $Q(t)=\widehat{Q}$ and $R(t)=\widehat{R}$, using Pontryagin’s maximum principle [41], we obtain the following adjoint equations:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\lambda }_1(t)}{\mathrm{d}t}=-\frac{\partial H}{\partial I_1(t)}=& {\lambda }_1(t)(\langle k\rangle {\alpha }_1{\widehat{S}}_1+{\mu }_1+d)-{\lambda }_3(t)\langle k\rangle {\alpha }_1{\widehat{S}}_1-{\lambda }_6(t){\mu }_1,\hfill \\ \hfill \frac{\mathrm{d}{\lambda }_2(t)}{\mathrm{d}t}=-\frac{\partial H}{\partial I_2(t)}=& {\lambda }_2(t)(\langle k\rangle {\alpha }_2{\widehat{S}}_2+{\mu }_2+d)-{\lambda }_4(t)\langle k\rangle {\alpha }_2{\widehat{S}}_2-{\lambda }_6(t){\mu }_2,\hfill \\ \hfill \frac{\mathrm{d}{\lambda }_3(t)}{\mathrm{d}t}=-\frac{\partial H}{\partial S_1(t)}=& -u_1+{\lambda }_1\langle k\rangle {\alpha }_1{\widehat{I}}_1-{\lambda }_3(t)[\langle k\rangle {\alpha }_1{\widehat{I}}_1-\frac{r_1(t){\sigma }_1}{{({\sigma }_1+{\widehat{S}}_1)}^2}-d-{\beta }_1-\rho ]\hfill \\ \hfill & -{\lambda }_4(t)\rho -{\lambda }_5(t)\frac{r_1(t){\sigma }_1}{{({\sigma }_1+{\widehat{S}}_1)}^2}-{\lambda }_6(t){\beta }_1,\hfill \\ \hfill \frac{\mathrm{d}{\lambda }_4(t)}{\mathrm{d}t}=-\frac{\partial H}{\partial S_2(t)}=& -u_2+{\lambda }_2\langle k\rangle {\alpha }_2{\widehat{I}}_2-{\lambda }_4(t)[\langle k\rangle {\alpha }_2{\widehat{I}}_2-\frac{r_2(t){\sigma }_2}{{({\sigma }_2+{\widehat{S}}_2)}^2}-{\beta }_2-d]\hfill \\ \hfill & -{\lambda }_5(t)\frac{r_2(t){\sigma }_2}{{({\sigma }_2+{\widehat{S}}_2)}^2}-{\lambda }_6(t){\beta }_2,\hfill \\ \hfill \frac{\mathrm{d}{\lambda }_5(t)}{\mathrm{d}t}=-\frac{\partial H}{\partial Q(t)}=& {\lambda }_5(t)(c+d)-{\lambda }_6(t)c,\hfill \\ \hfill \frac{\mathrm{d}{\lambda }_6(t)}{\mathrm{d}t}=-\frac{\partial H}{\partial R(t)}=& {\lambda }_6(t)d.\hfill \end{array}$$

By the optimality condition, the derivative of the Hamiltonian function with respect to $r_1(t)$ and $r_2(t)$ is as follows:

$$\begin{array}{cc}\hfill & \frac{\partial H}{\partial r_1(t)}{|}_{r_1(t)=r_1^{*}}=2u_3{r}_1(t)+\frac{({\lambda }_5(t)-{\lambda }_3(t)){\widehat{S}}_1}{{\sigma }_1+{\widehat{S}}_1}=0,\hfill \\ \hfill & \frac{\partial H}{\partial r_2(t)}{|}_{r_2(t)=r_2^{*}}=2u_4{r}_2(t)+\frac{({\lambda }_5(t)-{\lambda }_4(t)){\widehat{S}}_2}{{\sigma }_2+{\widehat{S}}_2}=0.\hfill \end{array}$$

Then, we obtain the optimal control

$$\begin{array}{c}\hfill r_1^{*}=\frac{({\lambda }_3(t)-{\lambda }_5(t)){\widehat{S}}_1}{2u_3({\sigma }_1+{\widehat{S}}_1)},r_2^{*}=\frac{({\lambda }_4(t)-{\lambda }_5(t)){\widehat{S}}_2}{2u_4({\sigma }_2+{\widehat{S}}_2)}.\end{array}$$

By combining the properties of bounded set U, the interval of $r_1^{*}$ and $r_2^{*}$ are shown in the following form:

$$\begin{array}{cc}\hfill & r_1^{*}=max\left\{min\{\frac{({\lambda }_3(t)-{\lambda }_5(t)){\widehat{S}}_1}{2u_3({\sigma }_1+{\widehat{S}}_1)},0\},r_1^{max}\right\},\hfill \\ \hfill & r_2^{*}=max\left\{min\{\frac{({\lambda }_4(t)-{\lambda }_5(t)){\widehat{S}}_2}{2u_4({\sigma }_2+{\widehat{S}}_2)},0\},r_2^{max}\right\}.\hfill \end{array}$$

□

Remark 4. 

In the optimal control, the control $r_1(t)$ and $r_2(t)$ are time-varying and associated with objective function (14). For different expected control time T and objective function, the control rates $r_1(t)$ and $r_2(t)$ may be different. This can be verified by the simulation in the later section.

4. 2I2SQR Rumor Model with Time-Delays

4.1. Model Formulation

In general, when the influence of rumor becomes greater, the social platform will take control measures to isolate the rumor disseminators. Because it takes a certain time to search for rumor disseminators on the network platform, the implementation of this measure sometimes lags behind the rumor-spreading process. Therefore, system (1) incorporating time-delays is formulated as follows:

$$\begin{array}{c}\hfill {\displaystyle \left\{\begin{array}{cc}\hfill & \frac{\mathrm{d}{I}_1(t)}{\mathrm{d}t}=B_1-\langle k\rangle {\alpha }_1{I}_1(t)S_1(t)-{\mu }_1{I}_1(t)-d{I}_1(t),\hfill \\ \hfill & \frac{\mathrm{d}{I}_2(t)}{\mathrm{d}t}=B_2-\langle k\rangle {\alpha }_2{I}_2(t)S_2(t)-{\mu }_2{I}_2(t)-d{I}_2(t),\hfill \\ \hfill & \frac{\mathrm{d}{S}_1(t)}{\mathrm{d}t}=\langle k\rangle {\alpha }_1{I}_1(t)S_1(t)-\frac{r_1{S}_1(t-{\tau }_1)}{{\sigma }_1+S_1(t-{\tau }_1)}-{\beta }_1{S}_1(t)-d{S}_1(t)-\rho S_1(t),\hfill \\ \hfill & \frac{\mathrm{d}{S}_2(t)}{\mathrm{d}t}=\langle k\rangle {\alpha }_2{I}_2(t)S_2(t)+\rho S_1(t)-\frac{r_2{S}_2(t-{\tau }_2)}{{\sigma }_2+S_2(t-{\tau }_2)}-{\beta }_2{S}_2(t)-d{S}_2(t),\hfill \\ \hfill & \frac{\mathrm{d}Q(t)}{\mathrm{d}t}=\frac{r_1{S}_1(t-{\tau }_1)}{{\sigma }_1+S_1(t-{\tau }_1)}+\frac{r_2{S}_2(t-{\tau }_2)}{{\sigma }_2+S_2(t-{\tau }_2)}-cQ(t)-dQ(t),\hfill \\ \hfill & \frac{\mathrm{d}R(t)}{\mathrm{d}t}={\beta }_1{S}_1(t)+{\beta }_2{S}_2(t)+{\mu }_1{I}_1(t)+cQ(t)+{\mu }_2{I}_2(t)-dR(t).\hfill \end{array}\right.}\end{array}$$

(17)

where ${\tau }_1,{\tau }_2>0$ represent the time delays, and other parameters and variables are consistent with the explanation in system (1). The initial conditions of system (17) are given by

$$\begin{array}{cc}\hfill & I_i(\theta )>0,S_i(\theta )>0,Q(\theta )\ge 0,R(\theta )\ge 0,\theta \in [-\tau ,0],i=1,2,\hfill \end{array}$$

where $\tau =max\{{\tau }_1,{\tau }_2\}$.

Similar to the previous analysis, we simplify system (17) as follows:

$$\begin{array}{c}\hfill {\displaystyle \left\{\begin{array}{cc}\hfill & \frac{\mathrm{d}{I}_1(t)}{\mathrm{d}t}=B_1-\langle k\rangle {\alpha }_1{I}_1(t)S_1(t)-{\mu }_1{I}_1(t)-d{I}_1(t),\hfill \\ \hfill & \frac{\mathrm{d}{I}_2(t)}{\mathrm{d}t}=B_2-\langle k\rangle {\alpha }_2{I}_2(t)S_2(t)-{\mu }_2{I}_2(t)-d{I}_2(t),\hfill \\ \hfill & \frac{\mathrm{d}{S}_1(t)}{\mathrm{d}t}=\langle k\rangle {\alpha }_1{I}_1(t)S_1(t)-\frac{r_1{S}_1(t-{\tau }_1)}{{\sigma }_1+S_1(t-{\tau }_1)}-{\beta }_1{S}_1(t)-d{S}_1(t)-\rho S_1(t),\hfill \\ \hfill & \frac{\mathrm{d}{S}_2(t)}{\mathrm{d}t}=\langle k\rangle {\alpha }_2{I}_2(t)S_2(t)+\rho S_1(t)-\frac{r_2{S}_2(t-{\tau }_2)}{{\sigma }_2+S_2(t-{\tau }_2)}-{\beta }_2{S}_2(t)-d{S}_2(t).\hfill \end{array}\right.}\end{array}$$

(18)

Because the equilibria of system (18) are the same as that of system (2), which does not depend on the time delays, then Theorem 1 is still established in this section. Hence, we omit the theoretical analysis of the existence of equilibria.

4.2. Stability and Hopf Bifurcation of Equilibria

In this subsection, we will first discuss the stability and Hopf bifurcation of rumor-free equilibrium $E_0$, and then analyze the dynamics of the rumor equilibrium.

The Jacobian matrix of system (18) at $E_0$ is expressed as follows:

$$\begin{array}{c}\hfill \mathcal{J}(E_0)=\left(\begin{array}{cccc}-({\mu }_1+d)& 0& -\frac{{\alpha }_1{B}_1\langle k\rangle }{d+{\mu }_1}& 0\\ 0& -({\mu }_2+d)& 0& -\frac{{\alpha }_2{B}_2\langle k\rangle }{d+{\mu }_2}\\ 0& 0& -\frac{r_1{e}^{-\lambda {\tau }_1}}{{\sigma }_1}-A_1& 0\\ 0& 0& \rho & -\frac{r_2{e}^{-\lambda {\tau }_2}}{{\sigma }_2}-A_2\end{array}\right),\end{array}$$

where $A_1=({\beta }_1+d+\rho )-\frac{{\alpha }_1{B}_1\langle k\rangle }{d+{\mu }_1}$, $A_2={\beta }_2+d-\frac{{\alpha }_2{B}_2\langle k\rangle }{d+{\mu }_2}$. Accordingly, the characteristic equation is

$$\begin{array}{c}\hfill (\lambda +{\mu }_1+d)(\lambda +{\mu }_2+d)(\lambda +A_1+\frac{r_1}{{\sigma }_1}{e}^{-\lambda {\tau }_1})(\lambda +A_2+\frac{r_2}{{\sigma }_2}{e}^{-\lambda {\tau }_2})=0.\end{array}$$

(19)

It is clear that ${\lambda }_1=-({\mu }_1+d)<0$, ${\lambda }_2=-({\mu }_2+d)<0$, and other eigenvalues are determined by the following equations:

$$\begin{array}{c}\hfill \lambda +A_1+\frac{r_1}{{\sigma }_1}{e}^{-\lambda {\tau }_1}=0,\lambda +A_2+\frac{r_2}{{\sigma }_2}{e}^{-\lambda {\tau }_2}=0.\end{array}$$

Let

$$\begin{array}{cc}\hfill & f(\lambda )=\lambda +A_1+\frac{r_1}{{\sigma }_1}{e}^{-\lambda {\tau }_1}=0,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & g(\lambda )=\lambda +A_2+\frac{r_2}{{\sigma }_2}{e}^{-\lambda {\tau }_2}=0.\hfill \end{array}$$

(21)

By Theorem 2, we know that all solutions of Equation (19) are negative when ${\tau }_1={\tau }_2=0$ and ${\Re }_0<1$. When ${\tau }_1,{\tau }_2>0$, assume that $\lambda =i\omega (\omega >0)$ is a solution of (19). Substituting $\lambda =i\omega$ into (19), the following two cases are discussed:

Case (1). If $f(\lambda )=0$, then it yields

$$\begin{array}{c}\hfill i\omega +A_1+\frac{r_1}{{\sigma }_1}(cos\omega {\tau }_1-isin\omega {\tau }_1)=0.\end{array}$$

Separating the real parts and imaginary parts, one has

$$\begin{array}{c}\hfill A_1+\frac{r_1}{{\sigma }_1}cos\omega {\tau }_1=0,\omega -\frac{r_1}{{\sigma }_1}sin\omega {\tau }_1=0.\end{array}$$

Squaring and adding the above equations gives

$$\begin{array}{c}\hfill {\omega }^2={(\frac{r_1}{{\sigma }_1})}^2-A_1^2=\frac{r_1+{\sigma }_1({\beta }_1+d+\rho )}{{\sigma }_1}(1-{\Re }_{01})\left(\frac{r_1}{{\sigma }_1}+\frac{\langle k\rangle {\alpha }_1{B}_1}{d+{\mu }_1}-({\beta }_1+d+\rho )\right).\end{array}$$

(22)

Assume that $\frac{r_1}{{\sigma }_1}+\frac{\langle k\rangle {\alpha }_1{B}_1}{d+{\mu }_1}>{\beta }_1+d+\rho$, then there exist a positive real number ${\omega }_{10}$ such that Equation (20) has purely imaginary roots $\lambda =\pm i{\omega }_{10}({\omega }_{10}>0)$. Then, according to the above analysis, we have

$$\begin{array}{c}\hfill cos({\omega }_{10}{\tau }_1)=-\frac{A_1{\sigma }_1}{r_1}=\frac{{\sigma }_1}{r_1}\left[\frac{\langle k\rangle {\alpha }_1{B}_1}{d+{\mu }_1}-({\beta }_1+d+\rho )\right],\end{array}$$

and ${\tau }_{1i}=\frac{1}{{\omega }_{10}}arccos(\frac{{\sigma }_1}{r_1}[\frac{\langle k\rangle {\alpha }_1{B}_1}{d+{\mu }_1}-({\beta }_1+d+\rho )])+\frac{2\pi i}{{\omega }_{10}}$, $i=0,1,2,\cdots .$

Additionally, the derivative of Equation (20) with respect to ${\tau }_1$ yields

$$\begin{array}{c}\hfill {(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_1})}^{-1}=\frac{1}{\lambda \frac{r_1}{{\sigma }_1}{e}^{-\lambda {\tau }_1}}-\frac{{\tau }_1}{\lambda }=\frac{1}{-\lambda (\lambda +A_1)}-\frac{{\tau }_1}{\lambda }.\end{array}$$

(23)

Taking $\lambda =i{\omega }_{10}$, ${\tau }_1={\tau }_{10}$ into (23), and separating the real part and the imaginary part, one has

$$\begin{array}{c}\hfill Re{(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_1})}^{-1}{|}_{\lambda =i{\omega }_{10},{\tau }_1={\tau }_{10}}=\frac{1}{A_1^2+{\omega }_{10}^2}>0.\end{array}$$

Case (2). If $g(\lambda )=0$, the analysis process is similar to Case (1). Hence, we obtain that

$$\begin{array}{c}\hfill {\omega }^2={(\frac{r_2}{{\sigma }_2})}^2-A_2^2=\frac{r_2+{\sigma }_2({\beta }_2+d)}{{\sigma }_2}(1-{\Re }_{02})\left(\frac{r_2}{{\sigma }_2}+\frac{\langle k\rangle {\alpha }_2{B}_2}{d+{\mu }_2}-({\beta }_2+d)\right).\end{array}$$

(24)

When $\frac{r_2}{{\sigma }_2}+\frac{\langle k\rangle {\alpha }_2{B}_2}{d+{\mu }_2}>{\beta }_2+d$, based on (24), there exists a positive real number equation ${\omega }_{20}$ such that Equation (21) has a purely imaginary root $\lambda =i{\omega }_{20}$. Through the analysis, it has

$$\begin{array}{c}\hfill {\tau }_{2j}=\frac{1}{{\omega }_{20}}arccos\left(\frac{{\sigma }_2}{r_2}\left[\frac{\langle k\rangle {\alpha }_2{B}_2}{d+{\mu }_2}-({\beta }_2+d)\right]\right)+\frac{2\pi j}{{\omega }_{20}},j=0,1,2,\cdots .\end{array}$$

(25)

and

$$\begin{array}{c}\hfill Re{(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_2})}^{-1}{|}_{\lambda =i{\omega }_{20},{\tau }_2={\tau }_{20}}=\frac{1}{A_2^2+{\omega }_{20}^2}>0.\end{array}$$

In summary, the following theorem is given.

Theorem 7. 

For system (18), if ${\Re }_0<1$, $\frac{r_1}{{\sigma }_1}+\frac{\langle k\rangle {\alpha }_1{B}_1}{d+{\mu }_1}>{\beta }_1+d+\rho$ and $\frac{r_2}{{\sigma }_2}+\frac{\langle k\rangle {\alpha }_2{B}_2}{d+{\mu }_2}>{\beta }_2+d$ are satisfied, the following statements hold.

(1). When ${\tau }_1\in [0,{\tau }_{10})$ and ${\tau }_2\in [0,{\tau }_{20})$, the rumor-free equilibrium $E_0$ is locally asymptotically stable.

(2). When ${\tau }_1>{\tau }_{10}$ or ${\tau }_2>{\tau }_{20}$, the rumor-free equilibrium $E_0$ is unstable. Moreover, system (16) has a Hopf bifurcation at $E_0$ when ${\tau }_1={\tau }_{10}$ or ${\tau }_2={\tau }_{20}$.

Next, we will discuss the stability and Hopf bifurcation of the rumor equilibrium $E_1^{*}$. The Jacobian matrix of system (18) at $E_1^{*}$ is given as follows:

$$\begin{array}{c}\hfill \mathcal{J}(E_1^{*})=\left(\begin{array}{cccc}-{\Omega }_1& 0& -{\alpha }_1{I}_1^{*}\langle k\rangle & 0\\ 0& -{\Omega }_2& 0& -{\alpha }_2{I}_2^{*}\langle k\rangle \\ {\alpha }_1{S}_1^{*}\langle k\rangle & 0& -{\Omega }_3-\frac{r_1{\sigma }_1{e}^{-\lambda {\tau }_1}}{{({\sigma }_1+S_1^{*})}^2}& 0\\ 0& {\alpha }_2{S}_2^{*}\langle k\rangle & \rho & -{\Omega }_4-\frac{r_2{\sigma }_2{e}^{-\lambda {\tau }_2}}{{({\sigma }_2+S_2^{*})}^2}\end{array}\right),\end{array}$$

where ${\Omega }_1={\alpha }_1{S}_1^{*}\langle k\rangle +{\mu }_1+d$, ${\Omega }_2={\alpha }_2{S}_2^{*}\langle k\rangle +{\mu }_2+d$, ${\Omega }_3=-{\alpha }_1{I}_1^{*}\langle k\rangle +({\beta }_1+d+\rho )$, and ${\Omega }_4=-{\alpha }_2{I}_2^{*}\langle k\rangle +({\beta }_2+d)$.

The characteristic equation for matrix $\mathcal{J}(E_1^{*})$ is expressed as

$$\begin{array}{cc}\hfill & \left({\lambda }^2+({\Omega }_1+{\Omega }_3)\lambda +C+(\lambda +{\Omega }_1)U_1{e}^{-\lambda {\tau }_1}\right)\hfill \\ \hfill \times & \left({\lambda }^2+({\Omega }_2+{\Omega }_4)\lambda +D+(\lambda +{\Omega }_2)U_2{e}^{-\lambda {\tau }_2}\right)=0,\hfill \end{array}$$

where $C=({\alpha }_1{S}_1^{*}\langle k\rangle +{\mu }_1+d)({\beta }_1+d+\rho )-{\alpha }_1{I}_1^{*}\langle k\rangle ({\mu }_1+d)$, $D=({\alpha }_2{S}_2^{*}\langle k\rangle +{\mu }_2+d)({\beta }_2+d)-{\alpha }_2{I}_2^{*}\langle k\rangle ({\mu }_2+d)$, $U_1=\frac{r_1{\sigma }_1}{{({\sigma }_1+S_1^{*})}^2}$, $U_2=\frac{r_2{\sigma }_2}{{({\sigma }_2+S_2^{*})}^2}$. The stability of the matrix $\mathcal{J}(E_1^{*})$ is determined by the solution of the following equations

$$\begin{array}{cc}\hfill & {\lambda }^2+({\Omega }_1+{\Omega }_3)\lambda +C+(\lambda +{\Omega }_1)U_1{e}^{-\lambda {\tau }_1}=0,\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & {\lambda }^2+({\Omega }_2+{\Omega }_4)\lambda +D+(\lambda +{\Omega }_2)U_2{e}^{-\lambda {\tau }_2}=0.\hfill \end{array}$$

(27)

For the case of ${\tau }_1={\tau }_2=0$, according to the proof of Theorem 4, the rumor equilibrium $E_1^{*}$ is locally asymptotically stable under the conditions of $H_i(S_i^{*})>0$ and $G_i(S_i^{*})>0.$

For the case of ${\tau }_1,{\tau }_2>0$, assuming that the solutions of Equations (26) and (27) are $\lambda =i{\omega }_1$ and $\lambda =i{\omega }_2$, respectively. Then, separating the real part and the imaginary part of Equations (26) and (27), the following equations can be obtained:

$$\begin{array}{c}\hfill {\displaystyle \left\{\begin{array}{cc}\hfill & {\omega }_1{U}_{1}sin{\omega }_1{\tau }_1+{\Omega }_1{U}_{1}cos{\omega }_1{\tau }_1={\omega }_1^2-C,\hfill \\ \hfill & {\omega }_1{U}_{1}cos{\omega }_1{\tau }_1-{\Omega }_1{U}_{1}sin{\omega }_1{\tau }_1=-({\Omega }_1+{\Omega }_3){\omega }_1,\hfill \end{array}\right.}\end{array}$$

(28)

$$\begin{array}{c}\hfill {\displaystyle \left\{\begin{array}{cc}\hfill & {\omega }_2{U}_{2}sin{\omega }_2{\tau }_2+{\Omega }_2{U}_{2}cos{\omega }_2{\tau }_2={\omega }_2^2-D,\hfill \\ \hfill & {\omega }_2{U}_{2}cos{\omega }_2{\tau }_2-{\Omega }_2{U}_{2}sin{\omega }_2{\tau }_2=-({\Omega }_2+{\Omega }_4){\omega }_2.\hfill \end{array}\right.}\end{array}$$

(29)

By adding the squares of the two equations in Equations (28) and (29), and letting $x={\omega }_1^2$ and $y={\omega }_2^2$, it has

$$\begin{array}{cc}\hfill & x^2+M_{1}x+M_2=0,\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill & y^2+L_{1}y+L_2=0,\hfill \end{array}$$

(31)

where $M_1={({\Omega }_1+{\Omega }_3)}^2-2C-U_1^2$, $M_2=C^2-{\Omega }_1^2{U}_1^2$, $L_1={({\Omega }_2+{\Omega }_4)}^2-2D-U_2^2$, $L_2=D^2-{\Omega }_2^2{U}_2^2$. Denote $G(x)=x^2+M_{1}x+M_2$, $W(y)=y^2+L_{1}y+L_2$.

For Equation (30), the following conclusions are true:

$({\mathcal{H}}_1)$. If $M_2<0$, Equation (30) has a positive root.

$({\mathcal{H}}_2)$. If $M_2\ge 0$, and $M_1\ge 0$, Equation (30) does not have any positive roots.

$({\mathcal{H}}_3)$. If $M_2>0$, $M_1<0$, and $M_1^2-4M_2\ge 0$, Equation (30) has two positive roots.

$({\mathcal{H}}_4)$. If $M_2>0$, $M_1<0$, and $M_1^2-4M_2<0$, Equation (30) does not have any positive roots.

Similar to Equation (31), the following conclusions are true:

$({\mathcal{K}}_1)$. If $L_2<0$, Equation (31) has a positive root.

$({\mathcal{K}}_2)$. If $L_2>0$, and $L_1<0$, and $L_1^2-4L_2<0$, Equation (31) does not have any positive roots.

$({\mathcal{K}}_3)$. If $L_2>0$, $L_1<0$, and $L_1^2-4L_2\ge 0$, Equation (31) has two positive roots.

$({\mathcal{K}}_4)$. If $L_2\ge 0$, and $L_1\ge 0$, Equation (31) does not have any positive roots.

For the case of $({\mathcal{H}}_1)$, Equation (30) has a positive root $x_0={\omega }_{10}^2$. By Equation (28), we have

$$\begin{array}{c}\hfill cos({\omega }_1{\tau }_1)=\frac{{\omega }_1^2-C-({\Omega }_1+{\Omega }_3){\omega }_1}{({\omega }_1{U}_1+{\Omega }_1{U}_1)},\end{array}$$

then

$$\begin{array}{c}\hfill {\tau }_{{\mathcal{H}}_{10}}^j=\frac{1}{{\omega }_{10}}arccos(\frac{{\omega }_{10}^2-C-({\Omega }_1+{\Omega }_3){\omega }_{10}}{{\omega }_{10}{U}_1+{\Omega }_1{U}_1})+\frac{2\pi j}{{\omega }_{10}},\end{array}$$

(32)

where $j=0,1,2,\dots$, $\pm {\omega }_{10}$ is a pair of pure imaginary roots of Equation (30). Further, taking the derivative of Equation (26) with respect to ${\tau }_1$, we have

$$\begin{array}{c}\hfill {(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_1})}^{-1}=\frac{2\lambda +{\Omega }_1+{\Omega }_3}{-[{\lambda }^2+({\Omega }_1+{\Omega }_3)\lambda +C]\lambda }+\frac{1}{\lambda (\lambda +{\Omega }_1)}-\frac{{\tau }_1}{\lambda }.\end{array}$$

Then

$$\begin{array}{c}\hfill Re{(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_1})}^{-1}{|}_{\lambda =i{\omega }_{10},{\tau }_1={\tau }_{{\mathcal{H}}_{10}}^0}=\frac{2{({\omega }_{10})}^2-2C+{({\Omega }_1+{\Omega }_3)}^2-U_1^2}{({\omega }_{10}^2+{\Omega }_1^2)U_1^2}=\frac{\mathcal{G}({\omega }_{10})}{({\omega }_{10}^2+{\Omega }_1^2)U_1^2}.\end{array}$$

(33)

Since $M_2<0$, $x_0={\omega }_{10}^2>0$ is a positive root of (30), it has $\mathcal{G}({\omega }_{10})>0$. Hence

$$\begin{array}{c}\hfill Re{(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_1})}^{-1}{|}_{\lambda =i{\omega }_{10},{\tau }_1={\tau }_{{\mathcal{H}}_{10}}^0}=\frac{\mathcal{G}({\omega }_{10})}{({\omega }_{10}^2+{\Omega }_1^2)U_1^2}>0.\end{array}$$

Similarly, for the case of $({\mathcal{K}}_1)$, if $L_2<0$, Equation (31) has a positive root $y_0={\omega }_{20}^2$. Then, we can obtain

$$\begin{array}{c}\hfill {\tau }_{{\mathcal{K}}_{10}}^j=\frac{1}{{\omega }_{20}}arccos(\frac{{\omega }_{20}^2-D-({\Omega }_2+{\Omega }_4){\omega }_{20}}{{\omega }_{20}{U}_2+{\Omega }_2{U}_2})+\frac{2\pi j}{{\omega }_{20}},j=0,1,2,\cdots ,\end{array}$$

(34)

$$\begin{array}{c}\hfill Re{(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_2})}^{-1}{|}_{\lambda =i{\omega }_{20},{\tau }_2={\tau }_{{\mathcal{K}}_{10}}^0}=\frac{2{({\omega }_{20})}^2-2D+{({\Omega }_2+{\Omega }_4)}^2-U_2^2}{({({\omega }_{20})}^2+{\Omega }_2^2)U_2^2}=\frac{\mathcal{W}({\omega }_{20})}{({({\omega }_{20})}^2+{\Omega }_2^2)U_2^2}.\end{array}$$

(35)

Similar to the analysis of $({\mathcal{H}}_1)$, one has

$$\begin{array}{c}\hfill Re{(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_2})}^{-1}{|}_{\lambda =i{\omega }_{20},{\tau }_2={\tau }_{{\mathcal{K}}_{10}}^0}=\frac{\mathcal{W}({\omega }_{20})}{({({\omega }_{20})}^2+{\Omega }_2^2)U_2^2}>0.\end{array}$$

Therefore, the following theorem is given.

Theorem 8. 

For system (18), if $H_i(S_i^{*})>0$ and $G_i(S_i^{*})>0$, $i=1,2$, the following results are satisfied.

(1). If one of $({\mathcal{H}}_i)$ and one of $({\mathcal{K}}_i)$ for $i=2,4$ are satisfied, the rumor equilibrium point $E_1^{*}$ is locally asymptotically stable for all ${\tau }_1,{\tau }_2\ge 0$.

(2). If $({\mathcal{H}}_1)$ and $({\mathcal{K}}_1)$ are satisfied, the rumor equilibrium point $E_1^{*}$ is locally asymptotically stable when ${\tau }_1\in [0,{\tau }_{{\mathcal{H}}_{10}}^0)$ and ${\tau }_2\in [0,{\tau }_{{\mathcal{K}}_{10}}^0)$, and the rumor equilibrium point $E_1^{*}$ is unstable when ${\tau }_1>{\tau }_{{\mathcal{H}}_{10}}^0$ or ${\tau }_2>{\tau }_{{\mathcal{K}}_{10}}^0$. Therefore, system (18) has a Hopf bifurcation at $E_1^{*}$ when ${\tau }_1={\tau }_{{\mathcal{H}}_{10}}^0$ or ${\tau }_2={\tau }_{{\mathcal{K}}_{10}}^0$.

For the case $({\mathcal{H}}_3)$, when $L_2\ge 0$, $L_1\ge 0$, and $L_1^2-4L_2\ge 0$, Equation (30) has two positive roots, defined by $x_1={\omega }_{11}^2$, $x_2={\omega }_{12}^2$, and $x_1>x_2$. Choose ${\omega }_{11}>0$, ${\omega }_{12}>0$ then ${\omega }_{11}>{\omega }_{12}$, $\mathcal{G}({\omega }_{11})>0$, $\mathcal{G}({\omega }_{12})<0$. According to Equation (32), we obtain that

$$\begin{array}{c}\hfill {\tau }_{{\mathcal{H}}_{3}m}^{\nu }=\frac{1}{{\omega }_{1m}}arccos\left(\frac{{\omega }_{1m}^2-C-({\Omega }_1+{\Omega }_3){\omega }_{1m}}{{\omega }_{1m}{U}_1+{\Omega }_1{U}_1}\right)+\frac{2\pi \nu }{{\omega }_{1m}},\nu =0,1,\cdots ,m=1,2.\end{array}$$

(36)

Let ${\tau }_1^{*}={\tau }_{{\mathcal{H}}_{3m}}^0=min\{{\tau }_{{\mathcal{H}}_{3m}}^{\nu }:m=1,2,\nu =0,1,2,\cdots \}$ and assume that ${\lambda }_{11}=v_{11}+i{\omega }_{11}$, ${\lambda }_{12}=v_{12}+i{\omega }_{12}$, then

$$\begin{array}{c}\hfill Re{(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_1})}^{-1}{|}_{\lambda =i{\omega }_{11},{\tau }_1={\tau }_{{\mathcal{H}}_{31}}^{\nu }}>0,\nu =0,1,2,\cdots .\end{array}$$

$$\begin{array}{c}\hfill Re{(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_1})}^{-1}{|}_{\lambda =i{\omega }_{12},{\tau }_1={\tau }_{{\mathcal{H}}_{32}}^{\nu }}<0,\nu =0,1,2,\cdots .\end{array}$$

Similarly, for the case of $({\mathcal{K}}_3)$, defined $y_1={\omega }_{21}^2$, $y_2={\omega }_{22}^2$, and $y_1>y_2$, then ${\omega }_{21}>{\omega }_{22}$, $\mathcal{W}({\omega }_{21})>0$, $\mathcal{W}({\omega }_{22})<0$. According to (36), it has

$$\begin{array}{c}\hfill {\tau }_{{\mathcal{K}}_{3}m}^p=\frac{1}{{\omega }_{2m}}arccos\left(\frac{{\omega }_{2m}^2-D-({\Omega }_2+{\Omega }_4){\omega }_{2m}}{{\omega }_{2m}{U}_2+{\Omega }_2{U}_2}\right)+\frac{2\pi p}{{\omega }_{2m}},p=0,1,\cdots ,m=1,2\end{array}$$

(37)

and assume that ${\lambda }_{21}=v_{21}+i{\omega }_{21}$, ${\lambda }_{22}=v_{22}+i{\omega }_{22}$, then

$$\begin{array}{c}\hfill Re{(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_2})}^{-1}{|}_{\lambda =i{\omega }_{21},{\tau }_2={\tau }_{{\mathcal{K}}_{31}}^p}>0,p=0,1,2,\cdots .\end{array}$$

$$\begin{array}{c}\hfill Re{(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_2})}^{-1}{|}_{\lambda =i{\omega }_{22},{\tau }_2={\tau }_{{\mathcal{K}}_{32}}^p}<0,p=0,1,2,\cdots .\end{array}$$

Lemma 4. 

If $2{\omega }_1{B}_i+{\alpha }_i\langle k\rangle B_i{I}_i^{*}>{(B_i)}^2{I}_i^{*}+{\alpha }_i\langle k\rangle ({\mu }_i+d)$, and ${\omega }_{i1}>{\omega }_{i2}$ for $i=1,2$, then ${\tau }_{{\mathcal{H}}_{31}}^{\nu }<{\tau }_{{\mathcal{H}}_{32}}^{\nu }$ and ${\tau }_{{\mathcal{K}}_{31}}^{\nu }<{\tau }_{{\mathcal{K}}_{32}}^{\nu }$, $\nu =1,2,\cdots$.

Proof. 

Define

$$\begin{array}{c}\hfill g(\omega )=\frac{{\omega }_1^2-C-({\Omega }_1+{\Omega }_3){\omega }_1}{({\omega }_1{U}_1+{\Omega }_1{U}_1)},\overline{g}(\omega )=arccosg(\omega ),\overline{G}(\omega )=\frac{1}{\omega }\overline{g}(\omega ).\end{array}$$

Then, the derivative of $g(\omega )$ is

$$\begin{array}{cc}\hfill g^{{}^{\prime }}(\omega )& =\frac{({\omega }_1^2+2{\omega }_1{\Omega }_1-({\Omega }_1+{\Omega }_3){\Omega }_1+C)U_1}{{({\omega }_1{U}_1+{\Omega }_1{U}_1)}^2}\hfill \\ \hfill & =\frac{[{\omega }_1^2+(2{\omega }_1-B_1{I}_1^{*}+{\alpha }_1{I}_1^{*}\langle k\rangle )B_1{I}_1^{*}-{\alpha }_1{I}_1^{*}\langle k\rangle ({\mu }_1+d)]U_1}{{({\omega }_1{U}_1+{\Omega }_1{U}_1)}^2}\hfill \\ \hfill & \ge \frac{[2{\omega }_1{B}_1+{\alpha }_1{I}_1^{*}\langle k\rangle B_1-{(B_1)}^2{I}_1^{*}-{\alpha }_1\langle k\rangle ({\mu }_1+d)]I_1^{*}}{U_1{({\omega }_1+{\Omega }_1)}^2}\hfill \\ \hfill & >0.\hfill \end{array}$$

Since $g(\omega )$ is always an increasing function about $\omega$, it is obtained that $\overline{g}(\omega )=arccosg(\omega )$ is a decreasing function about $\omega$. Therefore, we can obtain that

$$\begin{array}{c}\hfill {\overline{G}}^{{}^{\prime }}(\omega )=-\frac{1}{{\omega }^2}\overline{g}(\omega )+\frac{1}{\omega }{\overline{g}}^{{}^{\prime }}(\omega )<0.\end{array}$$

By conclusion, ${\tau }_{{\mathcal{H}}_{31}}^{\nu }<{\tau }_{{\mathcal{H}}_{32}}^{\nu }$ when ${\omega }_{11}>{\omega }_{12}$. Similarly, we obtain that ${\tau }_{{\mathcal{K}}_{31}}^{\nu }<{\tau }_{{\mathcal{K}}_{32}}^{\nu }$ for ${\omega }_{21}>{\omega }_{22}$, where $\nu =1,2,\cdots$. Hence, the following theorem is given. □

Theorem 9. 

For system (18), if Lemma 4 and $H_i(S_i^{*})>0$, $G_i(S_i^{*})>0$ for $i=1,2$, are satisfied, the following results can be obtained:

(1). If both $({\mathcal{H}}_3)$ and $({\mathcal{K}}_3)$ are satisfied, the rumor equilibrium point $E_1^{*}$ is locally asymptotically stable for ${\tau }_1\in [0,{\tau }_{{\mathcal{H}}_{31}}^0)$ and ${\tau }_2\in [0,{\tau }_{{\mathcal{K}}_{31}}^0)$, and unstable for ${\tau }_1>{\tau }_{{\mathcal{H}}_{31}}^0$ or ${\tau }_2>{\tau }_{{\mathcal{K}}_{31}}^0$. Hence, system (18) has a Hopf bifurcation at $E_1^{*}$ when ${\tau }_1={\tau }_{{\mathcal{H}}_{31}}^0$ or ${\tau }_2={\tau }_{{\mathcal{K}}_{31}}^0$.

(2). If $({\mathcal{H}}_1)$ and $({\mathcal{K}}_3)$ are satisfied, the rumor equilibrium point $E_1^{*}$ is locally asymptotically stable for ${\tau }_1\in [0,{\tau }_{{\mathcal{H}}_{10}}^0)$ and ${\tau }_2\in [0,{\tau }_{{\mathcal{K}}_{31}}^0)$, and unstable for ${\tau }_1>{\tau }_{{\mathcal{H}}_{10}}^0$ or ${\tau }_2>{\tau }_{{\mathcal{K}}_{31}}^0$. Therefore, system (18) has a Hopf bifurcation at $E_1^{*}$ when ${\tau }_1={\tau }_{{\mathcal{H}}_{10}}^0$ or ${\tau }_2={\tau }_{{\mathcal{K}}_{31}}^0$.

(3). If $({\mathcal{H}}_3)$ and $({\mathcal{K}}_1)$ are satisfied, the rumor equilibrium point $E_1^{*}$ is locally asymptotically stable when ${\tau }_1\in [0,{\tau }_{{\mathcal{H}}_{31}}^0)$ and ${\tau }_2\in [0,min\{{\tau }_{{\mathcal{K}}_{10}}^0)$, and the rumor equilibrium point $E_1^{*}$ is unstable when ${\tau }_1>{\tau }_{{\mathcal{H}}_{31}}^0$ or ${\tau }_2>{\tau }_{{\mathcal{K}}_{10}}^0$. Therefore, system (18) has a Hopf bifurcation at $E_1^{*}$ when ${\tau }_1={\tau }_{{\mathcal{H}}_{31}}^0$ or ${\tau }_2={\tau }_{{\mathcal{K}}_{10}}^0$.

(4). If one of $({\mathcal{H}}_i)(i=1,3)$ and one of $({\mathcal{K}}_j)(j=2,4)$ are satisfied, the rumor equilibrium point $E_1^{*}$ is locally asymptotically stable for ${\tau }_1\in [0,{\tau }_{{\mathfrak{H}}_{i0}}^0)$, and unstable for ${\tau }_1>{\tau }_{{\mathcal{H}}_{i0}}^0$. Hence, system (18) has a Hopf bifurcation at $E_1^{*}$ when ${\tau }_1={\tau }_{{\mathcal{H}}_{i0}}^0$.

(5). If one of $({\mathcal{H}}_i)(i=2,4)$ and one of $({\mathcal{K}}_j)(j=1,3)$ are satisfied, the rumor equilibrium point $E_1^{*}$ is locally asymptotically stable for ${\tau }_2\in [0,{\tau }_{{\mathcal{K}}_{j0}}^0)$, and unstable for ${\tau }_2>{\tau }_{{\mathcal{K}}_{j0}}^0$. Hence, system (18) has a Hopf bifurcation at $E_1^{*}$ when ${\tau }_2={\tau }_{{\mathcal{K}}_{j0}}^0$.

Remark 5. 

From the proof of the above theorems, we can see that the stability analysis of rumor equilibrium in the model with time delay is more complex than that in the model without time delays. This shows that the dynamic behavior of the rumor equilibrium is affected by the time delay.

Remark 6. 

In existing works [36,37,38,39], some rumor-spreading models were proposed on OSNs in a multilingual environment. In these models, the isolation mechanism, immune mechanism and time delay were not considered. However, in this paper, we comprehensively consider these factors and propose a new rumor-spreading model on OSNs. The dynamic behavior is carefully studied. Moreover, the optimal control method is applied to suppress the spreading of the rumor.

5. Numerical Simulations

In this section, we verify the effectiveness of the theoretical results by selecting different parameters. In

Table 1. The set of parametric values for simulation.

Parameters	Set 1	Set 2	Set 3	Set 4	Set 5	Set 6	Set 7
$B_1$	0.1	0.11	0.11	0.11	0.09	0.03	0.1
$B_2$	0.1	0.1	0.1	0.1	0.13	0.18	0.1
${\alpha }_1$	0.27	0.25	0.25	0.23	0.2	0.32	0.23
${\alpha }_2$	0.26	0.14	0.25	0.15	0.25	0.25	0.15
${\mu }_1$	0.01	0.01	0.01	0.01	0.01	0.01	0.01
${\mu }_2$	0.3	0.01	0.01	0.01	0.01	0.3	0.01
${\beta }_1$	0.13	0.13	0.13	0.13	0.24	0.08	0.13
${\beta }_2$	0.12	0.15	0.15	0.14	0.15	0.08	0.14
c	0.1	0.1	0.1	0.1	0.1	0.1	0.1
${\sigma }_1$	0.003	0.52	0.52	0.41	0.5	0.58	0.5
${\sigma }_2$	0.0028	0.51	0.51	0.4	0.52	0.5	0.5
$r_1$	0.001	0.02	0.02	0.05	0.1	0.1	0.06
$r_2$	0.001	0.05	0.03	0.05	0.08	0.0005	0.06
d	0.2	0.21	0.21	0.21	0.21	0.21	0.2
$\rho$	0.08	0.03	0.03	0.08	0.05	0.08	0.08
$\langle k\rangle$	3.5	7	7	7	7	3.5	7

, based on the range of parameters in the model, we randomly select seven groups of parameters to simulate the evolution of the rumor under different threshold conditions, and verify the correctness of the theoretical results.

5.1. The Stability of Rumor-Free Equilibrium

Example 1. 

For system (1) with the parameters in Set 1, we obtain ${\Re }_{01}=0.61<1$ and ${\Re }_{02}=0.25<1$. Then, it has a unique rumor-free equilibrium $E_0=(0.48,0.2,0,0,0,0.32)$. According to Theorem 2, $E_0$ is locally asymptotically stable. For different initial values, the dynamics of system (1) is shown in Figure 3. In Figure 3a, we can see that $S_1(t)$ and $S_2(t)$ will converge asymptotically to 0 under different initial values. In Figure 3b, $Q(t)$ gradually tends to 0. $I_1(t)$, $I_2(t)$ and $R(t)$ converge to 0.48, 0.2 and 0.32, respectively. It can be verified that $E_0$ is locally asymptotically stable. In other words, the rumor will disappear automatically. According to the expression of $E_0$, one can find that the stable states of $I_1(t)$ and $I_2(t)$ are $\frac{B_1}{d+{\mu }_1}$ and $\frac{B_2}{d+{\mu }_2}$, respectively. In Set 1, we choose $B_1=B_2$, then the main factor affecting $I_1(t)$ and $I_2(t)$ are the parameters ${\mu }_1$ and ${\mu }_2$. Because ${\mu }_1=0.01$ and ${\mu }_2=0.3$, the stable value of $I_1(t)$ is greater than the stable value of $I_2(t)$ in the numerical simulation. Moreover, as ${\mu }_2>{\mu }_1$, the descending speed and amplitude of $I_2(t)$ are greater than $I_1(t)$. In addition, the evolutionary processes of $I_1(t)$ and $I_2(t)$ are different under different initial values.

For the case with time delays, we can verify that the conditions of Theorem 7 hold. By calculation, ${\tau }_{10}=4.38$ and ${\tau }_{20}=5.12$. As shown in Figure 4, we find that the time-delay affects the convergence rate and the stability of the rumor-free equilibrium point. When ${\tau }_1\in [0,{\tau }_{10})$ and ${\tau }_2\in [0,{\tau }_{20})$, $E_0$ is locally asymptotically stable. When ${\tau }_1>{\tau }_{10}$ or ${\tau }_2>{\tau }_{20}$, $E_0$ is unstable. This is consistent with the results of Theorem 7.

In the analysis of Theorem 7, based on (22), we find that the ratios of $r_i$ to ${\sigma }_i$, $i=1,2$ have an important influence on the threshold condition of the time delay. From Figure 5, one can find that with the increase in the ratio, the time-delay threshold decreases, and they are inversely proportional.

However, the Hopf bifurcation exists at $E_0$ when the time delay is greater than the threshold in theoretical. In Figure 4, it is easy to find that there exist $S_1(t)<0$ and $S_2(t)<0$. This is inconsistent with the reality. In real life, $S_i(t)\ge 0$ for $i=1,2$. Therefore, the Hopf bifurcation will not exist. We redrew the dynamic behavior of $S_i(t)$ in combination with $S_i(t)\ge 0$ for $i=1,2$. It can be seen from Figure 6 that the time delay will inhibit the spread of the rumor.

5.2. The Stability of Rumor Equilibrium

Example 2. 

Consider system (1) with parameters in Set 2 of Table 1. It can be obtained that ${\Re }_{01}=2.14>1$ and ${\Re }_{02}=0.97<1$, ${\alpha }_2{\sigma }_2\langle k\rangle =0.26>{\mu }_2+d=0.22$, then system (1) has a unique rumor equilibrium $E_1^{*}$. As shown in Figure 7, rumor equilibrium $E_1^{*}$ is locally asymptotically stable. By calculating the parameters in Set 3 of Table 1, the results of ${\Re }_{01}=2.14>1$ and ${\Re }_{02}=1.9>1$, ${\alpha }_2{\sigma }_2\langle k\rangle =0.26>{\mu }_2+d=0.22$ are gained. Figure 8 shows that the system (1) achieves local asymptotic stability. Similarly, Figure 9 shows the locally asymptotically stability of the equilibrium point $E_1^{*}$ drawn with the parameters in Set 4 of Table 1. By comparing with the results, we find that when ${\Re }_{01}>1$, the value of $S_2$ at the equilibrium point increases with the increase in ${\Re }_{02}$. Figure 10 describes the asymptotical stability of the equilibrium $E_1^{*}$. These simulation results show that when ${\Re }_{01}>1$, the rumor in the network will continue to prevail in both languages.

Furthermore, Figure 11 shows the states of $S_1(t)$ and $S_2(t)$ with parameters in Set 5 and Set 6 of Table 1, respectively. One can find that when ${\Re }_{01}<1$, the rumor will only continue to spread in the official language.

5.3. The Hopf Bifurcation of Rumor Equilibrium

Example 3. 

In order to verify the influence of time delays on rumor spreading, we select the parameters in Set 7 of Table 1 in model (15). By calculating, we can obtian that $({\mathcal{H}}_1)$ and $({\mathcal{K}}_1)$ are satisfied. In particular, if ${\tau }_1={\tau }_2=\tau$, it has $min\{{\tau }_{10}^0,{\tau }_{20}^0\}=11.2549$. As shown in Figure 12, $S_1(t)$ and $S_2(t)$ are locally asymptotically stable if $\tau <11.2549$, and Figure 13 shows that $S_1(t)$ and $S_2(t)$ are unstable when $\tau =11.264$. Figure 14 describes the phase diagram of the equilibrium $E_1^{*}$. Therefore, it can be seen that the equilibrium $E_1^{*}$ becomes unstable as the time delay increasing. When the time delay is greater than the threshold, the stability of the equilibrium will disappear. This means that the rumor disseminator will exist continuously, and the density of the disseminator will fluctuate within a certain range over time.

In model (15), when ${\tau }_1\ne {\tau }_2>0$, the parameters are selected as the ones in Set 7 of Table 1. Through verification, $({\mathcal{H}}_1)$ and $({\mathcal{K}}_1)$ are satisfied in Theorem 8, and ${\tau }_{{\mathcal{H}}_{10}}^0=11.5$, ${\tau }_{{\mathcal{K}}_{10}}^0=9.5$. Figure 15, Figure 16 and Figure 17 show the densities of $S_1(t)$ and $S_2(t)$ with different time delays. We observe that the system is stable at the rumor equilibrium if the time delays are less than their threshold values in Figure 15. However, the system undergoes Hopf bifurcation when the time delays are equal to their threshold values. This situation can be seen in Figure 16. In addition, the system is unstable at the rumor equilibrium if the time delays are great than their threshold values, shown in Figure 17.

5.4. Feasibility of Optimal Control

Example 4. 

For optimal control system (15), the parameters are shown in Set 3 of Table 1. In the objective function (16), the parameters are $u_1=3$, $u_2=2$, $u_3=5$ and $u_4=1$, respectively. In order to explore the influence of optimal control on the system (15), the trajectories of the system (15) with and without optimal control are simulated, as shown in Figure 18. Through observation, we find that the values of $S_i(t)$ $(i=1,2)$ and $R(t)$ will decrease under the optimal control, and the densities of $I_i(t)$ $(i=1,2)$ and $Q(t)$ will be more under optimal control.

In (16), we choose $T=10$, then the control strength $r_1(t)$ and $r_2(t)$ are shown in Figure 19a. It can be observed that the control $r_1(t)$ and $r_2(t)$ decreases with time and gradually tends to 0 at $t=10$. Figure 19b shows the performance of the objective function. It can be seen that $J(t)$ is monotone and nondecreasing. Obviously, in a certain period of time, the control intensity of the system gradually decreases to 0, while the control consumption gradually increases to the maximum.

In order to facilitate the control, the control time in the optimal control will be given in advance. In the specific simulation, we give the trajectories of $r_1(t)$ and $r_2(t)$ when $T=2$, $T=5$, $T=8$ and $T=10$, as shown in Figure 20a. Figure 20b shows the consumption at $T=2$, $T=5$, $T=8$ and $T=10$, respectively. Through observation, we find that the larger the given time T, the greater the consumption of $J(t)$. Figure 20c shows the trajectories of $J(t)$ under different control forces. We found that among all these controls, $r_1=1$ and $r_2=1$ have the largest consumption, while the consumption of optimal control is minimal. That means that the optimal control proposed in this paper not only can control the rumor within the specified time, but also consumes the least resources.

6. Conclusions

In this paper, we studied the modeling, dynamic analysis and control of rumor propagation with quarantine control on multilingual OSNs. Firstly, we proposed a 2I2SQR rumor-spreading model without considering propagation delay. We gave the rumor-free equilibrium and calculated the basic reproduction number by using the next generation matrix method. Based on backward bifurcation theory, we discussed the existence of a rumor equilibrium, which is somewhat more complicated. It was found that the number of rumor equilibria is different under different parameters, including no rumor equilibrium, only one rumor equilibrium and the case with two rumor equilibria. By using the Jacobian matrix method and a differential inequality lemma, the local and global asymptotic stabilities of rumor free-equilibrium were proved. The local stability and bifurcation of the rumor equilibrium were discussed, and some related theorems were given. Secondly, in order to suppress the spreading of rumor with the lowest cost in an expected period, a continuous optimal control measure was proposed to curb the propagation of rumor. According to Pontryagin’s maximum principle, an optimal isolation strength was obtained. Moreover, considering the time delays from rumor spreading to isolation control, we proposed a 2I2SQR rumor propagation model with time delays on multilingual OSNs. The local stability of rumor-free and rumor equilibria was investigated. It was found that, different from the model without delay, both the rumor-free equilibrium and the rumor equilibrium are unstable and have Hopf bifurcations when the propagation delay is greater than some certain threshold conditions. Finally, the theory results were verified by some numerical simulations. In this research, the uniform network structure was considered in the model, and a continuous control strategy was proposed. In our future work, we will explore a more general heterogeneous network structure and propose some new control strategies to suppress the spread of rumors.