Stochastic Bifurcation of an SIS Epidemic Model with Treatment and Immigration

Weipeng Zhang^1,2*, Dan Gu²

¹Department of Mathematics, Jinan University, Guangzhou, China.
²School of Mathematics and Statistics, Northeast Normal University, Changchun, China.

DOI: 10.4236/jamp.2024.126135                

Abstract

In this paper, we investigate an SIS model with treatment and immigration. Firstly, the two-dimensional model is simplified by using the stochastic averaging method. Then, we derive the local stability of the stochastic system by computing the Lyapunov exponent of the linearized system. Further, the global stability of the stochastic model is analyzed based on the singular boundary theory. Moreover, we prove that the model undergoes a Hopf bifurcation and a pitchfork bifurcation. Finally, several numerical examples are provided to illustrate the theoretical results.

Keywords

Epidemic Model, Stochastic Averaging Method, Singular Boundary Theory, Stochastic Bifurcation

1. Introduction

The study of infectious disease models plays a significant role in understanding and managing public health, socioeconomic factors, and maintaining social order. Researchers can perform mathematical analysis of epidemics based on classical models of infectious diseases by incorporating new factors and variables, providing a powerful tool for disease monitoring, prevention and control. The SIS epidemic model is one of the fundamental and critical dynamical models used to describe the spread of certain diseases, such as flu, the common cold, some venereal diseases, etc. It divides the total population into susceptible individuals S, and infected individuals I, for modeling. It suggests that individuals who have been infected and cured may still remain susceptible to the disease without gaining lifelong immunity. In other words, cured individuals can become reinfected with the same disease.

The transmission of infectious diseases is closely related to population immigration. The movement of population between different regions within a country has a significant impact on the spread of infectious diseases. Therefore, it is of crucial importance to incorporate the factor of immigration into epidemic models and investigate their influence on disease transmission. For instance, the authors have taken into account the factor of immigration in their study on epidemic models [1] [2] [3].

As medical technology advances, many diseases now have corresponding treatments. Flu, for example, can be treated with medication. Thus, in addition to immigration factors, treatment also affects the spread of the disease. To better reflect real-world situations, many researchers have incorporated treatment into epidemiological models [4] [5] [6]. For example, Buonomo et al. [4] considered an epidemiological model with treatment and vaccination, assuming that the cure rate of an infectious disease is proportional to the number of infected individuals. However, healthcare resources are actually limited, and it is important for each city to have appropriate treatment capacity to avoid unnecessary costs and the risk of disease outbreaks. If the treatment capacity is too large, the city will afford unnecessary expenses. Conversely, if the capacity is too small, it may not be able to effectively control disease outbreaks. Therefore, it is crucial to determine the treatment capacity for a specific disease in order to achieve rational allocation and efficient utilization of resources. Further, in [5] [6], researchers use

$T\left(I\right)=\{\begin{array}{ll}\rho I,\hfill & 0\le I\le I_k,\hfill \\ K=\rho I_k,\hfill & I>I_k,\hfill \end{array}$

as the treatment rate for infectious diseases. Specifically, when the population of infected individuals is high, the treatment rate remains constant. In turn, when the number of infected individuals is small, the treatment rate is proportional to the number of infected individuals.

Saikh et al. [7] considered the following SIS model with immigration and treatment

$\{\begin{array}{l}\text{d}S\left(t\right)=\left(R-\lambda SI-dS+\left(1-r\right)\Lambda +T\left(I\right)\right)\text{d}t\mathrm{,}\\ \text{d}I\left(t\right)=\left(\lambda SI-\left(d+\mu \right)I+r\Lambda -T\left(I\right)\right)\text{d}t\mathrm{,}\end{array}$

where $S\left(t\right)$ and $I\left(t\right)$ denote the susceptible and infectious individuals at time t, respectively. $T\left(I\right)$ represents a segmented treatment function, Λ is a constant flow of immigrants in unit time and r denotes the proportion of infected individuals among the immigrants.

In epidemiological models, the bilinear incidence rate $\beta SI$ is commonly used to simulate disease transmission. However, epidemiological models with nonlinear incidence rates exhibit more complex dynamical behavior. Liu et al. [8] [9] investigated the dynamics of SEIRS and SIRS models with a general form of nonlinear incidence rate $\beta I^p{S}^q$ . Based on this, Liu [10] studied a model with nonlinear incidence rate $\beta S{I}^2$ and analyzed its dynamical behavior. Motivated by the above literature and [11], we consider an SIS model with immigration and treatment and assume that the population recruitment rate is related to the total population size. Therefore, we obtain that when $0\le I\le I_k$ ,

$\{\begin{array}{l}\text{d}S\left(t\right)=\left(b\left(S+I\right)-\lambda S{I}^2-dS+\left(1-r\right)\Lambda +\rho I\right)\text{d}t\mathrm{,}\\ \text{d}I\left(t\right)=\left(\lambda S{I}^2-\left(d+\mu \right)I+r\Lambda -\rho I\right)\text{d}t\mathrm{,}\end{array}$ (1)

and when $I>I_k$ ,

$\{\begin{array}{l}\text{d}S\left(t\right)=\left(b\left(S+I\right)-\lambda S{I}^2-dS+\left(1-r\right)\Lambda +K\right)\text{d}t\mathrm{,}\\ \text{d}I\left(t\right)=\left(\lambda S{I}^2-\left(d+\mu \right)I+r\Lambda -K\right)\text{d}t\mathrm{.}\end{array}$ (2)

Based on the biological significance, all the parameters of model (1) and (2) are positive constants. Furthermore, the significance of parameters in model (1) and (2) is shown in Table 1.

Table 1. The significance of parameters in model (1) and (2).

Parameter	Significance
b	Birth rate
λ	Force of infection or disease incidence
d	Natural death rate
μ	Disease related death rate
ρ	Treatment rate
Λ	Constant flow of immigrants in unit time

By computation, we can obtain that if $r=0$ and $d>b$ , model (1) has a disease-free equilibrium $E_0\left(\Lambda /\left(d-b\right)\mathrm{,0}\right)$ . If $B_1=2\sqrt{A_1{C}_1}$ , ${\left(B_1/A_1-B/A\right)}^2>\left(b-d\right)/\lambda >0$ and $\mu >b-d$ , system (1) has a positive equilibrium point $E_1\left(S_1\mathrm{,}{I}_1\right)$ , where

$\begin{array}{l}A=\lambda \left(b-d-\mu \right)\mathrm{,}B=\lambda \Lambda \mathrm{,}\\ A_1={\lambda }^2{\Lambda }^2-3\lambda \left(b-d-\mu \right)\left(b-d\right)\left(d+\mu +\rho \right)\mathrm{,}\\ C_1={\left(b-d\right)}^2{\left(d+\mu +\rho \right)}^2+3\lambda {\Lambda }^{2}r\left(b-d\right)\mathrm{,}\\ S_1\mathrm{<mo>=</mo>}\frac{\left(b+\rho \right)I_1+\left(1-r\right)\Lambda }{\lambda I_1^2-b+d}\mathrm{,}{I}_1=\frac{B_1}{A_1}-\frac{B}{A}\mathrm{.}\end{array}$

Similarly, if $r=K/\Lambda$ and $d>b$ , model (2) also has a disease-free equilibrium point $E_0$ . If $B_2=2\sqrt{A_2{C}_2}$ , ${\left(B_2/A_2-B/A\right)}^2>\left(b-d\right)/\lambda >0$ , $\mu >b-d$ and $r>K/\Lambda$ , then model (2) has a positive equilibrium point $E_2\left(S_2\mathrm{,}{I}_2\right)$ , where

$\begin{array}{l}{A}_2={\lambda }^2{\Lambda }^2-3\lambda \left(b-d-\mu \right)\left(b-d\right)\left(d+\mu \right),\\ C_2={\left(b-d\right)}^2{\left(d+\mu \right)}^2-3\lambda \Lambda \left(b-d\right)\left(K-r\Lambda \right),\\ S_2=\frac{b{I}_2+\left(1-r\right)\Lambda +K}{\lambda I_2^2-b+d},I_2=\frac{B_2}{A_2}-\frac{B}{A}.\end{array}$

Actually, many small independent random fluctuations, such as temperature, humidity, etc., often affect population sizes through their small variations [12]. Therefore, it is necessary to take into account the effect of random factors in the model. We introduce random perturbations to the parameters b and μ, that is,

$b\text{d}t\to b\text{d}t+{\sigma }_1\text{d}{B}_1\left(t\right)\mathrm{,<mi>\; </mi>}-\mu \text{d}t\to -\mu \text{d}t+{\sigma }_2\text{d}{B}_2\left(t\right)\mathrm{,}$

where $B_i\left(t\right)$ $\left(i=\mathrm{1,2}\right)$ are independent one-dimensional Brownian motions, and ${\sigma }_i$ $\left(i=\mathrm{1,2}\right)$ are the noise intensities. Thus, system (1) takes the following form

$\{\begin{array}{l}\text{d}S\left(t\right)=\left[b\left(S+I\right)-\lambda S{I}^2-dS+\left(1-r\right)\Lambda +\rho I\right]\text{d}t+{\sigma }_1\left[\left(S-S_1\right)+\left(I-I_1\right)\right]\text{d}{B}_1\left(t\right)\mathrm{,}\\ \text{d}I\left(t\right)=\left[\lambda S{I}^2-\left(d+\mu +\rho \right)I+r\Lambda \right]\text{d}t+{\sigma }_2\left(I-I_1\right)\text{d}{B}_2\left(t\right)\mathrm{,}\end{array}$ (3)

and system (2) transforms into

$\{\begin{array}{l}\text{d}S\left(t\right)=\left[b\left(S+I\right)-\lambda S{I}^2-dS+\left(1-r\right)\Lambda +K\right]\text{d}t+{\sigma }_1\left[\left(S-S_2\right)+\left(I-I_2\right)\right]\text{d}{B}_1\left(t\right)\mathrm{,}\\ \text{d}I\left(t\right)=\left[\lambda S{I}^2-\left(d+\mu \right)I+r\Lambda -K\right]\text{d}t+{\sigma }_2\left(I-I_2\right)\text{d}{B}_2\left(t\right)\mathrm{.}\end{array}$ (4)

In recent decades, the problem of bifurcation and periodic oscillations of the solutions of epidemic models has received special attention, as researchers have found that the incidence of several infectious diseases such as rubella, mumps, measles and influenza is periodic. As a result, many scholars have studied the bifurcation of some deterministic systems [4] [5] [6] [13]. In real life, however, environmental noise is also an important part of the ecosystem [14], which has a profound impact on population and microbial studies. Deterministic models, which do not take into account environmental fluctuations, are unable to accurately reflect the actual situation. Since the transmission process of infectious diseases is stochastic in nature, the stochastic differential equation is more suitable to describe the spread of infectious diseases [15]. In previous studies, many scholars have studied the stochastic bifurcation of Marchuk model, HAB model, Logistic model and chemostat model [16] [17] [18] [19] [20] [21]. Most of them utilize stochastic averaging method to simplify stochastic models, analyze local stability by solving Lyapunov exponent, investigate global stability based on singular boundary theory, and analyze stochastic D-bifurcation and P-bifurcation of models. In addition, some scholars have studied pitchfork bifurcation and Hopf bifurcation of stochastic systems [17] [18]. To the best of our knowledge, there is no literature considering the stochastic bifurcation of model (3) and (4). Hence, this paper will focus on the bifurcation in stochastic system (3) and (4).

The rest of the paper is organized as follows. Section 2 provides some preliminaries for this paper. Section 3 simplifies model (3) and (4) by using stochastic averaging method. Section 4 investigates the stochastic stability of model (3) and (4). Section 5 obtains conditions for the occurrence of stochastic Hopf bifurcation and stochastic pitchfork bifurcation. Section 6 presents numerical examples to verify the theoretical results. Finally, Section 7 ends the paper with a conclusion.

2. Preliminaries

In this section, we present some preliminaries that will be used in the subsequent sections to establish stochastic stability and stochastic bifurcation. Before proving the main theorem, we give some definitions.

Definition 2.1. (D-bifurcation) [17] [22]. Dynamical bifurcation is concerned with a family random dynamical systems which is differential and has invariant measure ${\mu }_{\alpha }$ . If there exist a constant ${\alpha }_D$ satisfying in any neighborhood of ${\alpha }_D$ , there exist another constant $\alpha$ and the corresponding invariant measure ${\nu }_{\alpha }\ne {\mu }_{\alpha }$ satisfying ${\nu }_{\alpha }\to {\mu }_{\alpha }$ as $\alpha \to {\alpha }_D$ . Then, the constant ${\alpha }_D$ is a point of D-bifurcation.

Definition 2.2. (P-bifurcation) [17] [22]. Phenomenological bifurcation is concerned with the change in the shape of stationary probability density of a family random dynamical systems as the change of the parameter. If there exists a constant ${\alpha }_0$ satisfying in any neighborhood of ${\alpha }_P$ , there exist other two constant ${\alpha }_1$ , ${\alpha }_2$ and their corresponding invariant measures $p_{{\alpha }_1}$ , $p_{{\alpha }_2}$ satisfying $p_{{\alpha }_1}$ and $p_{{\alpha }_2}$ are not equivalent. Then the constant ${\alpha }_P$ is a point of phenomenological bifurcation.

Definition 2.3. (Stochastic pitchfork bifurcation) [17].

In the viewpoint of dynamical bifurcation: If there exists a constant ${\mu }_0$ satisfies the following conditions:

(i) When $\mu <{\mu }_0$ , the stochastic differential equation has only one invariant measure ${\nu }_0$ , moreover ${\nu }_0$ is stable, i.e., the maximum Lyapunov exponent is negative.

(ii) When $\mu ={\mu }_0$ , the invariant measures ${\nu }_0$ loses its stability and becomes unstable, i.e., the maximum Lyapunov exponent is positive.

(iii) When $\mu >{\mu }_0$ , the stochastic differential equation has three invariant measures ${\nu }_0\mathrm{,<mi>\; </mi>}{\nu }_1$ and ${\nu }_2$ , both ${\nu }_1$ and ${\nu }_2$ are stable, i.e., their maximum Lyapunov exponents are both negative.

If the bifurcation of a stochastic differential equation has the above characteristics, then the stochastic differential equation undergoes a stochastic pitchfork bifurcation at $\mu ={\mu }_0$ .

Definition 2.4. (Stochastic Hopf bifurcation) [17] In the viewpoint of phenomenological bifurcation: The stationary solution of the FPK equation which is corresponded with the stochastic differential equation changes from peak to crater.

3. Model Transformation

3.1. Analysis of Model (3)

Let $x_1=S-S_1$ , $y_1=I-I_1$ , and substitute them into model (3),

$\{\begin{array}{l}\text{d}{x}_1\left(t\right)=[\left(b-d-\lambda I_1^2\right)x_1\left(t\right)+\left(b+\rho -2\lambda S_1{I}_1\right)y_1\left(t\right)-2\lambda I_1{x}_1\left(t\right)y_1\left(t\right)\\ -\lambda S_1{y}_1^2\left(t\right)-\lambda x_1\left(t\right)y_1^2\left(t\right)]\text{d}t+{\sigma }_1\left(x_1\left(t\right)+y_1\left(t\right)\right)\text{d}{B}_1\left(t\right)\mathrm{,}\\ \text{d}{y}_1\left(t\right)\mathrm{<mo>=</mo>}[\lambda I_1^2{x}_1\left(t\right)+\left(2\lambda S_1{I}_1-d-\mu -\rho \right)y_1\left(t\right)+2\lambda I_1{x}_1\left(t\right)y_1\left(t\right)\\ +\lambda S_1{y}_1^2\left(t\right)+\lambda x_1\left(t\right)y_1^2\left(t\right)]\text{d}t+{\sigma }_2{y}_1\left(t\right)\text{d}{B}_2\left(t\right)\mathrm{,}\end{array}$ (5)

where

$\begin{array}{l}{f}_1\left(x_1,y_1,t\right)=\left(b-d-\lambda I_1^2\right)x_1\left(t\right)+\left(b+\rho -2\lambda S_1{I}_1\right)y_1\left(t\right)-2\lambda I_1{x}_1\left(t\right)y_1\left(t\right)\\ -\lambda S_1{y}_1^2\left(t\right)-\lambda x_1\left(t\right)y_1^2\left(t\right),\\ f_2\left(x_1,y_1,t\right)=\lambda I_1^2{x}_1\left(t\right)+\left(2\lambda S_1{I}_1-d-\mu -\rho \right)y_1\left(t\right)+2\lambda I_1{x}_1\left(t\right)y_1\left(t\right)\\ +\lambda S_1{y}_1^2\left(t\right)+\lambda x_1\left(t\right)y_1^2\left(t\right),\\ g_1\left(x_1,y_1,t\right)={\sigma }_1\left(x_1\left(t\right)+y_1\left(t\right)\right),\\ g_2\left(x_1,y_1,t\right)={\sigma }_2{y}_1\left(t\right).\end{array}$

Assume that the coefficient of $x_1^j{y}_1^k$ in $f_i\left(x_1\mathrm{,}{y}_1\mathrm{,}t\right)$ is ${\xi }_{ijk}$ , and the coefficient of $x_1^j{y}_1^k$ in $g_i\left(x_1\mathrm{,}{y}_1\mathrm{,}t\right)$ is ${\eta }_{ijk}$ . Introducing standard rescalings [23],

$x=\overline{x},\mathrm{<mi>\; </mi>}y=\overline{y},\mathrm{<mi>\; </mi>}t=\overline{t},\mathrm{<mi>\; </mi>}{\xi }_{ijk}=\epsilon {\overline{\xi }}_{ijk},\mathrm{<mi>\; </mi>}{\eta }_{ijk}=\sqrt{\epsilon }{\overline{\eta }}_{ijk},$

where $\epsilon \ll 1$ is a sufficiently small, $i\in \left\{\mathrm{1,2}\right\}$ , $j,k\in \left\{\mathrm{0,1,2}\right\}$ . To simplify the notation, we still use $x,y,t,{\xi }_{ijk},{\eta }_{ijk}$ instead of $\overline{x},\overline{y},\overline{t},{\overline{\xi }}_{ijk},{\overline{\eta }}_{ijk}$ , then

$\{\begin{array}{l}\text{d}{x}_1\left(t\right)=\epsilon [\left(b-d-\lambda I_1^2\right)x_1\left(t\right)+\left(b+\rho -2\lambda S_1{I}_1\right)y_1\left(t\right)-2\lambda I_1{x}_1\left(t\right)y_1\left(t\right)\\ -\lambda S_1{y}_1^2\left(t\right)-\lambda x_1\left(t\right)y_1^2\left(t\right)]\text{d}t+\sqrt{\epsilon }{\sigma }_1\left(x_1\left(t\right)+y_1\left(t\right)\right)\text{d}{B}_1\left(t\right)\mathrm{,}\\ \text{d}{y}_1\left(t\right)\mathrm{<mo>=</mo>}\epsilon [\lambda I_1^2{x}_1\left(t\right)+\left(2\lambda S_1{I}_1-d-\mu -\rho \right)y_1\left(t\right)+2\lambda I_1{x}_1\left(t\right)y_1\left(t\right)\\ +\lambda S_1{y}_1^2\left(t\right)+\lambda x_1\left(t\right)y_1^2\left(t\right)]\text{d}t+\sqrt{\epsilon }{\sigma }_2{y}_1\left(t\right)\text{d}{B}_2\left(t\right).\end{array}$ (6)

Next, by applying a polar coordinate transformation to $x_1\left(t\right)$ , $y_1\left(t\right)$ , i.e., $x_1=a\cos \theta$ , $y_1=a\sin \theta$ , and combining it with the Itô formula, we have

$\{\begin{array}{l}\text{d}a=\epsilon \{[\left(b-d-\lambda I_1^2\right){\cos }^2\theta +\left(b+\rho -2\lambda S_1{I}_1+\lambda I_1^2\right)\cos \theta \sin \theta \stackrel{}{}\\ +\left(2\lambda S_1{I}_1-d-\mu -\rho \right){\sin }^2\theta +\frac{1}{2}\left({\sigma }_1^2+{\sigma }_2^2\right){\sin }^2\theta {\cos }^2\theta +{\sigma }_1^2{\sin }^3\theta \cos \theta \\ +\frac{1}{2}{\sigma }_1^2{\sin }^4\theta ]a+[-2\lambda I_1\sin \theta {\cos }^2\theta +\left(2\lambda I_1-\lambda S_1\right){\sin }^2\theta \cos \theta \\ \stackrel{}{}+\lambda S_1{\sin }^3\theta ]a^2+\lambda \left[{\sin }^3\theta \cos \theta -{\sin }^2\theta {\cos }^2\theta \right]a^3\}\text{d}t\\ +\sqrt{\epsilon }{\sigma }_{1}a\left({\cos }^2\theta +\sin \theta \cos \theta \right)\text{d}{B}_1\left(t\right)+\sqrt{\epsilon }{\sigma }_{2}a{\sin }^2\theta \text{d}{B}_2\left(t\right)\mathrm{,}\\ \text{d}\theta =\epsilon \{[\left(2\lambda S_1{I}_1-b-\rho \right){\sin }^2\theta +\lambda I_1^2{\cos }^2\theta +\left(\lambda I_1^2+2\lambda S_1{I}_1-b-\mu -\rho \right)\sin \theta \cos \theta \\ +{\sigma }_1^2\sin \theta {\cos }^3\theta +2{\sigma }_1^2{\sin }^2\theta {\cos }^2\theta +\left({\sigma }_1^2-{\sigma }_2^2\right){\sin }^3\theta \cos \theta ]\\ \left[\lambda S_1{\sin }^3\theta +\lambda \left(2I_1+S_1\right){\sin }^2\theta \cos \theta +2\lambda I_1\sin \theta {\cos }^2\theta \right]a\\ +\lambda \left({\sin }^2\theta {\cos }^2\theta +{\sin }^3\theta \cos \theta \right)a^2\}\text{d}t\\ -\sqrt{\epsilon }{\sigma }_1\left(\sin \theta \cos \theta +{\sin }^2\theta \right)\text{d}{B}_1\left(t\right)+\sqrt{\epsilon }{\sigma }_2\sin \theta \cos \theta \text{d}{B}_2\left(t\right)\mathrm{.}\end{array}$

Based on Khasminskii limiting theorem [24], the response process $\left\{a\left(t\right)\mathrm{,}\theta \left(t\right)\right\}$ of system (3) weakly converges to a two-dimensional Markov diffusion process. Thus, by stochastic averaging method [24], we have the Itô stochastic differential equation as follows.

$\{\begin{array}{l}\text{d}a=m_1\left(a\right)\text{d}t+{\kappa }_1\left(a\right)\text{d}{B}_a\left(t\right)\mathrm{,}\\ \text{d}\theta =m_2\left(a\right)\text{d}t+{\kappa }_2\left(a\right)\text{d}{B}_{\theta }\left(t\right)\mathrm{,}\end{array}$ (7)

where $B_a\left(t\right)$ and $B_{\theta }\left(t\right)$ are independent standard Brownian motions, the drift coefficient $m_i\left(a\right)\left(i=1,2\right)$ and the square of diffusion coefficient ${\kappa }_i^2\left(a\right)\left(i=1,2\right)$ are respectively

$\begin{array}{l}{m}_1\left(a\right)=\epsilon \left[\left(\lambda \left(S_1{I}_1-\frac{1}{2}{I}_1^2\right)+\frac{1}{2}\left(b-2d-\mu -\rho \right)+\frac{1}{16}\left(20{\sigma }_1^2+11{\sigma }_2^2\right)\right)a-\frac{1}{8}\lambda a^3\right]\mathrm{,}\\ m_2\left(a\right)=\epsilon \left[-\lambda \left(S_1{I}_1+\frac{1}{2}{I}_1^2\right)-\frac{1}{4}\left({\sigma }_1^2-2b-2\rho \right)-\frac{1}{8}\lambda a^2\right]\mathrm{,}\\ {\kappa }_1^2\left(a\right)=\frac{1}{8}\epsilon \left(4{\sigma }_1^2+3{\sigma }_2^2\right)a^2\mathrm{,}\\ {\kappa }_2^2\left(a\right)=\frac{1}{8}\epsilon \left(4{\sigma }_1^2+{\sigma }_2^2\right)\mathrm{.}\end{array}$

For convenience, we can rewrite system (7) as

$\{\begin{array}{l}\text{d}a=\left[\left({\nu }_1+\frac{{\nu }_2}{16}\right)a+\frac{{\nu }_3}{8}{a}^3\right]\text{d}t+{\left(\frac{{\nu }_4}{8}{a}^2\right)}^{\frac{1}{2}}\text{d}{B}_a\left(t\right)\mathrm{,}\\ \text{d}\theta \mathrm{<mo>=</mo>}\left({\nu }_5+{\nu }_6{a}^2\right)\text{d}t+{\left(\frac{{\nu }_7}{8}\right)}^{\frac{1}{2}}\text{d}{B}_{\theta }\left(t\right)\mathrm{,}\end{array}$ (8)

where

$\begin{array}{l}{\nu }_1=\lambda \epsilon \left(S_1{I}_1-\frac{1}{2}{I}_1^2\right)-\frac{1}{2}\epsilon \left(\mu +\rho \right),\\ {\nu }_2=\epsilon \left(20{\sigma }_1^2+11{\sigma }_2^2+8b-16d\right),\\ {\nu }_3=-\lambda \epsilon ,\\ {\nu }_4=\epsilon \left(4{\sigma }_1^2+3{\sigma }_2^2\right),\end{array}$

$\begin{array}{l}{\nu }_5=-\lambda \epsilon \left(S_1{I}_1+\frac{1}{2}{I}_1^2\right)-\frac{1}{4}\epsilon \left({\sigma }_1^2-2b-2\rho \right),\\ {\nu }_6=-\frac{1}{8}\lambda \epsilon ,\\ {\nu }_7=\epsilon \left(4{\sigma }_1^2+{\sigma }_2^2\right).\end{array}$

Obviously, the two equations of system (8) are not coupled. Hence, to analyze the stability and bifurcation of the stochastic system (8), we only focus on the following equation

$\text{d}a=\left[\left({\nu }_1+\frac{{\nu }_2}{16}\right)a+\frac{{\nu }_3}{8}{a}^3\right]\text{d}t+{\left(\frac{{\nu }_4}{8}{a}^2\right)}^{\frac{1}{2}}\text{d}{B}_a\left(t\right)\mathrm{.}$ (9)

3.2. Analysis of Model (4)

Let $x_2=S-S_2$ , $y_2=I-I_2$ and substitute them into (4), we can obtain the following stochastic differential equation

$\{\begin{array}{l}\text{d}{x}_2\left(t\right)=[\left(b-d-\lambda I_2^2\right)x_2\left(t\right)+\left(b-2\lambda S_2{I}_2\right)y_2\left(t\right)-2\lambda I_2{x}_2\left(t\right)y_2\left(t\right)\\ -\lambda S_2{y}_2^2\left(t\right)-\lambda x_2\left(t\right)y_2^2\left(t\right)]\text{d}t+{\sigma }_1\left(x_2\left(t\right)+y_2\left(t\right)\right)\text{d}{B}_1\left(t\right)\mathrm{,}\\ \text{d}{y}_2\left(t\right)=[\lambda I_2^2{x}_2\left(t\right)+\left(2\lambda S_2{I}_2-d-\mu \right)y_2\left(t\right)+2\lambda I_2{x}_2\left(t\right)y_2\left(t\right)\\ +\lambda S_2{y}_2^2\left(t\right)+\lambda x_2\left(t\right)y_2^2\left(t\right)]\text{d}t+{\sigma }_2{y}_2\left(t\right)\text{d}{B}_2\left(t\right)\mathrm{.}\end{array}$ (10)

After introducing standard rescalings [23], we have

$\{\begin{array}{l}\text{d}{x}_2\left(t\right)=\epsilon [\left(b-d-\lambda I_2^2\right)x_2\left(t\right)+\left(b-2\lambda S_2{I}_2\right)y_2\left(t\right)-2\lambda I_2{x}_2\left(t\right)y_2\left(t\right)\\ -\lambda S_2{y}_2^2\left(t\right)-\lambda x_2\left(t\right)y_2^2\left(t\right)]\text{d}t+\sqrt{\epsilon }{\sigma }_1\left(x_2\left(t\right)+y_2\left(t\right)\right)\text{d}{B}_1\left(t\right)\mathrm{,}\\ \text{d}{y}_2\left(t\right)=\epsilon [\lambda I_2^2{x}_2\left(t\right)+\left(2\lambda S_2{I}_2-d-\mu \right)y_2\left(t\right)+2\lambda I_2{x}_2\left(t\right)y_2\left(t\right)\\ +\lambda S_2{y}_2^2\left(t\right)+\lambda x_2\left(t\right)y_2^2\left(t\right)]\text{d}t+\sqrt{\epsilon }{\sigma }_2{y}_2\left(t\right)\text{d}{B}_2\left(t\right)\mathrm{.}\end{array}$ (11)

Let $x_2=\tilde{a}\cos \tilde{\theta }$ , $y_2=\tilde{a}\sin \tilde{\theta }$ and use the Itô formula for model (11). For convenience, we still use a and $\theta$ instead of $\tilde{a}$ and $\tilde{\theta }$ , and obtain

$\{\begin{array}{l}\text{d}a=\epsilon \{[\left(b-d-\lambda I_2^2\right){\cos }^2\theta +\left(b-2\lambda S_2{I}_2+\lambda I_2^2\right)\cos \theta \sin \theta \stackrel{}{}\\ +\left(2\lambda S_2{I}_2-d-\mu \right){\sin }^2\theta +\frac{1}{2}\left({\sigma }_1^2+{\sigma }_2^2\right){\sin }^2\theta {\cos }^2\theta +{\sigma }_1^2{\sin }^3\theta \cos \theta \\ +\frac{1}{2}{\sigma }^2{\sin }^4\theta ]a+[-2\lambda I_2\sin \theta {\cos }^2\theta +\left(2\lambda I_2-\lambda S_2\right){\sin }^2\theta \cos \theta \\ \stackrel{}{}+\lambda S_2{\sin }^3\theta ]a^2+\lambda \left[{\sin }^3\theta \cos \theta -{\sin }^2\theta {\cos }^2\theta \right]a^3\}\text{d}t\\ +\sqrt{\epsilon }{\sigma }_{1}a\left({\cos }^2\theta +\sin \theta \cos \theta \right)\text{d}{B}_1\left(t\right)+\sqrt{\epsilon }{\sigma }_{2}a{\sin }^2\theta \text{d}{B}_2\left(t\right)\mathrm{,}\\ \text{d}\theta =\epsilon \{[\left(2\lambda S_2{I}_2-b\right){\sin }^2\theta +\lambda I_2^2{\cos }^2\theta +\left(\lambda I_2^2+2\lambda S_2{I}_2-b-\mu \right)\sin \theta \cos \theta \\ +{\sigma }_1^2\sin \theta {\cos }^3\theta +2{\sigma }_1^2{\sin }^2\theta {\cos }^2\theta +\left({\sigma }_1^2-{\sigma }_2^2\right){\sin }^3\theta \cos \theta ]\\ \left[\lambda S_2{\sin }^3\theta +\lambda \left(2I_2+S_2\right){\sin }^2\theta \cos \theta +2\lambda I_2\sin \theta {\cos }^2\theta \right]a\\ +\lambda \left({\sin }^2\theta {\cos }^2\theta +{\sin }^3\theta \cos \theta \right)a^2\}\text{d}t\\ -\sqrt{\epsilon }{\sigma }_1\left(\sin \theta \cos \theta +{\sin }^2\theta \right)\text{d}{B}_1\left(t\right)+\sqrt{\epsilon }{\sigma }_2\sin \theta \cos \theta \text{d}{B}_2\left(t\right)\mathrm{.}\end{array}$(12)

Then according to the stochastic averaging method, we have

$\{\begin{array}{l}\text{d}a=m_3\left(a\right)\text{d}t+{\kappa }_1\left(a\right)\text{d}{B}_a\left(t\right)\mathrm{,}\hfill \\ \text{d}\theta =m_4\left(a\right)\text{d}t+{\kappa }_2\left(a\right)\text{d}{B}_{\theta }\left(t\right)\mathrm{,}\hfill \end{array}$ (13)

where $B_a\left(t\right)$ and $B_{\theta }\left(t\right)$ are independent standard Brownian motions. Moreover, drift coefficient $m_i\left(a\right)\left(i=3,4\right)$ and the square of diffusion coefficient ${\kappa }_i^2\left(a\right)\left(i=1,2\right)$ are respectively

$\begin{array}{l}{m}_3\left(a\right)=\epsilon \left[\left(\lambda \left(S_2{I}_2-\frac{1}{2}{I}_2^2\right)+\frac{1}{2}\left(b-2d-\mu \right)+\frac{1}{16}\left(20{\sigma }_1^2+11{\sigma }_2^2\right)\right)a-\frac{1}{8}\lambda a^3\right],\\ m_4\left(a\right)=\epsilon \left[-\lambda \left(S_2{I}_2+\frac{1}{2}{I}_2^2\right)-\frac{1}{4}\left({\sigma }_1^2-2b\right)-\frac{1}{8}\lambda a^2\right],\end{array}$

$\begin{array}{l}{\kappa }_1^2\left(a\right)=\frac{1}{8}\epsilon \left(4{\sigma }_1^2+3{\sigma }_2^2\right)a^2,\\ {\kappa }_2^2\left(a\right)=\frac{1}{8}\epsilon \left(4{\sigma }_1^2+{\sigma }_2^2\right).\end{array}$

For convenience, we can rewrite model (13) as follows

$\{\begin{array}{l}\text{d}a=\left[\left({\delta }_1+\frac{{\nu }_2}{16}\right)a+\frac{{\nu }_3}{8}{a}^3\right]\text{d}t+{\left(\frac{{\nu }_4}{8}{a}^2\right)}^{\frac{1}{2}}\text{d}{B}_a\left(t\right)\mathrm{,}\\ \text{d}\theta =\left({\delta }_2+{\nu }_6{a}^2\right)\text{d}t+{\left(\frac{{\nu }_7}{8}\right)}^{\frac{1}{2}}\text{d}{B}_{\theta }\left(t\right)\mathrm{,}\end{array}$ (14)

where

$\begin{array}{l}{\delta }_1=\lambda \epsilon \left(S_2{I}_2-\frac{1}{2}{I}_2^2\right)-\frac{1}{2}\epsilon \mu \mathrm{,}\\ {\delta }_2=-\lambda \epsilon \left(S_2{I}_2+\frac{1}{2}{I}_2^2\right)-\frac{1}{4}\epsilon \left({\sigma }_1^2-2b\right)\mathrm{,}\end{array}$

the other parameters are defined the same as model (8). Similarly, we study the stability and bifurcation of the stochastic model (14) and only need to discuss the following equation

$\text{d}a=\left[\left({\delta }_1+\frac{{\nu }_2}{16}\right)a+\frac{{\nu }_3}{8}{a}^3\right]\text{d}t+{\left(\frac{{\nu }_4}{8}{a}^2\right)}^{\frac{1}{2}}\text{d}{B}_a\left(t\right)\mathrm{.}$ (15)

4. Stochastic Stability

From the discussion in the previous section, we know that when $0\le I\le I_k$ ($I>I_k$ ), the stochastic stability of system (3) (or system (4)) at the positive equilibrium point $E_1\left(S_1\mathrm{,}{I}_1\right)\left(E_2\left(S_2\mathrm{,}{I}_2\right)\right)$ is equivalent to the stochastic stability of system (9) (or system (15)) at the equilibrium point $a=0$ .

4.1. Stochastic Local Stability

To describe the local stability of positive equilibrium points of system (3) (or system (4)), we investigate the stability of the average amplitude Equation (9) (or Equation (15)) at equilibrium point $a=0$ by computing the maximum Lyapunov exponent of the linearized system.

Theorem 4.1. The following conclusions hold.

(i) When ${\nu }_4>16{\nu }_1+{\nu }_2$ , the trivial solution of system (9) is asymptotically stable with probability one. Then the stochastic system (3) (or system (5)) is stable at the equilibrium point E₁ (or (0, 0)).

(ii) When ${\nu }_4<16{\nu }_1+{\nu }_2$ , the trivial solution of system (9) is unstable with probability one. Then the stochastic system (3) (or system (5)) is unstable at the equilibrium point E₁ (or (0, 0)).

Proof. Firstly, the linearized equation of (9) at $a=0$ is as follows.

$\text{d}a=\left({\nu }_1+\frac{{\nu }_2}{16}\right)a\text{d}t+{\left(\frac{{\nu }_4}{8}{a}^2\right)}^{\frac{1}{2}}\text{d}{B}_a\left(t\right)\mathrm{.}$ (16)

Then it follows that the solution of system (16) is

$a\left(t\right)=a\left(0\right)\exp \left[{\displaystyle {\int }_0^t}\left({\nu }_1+\frac{{\nu }_2}{16}-\frac{{\nu }_4}{16}\right)\text{d}s+{\displaystyle {\int }_0^t}{\left(\frac{{\nu }_4}{8}{a}^2\right)}^{\frac{1}{2}}\text{d}{B}_a\left(s\right)\right]\mathrm{.}$ (17)

Consequently, the maximum Lyapunov exponent of system (16) is

$\kappa =\underset{t\to +\infty }{lim}\frac{1}{t}\ln \Vert a\left(t\right)\Vert ={\nu }_1+\frac{{\nu }_2}{16}-\frac{{\nu }_4}{16}\mathrm{.}$

Further, we have

(i) When ${\nu }_4>16{\nu }_1+{\nu }_2$ , i.e., $\kappa <0$ , utilizing the Oseledec’s multiplicative ergodic theorem [25], the trivial solution of system (16) is asymptotically stable with probability one. Since Equation (16) exhibits robustness, we conclude that the stochastic system (9) is stable at $a=0$ .

(ii) When ${\nu }_4<16{\nu }_1+{\nu }_2$ , i.e., $\kappa >0$ , thus the trivial solution of the linear Itô stochastic differential Equation (16) is unstable with probability one, i.e., the stochastic system (9) is unstable at $a=0$ .

The proof is now complete.

Theorem 4.2. The following conclusions hold.

(i) When ${\nu }_4>16{\delta }_1+{\nu }_2$ , the trivial solution of system (15) is asymptotically stable with probability one. Then the stochastic system (4) (or system (10)) is stable at the equilibrium point E₂ (or (0, 0)).

(ii) When ${\nu }_4<16{\delta }_1+{\nu }_2$ , the trivial solution of system (15) is unstable with probability one. Then the stochastic system (4) (or system (10)) is unstable at the equilibrium point E₂ (or (0, 0)).

Proof. Since the proof is similar to that of Theorem 4.1, we omit it.

4.2. Stochastic Global Stability

In this subsection, we will apply the singular boundary theory to investigate the global stability of the stochastic system (3) and (4).

Theorem 4.3. When $16{\nu }_1+{\nu }_2<{\nu }_4$ , system (3) is globally stable at the equilibrium point E₁.

Proof. As $a\to 0^{+}$ , we have

$\begin{array}{l}{m}_1\left(a\right)=\left({\nu }_1+\frac{{\nu }_2}{16}\right)a+o\left(a\right)\mathrm{,}\\ {\kappa }_1^2\left(a\right)=\frac{{\nu }_4}{8}{a}^2\mathrm{.}\end{array}$

When $a=0$ , we can find ${\kappa }_1^2\left(a\right)=0$ . Meanwhile, it’s easy to notice that $m_1\left(a\right)=+\infty$ if $a=+\infty$ . Based on the singular boundary theory, $a=0$ and $a=+\infty$ are the first and second kind of singular boundaries of system (9), respectively. Therefore, we can calculate the diffusion exponent ${\alpha }_l$ , drift exponent ${\beta }_l$ and character value $c_l$ of boundary $a=0$ and the results are as follows.

${\alpha }_l=2,\mathrm{<mi>\; </mi>}{\beta }_l=1,\mathrm{<mi>\; </mi>}{c}_l=\underset{a\to 0^{+}}{\lim }\frac{2m_1\left(a\right)a^{{\alpha }_l-{\beta }_l}}{{\kappa }_1^2\left(a\right)}=\frac{16{\nu }_1+{\nu }_2}{{\nu }_4}.$

Thus, we have that

(i) if $16{\nu }_1+{\nu }_2<{\nu }_4$ , i.e., $c_l<1$ , the left boundary $a=0$ is attractively natural;

(ii) if $16{\nu }_1+{\nu }_2={\nu }_4$ , i.e., $c_l=1$ , the left boundary $a=0$ is strictly natural;

(iii) if $16{\nu }_1+{\nu }_2>{\nu }_4$ , i.e., $c_l>1$ , the left boundary $a=0$ is repulsively natural.

As $a\to +\infty$ , it implies that $m_1\left(a\right)$ and ${\kappa }_1^2\left(a\right)$ asymptotically converge to

$\begin{array}{l}{m}_1\left(a\right)=\frac{{\nu }_3}{8}{a}^3\mathrm{,}\\ {\kappa }_1^2\left(a\right)=\frac{{\nu }_4}{8}{a}^2\mathrm{,}\end{array}$

respectively. We can also calculate the diffusion exponent ${\alpha }_r$ , drift exponent ${\beta }_r$ and character value $c_r$ of boundary $a=+\infty$ , and the results are as follows

${\alpha }_r=2,\mathrm{<mi>\; </mi>}{\beta }_r=3,c_r=-\underset{a\to +\infty }{\lim }\frac{2m_1\left(a\right)a^{{\alpha }_r-{\beta }_r}}{{\kappa }_1^2\left(a\right)}=-\frac{2{\nu }_3}{{\nu }_4}.$

Due to ${\nu }_3=-\lambda \epsilon <0$ , it follows that $m_1\left(+\infty \right)=-\infty <0$ . Therefore, $a=+\infty$ is an entrance boundary.

Hence, when the boundary $a=0$ is attractively natural boundary and the boundary $a=+\infty$ is an entrance boundary, all solution curves enter the inner system from the right boundary and are subsequently attracted by the left boundary, which suggests that the trivial solution $a=0$ of system (9) is globally stable [16]. Thus, the positive equilibrium point E₁ of system (3) is globally stable if $16{\nu }_1+{\nu }_2<{\nu }_4$ . The proof is now complete.

Theorem 4.4. When $16{\delta }_1+{\nu }_2<{\nu }_4$ , system (4) is globally stable at the equilibrium point E₂.

Proof. Since the proof is similar to that of Theorem 4.3, we omit it.

5. Stochastic Bifurcation

5.1. Stochastic Hopf Bifurcation

Next, we will investigate the change of steady-state probability density function to study the Hopf bifurcation of system (9) and (15).

Firstly, from Definition 2.4, we study the P-bifurcation of the stochastic system (9). We get its Fokker-Planck equation

$\frac{\partial p\left(a\right)}{\partial t}=-\frac{\partial }{\partial a}\left[\left(\left({\nu }_1+\frac{{\nu }_2}{16}\right)a+\frac{{\nu }_3}{8}{a}^3\right)p\left(a\right)\right]+\frac{1}{2}\frac{{\partial }^2}{\partial a^2}\left(\frac{{\nu }_4}{8}{a}^{2}p\left(a\right)\right)\mathrm{.}$

And the stationary probability density $p\left(a\right)$ satisfies the following degenerate equation

$-\frac{\partial }{\partial a}\left[\left(\left({\nu }_1+\frac{{\nu }_2}{16}\right)a+\frac{{\nu }_3}{8}{a}^3\right)p\left(a\right)\right]+\frac{1}{2}\frac{{\partial }^2}{\partial a^2}\left(\frac{{\nu }_4}{8}{a}^{2}p\left(a\right)\right)=0.$

Then we obtain

$p\left(a\right)=\{\begin{array}{ll}\delta \left(a\right)\mathrm{,}\hfill & {\nu }_4\ge 16{\nu }_1+{\nu }_2\mathrm{,}\hfill \\ \frac{a^{\frac{16{\nu }_1+{\nu }_2-2{\nu }_4}{{\nu }_4}}\exp \left(\frac{{\nu }_3}{{\nu }_4}{a}^2\right)}{\Gamma \left(\frac{16{\nu }_1+{\nu }_2-{\nu }_4}{2{\nu }_4}\right){\left(-\frac{{\nu }_4}{{\nu }_3}\right)}^{\frac{16{\nu }_1+{\nu }_2-{\nu }_4}{2{\nu }_4}}}\mathrm{,}\hfill & {\nu }_4<16{\nu }_1+{\nu }_2\mathrm{.}\hfill \end{array}$ (18)

According to Namachchivaya’s theory [26], if $p\left(a\right)$ reaches its maximum value at point $a^{\mathrm{<mo>*</mo>}}$ , the sample trajectory will tend to remain in the neighborhood of $a^{\mathrm{<mo>*</mo>}}$ for a longer time with a bigger probability, which implies that $a^{\mathrm{<mo>*</mo>}}$ is stable in terms of probability. If $p\left(a\right)$ has a minimum value, the situation just be the opposite. If we have

${\frac{\text{d}p\left(a\right)}{\text{d}a}|}_{a=a^{\mathrm{<mo>*</mo>}}}=\mathrm{0,<mi>\; </mi>}{\frac{{\text{d}}^{2}p\left(a\right)}{\text{d}{a}^2}|}_{a=a^{\mathrm{<mo>*</mo>}}}<\mathrm{0,}$ (19)

then $p\left(a\right)$ has a maximum value at $a^{\mathrm{<mo>*</mo>}}$ . To obtain the extreme value point of the probability density $p\left(a\right)$ , we need to solve $\text{d}p\left(a\right)/\text{d}a=0$ , that is

$\frac{\left(\frac{16{\nu }_1+{\nu }_2-2{\nu }_4}{{\nu }_4}+\frac{2{\nu }_3}{{\nu }_4}{a}^2\right)a^{\frac{16{\nu }_1+{\nu }_2-3{\nu }_4}{{\nu }_4}}\exp \left(\frac{{\nu }_3}{{\nu }_4}{a}^2\right)}{\Gamma \left(\frac{16{\nu }_1+{\nu }_2-{\nu }_4}{2{\nu }_4}\right){\left(-\frac{{\nu }_4}{{\nu }_3}\right)}^{\frac{16{\nu }_1+{\nu }_2-{\nu }_4}{2{\nu }_4}}}=0.$ (20)

Thus we have $a=0$ or

$a=a^{\mathrm{<mo>*</mo>}}=\sqrt{\frac{2{\nu }_4-16{\nu }_1-{\nu }_2}{2{\nu }_3}}\mathrm{,<mi>\; </mi>}\text{for}\mathrm{<mi>\; </mi>}{\nu }_4<8{\nu }_1+\frac{1}{2}{\nu }_2\mathrm{.}$ (21)

Through an analysis of the value range of ${\nu }_4$ , we get the following three cases.

(i) When $8{\nu }_1+{\nu }_2/2\le {\nu }_4<16{\nu }_1+{\nu }_2$ , the probability density function $p\left(a\right)$ tends to infinity as $a\to 0^{+}$ . Then, the solution trajectories of the stochastic system (9) are concentrated in a neighborhood of point $a=0$ .

(ii) When $\left(16{\nu }_1+{\nu }_2\right)/3\le {\nu }_4<8{\nu }_1+{\nu }_2/2$ , the probability function $p\left(a\right)$ reaches a maximum value at point $a^{\mathrm{<mo>*</mo>}}=\sqrt{\frac{2{\nu }_4-16{\nu }_1-{\nu }_2}{2{\nu }_3}}$ and a minimum value at point $a=0$ , but the derivative of $p\left(a\right)$ does not exist at point $a=0$ . In this situation, the solution trajectories of the stochastic system (9) are concentrated in a neighborhood of point $a^{\mathrm{<mo>*</mo>}}$ .

(iii) When ${\nu }_4<\left(16{\nu }_1+{\nu }_2\right)/3$ , the probability function $p\left(a\right)$ reaches a maximum value at point $a^{\mathrm{<mo>*</mo>}}=\sqrt{\frac{2{\nu }_4-16{\nu }_1-{\nu }_2}{2{\nu }_3}}$ and still has a minimum value at point $a=0$ . At the same time, the derivative of the probability function $p\left(a\right)$ exists at point $a=0$ .

In summary, we obtain the following result.

Theorem 5.1. System (9) undergoes a stochastic P-bifurcation as the parameter ${\nu }_4$ passes through the values of $8{\nu }_1+{\nu }_2/2$ and $\left(16{\nu }_1+{\nu }_2\right)/3$ .

Next, by $\rho \left(x\mathrm{,}y\right)=\left|J\right|\rho \left(a\mathrm{,}\theta \right)$ and $p\left(a\right)={\displaystyle {\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}}\rho \left(a\mathrm{,}\theta \right)\text{d}\theta$ [27], where the determinant of the Jacobian matrix J of the nonlinear transformation is given by $\left|J\right|=1/a$ and $\rho$ is the joint probability density to a and $\theta$ , we have

$\rho \left(x\mathrm{,}y\right)=\frac{{\left(x^2+y^2\right)}^{\frac{16{\nu }_1+{\nu }_2-3{\nu }_4}{2{\nu }_4}}\exp \left(\frac{{\nu }_3}{{\nu }_4}\left(x^2+y^2\right)\right)}{\pi \Gamma \left(\frac{16{\nu }_1+{\nu }_2-{\nu }_4}{2{\nu }_4}\right){\left(-\frac{{\nu }_4}{{\nu }_3}\right)}^{\frac{16{\nu }_1+{\nu }_2-{\nu }_4}{2{\nu }_4}}}\mathrm{.}$

Let $\nabla$ denotes the gradient in ${ℝ}^2$ . Computing $\nabla \rho \left(x\mathrm{,}y\right)=0$ yields critical point of $\rho \left(x\mathrm{,}y\right)$ . Therefore, we have the following results.

(i) When ${\nu }_4\ge \left(16{\nu }_1+{\nu }_2\right)/3$ , the joint probability density $\rho \left(x\mathrm{,}y\right)$ approaches infinity as $x\to 0$ and $y\to 0$ .

(ii) When $\left(16{\nu }_1+{\nu }_2\right)/5\le {\nu }_4<\left(16{\nu }_1+{\nu }_2\right)/3$ , a maximum value of $\rho \left(x\mathrm{,}y\right)$ is reached at the points of cycle $x^2+y^2=\left(3{\nu }_4-16{\nu }_1-{\nu }_2\right)/2{\nu }_3$ , and a minimum value appears at point $O\left(\mathrm{0,0}\right)$ . Meanwhile, it follows that the partial derivatives of the joint probability density $\rho \left(x\mathrm{,}y\right)$ do not exist at point $O\left(\mathrm{0,0}\right)$ .

(iii) When ${\nu }_4<\left(16{\nu }_1+{\nu }_2\right)/5$ , the joint probability density $\rho \left(x\mathrm{,}y\right)$ reaches a maximum value at the points of cycle $x^2+y^2=\left(3{\nu }_4-16{\nu }_1-{\nu }_2\right)/2{\nu }_3$ , and a minimum value at point $O\left(\mathrm{0,0}\right)$ , where the partial derivatives of the joint probability density $\rho \left(x\mathrm{,}y\right)$ exist.

Thus, we obtain the following conclusion.

Theorem 5.2. The stochastic system (9) exhibits a P-bifurcation as the parameter ${\nu }_4$ passes through the critical values of $\left(16{\nu }_1+{\nu }_2\right)/3$ and $\left(16{\nu }_1+{\nu }_2\right)/5$ .

Finally, we study P-bifurcation of the stochastic system (15). Similar to the analysis of system (9), we can obtain its Fokker-Planck equation

$\frac{\partial \zeta \left(a\right)}{\partial t}=-\frac{\partial }{\partial a}\left[\left(\left({\delta }_1+\frac{{\nu }_2}{16}\right)a+\frac{{\nu }_3}{8}{a}^3\right)\zeta \left(a\right)\right]+\frac{1}{2}\frac{{\partial }^2}{\partial a^2}\left(\frac{{\nu }_4}{8}{a}^2\zeta \left(a\right)\right)\mathrm{.}$

Then $\zeta \left(a\right)$ is governed by the following degenerate equation

$-\frac{\partial }{\partial a}\left[\left(\left({\delta }_1+\frac{{\nu }_2}{16}\right)a+\frac{{\nu }_3}{8}{a}^3\right)\zeta \left(a\right)\right]+\frac{1}{2}\frac{{\partial }^2}{\partial a^2}\left(\frac{{\nu }_4}{8}{a}^2\zeta \left(a\right)\right)=0.$

Therefore, $\zeta \left(a\right)$ can be written as

$\zeta \left(a\right)=\{\begin{array}{ll}\delta \left(a\right),\hfill & {\nu }_4\ge 16{\delta }_1+{\nu }_2,\hfill \\ \frac{a^{\frac{16{\delta }_1+{\nu }_2-2{\nu }_4}{{\nu }_4}}\exp \left(\frac{{\nu }_3}{{\nu }_4}{a}^2\right)}{\Gamma \left(\frac{16{\delta }_1+{\nu }_2-{\nu }_4}{2{\nu }_4}\right){\left(-\frac{{\nu }_4}{{\nu }_3}\right)}^{\frac{16{\delta }_1+{\nu }_2-{\nu }_4}{2{\nu }_4}}},\hfill & {\nu }_4<16{\delta }_1+{\nu }_2.\hfill \end{array}$ (22)

Computing $\text{d}\zeta \left(a\right)/\text{d}a=0$ yields that

$\frac{\left[\frac{16{\delta }_1+{\nu }_2-2{\nu }_4}{{\nu }_4}+\frac{2{\nu }_3}{{\nu }_4}{a}^2\right]a^{\frac{16{\delta }_1+{\nu }_2-3{\nu }_4}{{\nu }_4}}\exp \left(\frac{{\nu }_3}{{\nu }_4}{a}^2\right)}{\Gamma \left(\frac{16{\delta }_1+{\nu }_2-{\nu }_4}{2{\nu }_4}\right){\left(-\frac{{\nu }_4}{{\nu }_3}\right)}^{\frac{16{\delta }_1+{\nu }_2-{\nu }_4}{2{\nu }_4}}}=0.$ (23)

Then we have $a=0$ or

$a=a^{\mathrm{<mo>*</mo>}}=\sqrt{\frac{2{\nu }_4-16{\delta }_1-{\nu }_2}{2{\nu }_3}}\mathrm{,<mi>\; </mi>}\text{for}\mathrm{<mi>\; </mi>}{\nu }_4<8{\delta }_1+\frac{1}{2}{\nu }_2\mathrm{.}$ (24)

Hence, we have the following three cases.

(i) When $8{\delta }_1+{\nu }_2/2\le {\nu }_4<16{\delta }_1+{\nu }_2$ , the probability density function $\zeta \left(a\right)$ tends to infinity as $a\to 0^{+}$ . This shows that the solution trajectories of system (15) are concentrated in a neighborhood of point $a=0$ .

(ii) When $\left(16{\delta }_1+{\nu }_2\right)/3\le {\nu }_4<8{\delta }_1+{\nu }_2/2$ , the probability function $\zeta \left(a\right)$ reaches a maximum value at $a^{\mathrm{<mo>*</mo>}}=\sqrt{\frac{2{\nu }_4-16{\delta }_1-{\nu }_2}{2{\nu }_3}}$ and a minimum value at $a=0$ , but the derivative of $\zeta \left(a\right)$ does not exist at $a=0$ . This indicates that the solution trajectories of system (15) are concentrated in a neighborhood of point $a^{\mathrm{<mo>*</mo>}}$ .

(iii) When ${\nu }_4<\left(16{\delta }_1+{\nu }_2\right)/3$ , the probability function $\zeta \left(a\right)$ has a maximum value at point $a^{\mathrm{<mo>*</mo>}}=\sqrt{\frac{2{\nu }_4-16{\delta }_1-{\nu }_2}{2{\nu }_3}}$ and still has a minimum value at point $a=0$ , but its derivative exists at point $a=0$ .

Consequently, we obtain the following result.

Theorem 5.3 System (15) undergoes a stochastic P-bifurcation as the parameter ${\nu }_4$ passes through the values of $8{\delta }_1+{\nu }_2/2$ and $\left(16{\delta }_1+{\nu }_2\right)/3$ .

In the following, we still focus on the main features of the joint probability density $\tau \left(u\mathrm{,}v\right)$ . Utilizing the relation $\tau \left(u\mathrm{,}v\right)=\left|J\right|\tilde{\rho }\left(a\mathrm{,}\theta \right)$ and $\zeta \left(a\right)={\displaystyle {\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}}\tilde{\rho }\left(a\mathrm{,}\theta \right)\text{d}\theta$, where $\tilde{\rho }$ is the joint probability density to a and $\theta$ [27], we have

$\tau \left(u\mathrm{,}v\right)=\frac{{\left(u^2+v^2\right)}^{\frac{16{\delta }_1+{\nu }_2-3{\nu }_4}{2{\nu }_4}}\exp \left(\frac{{\nu }_3}{{\nu }_4}\left(u^2+v^2\right)\right)}{\pi \Gamma \left(\frac{16{\delta }_1+{\nu }_2-{\nu }_4}{2{\nu }_4}\right){\left(-\frac{{\nu }_4}{{\nu }_3}\right)}^{\frac{16{\delta }_1+{\nu }_2-{\nu }_4}{2{\nu }_4}}}\mathrm{.}$

Through the analysis, we get the following results.

(i) If ${\nu }_4\ge \left(16{\delta }_1+{\nu }_2\right)/3$ , the joint probability density $\tau \left(u\mathrm{,}v\right)$ tends to infinity as $u\to 0$ and $v\to 0$ .

(ii) If $\left(16{\delta }_1+{\nu }_2\right)/5\le {\nu }_4<\left(16{\delta }_1+{\nu }_2\right)/3$ , a maximum value of $\tau \left(u\mathrm{,}v\right)$ is reached at the points of cycle $u^2+v^2=\left(3{\nu }_4-16{\delta }_1-{\nu }_2\right)/2{\nu }_3$ , and a minimum value appears at $O\left(\mathrm{0,0}\right)$ . Meanwhile, it follows that the partial derivatives of the joint probability density $\tau \left(u\mathrm{,}v\right)$ do not exist at the original point $O\left(\mathrm{0,0}\right)$ .

(iii) If ${\nu }_4<\left(16{\delta }_1+{\nu }_2\right)/5$ , the joint probability density $\tau \left(u\mathrm{,}v\right)$ has a maximum value at the points of cycle $u^2+v^2=\left(3{\nu }_4-16{\delta }_1-{\nu }_2\right)/2{\nu }_3$ , and a minimum value occurs at point $O\left(\mathrm{0,0}\right)$ , where the partial derivatives of the joint probability density $\tau \left(u\mathrm{,}v\right)$ exist.

Thus, we obtain the following conclusion.

Theorem 5.4. The stochastic system (15) exhibits a P-bifurcation as the parameter ${\nu }_4$ passes through the values of $\left(16{\delta }_1+{\nu }_2\right)/3$ and $\left(16{\delta }_1+{\nu }_2\right)/5$ .

5.2. Stochastic Pitchfork Bifurcation

Theorem 5.5. Let $\alpha ={\nu }_1+{\nu }_2/16-{\nu }_4/16$ . If the following three conditions hold,

(i) for $\alpha \in ℝ$ , the integrability conditions (IC) for ${\mu }_{1,\omega }^{\alpha }={\delta }_0$ is trivially satisfied and $\lambda \left({\mu }_1^{\alpha }\right)=\alpha$ ;

(ii) for $\alpha >0$ , ${\mu }_{\mathrm{2,}\omega }^{\alpha }={\delta }_{-d_{\alpha }^{-}\left(\omega \right)}$ is ${ℱ}_{-\infty }^0$ measurable, $\lambda \left({\mu }_2^{\alpha }\right)=-2\alpha <0$ ;

(iii) for $\alpha >0$ , ${\mu }_{\mathrm{3,}\omega }^{\alpha }={\delta }_{d_{\alpha }^{-}\left(\omega \right)}$ is ${ℱ}_{-\infty }^0$ measurable, $\lambda \left({\mu }_3^{\alpha }\right)=-2\alpha <0$ .

Then system (9) undergoes a stochastic pitchfork bifurcation.

Proof. Let ${\eta }_t={\left(-\frac{{\nu }_3}{8}\right)}^{\frac{1}{2}}a$ . Then system (9) can be translated into

$\text{d}{\eta }_t=\left[\left({\nu }_1+\frac{{\nu }_2}{16}\right){\eta }_t-{\eta }_t^3\right]\text{d}t+{\left(\frac{{\nu }_4}{8}\right)}^{\frac{1}{2}}{\eta }_t\text{d}{B}_a\left(t\right)\mathrm{,}$

which is equivalent to the following Stratonovich stochastic differential equation

$\text{d}{\eta }_t=\left[\left({\nu }_1+\frac{{\nu }_2}{16}-\frac{{\nu }_4}{16}\right){\eta }_t-{\eta }_t^3\right]\text{d}t+{\left(\frac{{\nu }_4}{8}\right)}^{\frac{1}{2}}{\eta }_t\circ \text{d}{B}_a\left(t\right)\mathrm{.}$ (25)

Let $\alpha ={\nu }_1+{\nu }_2/16-{\nu }_4/16$ and $\sigma ={\left({\nu }_4/8\right)}^{\frac{1}{2}}$ . Then system (25) can be written as

$\text{d}{\eta }_t=\left(\alpha {\eta }_t-{\eta }_t^3\right)\text{d}t+\sigma {\eta }_t\circ \text{d}{B}_a\left(t\right)\mathrm{.}$ (26)

Then, a local random dynamical system is generated as follows.

${\phi }_{\alpha }\left(t\mathrm{,}\omega \right)\eta =\frac{\eta \exp \left(\alpha t+\sigma B_a\left(t\right)\right)}{{\left(1+2{\eta }^2{\displaystyle {\int }_0^t}\exp \left(2\alpha s+2\sigma B_a\left(s\right)\right)\text{d}s\right)}^{\frac{1}{2}}}\mathrm{,}$ (27)

where $\eta$ is the initial value of ${\eta }_t$ . The domain $D_{\alpha }\left(t,\omega \right)$ and range $R_{\alpha }\left(t\mathrm{,}\omega \right)$ of ${\phi }_{\alpha }\left(t\mathrm{,}\omega \right):D_{\alpha }\left(t\mathrm{,}\omega \right)\to R_{\alpha }\left(t\mathrm{,}\omega \right)$ can be determined as follows

$D_{\alpha }\left(t,\omega \right)=\{\begin{array}{ll}ℝ,\hfill & t\ge 0,\hfill \\ \left(-d_{\alpha }\left(t,\omega \right),d_{\alpha }\left(t,\omega \right)\right),\hfill & t<0,\hfill \end{array}$ (28)

and

$R_{\alpha }\left(t,\omega \right)=D_{\alpha }\left(-t,\vartheta \left(t\right)\omega \right)=\{\begin{array}{ll}\left(-r_{\alpha }\left(t,\omega \right),r_{\alpha }\left(t,\omega \right)\right),\hfill & t>0,\hfill \\ ℝ,\hfill & t\le 0,\hfill \end{array}$ (29)

where $\vartheta (\cdot )$ means a flow of Ω and

$d_{\alpha }\left(t\mathrm{,}\omega \right)=\frac{1}{{\left(2\left|{\displaystyle {\int }_0^t}\exp \left(2\alpha s+2\sigma B_a\left(s\right)\right)\text{d}s\right|\right)}^{\frac{1}{2}}}>\mathrm{0,}$

and

$r_{\alpha }\left(t\mathrm{,}\omega \right)=d_{\alpha }\left(-t\mathrm{,}\vartheta \left(t\right)\omega \right)=\frac{\exp \left(\alpha t+\sigma B_a\left(t\right)\right)}{{\left(2\left|{\displaystyle {\int }_0^t}\exp \left(2\alpha s+2\sigma B_a\left(s\right)\right)\text{d}s\right|\right)}^{\frac{1}{2}}}>0.$

Define $E_{\alpha }\left(\omega \right)\mathrm{<mo>:</mo>}={\displaystyle {\cap }_{t\in ℝ}{D}_{\alpha }\left(t\mathrm{,}\omega \right)}$ , which is the collection of initial values $\eta$ ensuring the non-explosion of ${\phi }_{\alpha }\left(t\mathrm{,}\omega \right)\eta$ , i.e.,

$E_{\alpha }\left(\omega \right)=\{\begin{array}{ll}\left(-d_{\alpha }^{-}\left(t\mathrm{,}\omega \right)\mathrm{,}{d}_{\alpha }^{-}\left(t\mathrm{,}\omega \right)\right)\mathrm{,}\hfill & \alpha >\mathrm{0,}\hfill \\ \left\{0\right\}\mathrm{,}\hfill & \alpha \le \mathrm{0,}\hfill \end{array}$ (30)

where

$0<d_{\alpha }^{-}\left(t\mathrm{,}\omega \right)=\frac{1}{{\left(2\left|{\displaystyle {\int }_0^{\infty }}\exp \left(2\alpha s+2\sigma B_a\left(s\right)\right)\text{d}s\right|\right)}^{\frac{1}{2}}}<\infty \mathrm{.}$

The linearized stochastic dynamical system ${\vartheta }_t=D{\phi }_{\alpha }\left(t\mathrm{,}\omega \mathrm{,}\eta \right)\vartheta$ satisfies

$\text{d}{\vartheta }_t=\left(\alpha -3{\left({\phi }_{\alpha }\left(t\mathrm{,}\omega \mathrm{,}\eta \right)\right)}^2\right){\vartheta }_t\text{d}t+\sigma {\vartheta }_t\circ \text{d}{B}_a\left(t\right)\mathrm{,}$ (31)

then

$D{\phi }_{\alpha }\left(t\mathrm{,}\omega \mathrm{,}\eta \right)\vartheta =\vartheta \exp \left(\alpha t+\sigma B_a\left(t\right)-3{\displaystyle {\int }_0^t}{\left({\phi }_{\alpha }\left(s\mathrm{,}\omega \mathrm{,}\eta \right)\right)}^2\text{d}s\right)\mathrm{.}$

If ${\mu }_{i,\omega }^{\alpha }={\delta }_{{\eta }_i\left(\omega \right)}\left(i=1,2,3\right)$ is a $\phi$ -invariant measure, its Lyapunov exponent is

$\begin{array}{c}\lambda \left({\mu }_{i\mathrm{,}\omega }^{\alpha }\right)=\underset{t\to \infty }{lim}\frac{1}{t}\ln \Vert D{\phi }_{\alpha }\left(t\mathrm{,}\omega \mathrm{,}\eta \right)\vartheta \Vert \\ =\alpha -3\underset{t\to \infty }{lim}\frac{1}{t}{\displaystyle {\int }_0^t}{\left({\phi }_{\alpha }\left(s\mathrm{,}\omega \mathrm{,}\eta \right)\right)}^2\text{d}s\\ =\alpha -3\mathbb{E}{\eta }_i^2\mathrm{,}\end{array}$

provided that the IC ${\eta }_i^2\in {\text{L}}^1\left(\mathbb{P}\right)$ holds. Thus, we obtain the following results.

(i) For $\alpha \in ℝ$ , the IC for ${\mu }_{\mathrm{1,}\omega }^{\alpha }={\delta }_0$ is trivially satisfied and we have

$\lambda \left({\mu }_{\mathrm{1,}\omega }^{\alpha }\right)=\alpha ={\nu }_1+\frac{{\nu }_2}{16}-\frac{{\nu }_4}{16}\mathrm{.}$

It is easy to notice that ${\mu }_{\mathrm{1,}\omega }^{\alpha }$ is stable for $\alpha <0$ and unstable for $\alpha >0$ .

(ii) For $\alpha >0$ , ${\mu }_{\mathrm{2,}\omega }^{\alpha }={\delta }_{-d_{\alpha }^{-}\left(\omega \right)}$ is ${ℱ}_{-\infty }^0$ measurable. Therefore, the probability density function $p_2^{\alpha }\left(\eta \right)$ of ${\mu }_{\mathrm{2,}\omega }^{\alpha }$ satisfies the following Fokker-Planck equation:

$\frac{\partial p_2^{\alpha }\left(\eta \right)}{\partial t}=-\frac{\partial }{\partial \eta }\left[\left(\left(\alpha +\frac{{\sigma }^2}{2}\right)\eta -{\eta }^3\right)p_2^{\alpha }\left(\eta \right)\right]+\frac{1}{2}\frac{{\partial }^2}{\partial {\eta }^2}\left({\sigma }^2{\eta }^2{p}_2^{\alpha }\left(\eta \right)\right)\mathrm{,}$

which has the unique probability density solution

$p_2^{\alpha }\left(\eta \right)=\{\begin{array}{ll}{C}_{\alpha }{\eta }^{\frac{2\alpha }{{\sigma }^2}-1}\exp \left(-\frac{{\eta }^2}{{\sigma }^2}\right)\mathrm{,}\hfill & \eta >\mathrm{0,}\hfill \\ \mathrm{0,}\hfill & \eta \le \mathrm{0,}\hfill \end{array}$

with the normalizing parameter $C_{\alpha }$ satisfying $C_{\alpha }^{-1}={\sigma }^{\frac{2\alpha }{{\sigma }^2}}\Gamma \left(\frac{\alpha }{{\sigma }^2}\right)$ . Since

${\mathbb{E}}_{{\mu }_{\mathrm{2,}\omega }^{\alpha }}{\eta }^2=\mathbb{E}{\left(d_{\alpha }^{-}\right)}^2={\displaystyle {\int }_0^{\infty }}{\eta }^2{p}_2^{\alpha }\left(\eta \right)\text{d}\eta <\infty \mathrm{,}$

the IC is satisfied. Furthermore,

${\left(d_{\alpha }^{-}\left(\vartheta \left(t\right)\omega \right)\right)}^2=\frac{\exp \left(2\alpha t+2\sigma B_a\left(t\right)\right)}{2{\displaystyle {\int }_{-\infty }^0}\exp \left(2\alpha s+2\sigma B_a\left(s\right)\right)\text{d}s}=\frac{{\psi }^{\prime }\left(t\right)}{2\psi \left(t\right)}\mathrm{,}$

where $\psi \left(t\right)={\displaystyle {\int }_{-\infty }^t}\exp \left(2\alpha s+2\sigma B_a\left(s\right)\right)\text{d}s$ . According to the ergodic theory, we derive

$\mathbb{E}{\left(d_{\alpha }^{-}\left(\omega \right)\right)}^2=\underset{t\to \infty }{lim}\frac{1}{t}{\displaystyle {\int }_0^t}{\left(d_{\alpha }^{-}\left(\vartheta \left(t\right)\omega \right)\right)}^2\text{d}s=\frac{1}{2}\underset{t\to \infty }{lim}\frac{1}{t}\ln \psi \left(t\right)=\alpha \mathrm{.}$

Finally, we obtain $\lambda \left({\mu }_{2,\omega }^{\alpha }\right)=-2\alpha <0$ , which implies the invariant measure ${\mu }_{\mathrm{2,}\omega }^{\alpha }$ is stable for $\alpha >0$ .

(iii) For $\alpha >0$ , the invariant measure ${\mu }_{3,\omega }^{\alpha }={\delta }_{d_{\alpha }^{-}\left(\omega \right)}$ is also ${ℱ}_{-\infty }^0$ measurable, and its probability density $p_3^{\alpha }\left(\eta \right)$ is equal to $p_2^{\alpha }\left(-\eta \right)$ . In addition, we have

$\mathbb{E}{\left(-d_{\alpha }^{-}\left(\omega \right)\right)}^2=\mathbb{E}{\left(d_{\alpha }^{-}\left(\omega \right)\right)}^2=\alpha$

and $\lambda \left({\mu }_{3,\omega }^{\alpha }\right)=-2\alpha <0$ . Thus, we can conclude that the invariant measure ${\mu }_{\mathrm{3,}\omega }^{\alpha }$ is stable for $\alpha >0$ .

Hence system (9) undergoes a D-bifurcation. Similar to the analysis of (18), we find that system (9) undergoes a P-bifurcation at $\alpha ={\nu }_4/16$ . In summary, system (9) exhibits a D-bifurcation at $\alpha =0$ and a P-bifurcation at $\alpha ={\nu }_4/16$ . Consequently, system (9) undergoes a stochastic pitchfork bifurcation as shown in Figure 1. This proof is complete.

Theorem 5.6. Let $\beta ={\delta }_1+{\nu }_2/16-{\nu }_4/16$ . If the following three conditions hold,

(i) for $\alpha \in ℝ$ , the IC for ${\mu }_{1,\omega }^{\beta }={\delta }_0$ is trivially satisfied and $\lambda \left({\mu }_1^{\beta }\right)=\beta$ ;

(ii) for $\alpha >0$ , ${\mu }_{2,\omega }^{\beta }={\delta }_{-d_{\beta }^{-}\left(\omega \right)}$ is ${ℱ}_{-\infty }^0$ measurable, $\lambda \left({\mu }_2^{\beta }\right)=-2\beta <0$ ;

(iii) for $\alpha >0$ , ${\mu }_{3,\omega }^{\beta }={\delta }_{d_{\beta }^{-}\left(\omega \right)}$ is ${ℱ}_{-\infty }^0$ measurable, $\lambda \left({\mu }_3^{\beta }\right)=-2\beta <0$ .

Then system (15) undergoes a stochastic pitchfork bifurcation.

Proof. Since the proof is similar to that of Theorem 5.5, we omit it.

6. Numerical Simulation

In this section, we employed MATLAB software for numerical simulation to illustrate the theoretical results.

Example 6.1 Let $\Lambda =1$ , $b=0.6$ , $d=0.52$ , $\lambda =0.02$ , $\rho =0.4$ , $K=0.935$ , ${\sigma }_1=0.1$ , ${\sigma }_2=0.1$ , then $I_k=K/\rho =2.3375$ . For system (5) and (10), we will analyze the stochastic stability at point $\left(\mathrm{0,0}\right)$ as follows.

Figure 1. The stochastic pitchfork bifurcation.

Case (I) When $0\le I\le I_k$ , we choose $r=0.1$ and the initial value $x_1\left(0\right)=0.5$ , $y_1\left(0\right)=0.7$ . When $\mu =0.35$ , then ${\nu }_4=0.07$ and $16{\nu }_1+{\nu }_2=-1.402$ , i.e., ${\nu }_4>16{\nu }_1+{\nu }_2$ . We can see that the origin $\left(\mathrm{0,0}\right)$ is stable, as shown in Figure 2(a). Similarly, when $\mu =0.4$ , we can derive that ${\nu }_4=0.07$ and $16{\nu }_1+{\nu }_2=0.673$ , i.e., ${\nu }_4<16{\nu }_1+{\nu }_2$ , then the origin $\left(\mathrm{0,0}\right)$ is unstable, as shown in Figure 2(b). Consequently, the results of Theorem 4.1 are verified.

Case (II) When $I>I_k$ , we choose $r=0.95$ and the initial value $x_2\left(0\right)=0.5$ , $y_2\left(0\right)=0.7$ . When $\mu =0.38$ , then ${\nu }_4=0.07$ and $16{\delta }_1+{\nu }_2=-0.76$ , i.e., ${\nu }_4>16{\delta }_1+{\nu }_2$ . Hence the origin $\left(\mathrm{0,0}\right)$ is stable and see Figure 3(a). When $\mu =0.6$ , we can obtain that ${\nu }_4=0.07$ and $16{\delta }_1+{\nu }_2=2.8$ , i.e., ${\nu }_4<16{\delta }_1+{\nu }_2$ . Then, the origin $\left(\mathrm{0,0}\right)$ is not stable, as shown in Figure 3(b). Thus, the results of Theorem 4.2 are verified.

Example 6.2. Let $\Lambda =1$ , $b=0.6$ , $d=0.52$ , $\lambda =0.02$ , $\rho =0.4$ , $K=0.935$ , ${\sigma }_1=0.18$ , ${\sigma }_2=0.2$ , then $I_k=K/\rho =2.3375$ . We will investigate the P-bifurcation of the stochastic system (9) and (15) by varying parameter μ, respectively.

Case (I) When $0\le I\le I_k$ and $r=0.1$ , varying parameter μ which satisfies the conditions of Theorem 5.1 enables us to observe the qualitative changes of density function $p\left(a\right)$ as follows.

(i) If $8{\nu }_1+{\nu }_2/2\le {\nu }_4<16{\nu }_1+{\nu }_2$ , then the probability density function $p\left(a\right)$ tends to infinity as $a\to 0^{+}$ . Furthermore, the solution trajectories of system (9) are concentrated in a neighborhood of point $a=0$ . The graph of $p\left(a\right)$ is shown in Figure 4(a).

(ii) If $\left(16{\nu }_1+{\nu }_2\right)/3\le {\nu }_4<8{\nu }_1+{\nu }_2/2$ , then the probability density function

Figure 2. Phase portrait of system (5). (a) $\mu =0.35$ ; (b) $\mu =0.4$ .

Figure 3. Phase portrait of system (10). (a) $\mu =0.38$ ; (b) $\mu =0.6$ .

Figure 4. Probability density function $p\left(a\right)$ of system (9). (a) $\mu =0.37$ ; (b) $\mu =0.378$ ; (c) $\mu =0.39$ .

$p\left(a\right)$ reaches a maximum value at point $a^{\mathrm{<mo>*</mo>}}\approx 1.86$ and a minimum value at point $a=0$ . The graph of $p\left(a\right)$ is shown in Figure 4(b). We can see the derivative of $p\left(a\right)$ does not exist at point $a=0$ .

(iii) If ${\nu }_4<\left(16{\nu }_1+{\nu }_2\right)/3$ , then the probability density function $p\left(a\right)$ reaches a maximum value at point $a^{\mathrm{<mo>*</mo>}}\approx 3.87$ and a minimum value at point $a=0$ , as shown in Figure 4(c). We can see the derivative of the probability density function $p\left(a\right)$ exists at point $a=0$ .

Case (II) When $I>I_k$ and $r=0.95$ , varying parameter μ which satisfies the conditions of Theorem 5.3, we can observe the qualitative changes of density function $\zeta \left(a\right)$ as follows.

(i) If $8{\delta }_1+{\nu }_2/2\le {\nu }_4<16{\delta }_1+{\nu }_2$ , then the probability density function $\zeta \left(a\right)$ tends to infinity as $a\to 0^{+}$ (see Figure 5(a)). This indicates that the solution trajectories of system (15) are concentrated in a neighborhood of point $a=0$ .

(ii) If $\left(16{\delta }_1+{\nu }_2\right)/3\le {\nu }_4<8{\delta }_1+{\nu }_2/2$ , then the probability density function $\zeta \left(a\right)$ reaches a maximum value at point $a^{\mathrm{<mo>*</mo>}}\approx 1.46$ and a minimum value at point $a=0$ , as shown in Figure 5(b). Moreover, we can see the derivative of $\zeta \left(a\right)$ does not exist at point $a=0$ .

(iii) If ${\nu }_4<\left(16{\delta }_1+{\nu }_2\right)/3$ , then the probability density function $\zeta \left(a\right)$ reaches a maximum value at point $a^{\mathrm{<mo>*</mo>}}\approx 3.79$ and a minimum value at point $a=0$ , as shown in Figure 5(c). Furthermore, the derivative of the probability density function $\zeta \left(a\right)$ exists at $a=0$ .

Figure 5. Probability density function $\zeta \left(a\right)$ of system (15). (a) $\mu =0.39$ ; (b) $\mu =0.4$ ; (c) $\mu =0.42$ .

Next, we still change parameter μ and observe the qualitative changes of the joint density function $\rho \left(x\mathrm{,}y\right)$ and $\tau \left(u\mathrm{,}v\right)$ .

Case (III) When $0\le I\le I_k$ and $r=0.1$ , varying parameter μ which satisfies the conditions of Theorem 5.2, we obtain the following results.

(i) If ${\nu }_4\ge \left(16{\nu }_1+{\nu }_2\right)/3$ , the joint probability density function $\rho \left(x\mathrm{,}y\right)$ tends to infinity as $x\to 0$ , $y\to 0$ . The graph of $\rho \left(x\mathrm{,}y\right)$ is shown in Figure 6.

(ii) If $\left(16{\nu }_1+{\nu }_2\right)/5\le {\nu }_4<\left(16{\nu }_1+{\nu }_2\right)/3$ , a maximum value arises at the points of cycle $x^2+y^2\approx 4.07$ , and a minimum value appears at point $\left(\mathrm{0,0}\right)$ . At the same time, we find that the partial derivatives of joint probability density $\rho \left(x\mathrm{,}y\right)$ do not exist at point $\left(\mathrm{0,0}\right)$ , as shown in Figure 7.

(iii) If ${\nu }_4<\left(16{\nu }_1+{\nu }_2\right)/5$ , a maximum value arises at the points of cycle $x^2+y^2\approx 13.23$ , and a minimum value appears at point $\left(\mathrm{0,0}\right)$ . Meanwhile, we find that the partial derivatives of joint probability density $\rho \left(x\mathrm{,}y\right)$ exist at point $\left(\mathrm{0,0}\right)$ , as shown in Figure 8.

Figure 6. In system (9), for $\mu =0.378$ , (a) joint density function $\rho \left(x\mathrm{,}y\right)$ , (b) projection of joint density function $\rho \left(x\mathrm{,}y\right)$ and (c) joint density function $\rho \left(x\mathrm{,0}\right)$ .

Figure 7. In system (9), for $\mu =0.385$ , (a) joint density function $\rho \left(x\mathrm{,}y\right)$ , (b) projection of joint density function $\rho \left(x\mathrm{,}y\right)$ and (c) joint density function $\rho \left(x\mathrm{,0}\right)$ .

Figure 8. In system (9), for $\mu =0.395$ , (a) joint density function $\rho \left(x\mathrm{,}y\right)$ , (b) projection of joint density function $\rho \left(x\mathrm{,}y\right)$ and (c) joint density function $\rho \left(x\mathrm{,0}\right)$ .

Case (IV) When $I>I_k$ and $r=0.95$ , varying parameter μ which satisfies the conditions of Theorem 5.4, we obtain the following conclusions.

(i) If ${\nu }_4\ge \left(16{\delta }_1+{\nu }_2\right)/3$ , the joint probability density function $\tau \left(u\mathrm{,}v\right)$ tends to infinity as $u\to 0$ , $v\to 0$ . The graph of $\tau \left(u\mathrm{,}v\right)$ is shown in Figure 9.

(ii) If $\left(16{\delta }_1+{\nu }_2\right)/5\le {\nu }_4<\left(16{\delta }_1+{\nu }_2\right)/3$ , a maximum value arises at the points of cycle $x^2+y^2\approx 2.21$ , and a minimum value appears at $\left(\mathrm{0,0}\right)$ . Further, we find that the partial derivatives of joint probability density $\tau \left(u\mathrm{,}v\right)$ do not exist at point $\left(\mathrm{0,0}\right)$ , as shown in Figure 10.

(iii) If ${\nu }_4<\left(16{\delta }_1+{\nu }_2\right)/5$ , a maximum value arises at the points of cycle $x^2+y^2\approx 18.74$ , and a minimum value appears at the point $\left(\mathrm{0,0}\right)$ . At the same time, we find that the partial derivatives of joint probability density $\tau \left(u\mathrm{,}v\right)$ exist at point $\left(\mathrm{0,0}\right)$ , as shown in Figure 11.

Figure 9. In system (15), for $\mu =0.4$ , (a) joint density function $\tau \left(u\mathrm{,}v\right)$ , (b) projection of joint density function $\tau \left(u\mathrm{,}v\right)$ and (c) joint density function $\tau \left(u\mathrm{,0}\right)$ .

Figure 10. In system (15), for $\mu =0.41$ , (a) joint density function $\tau \left(u\mathrm{,}v\right)$ , (b) projection of joint density function $\tau \left(u\mathrm{,}v\right)$ and (c) joint density function $\tau \left(u\mathrm{,0}\right)$ .

Remark 6.1. From Definition 2.2, we can see that when model (1) and (2) are perturbed by environmental noise, the stationary probability density function $p\left(a\right)$ and $\zeta \left(a\right)$ change from a decreasing shape to a single peak as the disease-related death rate μ changes, indicating that model (3) and (4) undergo a P-bifurcation. The joint probability density function $\rho \left(x\mathrm{,}y\right)$ and $\tau \left(u\mathrm{,}v\right)$

Figure 11. In system (15), for $\mu =0.44$ , (a) joint density function $\tau \left(u\mathrm{,}v\right)$ , (b) projection of joint density function $\tau \left(u\mathrm{,}v\right)$ and (c) joint density function $\tau \left(u\mathrm{,0}\right)$ .

change from peak to crater-like, so that model (3) and (4) undergo a Hopf bifurcation which are shown in Definition 2.4.

Example 6.3. Let $\Lambda =1$ , $b=0.6$ , $d=0.52$ , $\lambda =0.02$ , $K=0.935$ , $\rho =0.4$ , then $I_k=2.3375$ . We will investigate the effect of noise on the phase portraits of the stochastic system (5) and (10), respectively.

Case (I) When $0\le I\le I_k$ , $\mu =0.395$ and ${\sigma }_1={\sigma }_2=0$ , then the stochastic system (5) degenerates into the deterministic system. Figure 12(a) shows the limit cycle of the deterministic system. We find that the limit cycle of the deterministic system disappears with increasing noise intensity, as shown in Figure 13.

Case (II) When $I>I_k$ and $\mu =0.43$ , similarly, we obtain Figure 12(b) and Figure 14 of the stochastic system (10).

Example 6.4 Let $\Lambda =1$ , $b=0.6$ , $d=0.52$ , $\lambda =0.02$ , ${\sigma }_1=0.18$ , ${\sigma }_2=0.2$ ,

Figure 12. Phase portrait of system (5) and (10) when ${\sigma }_1={\sigma }_2=0$ . (a) Phase portrait of system (5); (b) Phase portrait of system (10).

Figure 13. Phase portrait of system (5) as the noise intensity increases. (a) ${\sigma }_1={\sigma }_2=0.03$ ; (b) ${\sigma }_1=0.1$ ,${\sigma }_2=0.2$ ; (c) ${\sigma }_1={\sigma }_2=0.2$ .

Figure 14. Phase portrait of system (10) as the noise intensity increases. (a) ${\sigma }_1={\sigma }_2=0.03$ ; (b) ${\sigma }_1=0.1$ ; (c) ${\sigma }_1={\sigma }_2=0.2$ .

Figure 15. The sample paths of system (5). (a) $\mu =0.378$ ; (b) $\mu =0.385$ ; (c) $\mu =0.395$ .

$K=0.935$ , $\rho =0.4$ , then $I_k=2.3375$ . We obtain the sample paths of system (5) and (10), respectively.

Case (I) When $0\le I\le I_k$ and $r=0.1$ , by changing the parameter μ, then we obtain the sample paths of $x_1\left(t\right)$ and $y_1\left(t\right)$ for the stochastic system (5) shown in Figure 15.

Case (II) When $I>I_k$ and $r=0.95$ , by changing the parameter μ, then we obtain the sample paths of $x_2\left(t\right)$ and $y_2\left(t\right)$ for the stochastic system (10) shown in Figure 16.

Figure 16. The sample paths of system (10). (a) $\mu =0.4$ ; (b) $\mu =0.41$ ; (c) $\mu =0.44$ .

7. Conclusions

In this paper, we devote our main attention to studying both the stochastic model (3) and model (4) which include treatment and immigration. The disturbance of random noise was found to play an important role. By utilizing the stochastic averaging method, we simplified the stochastic model (5) and model (10). Based on Oseledec’s multiplicative ergodic theorem and singular boundary theory, we investigate the stochastic stability of the positive equilibrium point E₁ and E₂. It can be seen that varying the parameter μ can affect the stochastic stability of the positive equilibrium point, which may be useful for the subsequent study of disease transmission. In Theorem 5.2 and Theorem 5.4, the stochastic Hopf bifurcation of system (3) and system (4) are analyzed. Furthermore, in Theorem 5.5 and Theorem 5.6, the stochastic pitchfork bifurcation of system (3) and system (4) are investigated, resulting in a complete pitchfork bifurcation diagram. In the numerical simulation, we analyze the Hopf bifurcation of the system by changing the parameter μ. This indicates that we can avoid the occurrence of stochastic bifurcation by adjusting disease-induced mortality μ, offering a new approach for formulating infectious disease prevention and control strategies. In addition, we observe that as the noise intensity gradually increases from zero, the limit cycle gradually breaks and disappears.

Time delays are common in ecological systems, for instance, the transition of susceptible individuals to infected ones takes a certain period. To better reflect real-world, one can consider introducing time delay into the stochastic SIS model and then study its bifurcations in future work.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Ma, Z.E. and Zhou, Y.C. (2009) Global Stability of a Class of Discrete Age-Structured SIS Models with Immigration. Mathematical Biosciences & Engineering, 6, 409-425.
[2]	Brauer, F. and van den Driessche, P. (2001) Models for Transmission of Disease with Immigration of Infectives. Mathematical Biosciences, 171, 143-154. https://doi.org/10.1016/s0025-5564(01)00057-8
[3]	Cui, Q., Du, Q. and Wang, L. (2020) Global Dynamics of a Generalized SIRS Epidemic Model with Constant Immigration. Mathematical Problems in Engineering, 2020, Article ID: 7845390. https://doi.org/10.1155/2020/7845390
[4]	Buonomo, B. and Lacitignola, D. (2012) Forces of Infection Allowing for Backward Bifurcation in an Epidemic Model with Vaccination and Treatment. Acta Applicandae Mathematicae, 122, 283-293. https://doi.org/10.1007/s10440-012-9743-x
[5]	Wang, W. (2006) Backward Bifurcation of an Epidemic Model with Treatment. Mathematical Biosciences, 201, 58-71. https://doi.org/10.1016/j.mbs.2005.12.022
[6]	Wang, Z. (2008) Backward Bifurcation in Simple SIS Model. Acta Mathematicae Applicatae Sinica, English Series, 25, 127-136. https://doi.org/10.1007/s10255-006-6160-9
[7]	Saikh, A. and Gazi, N.H. (2021) The Effect of the Force of Infection and Treatment on the Disease Dynamics of an SIS Epidemic Model with Immigrants. Results in Control and Optimization, 2, Article ID: 100007. https://doi.org/10.1016/j.rico.2021.100007
[8]	Liu, W., Hethcote, H.W. and Levin, S.A. (1987) Dynamical Behavior of Epidemiological Models with Nonlinear Incidence Rates. Journal of Mathematical Biology, 25, 359-380. https://doi.org/10.1007/bf00277162
[9]	Liu, W., Levin, S.A. and Iwasa, Y. (1986) Influence of Nonlinear Incidence Rates Upon the Behavior of SIRS Epidemiological Models. Journal of Mathematical Biology, 23, 187-204. https://doi.org/10.1007/bf00276956
[10]	Liu, X. (2011) Bifurcation of an Eco-Epidemiological Model with a Nonlinear Incidence Rate. Applied Mathematics and Computation, 218, 2300-2309. https://doi.org/10.1016/j.amc.2011.07.050
[11]	Gray, A., Greenhalgh, D., Hu, L., Mao, X. and Pan, J. (2011) A Stochastic Differential Equation SIS Epidemic Model. SIAM Journal on Applied Mathematics, 71, 876-902. https://doi.org/10.1137/10081856x
[12]	Li, D., Wei, F. and Mao, X. (2022) Stationary Distribution and Density Function of a Stochastic SVIR Epidemic Model. Journal of the Franklin Institute, 359, 9422-9449. https://doi.org/10.1016/j.jfranklin.2022.09.026
[13]	Huang, C., Zhang, H., Cao, J. and Hu, H. (2019) Stability and Hopf Bifurcation of a Delayed Prey-Predator Model with Disease in the Predator. International Journal of Bifurcation and Chaos, 29, Article ID: 1950091. https://doi.org/10.1142/s0218127419500913
[14]	Gard, T.C. (1986) Stability for Multispecies Population Models in Random Environments. Nonlinear Analysis: Theory, Methods & Applications, 10, 1411-1419. https://doi.org/10.1016/0362-546x(86)90111-2
[15]	Farnoosh, R. and Parsamanesh, M. (2017) Stochastic Differential Equation Systems for an SIS Epidemic Model with Vaccination and Immigration. Communications in Statistics—Theory and Methods, 46, 8723-8736. https://doi.org/10.1080/03610926.2016.1189571
[16]	Luo, C. and Guo, S. (2016) Stability and Bifurcation of Two-Dimensional Stochastic Differential Equations with Multiplicative Excitations. Bulletin of the Malaysian Mathematical Sciences Society, 40, 795-817. https://doi.org/10.1007/s40840-016-0313-7
[17]	Huang, Z., Yang, Q. and Cao, J. (2011) Stochastic Stability and Bifurcation for the Chronic State in Marchuk’s Model with Noise. Applied Mathematical Modelling, 35, 5842-5855. https://doi.org/10.1016/j.apm.2011.05.027
[18]	Huang, D., Wang, H., Feng, J. and Zhu, Z. (2006) Hopf Bifurcation of the Stochastic Model on HAB Nonlinear Stochastic Dynamics. Chaos, Solitons & Fractals, 27, 1072-1079. https://doi.org/10.1016/j.chaos.2005.04.086
[19]	Xu, C. (2020) Phenomenological Bifurcation in a Stochastic Logistic Model with Correlated Colored Noises. Applied Mathematics Letters, 101, Article ID: 106064. https://doi.org/10.1016/j.aml.2019.106064
[20]	Zhao, D. and Yuan, S. (2020) Noise-Induced Bifurcations in the Stochastic Chemostat Model with General Nutrient Uptake Functions. Applied Mathematics Letters, 103, Article ID: 106180. https://doi.org/10.1016/j.aml.2019.106180
[21]	Zhang, X. and Yuan, R. (2022) Stochastic Bifurcation and Density Function Analysis of a Stochastic Logistic Equation with Distributed Delay and Weak Kernel. Mathematics and Computers in Simulation, 195, 56-70. https://doi.org/10.1016/j.matcom.2021.12.023
[22]	Arnold, L. (1998) Random Dynamical Systems. Springer-Verlag.
[23]	Fatehi Nia, M. and Hossein Akrami, M. (2019) Stability and Bifurcation in a Stochastic Vocal Folds Model. Communications in Nonlinear Science and Numerical Simulation, 79, Article ID: 104898. https://doi.org/10.1016/j.cnsns.2019.104898
[24]	Zhu, W. (2003) Nonlinear Stochastic Dynamics and Control in Hamiltonian Formulation. Science Press.
[25]	Oseledets, V.I. (1968) A Multiplicative Ergodic Theorem. Lyapunov Characteristic Numbers for Dynamical Systems. Transactions of the Moscow Mathematical Society, 19, 197-231.
[26]	Sri Namachchivaya, N. (1990) Stochastic Bifurcation. Applied Mathematics and Computation, 39, 37s-95s. https://doi.org/10.1016/0096-3003(90)90003-l
[27]	Wagner, U. V. and Wedig, W.V. (2000) On the Calculation of Stationary Solutions of Multi-Dimensional Fokker-Planck Equations by Orthogonal Functions. Nonlinear Dynamics, 21, 289-306. https://doi.org/10.1023/a:1008389909132