Establishing Wolbachia in the wild mosquito population: The effects of wind and critical patch size

1.   Introduction

Mosquito-borne diseases such as malaria, dengue fever, West Nile virus, chikungunya, and Zika virus have been a great threat to public health. For dengue virus alone, it has been estimated that 3.9 billion people, in 128 countries, are at risk of infection ^[1]. In the year of 2014, more than 43,000 cases with locally acquired denguelike illness were reported in Guangdong province, China ^[2]. The human viruses including dengue, Zika, chikungunya and yellow fever are transmitted primarily by Aedes aegypti mosquitoes. Due to the lack of vaccines and efficient clinical cures ^[3], the only effective control strategy seems to be controlling the population of mosquitoes that transmit human viruses. Since massive spraying of insecticides and elimination of mosquito breeding sites are not sustainable to reduce mosquito density and might also lead to serious environmental problems, a promising strategy is the Wolbachia approach: releasing male and female Aedes aegypti mosquitoes with Wolbachia so that these mosquitoes can breed with the wild mosquito population and pass Wolbachia to the entire mosquito population. On one hand, the ability to transmit viruses to human for mosquitoes with Wolbachia is greatly reduced ^[4,5]. On the other hand, since the Wolbachia infection often induces cytoplasmic incompatibility (CI), which leads to early embryonic death when Wolbachia-infected males mate with uninfected females ^[6,7], the Wolbachia approach would greatly reduce the density of the mosquito population and can thus potentially eliminate the mosquito population and thus eradicate the mosquito-born infectious diseases.

To understand the Wolbachia infection dynamics, Zheng, Tang and Yu ^[8] proposed a delay differential equation model. To be self-contained, we briefly introduce their idea here. We denote by $r_{f}$ and $r_{m}$ the numbers of released female mosquitoes and released males carrying Wolbachia, respectively. Due to strong competition between adults, $r_{f}$ and $r_{m}$ satisfy

Here $T(t) = r_{f}+r_{m}+I_{f}+I_{m}+U_{f}+U_{m}$ denotes the total population size, with $U_{f},$ $U_{m}$, $I_{f}$ and $I_{m}$ being the numbers of uninfected reproductive females, uninfected reproductive males, and infected reproductive females and males other than those from releasing, respectively. Let $b_{I}$ (resp. $b_{U}$) be the natural birth rate of the infected (resp. uninfected) mosquitos and $0\leq\delta\leq 1$ be the proportion of mosquitos born female. Then the proportion of mosquitos born male is $1-\delta$. With strong CI and high maternal transmission, if the average waiting time from parent mating to the emergence of reproductive progenies for both infected and uninfected mosquitoes is negligible, then we have

Since both $r_{f}$ and $r_{m}$ approach $0$ as $t\rightarrow+\infty$. Let

Assuming equal determination case, which means that $\delta = {1}/{2}$, $I_{f} = I_{m}$ and $U_{f} = U_{m},$ setting $b_{1} = { b_{I}}/{2}$ and $b_{2} = { b_{U}}/{2}$ and considering the spatiotemporal factor, Huang et al. ^[9] came up with the following reaction-diffusion system:

In (1.4), $d_{1}$ and $d_{2}$ are the diffusion rates, $\Delta$ denotes the Laplace operator in the spatial variable $x$, and $\nu$ denotes the unit outward normal vector to the boundary of $\Omega$.

Noticing that the spread of mosquitoes can be greatly affected by the wind speed ^[10], therefore, the advection due to wind should be incorporated into modeling. Note also practical experience of releasing mosquitoes with Wolbachia to the wild suggests that a minimum release area is needed in order to achieve a stable local establishment and spread in continuous habitats ^[11]. Motivated by these two aspects, in this paper, we extend the model considered in ^[9], i.e., system (1.4), to a reaction-diffusion-advection model with spatially heterogeneous environment and flexible boundary conditions. More specifically, our model is described by the following system

where $u(x, t)$ and $v(x, t)$ represent the population densities of infected and uninfected mosquitoes at location $x$ and time $t$, respectively. The parameter $d$ denotes the random diffusion rate of $u$ and $v$. The functions $b_{1}(x)$ and $b_{2}(x)$ denote the halves of the birth rates of the infected and uninfected mosquitoes, respectively (^[8]). The parameter $\delta_{1}$ (or $\delta_{2}$) denotes the density dependent death rate for the infected (or uninfected) mosquito species. The advection constant $\alpha$ measures the result of wind transportation. The parameter $\beta < \infty$ measures the relative rate of population loss at the downstream due to wind flow and replenishment at the upstream.

We suppose that $b_{i}(x)$ satisfies the following hypothesis:

Here $\underline{b}$ and $\bar{b}$ are two positive constants.

We regard $\beta = +\infty$ as the Dirichlet boundary case, and system (1.5) becomes

Ideally, if the entire mosquito population is replaced by mosquitoes with Wolbachia, then the Wolbachia establishment is called successful. This is achieved if the solutions of (1.5) and (1.7) approach a semi-trivial steady state, $(\widetilde{u}, 0)$, where $\widetilde{u}$ satisfies the following equations,

for $\beta < \infty$ or

for the case with $\beta = \infty$.

We point out that there has been several mathematical models formulated to describe the Wolbachia spreading dynamics. These models include differential equations with/without time delays ^[8,12,13,14], reaction-diffusion equations ^[9,15], and stochastic equations ^[16]. These models focused on studying the subtle relation between the threshold releasing level for Wolbachia-infected mosquitoes and several important parameters including the CI intensity and the fecundity cost of Wolbachia infection. We should also point out that besides releasing mosquitoes with Wolbachia, an alternative control strategy is releasing the sterile mosquitoes to the wild mosquito population ^[17].

We organize the rest of this paper as follows. In Section 2, we give some useful lemmas to establish the relation between two related principal eigenvalues and the domain size, the advection rate as well as the diffusion rate. Our main results are presented in Section 3. To illustrate our results, we also present some numerical simulations in Section 4. We conclude this paper in the last section.

2.   Useful lemmas

In order to discuss the existence of $(\widetilde{u}, 0)$, corresponding to systems (1.5) and (1.7), we need to study the following subsystems:

and

respectively.

This leads to the study of the following linear eigenvalue problem

In the case that $\beta = +\infty$, the corresponding eigenvalue problem reads as

Throughout this paper, we denote the principal eigenvalue of (2.3) or (2.4) by $\sigma_{1}(\alpha, d, L)$.

2.1.   Dependence of $\sigma_{1}(\alpha, d, L)$ on $\alpha$, $d$ and $L$

First, we shall investigate how $\sigma_{1}(\alpha, d, L)$ depends on $\alpha$.

Lemma 2.1. Suppose that (1.6) is satisfied, then

(a) When $\beta\in(0, \frac{1}{2})$, if $b'_1(x)\leq 0$, then $\sigma_{1}(\alpha, d, L)$ is a strictly monotone decreasing function of $\alpha$; When $\beta\in[\frac{1}{2}, +\infty]$, then $\sigma_{1}(\alpha, d, L)$ is a strictly monotone decreasing function of $\alpha$.

(b) $0 < \underline{b}\leq\lim_{\alpha\rightarrow0}\sigma_{1}(\alpha, d, L)\leq\bar{b}$ provided that $\beta < +\infty$.

(c) $\lim_{\alpha\rightarrow+\infty}\sigma_{1}(\alpha, d, L) = -\infty$.

Proof. Set $\varphi = e^{-\frac{\alpha}{d}x}\phi$, then (2.3) becomes

and (2.4) becomes

Let us denote by $\varphi_{\alpha}$ the derivative of $\varphi$ with respect to $\alpha$. Multiplying the first equation of (2.5) by $e^{\frac{\alpha}{d}x}\varphi_{\alpha}$, we obtain

Integrating (2.7) over $(0, L)$ yields

On the other hand, it follows from (2.5) that

Multiplying (2.9) by $e^{\frac{\alpha}{d}x}\varphi(x)$ and integrating it over $(0, L)$, we have

Using integration by parts and (2.8), we obtain

If $\beta = +\infty$, (2.11) becomes

In general, when $\beta\in(0, \frac{1}{2})$, $\sigma_{1}(\alpha, d, L)$ may not be monotone in $\alpha$. To ensure $\sigma_{\alpha} < 0$, we need to show $\varphi_{x} < 0$ in $(0, L)$. To this end, we set $W = \frac{\phi_{x}}{\phi}$, and we can rewrite (2.3) as the following equation

Since $b'_1(x)\leq 0$, then by the maximum principle we have $W < \frac{\alpha}{d}$, which implies that $\varphi_{x} < 0$ in $(0, L)$ (^[18]). This implies that $\sigma_{1}(\alpha, d, L)$ is a strictly monotone decreasing function of $\alpha$.

For the case $\beta\in[\frac{1}{2}, +\infty)$, from (2.11), we have $\sigma_{\alpha} = \frac{-e^{\frac{\alpha}{d}L}\beta\varphi^{2}(L)+\int^{L}_{0}e^{\frac{\alpha}{d}x} \varphi_{x}\varphi}{\int^{L}_{0}\varphi^{2}e^{\frac{\alpha}{d}x}} = \frac{-e^{\frac{\alpha}{d}L}\beta\varphi^{2}(L)+\int^{L}_{0}e^{\frac{\alpha}{d}x} d\frac{\varphi^{2}}{2}}{\int^{L}_{0}\varphi^{2}e^{\frac{\alpha}{d}x}} = \frac{-e^{\frac{\alpha}{d}L}\beta\varphi^{2}(L)+e^{\frac{\alpha}{d}x}\frac{\varphi^{2}}{2}\mid^L_{0}-\frac{\alpha}{d}\int^{L}_{0}e^{\frac{\alpha}{d}x} \frac{\varphi^{2}}{2}}{\int^{L}_{0}\varphi^{2}e^{\frac{\alpha}{d}x}} = \frac{-(\beta-\frac{1}{2})e^{\frac{\alpha}{d}L}\varphi^{2}(L)-\frac{\varphi^{2}(0)}{2}-\frac{\alpha}{d}\int^{L}_{0}e^{\frac{\alpha}{d}x} \frac{\varphi^{2}}{2}}{\int^{L}_{0}\varphi^{2}e^{\frac{\alpha}{d}x}}$.

It is easy to see $\beta\in[\frac{1}{2}, +\infty)$ implies that $\sigma_{\alpha} < 0$. When $\beta = +\infty$, Using the same calculation as above and combining with (2.12), we can obtain that $\sigma_{\alpha} < 0$.

If $\alpha = 0$, then the boundary condition of (2.5) is simply a Neumann boundary condition. Furthermore, combining with (1.6), we have

Especially, $\lim_{\alpha\rightarrow0}\sigma_{1}(\alpha, d, L) = b_{1}$, when $b_{1}(x) = b_{1}$. It follows from Proposition 2.1 of ^[19] and (1.6) that

where $0 < \gamma < \min\{1, \beta\}$. Consequently, $\lim_{\alpha\rightarrow+\infty}\sigma_{1}(\alpha, d, L) = -\infty$.

Next we investigate the relationship between $\sigma_{1}(\alpha, d, L)$ and $L$.

Lemma 2.2. $\sigma_{1}(\alpha, d, L)$ is a strictly monotone increasing function of $L$.

Proof. First we consider system (2.5) by fixing the parameters $\alpha$ and $d$ and varying $L$. For any $0 < L_{1} < L_{2}$, we show that $\sigma_{1}(\alpha, d, L_{1}) < \sigma_{1}(\alpha, d, L_{2})$. Let $\varphi_{1}$ be the positive eigenfunction corresponding to $\sigma_{1}(\alpha, d, L_{1})$, and define

Since

Due to the strict positivity of the eigenfunction corresponding to $\sigma_{1}(\alpha, d, L_{2})$ in $[0, L_{2}]$, the above inequality should be strict, which implies $\sigma_{1}(\alpha, d, L_{1}) < \sigma_{1}(\alpha, d, L_{2})$. The case with $\beta = +\infty$ can be dealt with in a similar fashion and we skip the details here.

Lemma 2.3. Assume (1.6) holds. If $\beta = +\infty$, then $\sigma_{1}(\alpha, d, L)\to-\infty$ as $L\to 0$.

Proof. We rewrite $\sigma_{1}(\alpha, d, L)$ as $\sigma_{1}(\alpha, d, L, b_{1}(x))$ to emphasize the dependence of $\sigma_{1}(\alpha, d, L)$ on the function $b_{1}(x)$. By Corollary 2.2 of ^[20], one easily sees that $\sigma_{1}(\alpha, d, L, b_{1}(x))\leq\sigma_{1}(\alpha, d, L, \bar{b})$.

Next we consider the following boundary value problem (BVP)

Set $\psi = e^{-\frac{\alpha x}{2d}}\phi$, then (2.14) yields

From the boundary conditions we find

Since by Lemma 2.2, we have

Assume that $\sigma^{*}$ is finite, then from (2.16), we obtain $\sigma^{*} = \frac{4d\bar{b}-\alpha^{2}}{4d}$, and hence $\psi(x) = Ax+B$. The boundary conditions imply that $A = B = 0$. This is impossible. Thus $\lim_{L\rightarrow0^{+}}\sigma_{1}(\bar{b}) = -\infty$ and hence $\lim_{L\rightarrow0^{+}}\sigma_{1}(b_{1}(x)) = -\infty.$

Lemma 2.4. If $\beta = +\infty$ and $\alpha = 0$, then $\sigma_{1}(\alpha, d, L)$ is a strictly monotone decreasing function of $d$.

Proof. Let us denote by $\varphi_{d}$ the derivative of $\varphi$ with respect to $d$. Multiplying the first equation of (2.6) by $e^{\frac{\alpha}{d}x}\varphi_{d}$, it holds that

Integrating (2.17) over $(0, L)$, we obtain

On the other hand, it follows from (2.6) that

Multiplying (2.19) by $e^{\frac{\alpha}{d}x}\varphi(x)$ and integrating it over $(0, L)$, we have

Using integration by parts and (2.18), we obtain

We find $\sigma_{d} < 0$, when $\alpha = 0$.

If $\beta = +\infty$, we also have the following remark (See Theorem 3.1 of ^[21]).

Remark 2.5. (i) $\sigma_{1}(0, d, L)\rightarrow \max_{x\in[0, L]}b_1(x) > 0$ as $d\rightarrow0$;

(ii) $\sigma_{1}(0, d, L)\rightarrow -\infty$ as $d\rightarrow+\infty$.

3.   Main results

Our first result is the following theorem.

Theorem 3.1. If system (1.5) admits a semi-trivial steady state $(\widetilde{u}, 0)$, then the semi-trivial steady state $(\widetilde{u}, 0)$ is locally asymptotically stable.

Proof. Linearizing the second equation of system (1.5) at $(\widetilde{u}, 0)$, we obtain the following eigenvalue problem

Let $\widetilde{\phi} = e^{-\frac{\alpha}{d}x}\phi$, then problem (3.1) becomes

By the variational method, $\zeta_{1}$ can be characterized by

and when $\beta = +\infty$,

where $\widetilde{\phi_{1}}$ is the eigenfunction associated with $\zeta_{1}$. Thus $\zeta_{1} < 0$ and it follows from Proposition 3.1 of ^[20] that the semi-trivial steady state $(\widetilde{u}, 0)$ is locally asymptotically stable.

Theorem 3.2. If $\alpha = 0$, or $\alpha > 0$ and $\beta = 0$, then system (1.5) admits a semi-trivial steady state, $(\widetilde{u}, 0)$, which is locally asymptotically stable, where $\widetilde{u}$ is the unique positive steady state of problem (2.1).

Proof. Under the assumptions, the linear eigenvalue problem (2.3) reads as

and

respectively. For (3.5), a standard eigenvalue analysis gives the principle eigenvalue

Moreover, it is easy to verify that $\lambda_{1} > 0$ by using the test function $\phi = 1$. For (3.6), we multiply the equations by $e^{-\frac{\alpha}{d}x}$ to get

Thus the associated principle eigenvalue is

Using the test function $\phi\equiv e^{\frac{\alpha}{d}x}$, we see that $\eta_{1} > 0$. The conclusion then follows from Propositions 3.2 and 3.3 of ^[20] and Theorem 3.1.

For the spatially homogeneous case, by Lemmas 2.1 and 2.2, Theorems 2.1 and 2.3 of ^[19], we have the following result.

Theorem 3.3. Assume that $b_{1}(x)\equiv b_{1} > 0$ is a constant.

(a) If $\beta\in(0, \frac{1}{2})$ and $0 < \alpha < \sqrt{\frac{db_{1}}{\beta(1-\beta)}}$, then system (1.5) admits a semi-trivial steady state, $(\widetilde{u}, 0)$, with $\widetilde{u}$ being the unique positive steady state of problem (2.1), if and only if $L > L_{1}^{*}$, where

Moreover, if $\alpha\geq\sqrt{\frac{db_{1}}{\beta(1-\beta)}}$, then problem (2.1) only has a globally asymptotically stable zero steady state and system (1.5) does not admit a semi-trivial steady state in the form of $(\widetilde{u}, 0)$ satisfying $\widetilde{u} > 0$.

(b) When $\beta = \frac{1}{2}$. If $0 < \alpha < \sqrt{4db_{1}}$, then system (1.5) admits a semi-trivial steady state, $(\widetilde{u}, 0)$, with $\widetilde{u}$ being the unique positive steady state of problem (2.1), if and only if $L > L_{1}^{*}$ with

If $\alpha\geq\sqrt{4db_{1}}$, then system (1.5) does not admit a semi-trivial steady state in the form of $(\widetilde{u}, 0)$ satisfying $\widetilde{u} > 0$.

(c) Suppose $\beta\in(\frac{1}{2}, +\infty)$. If $0 < \alpha < \sqrt{4db_{1}}$, then system (1.5) admits a semi-trivial steady state, $(\widetilde{u}, 0)$, with $\widetilde{u}$ being the unique positive steady state of problem (2.1), if and only if $L > L_{1}^{*}$ with

If $\alpha\geq\sqrt{4db_{1}}$, then system (1.5) does not admit a semi-trivial steady state in the form of $(\widetilde{u}, 0)$ satisfying $\widetilde{u} > 0$.

(d) Suppose that $\beta\in(0, \frac{1}{2})$. If $d > \frac{\alpha^{2}\beta(1-\beta)}{b_{1}}$, then system (1.5) admits a semi-trivial steady state, $(\widetilde{u}, 0)$, with $\widetilde{u}$ being the unique positive steady state of problem (2.1), if and only if $L > L_{1}^{*}$ with

(e) For $\beta = \frac{1}{2}$, if $d > \frac{\alpha^{2}}{4b_{1}}$, then system (1.5) admits a semi-trivial steady state, $(\widetilde{u}, 0)$, with $\widetilde{u}$ being the unique positive steady state of problem (2.1), if and only if $L > L_{1}^{*} = \frac{2d\arctan\frac{\alpha\beta\sqrt{4db_{1}-\alpha^{2}}}{2db_{1}-\alpha^{2}\beta}}{\sqrt{4db_{1}-\alpha^{2}}}.$

(f) Suppose that $\beta\in(\frac{1}{2}, +\infty)$, when $d > \frac{\alpha^{2}}{4b_{1}}$, then system (1.5) admits a semi-trivial steady state, $(\widetilde{u}, 0)$, with $\widetilde{u}$ being the unique positive steady state of problem (2.1), if and only if $L > L_{1}^{*}$, where $L_{1}^{*}$ is given by

(g) For any $\alpha > 0$, $\beta\in(0, 1)$, if $d\leq\tilde{d}$, then system (1.5) does not admit a semi-trivial steady state in the form of $(\widetilde{u}, 0)$ satisfying $\widetilde{u} > 0$, where

Similarly, for system (1.7), we have the following result.

Theorem 3.4. Suppose $b_{1}(x)\equiv b_{1} > 0$ is a constant. Then we have the following conclusions.

(i) For any $d > 0$, if $0 < \alpha < \sqrt{4b_{1}d}$, then there exists a critical number

such that for $L\geq \widehat{L}_{1}^{*}$, system (1.7) admits a stable semi-trivial steady state $(\widetilde{u}, 0)$, where $\widetilde{u} > 0$ is the unique positive steady state of system (2.2). Moreover, $\widehat{L}_{1}^{*}$ is an increasing function of the advection rate $\alpha$ for $0 < \alpha < \sqrt{4b_{1}d}$.

(ii) For any $\alpha$, if $d > \frac{\alpha^{2}}{4b_{1}}$, then there exists a critical number

such that for $L\geq \widehat{L}_{1}^{*}$, system (1.7) admits a stable semi-trivial steady state $(\widetilde{u}, 0)$, where $\widetilde{u} > 0$ is the unique positive steady state of system (2.2). Further, $\widehat{L}_{1}^{*}$ is an decreasing function of the diffusion rate $d$, for $\frac{\alpha^{2}}{4b_{1}} < d < \frac{\alpha^{2}}{2b_{1}}$.

(iv) If $\alpha\geq\sqrt{4b_{1}d}$, then system (1.7) does not admit a semi-trivial steady state in the form of $(\widetilde{u}, 0)$, where $\widetilde{u} > 0$.

If $b_1(x)$ is not a constant, we have the following result.

Theorem 3.5. Suppose that (1.6) is satisfied. Then the following conclusions hold.

(I) If $\alpha > 0$, when $\beta\in(0, \frac{1}{2})$, $b'_{1}(x)\leq0$, or $\beta\in[\frac{1}{2}, +\infty)$, then there exists $\alpha^{*} > 0$ such that for $0 < \alpha < \alpha^{*}$, system (1.5) admits a semi-trivial steady state, $(\widetilde{u}, 0)$, which is stable, and $\widetilde{u}$ is the unique positive steady state of problem (2.1).

(II) For the case with $\beta = \infty$, we have the following conclusions.

(II.1) If $\alpha = 0$, then there exists $d^{*}$ such that for $0 < d < d^{*}$, system (1.7) admits a semi-trivial steady state $(\widetilde{u}, 0)$, which is stable and $\widetilde{u}$ is the unique positive steady state of system (2.2).

(II.2) If $L > \frac{\pi}{2}\sqrt{\frac{d}{\underline{b}}}$, then there exists $\tilde{\alpha}^{*}$ such that for $0 < \alpha < \tilde{\alpha}^{*}$, system (1.7) admits a semi-trivial steady state $(\widetilde{u}, 0)$, which is stable and $\widetilde{u}$ is the unique positive steady state of system (2.2).

(II.3) If $0 < \alpha < \sqrt{4\underline{b}d}$, then there exists a $L^{*} > 0$ such that for $L > L^{*}$, system (1.7) admits a semi-trivial steady state $(\widetilde{u}, 0)$, which is stable and $\widetilde{u}$ is the unique positive steady state of system (2.2).

Proof. The case (I) can be proved by (a), (b) and (c) of Lemma 2.1 and the proof of case (II.1) can be obtained by Lemma 2.4 and Remark 2.5. To prove case (II.2), we just need to verify the following three facts.

$(1)$ $\sigma_{1}(\alpha, b_{1}(x))$ is a strictly monotone decreasing function of $\alpha$;

$(2)$ $\lim_{\alpha\rightarrow0}\sigma_{1}(\alpha, b_{1}(x)) > 0$;

$(3)$ $\lim_{\alpha\rightarrow+\infty}\sigma_{1}(\alpha, b_{1}(x)) = -\infty$.

Note that (1) and (3) follow immediately from Lemma 2.1. Furthermore, (2) follows from noting that $\sigma_{1}(0, b_{1}(x))\geq\sigma_{1}(0, \underline{b})$ for $b_{1}(x)\geq\underline{b}$ and $\sigma_{1}(0, b_1(x))\geq\sigma_{1}(0, \underline{b}) = \underline{b}-\frac{d\pi^{2}}{4L^{2}} > 0$, when $\beta = +\infty$ and $L > \frac{\pi}{2}\sqrt{\frac{d}{\underline{b}}}$.

Next we prove (II.3). Note that $\beta = \infty$, the associated eigenvalue problem is

where $\phi$ is corresponding eigenfunction of $\sigma_{1}$.

If $0 < \alpha < \sqrt{4\underline{b}d}$, then it follows from Theorem 3.4 that $\sigma_{1}(\underline{b}) > 0$, as $L\rightarrow +\infty$. Hence, it is easy to know that $\sigma_{1}(b_{1}(x)) > 0$, as $L\rightarrow +\infty$, if $0 < \alpha < \sqrt{4\underline{b}d}$. Therefore the proof of case (II.3) follows directly from Lemma 2.2 and Lemma 2.3.

4.   Numerical simulations

In this section, we use several numerical simulations to illustrate our theoretical results. First we verify Theorem 3.3. We take parameter values: $\beta = 0.4$, $b_{1}(x) = 0.1077$, $b_{2}(x) = 0.1988$, $\delta_{1} = \delta_2 = 8.5034\times10^{-6}$, $d = 1.25\times10^{-2}$, $\alpha = 1.25\times10^{-2}$. For this set of parameters, we find that $\alpha < \sqrt{\frac{db_{1}}{\beta(1-\beta)}}\approx 0.0749$ and $L_{1}^{*}\approx 0.0472$. Theorem 3.3 applies here: if we take $L = 0.04 < L_1^*$, then as shown in Figure 1 (Left), the solution of (1.5) approaches the trivial steady state $(0, 0)$, while if we take $L = 0.05 > L_1^*$, as shown in Figure 1 (Right), the solution of (1.5) approaches the semi-trivial steady state $(\widetilde{u}, 0)$ implying that all mosquitoes are infected with Wolbachia.

Now we take parameter values as: $d = 1.25\times10^{-2}$, $\delta_{1} = \delta_2 = 8.5034\times10^{-6}$, $b_{1}(x) = e^{-x}+1$, $b_2(x) = 1.5+\sin (x)$, $\beta = 0.4$, $L = 0.1$. For this set of parameters, we can numerically find $\alpha*\approx 0.261$ such that the solutions of (1.5) approach $(\widetilde{u}, 0)$ if $\alpha\in [0, \alpha^*)$, and approach $(0, 0)$ if $\alpha\geq \alpha^*$ (See Figure 2).

Now we take $\beta = \infty$, $\alpha = 0.1789$, $d = 1.25\times10^{-2}$, $\delta_1 = \delta_2 = 0.85034\times 10^{-6}$, $b_1(x) = e^{-x}+1$, $b_{2}(x) = e^{-x}+1$. Then $\underline{b} = 1$, $\alpha < \sqrt{4d\underline{b}} = 0.2236$. Numerically we find $\hat{L}^{*}\approx 0.21$. As shown in Figure 3, if $L > L_1^*$, then the solutions of (1.5) approach the semi-trivial steady state $(\widetilde{u}, 0)$ implying a full establishment Wolbachia is achieved, while the establishment of Wolbachia fails if $L < L_1^*$.

5.   Conclusions

In this paper, we have proposed a reaction-diffusion-advection model in one-dimensional spatially inhomogeneous environment with general boundary conditions. Our results (Theorems 3.2–3.5) show that in order to fully establish Wolbachia in the wild mosquito population, i.e., all mosquitoes eventually carry Wolbachia, the wind related parameter $\alpha$ and the patch size $L$ over which the mosquitoes with Wolbachia are spreading should satisfy certain requirements. Generally speaking, the wind cannot be too strong, and the minimum patch size cannot be too small. The critical values of the advection parameter $\alpha$ and the patch size $L$ for the establishment of Wolbachia depend on the model parameters including the diffusion rate $d$, the birth rates $b_1(x)$ and $b_2(x)$, and the death rates $\delta_1$ and $\delta_2$. For instance, Theorem 3.3 (a) indicates, if $b_1(x) = b_1$ is a positive constant and $\beta\in (0, 1/2)$, then the full establishment of Wolbachia will not be successful if $\alpha$ is too large such that $\alpha\geq\sqrt{\frac{db_{1}}{\beta(1-\beta)}}$ (i.e., the wind is too strong) or if $\alpha$ is in a proper range ($0 < \alpha < \sqrt{\frac{db_{1}}{\beta(1-\beta)}}$), but the patch size $L$ is too small satisfying $L < L_1^*$ (See Figure 1).

Acknowledgements

The authors wish to thank three anonymous reviewers for their very helpful comments. YL and ZG were supported by the National Science Foundation of China (No. 11771104), Program for Chang Jiang Scholars and Innovative Research Team in University (IRT-16R16). LW was partially supported by a Discovery grant from the Natural Sciences and Engineering Research Council of Canada.

Conflict of interest

The authors declare that they have no competing interests.