Threshold dynamics of a stochastic mathematical model for Wolbachia infections

Abstract

A stochastic mathematical model is proposed to study how environmental heterogeneity and the augmentation of mosquitoes with $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$ bacteria affect the outcomes of dengue disease. The existence and uniqueness of the positive solutions of the system are studied. Then the V-geometrically ergodicity and stochastic ultimate boundedness are investigated. Further, threshold conditions for successful population replacement are derived and the existence of a unique ergodic steady-state distribution of the system is explored. The results show that the ratio of infected to uninfected mosquitoes has a great influence on population replacement. Moreover, environmental noise plays a significant role in control of dengue fever.

Keywords:

• Dengue fever
• population replacement
• ergodicity
• Lyapunov function

1. Introduction

Every year many people are infected with mosquito-borne diseases including dengue fever and dengue haemorrhagic fever, caused by dengue viruses. Dengue infections are one of the main reasons for illness in the tropics and subtropics [21]. Due to the lack of licensed vaccines or drugs, the most effective way to treat dengue is to control its mosquito vectors, but traditional insecticide spraying not only pollutes the environment but also causes insecticide resistance [14]. Many studies have demonstrated that releases of mosquitoes carrying the endosymbiotic bacteria Wolbachia provide a novel method to control dengue [4, 19]. This succeeds because cytoplasmic incompatibility (CI) is induced in the mosquitoes with Wolbachia, which causes early embryonic death when uninfected females mate with infected males, but this does not affect infected females, resulting in widespread distribution of the bacteria in nature [6, 16].

Based on the mechanisms of CI, two strategies including population suppression and population replacement have been proposed to control mosquitoes [5], with population suppression being realized by releasing $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$-infected males inundatively and population replacement being achieved when $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$-infected mosquitoes are released inoculatively [22, 34]. Recently, researchers in many laboratories around the world are trying to control dengue virus by releasing $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$-infected mosquitoes with some success [17, 33]. Many different types of mathematical models have been proposed to study the transmission dynamics between wild and $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$-infected mosquitoes, including discrete time models [11, 13], continuous time models [10, 20] and impulsive differential equation models [37, 38]. For example, Haygood and Turelli established a discrete model to analyse the impact of host population subdivision on the evolution of CI-causing bacteria strains in specific host species, its results showed that in the subdivided host population with local density regulation, the strain evolution tended to be stronger CI [13]. Farkas and Hinow introduced and studied some differential equation models of population dynamics of $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$ infection, and they found that under the condition of mutual compatibility, strains with higher transmission efficiency or lower infection mortality are superior to competitors [10]. Zhang et al. proposed the birth-pulse model of $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$ transmission through mosquito population and reached the conclusion that population eradication can be achieved only when the parameters lie in a specific regions and the initial density of non-infection is low enough, regardless of the infection ratio [37].

However, wild mosquitoes and the released $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$-infected mosquitoes in nature are inevitably influenced by environmental fluctuations including temperature, wind, rainfall, oxygen and so on [18, 32]. These environmental factors affect the breeding, growth and development of mosquitoes: (1) the larval indices in the wet season is greater than that in the dry season [7, 31], (2) higher temperature will hasten mosquitoes development [8]. Relevant research indicates that the synergistic effect of between temperature and precipitation may have a significant impact on mosquitoes ecology and mosquito-borne diseases [1]. Indeed, the population dynamics of dengue vector mosquitoes are strongly linked with temperature and rainfall fluctuations [36]. The significant changes that this fluctuation causes to the mosquito population, stochastic differential equations with white noises provide a more realistic description.

Therefore, in this paper, we propose a mathematical model comprising a system of stochastic differential equations, governing the evolution of mosquitoes with white noises, then derive threshold conditions for population replacement and study the ergodic steady-state distribution of the system.

2. Model formation and preliminaries

Throughout this paper, the total population of mosquitoes is denoted by $N\left(t\right)$ and we assume that the mosquito population is infected by a single strain of $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$. The transmission can only be passed from infected females to their offspring, but transmission is imperfect with probability $\tau \in (0,1]$. $N\left(t\right)$ can be subdivided into four subpopulations, namely uninfected females, $F_U$, infected females, $F_I$, uninfected males, $M_U$, and infected males, $M_I$. b and d are the density dependent birth rate and death rate for the mosquito population, respectively. f is the proportion of females in the offspring. When an infected male mates with an uninfected female, zygotic death of offspring caused by CI usually occurs with a probability $q\in [0,1]$. Although CI has a beneficial effect on infected females, a fitness cost effect D of mosquitoes infected with $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$ is assumed to be non-zero with D>0 being fitness cost and D<0 being fitness benefit, and the sign of D depends on the mosquito species and $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$ strains. Then the model with overlapping generations can be described as follows: (1) $\{\begin{array}{l}{\displaystyle \frac{d{F}_I\left(t\right)}{dt}=f\tau b{F}_I-(d+D)N{F}_I,}\\ {\displaystyle \frac{d{F}_U\left(t\right)}{dt}=fb(1-\tau )F_I+fb{F}_U(1-\frac{q{M}_I}{M_U+M_I})-dN{F}_U,}\\ {\displaystyle \frac{d{M}_I\left(t\right)}{dt}=(1-f)\tau b{F}_I-(d+D)N{M}_I,}\\ {\displaystyle \frac{d{M}_U\left(t\right)}{dt}=(1-f)b(1-\tau )F_I+(1-f)b{F}_U(1-\frac{q{M}_I}{M_U+M_I})-dN{M}_U.}\end{array}$(1) In reality, the proportion of infected males that mate with infected females is usually the same as in uninfected populations, i.e. we have $M_I/F_I=M_U/F_U$, so after one or two generations, the ratio of males to females is identical in both cases. However, modifications may be considered if the sex ratios change [15]. Therefore, the entire infected ($I\left(t\right)$) and uninfected ($U\left(t\right)$) populations are introduced to simplify system (1(1) $\{\begin{array}{l}{\displaystyle \frac{d{F}_I\left(t\right)}{dt}=f\tau b{F}_I-(d+D)N{F}_I,}\\ {\displaystyle \frac{d{F}_U\left(t\right)}{dt}=fb(1-\tau )F_I+fb{F}_U(1-\frac{q{M}_I}{M_U+M_I})-dN{F}_U,}\\ {\displaystyle \frac{d{M}_I\left(t\right)}{dt}=(1-f)\tau b{F}_I-(d+D)N{M}_I,}\\ {\displaystyle \frac{d{M}_U\left(t\right)}{dt}=(1-f)b(1-\tau )F_I+(1-f)b{F}_U(1-\frac{q{M}_I}{M_U+M_I})-dN{M}_U.}\end{array}$(1) ) [10, 20], with suitable parameters the model is rescaled to (2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=\tau bI-(d+D)(I+U)I,}\\ {\displaystyle \frac{\mathrm{d}U}{\mathrm{d}t}=(1-\tau )bI+bU(1-\frac{qI}{U+I})-d(I+U)U.}\end{array}$(2) Many studies have focused on the model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=\tau bI-(d+D)(I+U)I,}\\ {\displaystyle \frac{\mathrm{d}U}{\mathrm{d}t}=(1-\tau )bI+bU(1-\frac{qI}{U+I})-d(I+U)U.}\end{array}$(2) ) [10, 20, 39], not only providing threshold conditions for the existence and stability of all possible equilibria but also discussing the biological significance regarding mosquito population replacement. But the effects of environmental heterogeneity on the dynamics of $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$ spread have been ignored. All organisms in natural habitats are constrained by fluctuations of many environmental factors such as temperature, nutrition, oxygen, pH and so on [32], and mosquitoes are no exception. Hu and co-authors constructed a mathematical model to study how random switches in birth rates affect the dynamics of $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$ spread [18]. Hence, we introduce white noise to study the influences of stochastic perturbations and assume that they are directly proportional to the entire infected $I\left(t\right)$ and uninfected $U\left(t\right)$ populations. Some authors (see [23, 24]) have pointed out that this assumption is reasonable and well justified biologically. Because the specific $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$ can not only be successfully transmitted in mosquito populations but also act like a vaccine to stop mosquitoes from replicating and spreading dengue virus, augmentation of mosquitoes with $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$ bacteria has been used to realize the aims of population replacement, and to prevent the occurrence of diseases, such as dengue disease [17, 34]. Therefore, we proposed a stochastic system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=\tau bI-(d+D)(I+U)I,}\\ {\displaystyle \frac{\mathrm{d}U}{\mathrm{d}t}=(1-\tau )bI+bU(1-\frac{qI}{U+I})-d(I+U)U.}\end{array}$(2) ) with control to model the effects of environmental fluctuations and the augmentation of mosquitoes, and system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=\tau bI-(d+D)(I+U)I,}\\ {\displaystyle \frac{\mathrm{d}U}{\mathrm{d}t}=(1-\tau )bI+bU(1-\frac{qI}{U+I})-d(I+U)U.}\end{array}$(2) ) is modified as (3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) where θ is the quantity of mosquitoes infected with $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$ being continuously released. $B_1\left(t\right)$ and $B_2\left(t\right)$ denote independent Brownian motion defined in a complete probability space $(\mathrm{\Omega },\mathcal{F},\{{\mathcal{F}}_t{\}}_{t\ge 0},\mathcal{P})$, ${\alpha }_1^2$ and ${\alpha }_2^2$ are the intensities of the noise on the entire infected $I\left(t\right)$ and uninfected $U\left(t\right)$ populations, respectively.

Throughout the paper, $(\mathrm{\Omega },\mathcal{F},\{{\mathcal{F}}_t{\}}_{t\ge 0},\mathcal{P})$ is denoted as a complete probability space with filtration $\{{\mathcal{F}}_t{\}}_{t\ge 0}$ and satisfies: (a) right continuous and (b) $\left\{{\mathcal{F}}_0\right\}$ involves all $\mathcal{P}$-null sets. The independent Brownian motion $B_i\left(t\right)$ is defined on this probability space. Assume that $X\left(t_0\right)=X_0(0\le t_0<T<\mathrm{\infty })$ is an $\left\{{\mathcal{F}}_0\right\}$-measurable ${\mathbf{R}}^2$-valued random variable, where ${\mathbf{R}}_{+}^2=\{x\in {\mathbf{R}}^2:x_i>0\text{for any}1\le i\le 2\}$. Define functions $f:{\mathbf{R}}^2\times [t_0,T]\to {\mathbf{R}}^2$ and $g:{\mathbf{R}}^2\times [t_0,T]\to {\mathbf{R}}^{2\times 2}$ such that they are Borel measurable. Consider the following Itô-type stochastic differential equation (4) $dX\left(t\right)=f\left(X\right(t),t)\mathrm{d}t+g\left(X\right(t),t)dB\left(t\right),X\left(0\right)=X_0,$(4) and the equivalent system of (4(4) $dX\left(t\right)=f\left(X\right(t),t)\mathrm{d}t+g\left(X\right(t),t)dB\left(t\right),X\left(0\right)=X_0,$(4) ) is (5) $X\left(t\right)=X_0+{\int }_{t_0}^{t}f\left(X\right(s),s)\mathrm{d}s+{\int }_{t_0}^{t}g\left(X\right(s),s)dB\left(s\right)\text{on}{t}_0\le t\le T.$(5) Let $p_t(X_0,A)$ be the transition probability, and $p_t(X_0,A)=\mathcal{P}\left(X\right(t)\in A\mid X(0)=X_0)$ for any $t\in {\mathbf{R}}_{+}$, $X_0\in {\mathbf{R}}_{+}^2$ and $A\in \mathcal{B}\left({\mathbf{R}}_{+}^2\right)$.

Then we give some important definitions as follows [12, 25, 26, 35].

Definition 2.1

[26]

Let $X\left(t\right)=\left(I\right(t),U(t){)}^T$ be a solution of (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ) provided

1. {X(t)} is continuous and $\left\{{\mathcal{F}}_t\right\}$-adapted;

2. $f\left(X\right(t),t)\in {\mathcal{L}}^1\left(\right[t_0,T];{\mathbf{R}}^2)$ and $g\left(X\right(t),t)\in {\mathcal{L}}^2\left(\right[t_0,T];{\mathbf{R}}^{2\times 2})$;

3. for any $t\in [t_0,T]$ (4(4) $dX\left(t\right)=f\left(X\right(t),t)\mathrm{d}t+g\left(X\right(t),t)dB\left(t\right),X\left(0\right)=X_0,$(4) ) holds with probability 1.

Definition 2.2

[12, 25]

Let $X\left(t\right)=\left(I\right(t),U(t){)}^T$ be a solution of SDE (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ):

1. if $\underset{t\to +\mathrm{\infty }}{lim}U\left(t\right)=0$, then $U\left(t\right)$ becomes extinctive;

2. if $\underset{t\to +\mathrm{\infty }}{lim sup}U\left(t\right)>0$, then $U\left(t\right)$ becomes weakly persistent;

3. if for any $ϵ\in (0,1)$, there are two constants $\text{β}>0$ and $\delta >0$ such that

$\underset{t\to +\mathrm{\infty }}{lim inf}\mathcal{P}\left\{U\right(t)\ge \text{β}\}\ge 1-ϵ,\underset{t\to +\mathrm{\infty }}{lim inf}\mathcal{P}\left\{U\right(t)\le \delta \}\ge 1-ϵ,$then $U\left(t\right)$ is called stochastically persistent.

Lemma 2.1

[35]

Let $f\left(t\right)\in C(\mathrm{\Omega }\times [0,+\mathrm{\infty }),R_{+}),$ if there are constants ${\zeta }_0,t_1$ and $\zeta \ge 0$ such that $f\left(t\right)$ satisfies $\ln f\left(t\right)\le \zeta t-{\zeta }_0{\int }_0^{t}f\left(s\right)\mathrm{d}s+\sum _{i=1}^n{\text{β}}_i{B}_i\left(t\right)$ for any $t\ge t_1$, ${\text{β}}_i$ is also a constant, then $\underset{t\to +\mathrm{\infty }}{lim}sup\frac{1}{t}{\int }_0^{t}f\left(s\right)\mathrm{d}s\le \frac{\zeta }{{\zeta }_0}$; if there are constants ${\zeta }_0,t_1$ and $\zeta \ge 0$ such that $f\left(t\right)$ satisfies $\ln f\left(t\right)\ge \zeta t-{\zeta }_0{\int }_0^{t}f\left(s\right)\mathrm{d}s+\sum _{i=1}^n{\text{β}}_i{B}_i\left(t\right)$ for any $t\ge t_1$, then $\underset{t\to +\mathrm{\infty }}{lim}sup\frac{1}{t}{\int }_0^{t}f\left(s\right)\mathrm{d}s\ge \frac{\zeta }{{\zeta }_0}$.

3. Main results

3.1. Properties of the solutions

For any given initial conditions, a unique global solution of system (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ) exists if its coefficients satisfy the linear growth condition and local Lipschitz condition [26]. Now, we show that the solution of system (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ) is positive and global by using methods of Lyapunov analysis [9].

Theorem 3.1

Let $X\left(t\right)=\left(I\right(t),U(t\left)\right)$ be a solution of SDE (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ) with initial condition $\left(I\right(0),U(0\left)\right)\in {\mathbf{R}}_{+}^2$, then $X\left(t\right)$ is unique for any $t\ge 0$ and it further remains in ${\mathbf{R}}_{+}^2$ with probability 1, namely $\left(I\right(t),U(t\left)\right)\in {\mathbf{R}}_{+}^2$ for any $t\ge 0$ almost surely.

Proof.

Let $x=\ln I$ and $y=\ln U$, applying Itô's formula to system (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ) yields (6) $\{\begin{array}{ll}dx\left(t\right)& =(\tau b-(d+D\left)\right(e^x+e^y)+\frac{\theta }{e^x}-\frac{1}{2}{\alpha }_1^2)\mathrm{d}t+{\alpha }_{1}d{B}_1\left(t\right),\\ dy\left(t\right)& =\left(\right(1-\tau )b\frac{e^x}{e^y}+b(1-\frac{q{e}^x}{e^y+e^x})-d(e^x+e^y)-\frac{1}{2}{\alpha }_2^2)\mathrm{d}t\\ & +{\alpha }_{2}d{B}_2\left(t\right),\end{array}$(6) where initial conditions are $x_0=\ln I_0$ and $y_0=\ln U_0$. Note that the coefficients of system (6(6) $\{\begin{array}{ll}dx\left(t\right)& =(\tau b-(d+D\left)\right(e^x+e^y)+\frac{\theta }{e^x}-\frac{1}{2}{\alpha }_1^2)\mathrm{d}t+{\alpha }_{1}d{B}_1\left(t\right),\\ dy\left(t\right)& =\left(\right(1-\tau )b\frac{e^x}{e^y}+b(1-\frac{q{e}^x}{e^y+e^x})-d(e^x+e^y)-\frac{1}{2}{\alpha }_2^2)\mathrm{d}t\\ & +{\alpha }_{2}d{B}_2\left(t\right),\end{array}$(6) ) are locally Lipschitz continuous, there is a blow up time ${\tau }_e$ such that system (6(6) $\{\begin{array}{ll}dx\left(t\right)& =(\tau b-(d+D\left)\right(e^x+e^y)+\frac{\theta }{e^x}-\frac{1}{2}{\alpha }_1^2)\mathrm{d}t+{\alpha }_{1}d{B}_1\left(t\right),\\ dy\left(t\right)& =\left(\right(1-\tau )b\frac{e^x}{e^y}+b(1-\frac{q{e}^x}{e^y+e^x})-d(e^x+e^y)-\frac{1}{2}{\alpha }_2^2)\mathrm{d}t\\ & +{\alpha }_{2}d{B}_2\left(t\right),\end{array}$(6) ) exists with a unique local solution $\left(x\right(t),y(t\left)\right)$ on $[0,{\tau }_e)$. It follows from Itô's formula that $\left(I\right(t),U(t\left)\right)=(e^{x\left(t\right)},e^{y\left(t\right)})$ is just the unique local solution of system (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ) on $[0,{\tau }_e)$. To show the solution $\left(I\right(t),U(t\left)\right)$ is global we only need to prove ${\tau }_e=\mathrm{\infty }$. Let $n_0>0$ big enough for $I_0$ and $U_0$ lying within the interval $[\frac{1}{n_0},n_0]$. For any integer $n\ge n_0$, define the stopping time as ${\tau }_n=inf\{t\in [0,{\tau }_e):min\{I\left(t\right),U\left(t\right)\}\le \frac{1}{n}\text{or}max\{I\left(t\right),U\left(t\right)\}\ge n\}.$Let $\mathrm{\varnothing }=$ be empty set, and set $inf\mathrm{\varnothing }=\mathrm{\infty }$. ${\tau }_n$ increases as $n\to \mathrm{\infty }$. Let ${\tau }_{\mathrm{\infty }}=lim\underset{n\to \mathrm{\infty }}{sup}{\tau }_n$, clearly, ${\tau }_{\mathrm{\infty }}\le {\tau }_e$. For any $t\ge 0$, if ${\tau }_{\mathrm{\infty }}=\mathrm{\infty }$ a.s., then we have ${\tau }_e=\mathrm{\infty }$ and $\left(I\right(t),U(t\left)\right)\in {\mathbf{R}}_{+}^2$. Therefore, what we only need to do is to show ${\tau }_{\mathrm{\infty }}=\mathrm{\infty }$. Otherwise, there exists with two constants T>0 and $ϵ\in (0,1)$ so that $\mathcal{P}\{{\tau }_{\mathrm{\infty }}\le T\}>ϵ$. Consequently, there is an integer $n_1\ge n_0$ such that (7) $\mathcal{P}\{{\tau }_n\le T\}>ϵ,n\ge n_1.$(7) Define a $C^2$-function $V(I,U):{\mathbf{R}}_{+}\times {\mathbf{R}}_{+}\to \mathbf{R}$ as $V(I,U)=(I-1-\ln I)+(U-1-\ln U),$and $V(I,U)$ is positive for $I\ge 0,U\ge 0$. Making use of Itô's formula, then $\begin{array}{rl}dV& =\{(1-\frac{1}{I})(\tau bI-(d+D\left)\right(I+U)I+\theta )+\frac{{\alpha }_1^2}{2}\}\mathrm{d}t+{\alpha }_1(I-1)d{B}_1\left(t\right)\\ & +\{(1-\frac{1}{U})\left(\right(1-\tau )bI+bU(1-\frac{qI}{U+I})-d(I+U\left)U\right)+\frac{{\alpha }_2^2}{2}\}\mathrm{d}t\\ & +{\alpha }_2(U-1)d{B}_2\left(t\right)\\ & \doteq LV\mathrm{d}t+{\alpha }_1(I-1)d{B}_1\left(t\right)+{\alpha }_2(U-1)d{B}_2\left(t\right),\end{array}$where $\begin{array}{rl}LV& =(1-\frac{1}{I})(\tau bI-(d+D\left)\right(I+U)I+\theta )\\ & +(1-\frac{1}{U})\left(\right(1-\tau )bI+bU(1-\frac{qI}{U+I})-d(I+U\left)U\right)+\frac{{\alpha }_1^2+{\alpha }_2^2}{2}\\ & \le -(d+D)I^2+(2d+D+b)I-d{U}^2+(2d+D+b)U\\ & -(d+D)UI-\frac{bqUI}{U+I}-dIU-\frac{(1-\tau )bI}{U}-\frac{\theta }{I}-b+\frac{{\alpha }_1^2+{\alpha }_2^2}{2}+bq+\theta \\ & \le \frac{(2d+D+b{)}^2}{4(d+D)}+\frac{(2d+D+b{)}^2}{4d}+\frac{{\alpha }_1^2+{\alpha }_2^2}{2}+bq+\theta \\ & \doteq C,\end{array}$where C is a positive constant and independent of I, U and t. Thus $dV(I,U)\le C\mathrm{d}t+{\alpha }_1(I-1)d{B}_1\left(t\right)+{\alpha }_2(U-1)d{B}_2\left(t\right).$Set ${\tau }_n\wedge T=min\{{\tau }_n,T\}$, then we obtain $\begin{array}{rl}{\int }_0^{{\tau }_n\wedge T}dV\left(I\right(t),U(t\left)\right)& \le {\int }_0^{{\tau }_n\wedge T}C\mathrm{d}t+{\int }_0^{{\tau }_n\wedge T}{\alpha }_1(I-1)d{B}_1\left(t\right)\\ & +{\int }_0^{{\tau }_n\wedge T}{\alpha }_2(U-1)d{B}_2\left(t\right),\end{array}$calculating the mathematical expectation of the above inequality yields (8) $\begin{array}{rl}EV\left(I\right({\tau }_n\wedge T),U({\tau }_n\wedge T\left)\right)& \le V\left(I\right(0),U(0\left)\right)+E{\int }_0^{{\tau }_n\wedge T}C\mathrm{d}t\\ & \le V\left(I\right(0),U(0\left)\right)+CT.\end{array}$(8) For any $n\le n_1$, let ${\mathrm{{\rm Y}}}_n={\tau }_n\le T$. From (7(7) $\mathcal{P}\{{\tau }_n\le T\}>ϵ,n\ge n_1.$(7) ) we have $\mathcal{P}\left({\mathrm{{\rm Y}}}_n\right)>ϵ$. For any $t^{\ast }\in {\mathrm{{\rm Y}}}_n$, there is at least one $I({\tau }_n,t^{\ast })$ or $U({\tau }_n,t^{\ast })$ equalling n or $\frac{1}{n}$, thereby, $V\left(I\right({\tau }_n,t^{\ast }),U({\tau }_n,t^{\ast }\left)\right)\ge (n-1-\ln n)\wedge (\frac{1}{n}-1-\ln \frac{1}{n}).$Combinations of (7(7) $\mathcal{P}\{{\tau }_n\le T\}>ϵ,n\ge n_1.$(7) ) and (8(8) $\begin{array}{rl}EV\left(I\right({\tau }_n\wedge T),U({\tau }_n\wedge T\left)\right)& \le V\left(I\right(0),U(0\left)\right)+E{\int }_0^{{\tau }_n\wedge T}C\mathrm{d}t\\ & \le V\left(I\right(0),U(0\left)\right)+CT.\end{array}$(8) ) lead to $\begin{array}{rl}V\left(I\right(0),U(0\left)\right)+CT& \ge E\left\{{\mathcal{L}}_{{\mathrm{{\rm Y}}}_n\left(t^{\ast }\right)}V\right(I\left({\tau }_n\right),U\left({\tau }_n\right)\left)\right\}\\ & \ge ϵ\left\{\right(n-1-\ln n)\wedge (\frac{1}{n}-1-\ln \frac{1}{n})\},\end{array}$with ${\mathcal{L}}_{{\mathrm{{\rm Y}}}_n\left(t^{\ast }\right)}$ denoting as the indicator function of ${\mathrm{{\rm Y}}}_n$. Once $n\to \mathrm{\infty }$, we obtain $\mathrm{\infty }=V\left(I\right(0),U(0\left)\right)+CT<\mathrm{\infty },$which is a contradiction. Therefore, ${\tau }_{\mathrm{\infty }}=\mathrm{\infty }$. This completes the proof.

Theorem 3.1 indicates that any solutions of SDE (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ) will finally remain in a compact set of ${\mathbf{R}}_{+}^2$. In the following, we prove that the solutions of SDE (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ) satisfy another important property, i.e. the Markov process $X\left(t\right)=\left(I\right(t),U(t\left)\right)$ is V-geometrically ergodic.

Theorem 3.2

Let initial condition $X_0\in {\mathbf{R}}_{+}^2$, if ${\alpha }_1>0$ and ${\alpha }_2>0$, then Markov process $X\left(t\right)=\left(I\right(t),U(t\left)\right)$ is V-geometrically ergodic.

Proof.

Note that N = I + U, define a function $V\left(X\right(t\left)\right)$ such that $V\left(X\right(t\left)\right)\to \mathrm{\infty }$ as $\mid X\left(t\right)\mid \to \mathrm{\infty }$ for $X\left(t\right)\in {\mathbf{R}}_{+}^2$, where (9) $V\left(X\right(t\left)\right)=N+\frac{1}{N}.$(9) Applying Itô's formula yields (10) $\begin{array}{rl}LV\left(X\right(t\left)\right)& =\tau bI-(d+D)NI+\theta +(1-\tau )bI+bU(1-\frac{qI}{N})-dNU\\ & -\frac{\tau bI-(d+D)NI+\theta +(1-\tau )bI+bU(1-\frac{qI}{N})-dNU}{N^2}\\ & +\frac{{\alpha }_1^2{I}^2+{\alpha }_2^2{U}^2}{N^3}\\ & \le -d{N}^2-DNI+\theta +bN-\frac{bqUI}{N}+\frac{(d+D)NI}{N^2}-\frac{\theta }{N^2}-\frac{b}{N}\\ & +\frac{bqUI}{N^3}+\frac{dU}{N}+\frac{{\alpha }_1^2{I}^2+{\alpha }_2^2{U}^2}{N^3}\\ & \le -bV\left(X\right)+(2d+D+\theta )-d{N}^2+2bN-\frac{\theta }{N^2}+\frac{bq+{\alpha }_1^2+{\alpha }_2^2}{N}\\ & \le 2d+D+\theta +\frac{b^2}{d}+\frac{(bq+{\alpha }_1^2+{\alpha }_2^2{)}^2}{4\theta }-bV\left(X\right)\\ & =C^{\ast }-bV\left(X\right),\end{array}$(10) where $C^{\ast }=2d+D+\theta +\frac{b^2}{d}+\frac{(bq+{\alpha }_1^2+{\alpha }_2^2{)}^2}{4\theta }.$Hence the Lyapunov condition holds, for details, see [28].

Note that when ${\alpha }_1>0$ and ${\alpha }_2>0$, then SDE (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ) is uniformly elliptic [3]. Athreya et al. pointed out that there is a jointly continuous function defined as $p:{\mathbf{R}}_{+}\times {\mathbf{R}}_{+}^2\times {\mathbf{R}}_{+}^2\to (0,\mathrm{\infty })$. For all $(t,X_0,Y)$, $p_t(X_0,Y)$ is strictly positive so that for all measure sets A we obtain $p_t(X_0,A)={\int }_A{p}_t(X_0,Y)dY.$For any w>0, there is a positive constant $c=c(w,t)>0$ such that $inf\left\{p_t\right(X_0,Y):X_0,Y\in {\mathbf{R}}_{+}^2,\mid X_0\mid ,\mid Y\mid \le w\}\ge c$. Thus for any measurable set A, we have $p_t(X,A)={\int }_A{p}_t(X_0,Y)dY\ge c\text{Leb}(A\cap {\mathcal{B}}_w(0\left)\right)=c\text{Leb}\left({\mathcal{B}}_w\right(0\left)\right)\upsilon \left(A\right),$where Leb is Lebesgue measure and $\upsilon \left(A\right)=\text{Leb}(A\cap {\mathcal{B}}_w(0\left)\right)/\text{Leb}\left({\mathcal{B}}_w\right(0\left)\right)$. So the Minirization condition holds. This completes the proof.

Theorem 3.3

For any initial value $X_0\in {\mathbf{R}}_{+}^2$, the solution $X\left(t\right)=\left(I\right(t),U(t\left)\right)$ of system (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ) is stochastically ultimately bounded and permanent.

Proof.

Clearly, N = I + U, define a Lyapunov function $V\left(t\right)=N+1/N$ and choose a small enough ξ such that $0<\xi \le b$. Making use of Itô's formula, it follows from (10(10) $\begin{array}{rl}LV\left(X\right(t\left)\right)& =\tau bI-(d+D)NI+\theta +(1-\tau )bI+bU(1-\frac{qI}{N})-dNU\\ & -\frac{\tau bI-(d+D)NI+\theta +(1-\tau )bI+bU(1-\frac{qI}{N})-dNU}{N^2}\\ & +\frac{{\alpha }_1^2{I}^2+{\alpha }_2^2{U}^2}{N^3}\\ & \le -d{N}^2-DNI+\theta +bN-\frac{bqUI}{N}+\frac{(d+D)NI}{N^2}-\frac{\theta }{N^2}-\frac{b}{N}\\ & +\frac{bqUI}{N^3}+\frac{dU}{N}+\frac{{\alpha }_1^2{I}^2+{\alpha }_2^2{U}^2}{N^3}\\ & \le -bV\left(X\right)+(2d+D+\theta )-d{N}^2+2bN-\frac{\theta }{N^2}+\frac{bq+{\alpha }_1^2+{\alpha }_2^2}{N}\\ & \le 2d+D+\theta +\frac{b^2}{d}+\frac{(bq+{\alpha }_1^2+{\alpha }_2^2{)}^2}{4\theta }-bV\left(X\right)\\ & =C^{\ast }-bV\left(X\right),\end{array}$(10) ) that we obtain $\begin{array}{rl}d{\exp }^{\xi t}V\left(t\right)& =\xi {\exp }^{\xi t}V\left(t\right)\mathrm{d}t+{\exp }^{\xi t}dV\left(t\right)=\xi {\exp }^{\xi t}V\left(t\right)\mathrm{d}t\\ & +{\exp }^{\xi t}\left\{\right(C^{\ast }-bV\left(t\right))\mathrm{d}t+(1-\frac{1}{N^2})({\alpha }_{1}Id{B}_1\left(t\right)+{\alpha }_{2}Ud{B}_2\left(t\right)\left)\right\}.\end{array}$Integrating the above equation from 0 to t and taking mathematical expectation, then $\begin{array}{rl}E[{\exp }^{\xi t}V(t\left)\right]& =E\left[V\right(0\left)\right]+E[{\int }_0^t{\exp }^{\xi s}(\xi dV\left(s\right)+LV\left(s\right)\left)\mathrm{d}s\right]\\ & \le E\left[V\right(0\left)\right]+C^{\ast }E[{\int }_0^t{\exp }^{\xi s}\mathrm{d}s]\\ & =E\left[V\right(0\left)\right]+\frac{C^{\ast }}{\xi }({\exp }^{\xi t}-1).\end{array}$Thus $\begin{array}{rl}E\left[V\right(t\left)\right]& \le {\exp }^{-\xi t}E\left[V\right(0\left)\right]+\frac{C^{\ast }}{\xi }(1-{\exp }^{-\xi t})\\ & \le E\left[V\right(0\left)\right]+\frac{C^{\ast }}{\xi }\doteq \mathrm{\Theta }.\end{array}$From Markov inequality, choose a positive constant ${\mathrm{\Theta }}^{\ast }$ large enough so that $\mathrm{\Theta }/{\mathrm{\Theta }}^{\ast }<1$, $\mathcal{P}\{N+\frac{1}{N}>{\mathrm{\Theta }}^{\ast }\}\le \frac{1}{{\mathrm{\Theta }}^{\ast }}E[N+\frac{1}{N}]\le \frac{\mathrm{\Theta }}{{\mathrm{\Theta }}^{\ast }}\doteq ϵ.$Consequently, $1-ϵ\le \mathcal{P}\{N+\frac{1}{N}\le {\mathrm{\Theta }}^{\ast }\}\le \mathcal{P}\{\frac{1}{{\mathrm{\Theta }}^{\ast }}\le N\le {\mathrm{\Theta }}^{\ast }\}.$In view of $N^2\le 3\mid X{\mid }^2\le 3N^2$, so $\mathcal{P}\{\frac{1}{\sqrt{3}{\mathrm{\Theta }}^{\ast }}\le \frac{N}{\sqrt{3}}\le \mid X\mid \le N\le {\mathrm{\Theta }}^{\ast }\}\ge 1-ϵ.$It follows from the definitions that system (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ) is stochastically ultimately bounded and permanent. This ends the proof.

3.2. Population replacement

One of the feasible measures to prevent dengue diseases from spreading is to achieve population replacement, by means of releasing mosquitoes with $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$. This section will focus on the conditions for population replacement. To this end, for simplicity we denote $g_{\ast }\left(t\right)=\underset{t\to +\mathrm{\infty }}{lim}inf\frac{1}{t}{\int }_0^{t}g\left(s\right)\mathrm{d}s\text{and}{g}^{\ast }\left(t\right)=\underset{t\to +\mathrm{\infty }}{lim}sup\frac{1}{t}{\int }_0^{t}g\left(s\right)\mathrm{d}s.$Note that $M_I/F_I=M_U/F_U$, then we have $F_I/F_U=M_I/M_U=k$, which means that $I/U=(M_I+F_I)/(M_U+F_U)=M_I/F_I=k$. Based on this fact, we have the following main results.

Theorem 3.4

Theorem 3.4 If $b-\frac{1}{2}{\alpha }_2^2+(1-\tau )bk-\frac{bqk}{1+k}<0$, then the total uninfected mosquitoes become extinct.

Proof.

According to system (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ), defining a Lyapunov function $V\left(t\right)=\ln U\left(t\right)$, applying Itô's formula to the second equation of system (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ) and integrating the above equation from 0 to t we obtain (11) $\frac{1}{t}\ln \frac{U\left(t\right)}{U\left(0\right)}=b-\frac{1}{2}{\alpha }_2^2+(1-\tau )bk-\frac{bqk}{1+k}-d\frac{1}{t}{\int }_0^t\left(I\right(s)+U(s\left)\right)\mathrm{d}s+\frac{M_2\left(t\right)}{t}.$(11) where $M_2\left(t\right)={\int }_0^t{\alpha }_{2}d{B}_2\left(t\right)$. Because of $<M_2\left(t\right),M_2\left(t\right)>={\int }_0^t{\delta }_2^2\mathrm{d}s$ and the strong law of large numbers for local martingales we obtain (12) $\underset{t\to +\mathrm{\infty }}{lim}\frac{M_2\left(t\right)}{t}=0.$(12) When $t\to +\mathrm{\infty }$, taking the superior limit of equation (11(11) $\frac{1}{t}\ln \frac{U\left(t\right)}{U\left(0\right)}=b-\frac{1}{2}{\alpha }_2^2+(1-\tau )bk-\frac{bqk}{1+k}-d\frac{1}{t}{\int }_0^t\left(I\right(s)+U(s\left)\right)\mathrm{d}s+\frac{M_2\left(t\right)}{t}.$(11) ) and by using L'Hospital's rule we have $\underset{t\to +\mathrm{\infty }}{lim}sup\frac{\ln U\left(t\right)}{t}\le b-\frac{1}{2}{\alpha }_2^2+(1-\tau )bk-\frac{bqk}{1+k}<0.$So $\underset{t\to +\mathrm{\infty }}{lim}U\left(t\right)=0$, which indicates that the total uninfected mosquitoes become extinct. This completes the proof.

Remark 3.1

Notice that the threshold conditions of Theorem 4 determine the outcomes of population replacement, it implies that the ratio of infected to uninfected mosquitoes and environmental noise play significant roles in control of dengue fever.

For simplicity, denote $R_0=b-\frac{1}{2}{\alpha }_2^2+(1-\tau )bk-\frac{bqk}{1+k}.$It can be seen from Figure 1(a) and (b) that the larger the random fluctuation, the earlier $U\left(t\right)$ tends to zero, i.e. the faster the uninfected mosquitoes die out. The larger θ, the faster $U\left(t\right)$ tends to zero, which means that increasing the number of continuously released $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$-infected mosquitoes will contribute to the extinction of uninfected mosquitoes.

Figure 1. Extinction of the uninfected mosquitoes. (a) $\theta =0.2$; (b) $\theta =1$. All other parameter values were fixed as: q = 0.2, $\tau =1$, b = 0.02, d = 0.02, D = 0.01, $\theta =0.2$ and initial value $\left(I\right(0),U(0\left)\right)=(0.5,0.5).$

Since $R_0$ is very sensitive to key parameters, it is critical to conduct sensitivity analysis. In Figure 2, two cases for ${\alpha }_2$ and k are discussed respectively. If fixing parameter values as shown in Figure 2(a), it can be found that $R_0$ increases with the increase of k. For ${\alpha }_2=0.5$, a small value of k can ensure that $R_0<0$, which indicates that uninfected mosquitoes are extinct. For ${\alpha }_2=0.01$, it can be observed that $R_0>0$, implying uninfected mosquitoes continue to exist. If fixing parameter values as shown in Figure 2(b), it is noted that $R_0$ decreases from $R_0>0$ to $R_0<0$ when ${\alpha }_2$ increases. Obviously, for k = 0.1, $R_0$ reaches $R_0<0$ faster. In general, the smaller k and the larger ${\alpha }_2$ will hasten to the extinction of uninfected mosquitoes.

Figure 2. These plots show that sensitivity of $k,{\alpha }_2$ on $R_0$. (a) We set ${\alpha }_2=0.5$ and ${\alpha }_2=0.01$ ; (b) We set k = 0.1 and k = 0.8, and all other parameter values were fixed as: $\tau =0.3$, b = 0.1, q = 0.2.

Through the above discussion, it is not difficult to find that increasing θ (i.e. the release of $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$-infected mosquitoes) can increase k (i.e. the ratio of $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$-infected mosquitoes and uninfected mosquitoes). Because a small k can ensure the rapid extinction of mosquitoes, it suggests that when using the ‘mosquito control’ strategy, it is necessary to choose an appropriate threshold for the release of $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$-infected mosquitoes.

3.3. Stationary distribution and ergodicity for the system

In this section, we explore the existence of a unique ergodic steady-state distribution of the system (6(6) $\{\begin{array}{ll}dx\left(t\right)& =(\tau b-(d+D\left)\right(e^x+e^y)+\frac{\theta }{e^x}-\frac{1}{2}{\alpha }_1^2)\mathrm{d}t+{\alpha }_{1}d{B}_1\left(t\right),\\ dy\left(t\right)& =\left(\right(1-\tau )b\frac{e^x}{e^y}+b(1-\frac{q{e}^x}{e^y+e^x})-d(e^x+e^y)-\frac{1}{2}{\alpha }_2^2)\mathrm{d}t\\ & +{\alpha }_{2}d{B}_2\left(t\right),\end{array}$(6) ). If f is a bounded function on ${\mathrm{R}}_{+}$, defining $f^u=\underset{t\in R_{+}}{sup}f\left(t\right)$, then we only need to show that the following two properties hold true [27, 40],

1. there exists a bounded domain $E\in \mathrm{I}\mathrm{n}\mathrm{t}{\mathrm{R}}_{+}^2$ with regular boundary Γ such that its closure $\overline{E}\subset \mathrm{I}\mathrm{n}\mathrm{t}{\mathrm{R}}_{+}^2$, and a non-negative $C^2-$ function $V\left(x\right)$ exists such that for any $x\in \mathrm{I}\mathrm{n}\mathrm{t}{\mathrm{R}}_{+}^2\setminus E$, LV is negative;

2. for any bounded domain $\hat{E}\in \mathrm{I}\mathrm{n}\mathrm{t}{\mathrm{R}}_{+}^2,$ there is a positive constant ζ such that the diffusion matrix for system (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ) given by $b\left(Z\right)=\left(\begin{array}{cc}{\alpha }_1^2{I}^2& 0\\ 0& {\alpha }_2^2{U}^2\end{array}\right)$satisfies $\sum _{i,j=1}^2{b}_{ij}\left(Z\right){\xi }_i{\xi }_j>\zeta \Vert \xi {\Vert }^2$ for all $Z=(I,U)\in \hat{E}$, and $\xi =$ $({\xi }_1,{\xi }_2)\in {\mathrm{R}}^2.$

Theorem 3.5

If (13) $\tau b-\frac{1}{2}{\alpha }_1^2>0,b-\frac{1}{2}{\alpha }_2^2>0\text{and}d+D>0,$(13) then system (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ) has a unique ergodic stationary distribution.

Proof.

Let $V(I,U)=I^p+U^p+\frac{1}{I^{\sigma }}+\frac{1}{U^{\sigma }},0<p,\sigma <1,$where $\sigma >0$ is a sufficiently small constant. Take $\sigma \in (0,1)$ such that $\tau b-\frac{\sigma +1}{2}{\alpha }_1^2>0,b-\frac{\sigma +1}{2}{\alpha }_2^2>0,$then by use of Itô's formula on V yields $dV=LVdt+(p{I}^p-\sigma I^{-\sigma }){\alpha }_{1}d{B}_1\left(t\right)+(p{U}^p-\sigma U^{-\sigma }){\alpha }_{2}d{B}_2\left(t\right),$where $\begin{array}{rl}LV(I,U)& =p{I}^p[\tau b-(d+D)I-(d+D)U+\frac{\theta }{I}+\frac{p-1}{2}{\alpha }_1^2]\\ & +p{U}^p[\frac{(1-\tau )bI}{U}+b-\frac{qbI}{U+I}-dI-dU+\frac{p-1}{2}{\alpha }_2^2]\\ & -\sigma I^{-\sigma }[\tau b-(d+D)I-(d+D)U+\frac{\theta }{I}-\frac{\sigma +1}{2}{\alpha }_1^2]\\ & -\sigma U^{-\sigma }[\frac{(1-\tau )bI}{U}+b-\frac{qbI}{U+I}-dI-dU-\frac{\sigma +1}{2}{\alpha }_2^2]\\ & \le -p(d+D)I^{p+1}+p(\tau b+\frac{p-1}{2}{\alpha }_1^2)I^p+p\theta I^{p-1}\\ & -d{U}^{p+1}+(b+\frac{p-1}{2}{\alpha }_2^2)U^p+p(1-\tau )bI{U}^{p-1}\\ & -\sigma I^{-\sigma }[\tau b-\frac{\sigma +1}{2}{\alpha }_1^2]+\sigma (d+D)I^{1-\sigma }+\sigma I^{-\sigma }(d+D)U\\ & -\sigma U^{-\sigma }[b-\frac{\sigma +1}{2}{\alpha }_2^2]+\sigma qbI{U}^{-\sigma }+\sigma d{U}^{1-\sigma }+\sigma U^{-\sigma }dI\\ & \le -p(d+D)I^{p+1}+p(\tau b+\frac{p-1}{2}{\alpha }_1^2)I^p+p\theta I^{p-1}\\ & -\sigma I^{-\sigma }[\tau b-\frac{\sigma +1}{2}{\alpha }_1^2]+\sigma (d+D)I^{1-\sigma }+\sigma qbI+\sigma dI\\ & -d{U}^{p+1}+(b+\frac{p-1}{2}{\alpha }_2^2)U^p+pb{M}^p\\ & -\sigma U^{-\sigma }[b-\frac{\sigma +1}{2}{\alpha }_2^2]+\sigma I^{-\sigma }(d+D)U+\sigma d{U}^{1-\sigma },\end{array}$where M is an upper bound which satisfies $M\ge I$ and $M\ge U$ due to the solution $X\left(t\right)=\left(I\right(t),U(t\left)\right)$ of system (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ) being stochastically ultimately bounded.

It is easy to obtain that $LV(I,U)\le {\phi }_1\left(I\right)+{\phi }_2\left(U\right),$where $\begin{array}{rl}{\phi }_1\left(I\right)& =-p(d+D)I^{p+1}+p(\tau b+\frac{p-1}{2}{\alpha }_1^2)I^p+p\theta I^{p-1}\\ & -\sigma I^{-\sigma }[\tau b-\frac{\sigma +1}{2}{\alpha }_1^2]+\sigma (d+D)I^{1-\sigma }+\sigma qbI+\sigma dI,\\ {\phi }_2\left(U\right)& =-d{U}^{p+1}+(b+\frac{p-1}{2}{\alpha }_2^2)U^p+pb{M}^p\\ & -\sigma U^{-\sigma }[b-\frac{\sigma +1}{2}{\alpha }_2^2]+\sigma I^{-\sigma }(d+D)U+\sigma d{U}^{1-\sigma }.\end{array}$Case 1. If $I\to 0^{+}$, then $LV={\phi }_1\left(I\right)+{\phi }_2\left(U\right)\le \phi \left(I\right)+{\phi }_2^u\to -\mathrm{\infty }.$If $U\to 0^{+}$, then $LV={\phi }_1\left(I\right)+{\phi }_2\left(U\right)\le {\phi }_1^u+\phi \left(U\right)\to -\mathrm{\infty }.$Case 2. If $I\to +\mathrm{\infty }$, then $LV={\phi }_1\left(I\right)+{\phi }_2\left(U\right)\le \phi \left(I\right)+{\phi }_2^u\to -\mathrm{\infty }.$If $U\to +\mathrm{\infty }$, then $LV={\phi }_1\left(I\right)+{\phi }_2\left(U\right)\le {\phi }_1^u+\phi \left(U\right)\to -\mathrm{\infty }.$In conclusion, when $I\to 0^{+}$ or $U\to 0^{+}$ or $I\to +\mathrm{\infty }$ or $U\to +\mathrm{\infty }$, we can get $LV\to -\mathrm{\infty }$. Hence, choose $\kappa >0$ sufficiently small and let $E:=[\kappa ,\frac{1}{\kappa }]\times [\kappa ,\frac{1}{\kappa }]$, then $LV(I,U)\le -1\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}(I,U)\in \mathrm{I}\mathrm{n}\mathrm{t}{\mathrm{R}}_{+}^2\setminus E,$which means that the condition (1) holds true. Furthermore, $\begin{array}{rl}\sum _{i,j=1}^2{b}_{ij}(I,U){\xi }_i{\xi }_j& ={\alpha }_1^2{I}^2{\xi }_1^2+{\alpha }_2^2{U}^2{\xi }_2^2\\ & \ge \underset{(I,U)\in E}{min}\{{\alpha }_1^2{I}^2,{\alpha }_2^2{U}^2\}\Vert \xi {\Vert }^2\\ & \mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}(I,U)\in E,({\xi }_1,{\xi }_2)\in {\mathrm{R}}^2,\end{array}$and thus condition (2) has also been satisfied. Therefore, the system (3(3) $\{\begin{array}{ll}dI\left(t\right)& =\left(\tau bI\right(t)-(d+D\left)\right(I\left(t\right)+U\left(t\right)\left)I\right(t)+\theta )\mathrm{d}t+{\alpha }_{1}I\left(t\right)d{B}_1\left(t\right),\\ dU\left(t\right)& =\left(\right(1-\tau \left)bI\right(t)+bU(t)(1-\frac{qI\left(t\right)}{U\left(t\right)+I\left(t\right)})-d(I\left(t\right)+U\left(t\right)\left)U\right(t\left)\right)\mathrm{d}t\\ & +{\alpha }_{2}U\left(t\right)d{B}_2\left(t\right),\end{array}$(3) ) has a unique ergodic stationary distribution (Figure 3). This completes the proof.

Figure 3. Stationary distribution of deterministic model and stochastic model: (a) we set initial values as $\left(I\right(0),U(0)=(0.5,0.5)$; (b) we set initial values as $\left(I\right(0),U(0)=(1,1)$. The initial values of the solution illustrated by the black line were fixed as $\left(X_1\right(0),X_2(0),Y(0\left)\right)=(10,10,0.5)$, and all other parameters were fixed as: $\tau =0.3$, b = 0.1, d = 0.1, D = 0.01, $\theta =0.2$, q = 0.2, q = 0.3.

4. Conclusion

Many researchers are designing programmes to release $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$-infected mosquitoes to control dengue virus, a promising strategy that has attracted the attention of many mathematical researchers. Given that mosquitoes in natural habitats are inevitably affected by environmental fluctuations [18, 32, 36], we developed a mathematical model comprising a system of stochastic differential equations, governing the evolution of mosquitoes with white noise.

We first prove the existence and uniqueness of the positive solutions of the proposed system. Then we study the V-geometrical ergodicity and stochastic ultimately boundedness of the system. Further, threshold conditions for successful population replacement are derived, and it is shown that the system has a unique ergodic steady-state distribution. The results show that the ratio of infected to uninfected mosquitoes has a great influence on population replacement. Moreover, environmental noise plays a significant role in control of dengue fever.

The highlights are listed as follows: (1) the proposed model considers not only the influence of white noise but also introduces the continuously released $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$-infected mosquitoes; (2) the threshold conditions for the extinction of uninfected mosquitoes and the system's stationary distribution are obtained; (3) biologically, by increasing the release of $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$-infected mosquitoes, resulting in an increase in the ratio of $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$-infected mosquitoes and uninfected mosquitoes, and then reaching a suitable value, which can promote the rapid extinction of uninfected mosquitoes.

In fact, some studies have shown that different $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$ strains have different biological characteristics, which are different in reducing the incidence rate of dengue fever [2, 29, 30]. For example, the $\mathit{w}$Au strain has a high virus blocking rate, while the $\mathit{w}$Mel strain has only a moderate level, but both strains have high maternal transmission rates. In addition, the loss of $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$ infection is lower in the $\mathit{w}$Au strain and higher in the $\mathit{w}$Mel strain. Finally, the $\mathit{w}$Au strain does not exhibit cytoplasmic incompatibility, while $\mathit{w}$Mel exhibits cytoplasmic incompatibility. Therefore, further research on the impact of environmental factors on the effectiveness of different $\mathit{W}\mathit{o}\mathit{l}\mathit{b}\mathit{a}\mathit{c}\mathit{h}\mathit{i}\mathit{a}$ strains may add new perspectives, which will also become our future research direction.

Acknowledgments

The author is very grateful to the anonymous referee for a careful reading, helpful suggestions and valuable comments which led to the improvement of the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the National Natural Science Foundation of China under grants (11961024(Y. Tan), 11801047(J. Yang)), and by the Joint Training Base Construction Project for Graduate Students in Chongqing (JDLHPYJD2021016), and by the Group Building Scientific Innovation Project for universities in Chongqing (CXQT21021).