Effectiveness of matrix measure in finding periodic solutions for nonlinear systems of differential and integro-differential equations with delays

1.   Introduction

Rumor spreading as a social contagion process is very similar to the epidemic diffusion, so most rumor spreading models have evolved from epidemic models, such as SI, SIR and SIS. A classic rumor model was DK model proposed by Daley and Kendall in 1964 ^[2], in which the population was divided into three classes: people who did not know the rumor, people who spread the rumor and people who know but will never spread the rumor. Maki and Thomson modified the DK model into the MK model ^[8]. Since then, a number of scholars have proposed various rumor spreading models by improving the traditional epidemic models ^{[1,3,14,15,16]}. In 2019, Tian and Ding ^[11] formulated an ordinary differential equation (ODE) compartmental model for rumor, where the population was divided into five disjoint classes, namely the ignorants, the latents, the rumor-spreaders, the debunkers and the stiflers. At time $t$, the numbers in each of these classes are denoted by $I(t), L(t), R(t), D(t)$ and $S(t)$, respectively. The rumor spreading model considered by Tian and Ding ^[11] can be described by a system of ODEs

where $\Lambda$ is the constant immigration rate of the population, $\mu$ is the rumor-contacting rate, $k$ is the debunker-contacting rate, $\rho$ is the rate at which all existing users exit from the five classes (i.e. emigration rate), $\alpha$ is the rumor-spreading rate, $\beta$ is the debunking rate, $\gamma$ is the silent rate, $\delta$ is the rumor-stifling, $\xi$ is the reversal rate and $\theta$ represents the debunking-stifling rate. All parameters are assumed to be independent of time $t$ and positive.

The model used in the above study to describe rumor propagation behavior is deterministic model, whereas the random models used to study rumor propagation are few^[3]. But in the real world, rumor models are often affected by environmental noise. Especially in emergency events, when rumors are widely spread, the propagation process is affected by many uncertain factors, which increase the volatility of the propagation process. Therefore, it would be necessary and interesting to reveal how the environmental noise affects the rumor spreading model. Following the idea of Jia et al.^[3], in this paper, we assume that the stochastic perturbations are of white noise type which are proportional to the system variable, respectively. Then we obtain a stochastic analogue of the deterministic model (1.1) as follows

where $B_i(t)(i = 1, 2, 3, 4, 5)$ are independent Brownian motions and $\sigma_i > 0 (i = 1, 2, 3, 4, 5)$ are their intensities. All the other parameters in system (1.2) have the same meaning as in system (1.1). Throughout this paper, unless otherwise specified, let $(\Omega, \mathscr{F}, \{\mathscr{F}_t\}_{t\geq0}, \mathbb{P})$ be a complete probability space with a filtration $\{\mathscr{F}_t\}_{t\geq0}$ satisfying the usual conditions, that is, it is rightly continuous and increasing while $\mathscr{F}_0$ contains all $\mathbb{P}-$null sets, and let $B_i(t)(i = 1, 2, 3, 4, 5)$ be scalar Brownian motions defined on the probability space. For the sake of simplicity, we introduce the following notations: $\mathbb{R}_+^5 = \{(x_1, x_2, x_3, x_4, x_5)\in\mathbb{R}^5:x_i > 0, i = 1, 2, 3, 4, 5\}, a\vee b = \max\{a, b\}, a\wedge b = \min\{a, b\}, \sigma = \sigma_1\vee\sigma_2\vee\sigma_3\vee\sigma_4\vee\sigma_5.$

With the help of the Lyapunov function methods and the inequality techniques, the existence and uniqueness of an ergodic stationary distribution of the positive solutions to system (1.2) are presented. The main difficulties lies in the construction of Lyapunov function and the construction of a bounded closed domain. The main contribution of this paper are highlighted as follows: (ⅰ) a stochastic rumor spreading model is proposed and investigated; (ⅱ) some sufficient conditions for the existence of an ergodic stationary distribution; and (ⅲ) the stationary distribution implies that the rumor can be persistent in the mean. The subsequent part of this paper is as follows: In Section 2, we prove the existence and uniqueness of a global positive solution to system (1.2) with any positive initial value. In Section 3, by constructing a suitable stochastic Lyapunov function, we establish sufficient conditions for the existence and uniqueness of an ergodic stationary distribution of the positive solutions to model (1.2). In Section 4, some numerical simulations are provided to illustrate our theoretical results. Finally, some concluding remarks are presented to end this paper.

2.   Existence and uniqueness of the global positive solution

To analyze the dynamical behavior of a rumor spreading model, the first concerning thing is whether the solution is global and positive. In this section, motivated by the methods in ^[9] and we show that there is a unique global positive solution of system (1.2). The key is to construct a Lyapunov function.

Theorem 2.1. For any given initial value $(I(0), L(0), R(0), D(0), S(0))\in \mathbb{R}_{+}^5$, system (1.2) admits a unique solution $(I(t), L(t), R(t), D(t), S(t))\in \mathbb{R}_{+}^5$ on $t\geq 0$ and the solution will remain in $\mathbb{R}_{+}^5$ with probability one, namely $(I(t), L(t), R(t), D(t), S(t))\in \mathbb{R}_{+}^5$ for all $t\geq 0$ almost surely (a.s.).

Proof. Since the coefficients of system (1.2) satisfy the local Lipschitz condition, we know that, for any initial value $(I(0), L(0), R(0), D(0), S(0))\in \mathbb{R}_{+}^5$, there is a unique local solution $(I(t), L(t), R(t), D(t), S(t))\in \mathbb{R}_{+}^5$ on $t\in [0, \tau _{e}]$, where $\tau _{e}$ is the explosion time ^[10]. Now we prove the solution is global, i.e. to prove $\tau _{e} = \infty$ a.s. To this end, let $m_0 > 0$ be sufficiently large such that each component of $(I(0), L(0), R(0), D(0), S(0))$ all lies in the interval $[\frac{1}{m_0}, m_0 ]$. For each integer $m\geq m_0$, define the following stopping time

Throughout this paper, we set $\inf \emptyset = \infty$ (as usual $\emptyset$ denotes the empty set). Obviously, $\tau _{m}$ is an increasing function as $m \rightarrow \infty$. Set $\tau _{\infty} = \lim_{m \rightarrow \infty} \tau _{m}$. Then $\tau _{\infty} \leq \tau _{e}$ a.s. If $\tau _{\infty} = \infty$ a.s. is true, then $\tau _{e} = \infty$ a.s. and $(I(t), L(t), R(t), D(t), S(t))\in \mathbb{R}_{+}^5$ a.s. for all $t\geq 0$. In other words, in order to show this assertion, we only need to prove $\tau _{\infty} = \infty$ a.s. If the assertion is false, then there is a pair of constants $T > 0$ and $\bar{\epsilon} \in (0, 1)$ such that $\mathbb{P}\{\tau _{m}\leq T\}\geq \bar{\epsilon}$ for each integer $m \geq m_0$. Define a $C^2$-function $V:\mathbb{R}_{+}^5\rightarrow\mathbb{R}_{+}\cup\{0\}$ by

where $a$ is a positive constant to be determined later. The nonnegativity of this function can be seen from $u-1-\ln u\geq 0, \forall u > 0$. According to the general Itô formula (see, for example, Theorem 4.2.1 of ^[10]), we have

where $\mathcal{L}V:\mathbb{R}_{+}^5 \rightarrow \mathbb{R}$ is defined by

Choose $a = min\{\frac{\rho}{\mu}, \frac{\rho}{k}\}$, then we obtain

and $K$ is a positive constant. Thus,

For any $m\geq m_0$, integrating (2.1) on both sides from 0 to $\tau_m \wedge T$ and then taking expectation yield

Let $\Omega_m = \{\omega\in\Omega:\tau_m = \tau_m(\omega)\leq T\}$ for $m\geq m_0$. Then we have $\mathbb{P}(\Omega_m)\geq \bar{\epsilon}$. Note that, for every $\omega \in \Omega_m$, there exists $I(\tau_m, \omega)$ or $L(\tau_m, \omega)$ or $R(\tau_m, \omega)$ or $D(\tau_m, \omega)$ or $S(\tau_m, \omega)$ equaling either $m$ or $\frac{1}{m}$. Thus $V(I(\tau_m, \omega), L(\tau_m, \omega)), R(\tau_m, \omega), D(\tau_m, \omega), S(\tau_m, \omega)$ is no less than either

$m-a-a\ln\frac{m}{a}$ or $\frac{1}{m}-a-a\ln\frac{1}{ma} = \frac{1}{m}-a+a\ln(ma)$ or $m-1-\ln m$ or $\frac{1}{m}-1-\ln\frac{1}{m} = \frac{1}{m}-1+\ln m.$

So we have

Consequently,

where $1_{\Omega_m}$ denotes the indicator function of $\Omega_m$. Letting $m\rightarrow\infty$ leads to the contradiction

This completes the proof.

3.   Existence of ergodic stationary distribution

When considering rumor propagation model, we are also interested to know when the rumor will persist and prevail in a population. In the deterministic models, it can be solved by proving the rumor-epidemic equilibrium of the corresponding model is a global attractor or globally asymptotically stable. But there is no rumor-epidemic equilibrium in system (1.2). In this section, based on the theory of Khasminskii(^[7]), we prove that there is a stationary distribution which reveals that the rumor will persist in the mean. Here we present some theory about the stationary distribution which is introduced in (^[7]).

Definition 3.1. (^[7]) The transition probability function $P(s, x, t, A)$ is said to be time-homogeneous (and the corresponding Markov process is called time-homogeneous) if the function $P(s, x, t+s, A)$ is independent of $s$, where $0\leq s\leq t$, $x\in\mathbb{R}^l$ and $A\in \mathcal{B}$ and $\mathcal{B}$ denotes the $\sigma-$ algebra of Borel sets in $\mathbb{R}^l$.

Definition 3.2. (^[7]) Let $X(t)$ be a regular time-homogeneous Markov process in $\mathbb{R}^l$ described by the stochastic differential equation:

The diffusion matrix of the process $X(t)$ is defined as follows:

Lemma 3.1. (^[7]) The Markov process $X(t)$ has a unique ergodic stationary distribution $\pi(\cdot)$ if there exists a bounded domain $D \subset \mathbb{R}^n$ with regular boundary $\Gamma$ and

A1: there is a positive number $M$ such that

A2: there exists a nonnegative $C^2$-function $V$ such that $\mathcal{L}V$ is negative for any $\mathbb{R}^n \setminus D$, where $\mathcal{L}$ denotes the differential operator defined by

Then

for all $x\in \mathbb{R}^n$, where $f(\cdot)$ is a function integrable with respect to the measure $\pi$.

Theorem 3.1. Assume

Then, for any initial value $(I(0), L(0), R(0), D(0), S(0)) \in \mathbb{R}_{+}^5$, system (1.2) has a unique stationary distribution $\pi(\cdot)$ and the ergodicity holds.

Proof. Theorem 2.1 tells us that for any initial value $(I(0), L(0), R(0), D(0), S(0)) \in \mathbb{R}_{+}^5$, there is a unique global solution $(I(t), L(t), R(t), D(t), S(t)) \in \mathbb{R}_{+}^5.$ In order to prove this Theorem, it suffices to validate $A1$ and $A2$ in Lemma 3.1. First, we verify $A1$. The diffusion matrix of system (1.2) is given by

It is easy to see that the matrix $A$ is positive definite for any compact subset of $\mathbb{R}_{+}^5$, so condition $A1$ in Lemma 3.1 is satisfied. Now we prove condition $A2$. Define a $C^2$-function

where $n_1, n_2, n_3, n_4, n_5, \varrho$ and $M$ are positive constants, which will be determined later. It is easy to check that

where $U_k = \Pi_{i = 1}^5\left(\frac{1}{k}, k\right)$. In addition, $\tilde{Q}(I, L, R, D, S)$ is a continuous function. Hence, $\tilde{Q}(I, L, R, D, S)$ must have a minimum point $(I_0, L_0, R_0, D_0, S_0)\in \mathbb{R}_{+}^5$. Therefore, we define a nonnegative $C^2$-function $Q:\mathbb{R}_{+}^5\rightarrow\mathbb{R}_{+}$

Applying the general Itô formula ^[10] to $Q_1$, one obtains the differential operator $\mathcal{L}$ of $Q_1$ as follows:

Let

Then it follows that

Similarly, one has

where we choose $\varrho$ sufficiently small such that $\phi = \rho-\frac{\varrho\sigma^2}{2} > 0$ and

Moreover, one has

Making use of (2.1)–(3.3), we then derive that

For the convenience, we define

and

Now we are in the position to construct a bounded closed domain $U_\epsilon$ such that the condition $A2$ in Lemma 3.1 holds. To this end, we define a compact set as follows

where $\epsilon > 0$ is a sufficiently small constant such that

and

For convenience, we can divide $\mathbb{R}_{+}^5\setminus U_\epsilon$ into ten domains, where

Obviously, $\mathbb{R}_{+}^5\setminus U_\epsilon = \bigcup_{i = 1}^{10}U_i.$ Next, we will prove that $\mathcal{L}{Q(I, L, R, D, S)}\leq-1$ for any $(I, L, R, D, S)\in \mathbb{R}_{+}^5\setminus U_\epsilon$, which is equivalent to proving it on the above ten domains, respectively.

Case 1. For any $(I, L, R, D, S)\in U_1$, then (3.4) implies that

Case 2. For any $(I, L, R, D, S)\in U_2$, using (3.4) one obtains

Case 3. For any $(I, L, R, D, S)\in U_3$, in view of (3.4), we get

Case 4. For any $(I, L, R, D, S)\in U_4$, it follows from (3.4) that

Case 5. For any $(I, L, R, D, S)\in U_5$, according to (3.4), we derive

Case 6. For any $(I, L, R, D, S)\in U_6$, by (3.4), we obtain

Case 7. For any $(I, L, R, D, S)\in U_7$, note from (3.4) that

Case 8. For any $(I, L, R, D, S)\in U_8$, making use of (3.4) one obtains that

Case 9. For any $(I, L, R, D, S)\in U_9$, we know from (3.4) that

Case 10. For any $(I, L, R, D, S)\in U_{10}$, using (3.4), we can show that

It follows from (3.5)–(3.16) that

which proves that condition $A2$ holds. Thus the conditions in Lemma 3.1 are verified, the proof is completed.

Remark 3.1. Comparing with the threshold parameter $R_0$ in ^[11], the parameter $R_0^s$ in stochastic system (1.2) is less than $R_0$, which reveals that the extinction of the rumor is much easier than in the corresponding deterministic model (1.1). Moreover, taking attention to the expression of $R_0^s$, we can control the rumor propagation by environmental white noise.

4.   Numerical simulations

In this section, we give some examples to illustrate the obtained theoretical results. We illustrate our findings by the Milstein's Higher Order Method developed in ^[5]. According to this method, we can get the following discretization equation of system (1.2):

where time increment $\Delta t > 0$, and $\xi_{1, k}, \xi_{2, k}, \xi_{3, k}, \xi_{4, k}, \xi_{5, k}$ are independent Gaussian random variables which follows $N(0, 1)$. However, we would like to point out that the values of parameters of system (1.2) and the initial values in the following numerical simulations are chosen for illustration purposes and are not taken from any real life data for any rumors. To this end, we set$(I(0), L(0), R(0), D(0), S(0)) = (0.3, 0.6, 0.2, 0.9, 0.6)$.

Example 1. Choose the parameters: $\Lambda = 5, \mu = 0.045, k = 0.2, \alpha = 0.05, \beta = 0.05, \gamma = 0.05, \rho = 0.01, \delta = 0.01, \xi = 0.01, \theta = 0.2.$

By simple calculation, we have $R_0^s = 1.213 > 1,$ which means that the conditions of Theorem 3.1 hold. Therefore, system (1.2) has a unique ergodic stationary distribution. Figures 1 and 2 illustrate this fact.

Example 2. Choose the parameters: $\Lambda = 0.6, \mu = 0.04, k = 0.8, \alpha = 0.5, \beta = 0.05, \gamma = 0.5, \rho = 0.01, \delta = 0.1, \xi = 0.2, \theta = 0.2.$

In what follows, we start comparing the stochastic system and deterministic system. In Figure 3, the rumor-spreaders in stochastic system are shown in red lines, compared with the rumor-spreaders in deterministic system which are shown in blue lines. It reveals that the environmental disturbance may help to curb the outbreak of rumors. Further, if we increase the environmental noise, the simulation results in Figure 4 suggests that the extinction of rumor-spreaders in stochastic system is much more easier than that in the corresponding deterministic system.

5.   Conclusions

In the current paper, we have studied a stochastic rumor spreading model. We have established sufficient conditions for the existence and uniqueness of an ergodic stationary distribution of the positive solutions to system (1.2) by using the stochastic Lyapunov function method. The ergodic property can help us better understand cycling phenomena of a rumor spreading model, and so describe the persistence of a rumor spreading model in practice. More precisely, we have obtained the following result:

● Assume

Then, for any initial value $(I(0), L(0), R(0), D(0), S(0)) \in \mathbb{R}_{+}^5$, system (1.2) has a unique stationary distribution $\pi(\cdot)$ and the ergodicity holds. The stationary distribution indicates that the rumor can become persistent in vivo. The theoretic work extended the results of the corresponding deterministic system. The results show that the rumors will maintain its persistence if the environmental noise is sufficiently small, while large stochastic noise can suppress the spread of rumors.

Some interesting topics deserve further consideration. As we all know, time-delay occurs frequently in many practical engineering systems, which is usually the source of oscillation, instability and poor performance of the systems. Now time-delay has been considered into many stochastic models (see, for example, ^[4,6,12,13]). We leave time-delay case for our future work.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (11701460), the Natural Science Foundation of Sichuan Provincial Education Department (Grant No. 18ZB0581), the Meritocracy Research Funds of China West Normal University (17YC373) and the research startup foundation of China West Normal University (Grant No. 18Q060).

Conflict of interest

The authors declare that there are no conflicts of interest.