Threshold dynamics of a stochastic general SIRS epidemic model with migration

1.   Introduction

Zika virus disease, referred to as Zika, is a mosquito-borne infectious disease induced by Zika virus originally found in rhesus monkeys in the jungle of Uganda in 1947 ^[1] and afterward isolated from humans in Uganda and Tanzania in 1952 ^[2]. In the following decades, only a few cases from Africa and Southeast Asia were reported sporadically ^[3] until 2007 when Zika broke out on Yap Island in Micronesia in the western Pacific Ocean ^[4]. In early 2015, researchers detected Zika virus infection in Brazil ^[5]. The virus spread rapidly to Northern Europe, Australia, the United States, Canada ^[6,7,8,9], and then to Japan, China, India, and other countries ^[10,11,12], causing great harm to human health. At present, Zika is still prevalent at a low level in Central and South America. From January 1 to April 30, 2022, a total of 6171 suspected cases of Zika had been reported in Brazil, of which 541 cases were confirmed ^[13]. Zika remains an essential global public health challenge.

The reason why Zika spreads fast and widely is mainly its multiple transmission channels. The virus is primarily propagated to mankind via the biting of infected Aedes aegypti and Aedes albopictus ^[14]. Meanwhile, it can be spread among humans through heterosexual or homosexual sexual contact ^[15,16]. In addition, it can also be spread from infected female mosquitoes to their descendants vertically ^[17] and from the water contaminated by the urine of the infected person to the mosquitoes in the aquatic stage ^[18]. The latency of the virus in the human body is generally 3–14 days ^[14]. The majority of infected people are asymptomatic, and only a quarter are believed to develop slight symptoms such as fever, erythema, conjunctivitis, and arthrodynia, with only a handful of documented fatalities ^[14]. Although the mortality of Zika virus disease is meager, it is believed that Zika infection during the gestational period is one of the causes of microcephaly and other congenital malformations in developing fetuses and newborns ^[19]. Zika infection is also a trigger factor for Guillain-Barre syndrome, myelitis, and neuropathy, especially in adults and older children ^[20]. Unfortunately, there is still no allowable vaccination or antivirus drug for the virus.

As we all know, mathematical modeling is an effective and indispensable tool for a better understanding of population dynamics and epidemics ^[21]. Using this tool, many scholars have conducted rich and detailed research on the transmission of Zika disease, see ^{[22,23,24,25,26,27]} and their references. For example, Gao et al. ^[22] proposed an ordinary differential equation (ODE) model to examine the influences of media transmission and sexual transmission on the propagation of Zika disease and carried out a sensitive analysis of basic reproduction number. Agusto et al. ^[23] established an ODE model of Zika virus, including human vertical transmission, the birth of babies with microcephaly, and asymptomatic infection, and studied the dynamic behavior of the model. Considering the limitation of medical resources during the outbreak of Zika, Zhao et al. built an ODE model of Zika to investigate the effect of medical resources on the spread of Zika ^[24]. In addition, due to the impact of spatial differentiation and spatial mobility of human and vector populations on the dynamics of vector-borne diseases, some reaction-diffusion models for describing the spatial transmission of Zika virus have been developed accordingly, see ^[28,29,30] and references therein. However, the transmission of Zika virus is also influenced by temperature, wind, rain, fire, and other random environmental factors, which are ignored by the deterministic model. Using a stochastic differential equation (SDE) model to describe the epidemic dynamics can better reflect the actual phenomenon to some extent. The extinction and persistence of the epidemic driven by random noise have been studied in some works of literature, see ^{[31,32,33,34]}. Nevertheless, to our knowledge, there are few documents on the dynamic analysis of infectious diseases considering both random factors and spatial diffusion. Therefore, this paper intends to explore the permanence and extinction of Zika disease described by a random model with spatial diffusion, which includes human-mosquito transmission, human-human sexual transmission, and vertical transmission of mosquitoes, to fill this gap.

The outbreak and prevalence of Zika have brought enormous economic burdens and health losses to the local people and government. Therefore, from the perspective of epidemiology and social economics, it is an essential and meaningful problem to formulate the optimal control strategy for Zika virus, that is, to achieve the greatest limitation of the disease with the least cost. There have been some studies on the optimal control problem of Zika. For instance, the literature ^[35] introduced vaccination as a control variable (although the vaccine has not been publicly available yet) and characterized the most economical and effective vaccination strategy in a reaction-diffusion model of Zika virus by utilizing the optimal control theory. The control variables (the prevention through mosquito nets, the treatment of infection patients, and the spraying of insecticides) were selected into the ODE model by authors in ^[36,37] to establish an optimal control problem of Zika virus. Their numerical simulation results suggested that it might be more beneficial to eliminate Zika virus infection if all three control measures are considered. This paper also plans to adopt three control variables, namely, personal protection, treatment of infected humans (here we use saturation treatment function due to limited medical resources), and reduction in the number of mosquitoes, and draw them into the SRDM to generate a stochastic control model of Zika virus.

The rest of this paper is arranged as follows. In Section 2, we present the model description of Zika virus and prove the existence, uniqueness, and boundedness of the global positive solution of this model. In Sections 3 and 4, we discuss the conditions of disease extinction as well as the existence conditions of the stationary distribution. A stochastic optimal control problem is proposed and the expressions of the optimal controls are acquired in Section 5. Some numerical simulations are performed in Section 6 to declare and supplement the theoretical contents. At last, a summary is made in Section 7.

2.   Model formulation and preliminaries

2.1.   Model formulation

According to the transmission mechanism of Zika virus, we plot the flow chart in Figure 1. In the flow chart, $S_h(t, x)$, $E_h(t, x)$, $I_h(t, x)$, and $R_h(t, x)$ are the number of susceptible, exposed, infected, and recovered human population at time $t$ and location $x$, respectively. The total number of humans is $N_h = S_h+E_h+I_h+R_h$. The number of susceptible, exposed, and infected mosquitoes at time $t$ and position $x$, are recorded as $S_m(t, x)$, $E_m(t, x)$, and $I_m(t, x)$, respectively. Thus, $N_m = S_m+E_m+I_m$ is the total amount of mosquitoes. Since only female mosquitoes suck blood and transmit diseases, the mosquitoes in this paper only refer to female mosquitoes.

A susceptible human may be infected with Zika virus via the bite of an exposed or infected mosquito at a rate $\lambda_{mh}(t, x) = b_m(t, x)\beta_h(k_hE_m+I_m)/N_h$ or through sexual contact with an exposed or infected partner at a rate $\lambda_{hh}(t, x) = \beta_{hh}(k_{hh}E_h+I_h)/N_h$, where $b_m(t, x) = b(t, x)/N_m$, $b(t, x) = \frac{\alpha_hN_h\alpha_mN_m}{\alpha_hN_h+\alpha_mN_m}$ is the total number of bites per day at position $x$ ^[38,39], and so $b_m(t, x)$ is the average number of bites per mosquito per day at position $x$. Susceptible mosquitoes move to the exposed class after biting exposed or infected humans at a rate $\lambda_{hm}(t, x) = b_h(t, x)\beta_m(k_mE_h+I_h)/N_m$, here $b_h(t, x) = b(t, x)/N_h$ is the average number of bites per day for an infectious person at position $x$. This paper considers the infectivity of humans and mosquitoes during the exposure period and the modification parameters $0 < k_h, k_{hh}, k_m < 1$ measure the reduction in transmissibility during the exposure period relative to the infection period.

Zika virus can also be spread vertically from infected mothers to newborns and from infected female mosquitoes to their offspring ^[17,40]. This paper only deals with vertical transmission in mosquitoes, neglecting vertical transmission in humans because Zika has a very short transmission period compared to the human lifespan ^[26]. We assume that $k\theta\mu_mE_m+\theta\mu_mI_m$ of the mosquito's offspring will be infected, and thus enter $E_m$ class, where $\mu_m$ is the average birth rate of mosquitoes, $\theta$ is the proportion of congenital infections in the progeny of infectious female mosquitoes, and $0 < k < 1$ is also a modification parameter.

Since the symptoms of Zika are slight and rarely fatal, we ignore the human mortality caused by the disease. And because of the short lifespan of mosquitoes, we assume that infected mosquitoes will not recover until natural death and that these mosquitoes will not die from Zika.

Based on the above description and the flow chart in Figure 1, and taking into account the move of humans and mosquitoes, we establish the following reaction-diffusion system for Zika virus

for $t > 0, x\in Q$, with the boundary conditions

and initial conditions

here $Q$ is a bounded region possessing smooth boundary $\partial Q$ and ${\bf n}$ is the outward unit normal vector on $\partial Q$. $d_1, d_2, d_3,$ and $d_4$ represent the diffusion coefficients of susceptible, exposed, infected, and recovered human population, respectively, and $d_5, d_6,$ and $d_7$ denote the diffusion coefficients of susceptible, exposed, and infected mosquitoes population, respectively, $d_i > 0, \; i = 1, 2, \cdots, 7$. The meanings of the remaining parameters of model (2.1) are explained in Table 1.

Because parameters in the infectious disease model are often subject to environmental noise and exhibit random fluctuations to a certain extent, this paper intends to build a stochastic Zika model by perturbing the natural death rates $d_h$ and $d_m$ for humans and mosquitoes with white noise. In other words, we will replace $d_h$ and $d_m$ in model (2.1) with $d_h-\sigma_1\dot{B_1}(t)$ and $d_m-\sigma_2\dot{B_2}(t)$, respectively, where $B_1(t)$ and $B_2(t)$ are independent standard Brownian motions in the complete probability space $(\Omega, \mathcal{F}, \mathbb{P})$ with a filtration $\{\mathcal{F}_t\}_{t\geq0}$, which is increasing, right continuous, and satisfies that $\mathcal{F}_0$ involves all $\mathbb{P}$-null sets. $\sigma_1 > 0$ and $\sigma_2 > 0$ are the intensities of the noise. Then the corresponding stochastic system of model (2.1) has the following form

for $t > 0, x\in Q$, with the boundary conditions

and initial conditions

2.2.   Preliminaries

Let $H = H^1(Q) = \{\varphi|\varphi\in L^2(Q), \frac{\partial\varphi}{\partial x_i} \in L^2(Q) \text{ is generalized partial derivative, i = 1, 2, 3}\}$. $H$ is a Sobolev space and $H\hookrightarrow L^2(Q)\hookrightarrow H'$, where $H' = H^{-1}(Q)$ is the dual space of $H$. $\|\cdot\|$ and $\|\cdot\|_\ast$ are the norms of $H$ and $L^2(Q)$, respectively. $\|\varphi\|^2 = \|\varphi\|_\ast^2+\|\nabla\varphi\|_\ast^2$, and there exists a positive constant $c$ such that $\|\varphi\|_\ast\leq c\|\varphi\|$. $\langle\cdot, \cdot\rangle$ indicates the dual product of $H$ and $H'$. The norm of Euclidean space is denoted by $|\cdot|$. $\mathcal{H} = H^7$. Denote $H^+ = \{\varphi|\varphi\in L^2(Q; (0, \infty)), \frac{\partial\varphi}{\partial x_i} \in L^2(Q), i = 1, 2, 3\}$, $\mathcal{H}^+ = (H^+)^7$. In addition, $\mathbb{R}_+^l = \{(x_1, x_2, \cdots, x_l)\in \mathbb{R}^l:x_i > 0, i = 1, 2, \cdots, l\}$, $\mathbb{R}_+ = [0, \infty)$. $L^2_\mathcal{F}([0, T]\times Q; \mathbb{R}^l)$ is a set of square integrable and $\mathcal{F}_t$-adapted stochastic processes. The indicative function of set $A$ is denoted by $\chi_A$. $a\wedge b = \min\{a, b\}$, $a\vee b = \max\{a, b\}$. $\varphi_x$ represents the partial derivative of $\varphi$ to $x$. $B(t, x)$ is sometimes abbreviated to $B$ for convenience without causing confusion.

Theorem 2.1. For any initial value $X(0, x) = (S_h^0(x), E_h^0(x), I_h^0(x), R_h^0(x), S_m^0(x), E_m^0(x), I_m^0(x))\in\mathcal{H}^+$, stochastic Zika system (2.2) has a unique global positive solution $X(t, x) = (S_h(t, x), E_h(t, x), I_h(t, x), R_h(t, x), S_m(t, x), E_m(t, x), I_m(t, x))\in\mathcal{H}^+$ on $t\geq0$. Moreover, there is a positive constant $C_0$ such that

Theorem 2.1 is an important fundamental theorem, which gives the existence, uniqueness, and boundness of the positive solution of system (2.2), and its proof is shown in Appendix A. The following theorem discusses the $p$th moment boundedness of system (2.2) and its proof can be found in Appendix B.

Theorem 2.2. For any $p > 0$, we have

where $C$ is a constant related to $p, T,$ and the original condition and the parameters of system (2.2).

3.   Extinction of the disease

In this section, we will discuss the conditions for almost surely exponential extinction of Zika disease.

In general, consider an $l$-dimensional stochastic reaction-diffusion system by

with boundary condition $\frac{\partial\upsilon(t, x)}{\partial {\bf n}} = 0 \; (t > t_0, x\in\partial Q)$ and initial condition $\upsilon(t_0, x) = \upsilon_0(x)\; (x\in Q).$

We give the definition of the almost surely exponential stability of system (3.1) ^[44].

Definition 3.1. The trivial solution of system (3.1) is said to be almost surely exponentially stable if

for all $\upsilon_0\in \mathbb{R}_+^l$, where $\big|\log|\upsilon(\cdot, x)|\big|_Q: = \int_Q\log|\upsilon(\cdot, x)|dx$.

Next, the almost surely exponential extinction of Zika disease will be given in the following theorem.

Theorem 3.2. For any starting value $X(0, x)\in\mathcal{H}^+$ of system (2.2), if

and

then

Proof. Let $V(t, x) = \log\big(E_h(t, x)+\frac{\xi_h+d_h}{\xi_h}I_h(t, x)+E_m(t, x)+\frac{\xi_m+d_m}{\xi_m}I_m(t, x)\big)$, then

Let $A(t, x) = \big(E_h(t, x)+\frac{\xi_h+d_h}{\xi_h}I_h(t, x)+E_m(t, x)+\frac{\xi_m+d_m}{\xi_m}I_m(t, x)\big)$. By $It\hat{o}$'s formula, we have

Integrating the above inequality from $0$ to $t$, we get

where

In addition, the quadratic variations of $M_1$ and $M_2$ are respectively

Therefore

Thus, martingale's strong law of large numbers yields

Together with (3.2)–(3.4), , we obtain

which shows that

This completes the proof. □

Conditions (3.2) and (3.3) of Theorem 3.2 suggest that Zika virus will become exponentially extinct when the diffusion coefficients of infected people and mosquitoes are very small, that is, they hardly move, and the intensities of environmental noise are relatively large. There is no doubt that such conditions are very harsh. In what follows we will talk about the stationary distribution of system (2.2), which means the persistence of Zika disease.

4.   Stationary distribution

First of all, we give the definition of the stationary distribution of system (2.2) ^[45]. Let $\mathcal{P}(\mathcal{H})$ represent the space of all probability measures on $(\mathcal{H}, \mathcal{B}(\mathcal{H}))$, here $\mathcal{B}(\mathcal{H})$ denotes the Borel $\sigma$-algebra on $\mathcal{H}$. $\mathcal{C}_b(\mathcal{H})$ is the set of all bounded and continuous real-valued functions on $\mathcal{H}$.

Definition 4.1. A stationary distribution of the positive solution $X(t, x)$, $t\geq0$, of system (2.2) is defined as a probability measure $\pi\in\mathcal{P}(\mathcal{H})$ which satisfies

here $\pi(g): = \int_{\mathcal{H}}g(\phi)\pi(d\phi)$, $\mathbb{P}_tg(\phi): = \mathbb{E}g(|X(t, x, \phi)|_Q)$, and $g\in\mathcal{C}_b(\mathcal{H})$.

For $\pi_1, \pi_2\in\mathcal{P}(\mathcal{H})$, a measure on $\mathcal{P}(\mathcal{H})$ is defined by

where $\mathcal{N}: = \{g:\mathcal{H}\rightarrow \mathbb{R}, |g(\phi_1)-g(\phi_2)|\leq\|\phi_1-\phi_2\|\text{ for any }\phi_1, \phi_2\in \mathcal{H}\text{ and }|g(\cdot)|\leq1\}$. $\mathcal{P}(\mathcal{H})$ is complete under the measure $d(\cdot, \cdot)$ by ^[46], and then we can get an important lemma as follows, which provides an assertion for the existence of stationary distribution ^[45].

Lemma 4.1. Assuming that for arbitrary bounded subset $\mathcal{O}$ of $\mathcal{H}^+$, $p > 1,$

$(i) \lim_{t\rightarrow \infty}\sup_{\phi_1, \phi_2\in\mathcal{O}} \mathbb{E}\|X(t, x, \phi_1)-X(t, x, \phi_2)\|^p = 0$;

$(ii)\sup_{t\geq 0}\sup_{\phi\in\mathcal{O}} \mathbb{E}\|X(t, x, \phi)\|^p < \infty$.

Then, $X(t, x, \phi), t\geq 0$, has a stationary distribution for initial data $\phi\in\mathcal{H}^+$.

Applying Lemma 4.1, we can obtain the conditions for the existence of the steady-state distribution of system (2.2).

Theorem 4.2. Assume there are constants $p > 1$, $\vartheta > 0$, and $0 < c_i < 1\; (i = 1, 2, \cdots, 7)$ such that

and

where $\hat{A_1} = 2\beta_m^p\vee\beta_h^p(1+k_h^p)$, $\hat{A_2} = \xi_h^p\vee\gamma^p$, $\hat{A_3} = 2\beta_h^p\vee\beta_m^p(1+k_m^p)$, $\check{B_1} = c_1d_1\wedge c_2d_2\wedge c_3d_3\wedge c_4d_4$, $\check{B_2} = c_5d_5\wedge c_6d_6\wedge c_7d_7$, then process $X(t, x), t\geq 0,$ of system (2.2) has a unique stationary distribution $\pi\in\mathcal{P}(\mathcal{H})$.

Proof. To illustrate the existence of steady-state distribution of system (2.2), we need to prove that the conditions $(i)$ and $(ii)$ in Lemma 4.1 hold. Since Theorem 2.2 implies that $(ii)$ is true, we only need to verify $(i)$. To this end, for $p > 1$, $\vartheta > 0$, let

where

Applying $It\hat{o}$'s formula, we deduce

For simplicity, we assume that $N_h(t, x, \phi_1) = N_h(t, x, \phi_2) = N_h(t, x)$ and $N_m(t, x, \phi_1) = N_m(t, x, \phi_2) = N_m(t, x)$. Thus denote $\lambda_{mh1}(t, x) = b_m(t, x)\beta_h(k_hE_{m1}+I_{m1})/N_h = \frac{\alpha_h\alpha_m}{\alpha_hN_h+\alpha_mN_m}\beta_h(k_hE_{m1}+I_{m1})$, $\lambda_{mh2}(t, x) = b_m(t, x)\beta_h(k_hE_{m2}+I_{m2})/N_h = \frac{\alpha_h\alpha_m}{\alpha_hN_h+\alpha_mN_m}\beta_h(k_hE_{m2}+I_{m2})$, $\lambda_{hh1}(t, x) = \beta_{hh}(k_{hh}E_{h1}+I_{h1})/N_h$, and $\lambda_{hh2}(t, x) = \beta_{hh}(k_{hh}E_{h2}+I_{h2})/N_h$. By $It\hat{o}$'s formula and embedding theorem,

here $c_1 < 1$ is a constant and the last inequality sign takes advantage of the Young inequality. Similarly,

Substituting (4.4)–(4.10) into (4.3), integrating the two sides of (4.3), and seeking mathematical expectation, then

where

and $\vartheta+C_7 > 0$. Next, we take the $\sup\limits_{\phi_1, \phi_2\in \mathcal{O}}$ and use the Gronwall inequality to get

i.e.,

According to (4.1) and (4.2), $C_7 < 0$. Therefore,

Thus, the condition $(i)$ of Lemma 4.1 is proved. Let us now explain the uniqueness of steady-state distribution of system (2.2).

Suppose $\pi'\in\mathcal{P}(\mathcal{H})$ is another steady-state distribution for $X(t, x), t\geq 0,$ of system (2.2). $\mathcal{C}_{lb}(\mathcal{H})$ is a bounded and Lipschitz continuous function family on $\mathcal{H}$. Then by the definition of stationary distribution, the Hölder inequality, and (4.11), for $g\in\mathcal{C}_{lb}(\mathcal{H})$, we can derive that

here $C_8 > 0$ is a constant. Whereupon, the uniqueness of stationary distribution can be obtained by setting $t\rightarrow \infty$ in (4.12) when $C_7 < 0$. The proof is completed. □

From (4.1) and (4.2) of Theorem 4.2, we find that Zika disease will be persistent when the intensities of environmental noise are low, while the diffusion coefficients of humans and mosquitoes are relatively large, which is the opposite of exponential extinction.

5.   Optimal control problem

The objective of this section is to illustrate that anti-Zika control strategies can be implemented while minimizing the cost of implementing these measures. So we formulate a stochastic optimal control problem by introducing three control variables into system (2.2). The control $u_1(t, x)$ denotes the level of personal protective efforts among the population, so the correlative infectivity is decreased by the factor $(1 -u_1(t, x))$. The control $u_2(t, x)$ represents the level of treatment for infected people. We choose saturated treatment rate function $\frac{cu_2I_h}{1+\alpha I_h}$ with treatment rate $c > 0$ and saturation coefficient $\alpha\geq 0$ due to the limited medical resources (medical staff, medicines, hospital beds, etc.), where $\frac{c}{\alpha}$ is the largest medical resource provided per unit of time. The control $u_3(t, x)$ indicates the level of insecticides used to kill mosquitoes in mosquito breeding grounds, which increases the mosquito mortality rate from $d_m$ to $d_m+c_0u_3$ with killing efficacy $c_0$. In this thesis, $0\leq u_i\leq 1\; (i = 1, 2, 3)$ means that there is no effort (i.e., no control) when the control is zero, and the maximum control is put when the control is one. Let $u = (u_1, u_2, u_3)$. Thus the stochastic control system for Zika disease will be written as

for $t > 0, x\in Q$, with the boundary conditions

and initial conditions

Our optimal control study aims at minimizing the number of exposed and infected people, the total number of mosquitoes, and the cost of executing the control in time interval $[0, T]$ and region $Q$. In order to realize this goal, an objective functional is defined as

with

where $a_1, a_2,$ and $a_3$ are positive coefficients of weight of the exposed, infected human and the total mosquito populations, respectively, $b_1, b_2,$ and $b_3$ are positive coefficients of weigh for the linear costs of personal protection, the treatment for infected people, and mosquito control, respectively, $c_1, c_2,$ and $c_3$ are positive coefficients of weight for the quadratic costs, respectively. A function of $X(t, x)$ at the terminal time $T$ is denoted by $\varphi(X(T, x))$. The next task is to find the optimal control $\bar{u} = (\bar{u}_1, \bar{u}_2, \bar{u}_3)$ such that

here $U$ is the admissible control set as follows

Similar to Theorem 2.1, the existence, uniqueness, and boundness of the positive solution of system (5.1) can also be verified. Further, we can obtain the boundedness and convexity of $f_i(t, x, X, u)(i = 1, 2, \cdots, 7)$ and the compactness of $U$, and then the existence of optimal control $\bar{u}$ can be shown according to Theorem 3.1 in ^[47].

Next, we will makes use of the Pontryagin maximum principle ^[48] to obtain the optimal control. Denote $\lambda(t, x) = (\lambda_1(t, x), \lambda_2(t, x), \cdots, \lambda_7(t, x))$, $\mu(t, x) = (\mu_1(t, x), \mu_2(t, x), \cdots, \mu_7(t, x))$. $\bar{X} = (\bar{X}_1, \bar{X}_2, \bar{X}_3, \bar{X}_4, \bar{X}_5, \bar{X}_6, \bar{X}_7) = (\bar{S}_h, \bar{E}_h, \bar{I}_h, \bar{R}_h, \bar{S}_m, \bar{E}_m, \bar{I}_m)$ is the optimal state variable of system (5.1) corresponding to the optimal control $\bar{u}$. Then by the stochastic maximum principle, there exists a pair of processes $(\lambda(t, x), \mu(t, x))\in L^2_\mathcal{F}([0, T]\times Q; \mathbb{R}^7)\times L^2_\mathcal{F}([0, T]\times Q; \mathbb{R}^7)$ that satisfy the following SDE

where

and $\bar{N}_h = \bar{S}_h+\bar{E}_h+\bar{I}_h+\bar{R}_h$, $\bar{N}_m = \bar{S}_m+\bar{E}_m+\bar{I}_m$, $\bar{\lambda}_{mh} = \frac{\alpha_m\alpha_h\beta_h(k_h\bar{E}_m+\bar{I}_m)}{\alpha_m\bar{N}_m+\alpha_h\bar{N}_h}$, $\bar{\lambda}_{hh} = \frac{\beta_{hh}(k_{hh}\bar{E}_h+\bar{I}_h)}{\bar{N}_h}$, $\bar{\lambda}_{hm} = \frac{\alpha_m\alpha_h\beta_m(k_m\bar{E}_h+\bar{I}_h)}{\alpha_m\bar{N}_m+\alpha_h\bar{N}_h}$.

Define the following Hamilton function:

for $(t, X, u, \lambda, \mu)\in [0, T]\times\mathcal{H^+}\times U\times \mathbb{R}^7\times \mathbb{R}^7$. So, according to the maximum condition of stochastic maximum principle, from $\frac{\partial H(t, X, u, \lambda, \mu)}{\partial u_j} = 0$ and $0\leq u_j\leq 1, j = 1, 2, 3$, we get the following conclusion:

Theorem 5.1. Under objective functional (5.2), the expressions for the optimal controls of system (5.1) are

6.   Numerical simulations

In this part, some numerical simulations will be conducted to illustrate our theoretical results more intuitively. We can write the discrete form (13.1) of the state Eq (2.2) shown in Appendix C using Milstein's method ^[49]. The initial conditions of model (2.2) and model (5.1) are selected as $S_h^0(x) = 750000+200\cos\frac{\pi x}{30}, \; E_h^0(x) = 150+10\cos\frac{\pi x}{30}, \; I_h^0(x) = 10+5\cos\frac{\pi x}{30}, \; R_h^0(x) = 0, \; S_m^0(x) = 80000+100\cos\frac{\pi x}{30}, E_m^0(x) = 100+10\cos\frac{\pi x}{30}, \; I_m^0(x) = 10+2\cos\frac{\pi x}{30}, \; x\in(0,100).$

6.1.   Disease extinction

We choose parameters $\beta_h = 0.1, \beta_{hh} = 0.001, \beta_m = 0.3, \alpha_h = 0.1, \alpha_m = 0.19, \xi_h = 0.14, \xi_m = 0.14, d_h = 0.000039, d_m = 0.03, \gamma = 0.14, \mu_m = 0.029, \theta = 0.001, d_1 = 0.012, d_2 = d_3 = 0.0002, d_4 = d_5 = 0.008, d_6 = d_7 = 0.0001, \sigma_1 = \sigma_2 = 0.59.$ Please refer to Table $1$ for other parameter values. Through calculation, conditions (3.2) and (3.3) for the almost surely exponential extinction of Zika disease in Theorem 3.2 are established at this time. Under these parameter values, we draw the evolutions of $E_h, I_h, E_m$, and $I_m$ as shown in Figure 2. Here, the left is the spatio-temporal graphs, the right is the relevant projection graphs, where the curves of different colors indicate the variations of the population in different regions over time. Obviously, Figure 2 verifies the conclusion of Theorem 3.2.

It can be seen from the above, to make Zika disease almost surely exponentially extinct, on the one hand, the intensities of noise should be higher; on the other hand, the diffusion coefficients of the infected people and mosquitoes should be very small, which can guide us on how to eliminate Zika disease faster.

6.2.   The existence of the stationary distribution

In this subsection, the existence of the stationary distribution of system (2.2) will be numerically simulated. We take $\beta_h = 0.23, \beta_{hh} = 0.10, \beta_m = 0.33, \alpha_h = 0.1, \alpha_m = 0.26, \xi_h = 0.2, \xi_m = 0.1, \sigma_1 = 0.1, \sigma_2 = 0.1, d_1 = 0.6, d_2 = 0.6, d_3 = 0.55, d_4 = 0.58, d_5 = 0.3, d_6 = 0.3, d_7 = 0.25, p = 1.01, c_i = 0.96 (i = 1, 2, \cdots, 7).$ Other parameter values are the same as those in Figure 2. At this time, conditions (4.1) and (4.2) in Theorem 4.2 describing the existence of steady-state distribution hold. The trajectories of the solution of system (2.2) are presented in Figures 3 and 4, whose left side is the spatio-temporal graphs and the right side is the corresponding projection graphs, which indicate that the system will achieve a steady state over time $t$. Figures 5 and 6 show the evolutions of the solution of the system when spatial variable $x = 10$ and their corresponding histograms, from which we can see that the system has a stationary distribution.

6.3.   The impacts of noise and diffusion coefficients on disease

1) The impact of noise on disease

We choose the same parameters as in Figure 2 and take the values of noise as $\sigma_1 = \sigma_2 = 0$, $\sigma_1 = \sigma_2 = 0.05$, $\sigma_1 = \sigma_2 = 0.3$, $\sigma_1 = \sigma_2 = 0.8$. Figure 7 describes the variations of $I_h(t, x)$ and $I_m(t, x)$ under different noise intensities when $x = 10$. We can observe that a smaller noise intensity has a slight fluctuation in the number of infected people and mosquitoes, however, as the noise intensity enhances, the number of infected persons and mosquitoes decreases significantly. Therefore, we can consider random noise as a control strategy, such as human treatment and mosquito repellent spraying, to achieve the control of Zika disease.

2) The impact of diffusion coefficients on disease

We set diffusion intensity 1 to $d_1 = 0.6, d_2 = d_3 = 0.0002, d_4 = 0.5, d_5 = 0.3, d_6 = d_7 = 0.0001$, diffusion intensity 2 to $d_1 = 0.6, d_2 = d_3 = 0.1, d_4 = 0.58, d_5 = 0.3, d_6 = 0.18, d_7 = 0.1$, diffusion intensity 3 to $d_1 = d_2 = 0.6, d_3 = 0.55, d_4 = 0.58, d_5 = d_6 = 0.3, d_7 = 0.25$, $\sigma_1 = \sigma_2 = 0.01$, and the other parameter values are the same as those in Figure 2. The impacts of different diffusion strengths on infected people and mosquitoes are given in Figure 8, from which we find that with the increasing movement of infected people and mosquitoes, the number of infected people and mosquitoes also increases, indicating that controlling the movement of infected people and mosquitoes can reduce the risk of Zika disease.

6.4.   Optimal controls

Some numerical results of optimal controls will be presented in this subsection. We choose $c = 0.5$, $\alpha = 0.1$, $c_0 = 0.5$ in (5.1), $a_1 = 0.25$, $a_2 = 0.35$, $a_3 = 0.33$, $b_1 = 0.1$, $b_2 = 0.2$, $b_3 = 0.2$, $c_1 = 44$, $c_2 = 40$, $c_3 = 50$ in (5.2), and $\mu_i(t, x) = 2.0 (i = 1, 2, \cdots, 7)$ in (5.4). The remaining parameter values are consistent with those in Figure 2.

Figure 9 shows the space-time diagrams of optimal controls. Figure 10 is the time evolutions of optimal controls when space variable $x = 10$. From these two figures, we can see that the levels of human control (individual protection $u_1$ and medical treatment of the infected people $u_2$) are very high in the early and middle stages of the disease, but they are low in the later stages; however, the control level of mosquitoes ($u_3$) has maintained the maximum for a long time, which means that the strength of control for mosquitoes has exceeded that for humans. Figure 11 demonstrates the trajectories of $I_h(t, x)$ and $I_m(t, x)$ for $x = 10$ under the four conditions of no control, only controlling humans, only controlling mosquitoes, and controlling both humans and mosquitoes. We observe that the effects on the disease of implementing three control variables and no control are very significant. In addition, the control variables $u_1$ and $u_2$ have great influence on the changes of infected people, but have little influence on infected mosquitoes. However, $u_3$ has a great impact on both people and mosquitoes. Therefore, to sum up, reducing the number of mosquitoes is the primary factor to control Zika disease and personal protection and treatment of the infected humans are also two indispensable measures.

7.   Conclusions

This paper presents a stochastic Zika disease model with spatial diffusion, which includes human-mosquito transmission, human-human sex transmission, and vertical transmission of mosquitoes, and studies the dynamic behavior and optimal control of the model. Firstly, we give the conditions for almost surely exponential extinction of Zika disease, and the result signifies that the Zika disease will disappear when the diffusion coefficients of infected people and mosquitoes are very small and the fluctuations of environmental noise are relatively large. Secondly, we prove the sufficient conditions for the existence and uniqueness of the steady-state distribution representing the persistence of the disease, and research suggests that when the strengths of environmental noise are low and the diffusion coefficients of humans and mosquitoes are relatively large, Zika disease will continue to exist, which is contrary to the situation of disease extinction. In addition, numerical simulations have shown that increasing the intensity of random noise or decreasing the movement of infected people and mosquitoes can lessen the occurrence of Zika disease. Finally, we take three control variables, namely, individual protection, medical treatment of the infected people, and insecticides for spraying mosquitoes, into the model, and derive the expressions of optimal controls according to the Pontryagin maximum principle. Numerical simulations show that individual protection and treatment of infected persons are very effective for human beings, but reducing the number of mosquitoes is still the most important measure to control Zika.

The experiments demonstrate that the growth, survival, propagation, biting rate, transmission, and infection probability of Aedes aegypti and Aedes albopictus are closely related to the temperature, which is an essential factor affecting the dynamics of the spread of mosquito-borne diseases ^[50]. Thus, incorporating seasonality, establishing a stochastic periodic system, and studying the dynamics and control of Zika disease are our next exploration directions.

Use of AI tools declaration

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

Acknowledgments

We would like to thank the editor and the reviewers for their help, which greatly improve this paper. The research is funded by the National Natural Science Foundation of China under grant number 11971013.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix A.   The proof of Theorem 2.1

Proof. Because of the local Lipschitz continuity of the coefficients of system (2.2), for any given initial function $(S_h^0(x), E_h^0(x), I_h^0(x), R_h^0(x), S_m^0(x), E_m^0(x), I_m^0(x))\in\mathcal{H}^+$, there is a unique local solution $(S_h(t, x), E_h(t, x), I_h(t, x), R_h(t, x), S_m(t, x), E_m(t, x), I_m(t, x))$ in $t\in[0, \tau_e), x\in Q$, here $\tau_e$ is the moment of explosion. To validate the local solution is also global, we only need to prove that $\tau_e = \infty\; a.s.$ Let $k_0$ be large enough such that every component of $(|S_h^0(x)|_Q, |E_h^0(x)|_Q, |I_h^0(x)|_Q, |R_h^0(x)|_Q, |S_m^0(x)|_Q, |E_m^0(x)|_Q, |I_m^0(x)|_Q)$ is in the interval $(\frac{1}{k_0}, k_0]$. Then for every integer $k\geq k_0$, define a stopping time

We set $\inf\emptyset = \infty$ ($\emptyset$ denotes the empty set usually) throughout this paper. Apparently, $\tau_k$ is increasing constantly as $k\rightarrow \infty$. Let $\tau_\infty = \lim_{k\rightarrow \infty}\tau_k$, accordingly $\tau_\infty\leq\tau_e\; a.s.$ As long as we can verify that $\tau_\infty = \infty\; a.s.$, then $\tau_e = \infty\; a.s.$ This means that $(S_h(t, x), E_h(t, x), I_h(t, x), R_h(t, x), S_m(t, x), E_m(t, x), I_m(t, x))\in\mathcal{H}^+\; a.s.$ for all $t\geq0$. Before showing that $\tau_\infty = \infty\; a.s.$, let us prove the boundedness of solution for every $k$ when $t\in[0, \tau_k)$.

Let

then

here $|Q|$ stands for the volume of $Q$, $\nu = d_h\wedge d_m$.

Consider the following SDE

By the constant variation method, the solution of equation (A.1) can be obtained as

where $M(t) = \int_{0}^{t}e^{-\nu(t-s)}\int_{Q}\sigma_1(S_h(s, x)+E_h(s, x)+I_h(s, x)+R_h(s, x))dxdB_1(s) +\int_{0}^{t}e^{-\nu(t-s)}\int_{Q}\sigma_2(S_m(s, x)+E_m(s, x)+I_m(s, x))dxdB_2(s)$ is a continuous local martingale with $M(0) = 0\; a.s.$ Combining the stochastic comparison theorem, we can get that there is a constant $C_0 > 0$ such that $N(t)\leq Z(t)\leq C_0\; \; a.s.$ That is, for each $k$, when $t\in[0, \tau_k)$,

Next, we continue to prove that $\tau_\infty = \infty\; a.s.$ For any $T > 0$, define $V(t) = \langle S_h, S_h\rangle+\langle E_h, E_h\rangle+\langle I_h, I_h\rangle+\langle R_h, R_h\rangle+\langle S_m, S_m\rangle +\langle E_m, E_m\rangle+\langle I_m, I_m\rangle$, $t\in[0, \tau_k\wedge T)$. Using the $It\hat{o}$ formula, we have

Integrating the two ends of the above equation from $0$ to $\tau_k\wedge T$, taking the expectation, and using (A.2) and the fundamental inequality, yield

Thereby,

where

By taking advantage of the Gronwall inequality, we get

Denote $\varrho_k = \inf_{\|X(t, x)\|\geq k, 0 < t < \infty}V(t)$ for $k\geq k_0$. Obviously, $\varrho_k\rightarrow \infty$ $(k\rightarrow \infty)$. Combine (A.3) to get $C_1e^{C_2T}\geq \mathbb{E}V(\tau_k\wedge T) = \mathbb{E}[V(\tau_k)\chi_{\{\tau_k\leq T\}}]\geq\varrho_k \mathbb{P}(\tau_k\leq T)$. Setting $k\rightarrow \infty$, then $\mathbb{P}(\tau_k\leq T) = 0$. Thence $\mathbb{P}(\tau_\infty > T) = 1$. The proof is completed. □

Appendix B.   The proof of Theorem 2.2

Proof. Define $V(t) = \|S_h(t, x)\|^p+\|E_h(t, x)\|^p+\|I_h(t, x)\|^p+\|R_h(t, x)\|^p+\|S_m(t, x)\|^p+\|E_m(t, x)\|^p+\|I_m(t, x)\|^p$. First, consider $p\geq2$. Making use of the $It\hat{o}$ formula, we obtain

Integrating the two sides of the above equation and taking the supremum and expectation, we can get

Using the Young inequality as well as the Burkholder-Davis-Gundy inequality, we can further see that

After organizing, we have

where

The Gronwall inequality implies

Next suppose $0 < p < 2$. By the Hölder inequality, we have

here the second inequality sign makes use of the fundamental inequality $|x_1+x_2+\cdots+x_7|^r\leq7^{r-1}(|x_1|^r+|x_2|^r+\cdots+|x_7|^r), \; \forall r > 1$. This completes the proof. □

Appendix C.   The discrete form of model (2.2)

The discrete version of model (2.2) is as follows

where $\zeta_{1i}, \zeta_{2i}(i = 1, 2, \cdots)$ are independent standard normal variables.