1. Introduction

Infectious diseases have been attracting considerable attention in recent years, and various epidemic models have been proposed and analyzed for prevention and control strategies, especially for vector borne diseases [15,39]. For example, a West Nile virus (WNv) is an arbovirus of the Flavivirus kind in the family Flaviviridae that causes the epidemics of febrile illness and sporadic encephalitis [7]. WNv is found in temperate and tropical regions of the world, it was first isolated and identified from the blood of a febrile Ugandan woman during research on yellow fever virus in 1937 [3].

Although WNv is widely distributed in Africa, the Middle East, Asia and southern Europe, in North America, the first infected case was detected in 1999 during an outbreak of encephalitis in New York city [3,21,26,39]. Since 1999 this virus has spread spatially and prevail in much of North America [8,21], it is evident that the spread of WNv comes from the interplay of disease dynamics and bird and mosquito movement.

To the best of our knowledge, currently there are no effective vaccine or medicine for WNv. To reduce the rates of WNv infection, anti-WNv efforts are primarily based on personal protective measures like insect repellent and protective clothing, and public heath measures [4].

Many mathematical models for WNv have been proposed and analyzed, however most of the models are focused on the non-spatial transmission dynamics [4,36,39]. In fact, the spatial spreading is an important factor to affect the persistence and eradication of WNv. In 2006, Lewis et al. [21] investigated the spatial spread of WNv to describe the movement of birds and mosquitoes. The reaction-diffusion model was extended from the non-spatial model for cross infection between birds and mosquitoes that was proposed and developed by Wonham et al. in [39].

To utilize the cooperative characteristic of cross-infection dynamics and estimate the spatial spread rate of infection, Lewis et al. in [21] proposed the following simplified WNv model

$\left\{ \begin{array}{ll} \frac{\partial I_b}{\partial t} = D_1\Delta I_b +\alpha_{b}\beta_{b}\frac{(N_b-I_b)}{N_b} I_m -\gamma_{b}I_b ,&(x,t)\in \Omega \times (0,+\infty), \\ \frac{\partial I_m }{\partial t} = D_2\Delta I_m +\alpha_{m}\beta_b \frac{(A_m-I_m)}{N_b} I_b -d_{m} I_m ,&(x,t)\in \Omega \times (0,+\infty),\\ \end{array} \right. \label{Aa1}$

(1)

where the positive constants $N_b$ and $A_m$ denote the total population of birds and adult mosquitos; $I_b(x,t)$ and $I_m(x,t)$ represent the populations of infected birds and mosquitos at the location $x$ in the habitat $\Omega\subset {\mathbb{R}}^N$ and at time $t\geq 0$, respectively, and $I_b(x,0)+I_m(x,0)>0$. The parameters in the above system are defined as follows:

● $\alpha_{m}$, $\alpha_{b}$ : WNv transmission probability per bite to mosquitoes and birds, respectively;

● $\beta_{b}$ : biting rate of mosquitoes on birds;

● $d_m$ : adult mosquitos death rate;

● $\gamma_b$ : bird recovery rate from WNv.

Here, the positive constants $D_1$ and $D_2$ are diffusion coefficients for birds and mosquitoes, respectively.

Considering the spatially-independent model

$\left\{ \begin{array}{ll} \frac{d I_b(t)}{d t} = -\gamma_{b}I_b(t)+\alpha_{b}\beta_{b} \frac{(N_b-I_b (t))}{N_b}I_m(t),&t>0, \\ \frac{d I_m(t)}{d t} = -d_{m}I_m(t)+\alpha_{m}\beta_{b} \frac{(A_m-I_m (t))}{N_b} I_b (t),&t>0, \end{array} \right. \label{aode1}$

(2)

one can see that if $R_0 (: = \sqrt{\frac{\alpha_{m}\alpha_{b}\beta_b ^2 A_m}{d_m\gamma_b N_b}}\ )< 1$, the virus will vanish eventually, while for $R_0 > 1$, a nontrivial epidemic level appears and is globally asymptotically stable in the positive quadrant [21].

For the diffusive model (1), Lewis et al. proved the existence of traveling wave and calculated the spatial spread rate of infection [21]. The corresponding free boundary problem describing the expanding process has been discussed in [34]. It is worth mentioning that WNv usually spreads from one area to another because of the diffusions of birds and mosquitoes, so that its transmission is affected not only by the characteristics of pathogens, but also by the spatial difference of environment in which birds or mosquitoes reside. Considering the complexity of diffusions and the heterogeneity of environment, the system (1) can be extended to the following strongly-coupled parabolic system

$\left\{ \begin{array}{ll} \frac{\partial I_b}{\partial t}-\Delta [(d_1+\alpha_{1}I_b+\frac{\beta_{1} I_b}{\gamma_1+I_m})I_b] = \\ \;\;\;\;\;\;\;\;\;\;\; \alpha_{b}(x)\beta_{b}(x)\frac{(N_b-I_b)}{N_b} I_{m} -\gamma_{b}(x)I_b ,&(x,t)\in \Omega \times (0,+\infty),\\ \frac{\partial I_m}{\partial t}-\Delta [(d_2+\frac{\beta_{2} I_m}{\gamma_2+I_b}+\alpha_{2}I_m)I_m] = \\ \;\;\;\;\;\;\;\;\;\;\; \alpha_{m}(x)\beta_b (x)\frac{(A_m-I_m)}{N_b} I_{b} -d_{m}(x) I_m ,&(x,t)\in \Omega \times (0,+\infty),\\ I_b(x,t) = I_m(x,t) = 0,&(x,t)\in \partial\Omega \times (0,+\infty),\\ 0\leq I_b(x,0)\leq N_b, \;\;\;0\leq I_m(x,0)\leq A_m,&x\in \overline\Omega, \end{array} \right. \label{a3b}$

(3)

and the corresponding elliptic problem with Dirichlet boundary conditions becomes

$\left\{ \begin{array}{ll} -\Delta [(d_1+\alpha_{1}I_b+\frac{\beta_{1} I_b}{\gamma_1+I_m})I_b] = f_1(x,I_b,I_m),& x\in \Omega , \\ -\Delta [(d_2+\frac{\beta_{2} I_m}{\gamma_2+I_b}+\alpha_{2}I_m)I_m] = f_2(x,I_b,I_m),&x\in \Omega,\\ I_b(x) = I_m(x) = 0,&x\in \partial\partial\Omega, \end{array} \right. \label{a3}$

(4)

where

$\left. \begin{array}{ll} f_1(x,I_b,I_m) = \alpha_{b}(x)\beta_{b}(x)\frac{(N_b-I_b)}{N_b} I_m -\gamma_{b}(x)I_b,\\ f_2(x,I_b,I_m) = \alpha_{m}(x)\beta_b(x) \frac{(A_m-I_m)}{N_b} I_b -d_{m}(x) I_m, \end{array} \right.$

and the parameters $\alpha_{b}(x)$, $\beta_{b}(x)$, $\gamma_b(x)$, $\alpha_{m}(x)$ and $d_{m}(x)$ are all sufficiently smooth and strictly positive functions defined on $\overline\Omega$. $d_{i}$ $(i = 1,2)$ is positive constants represent the free-diffusion coefficients of population $I_b(x)$ and $I_m(x)$, respectively. $\alpha_{i}$, $\beta_{i}$ and $\gamma_i$ $(i = 1,2)$ are nonnegative constants; $\alpha_{i}$, $\beta_{i}$ and $\gamma_i$ are the self-diffusion coefficients, the cross-diffusion rates and cross-diffusion pressures, respectively. The homogeneous Dirichlet boundary condition in (4) means that there is no infection on the boundary and outside of the domain $\Omega$. More specifically, the diffusion terms can be written as

$\textrm{div} \{(d_1+2\alpha_{1}I_b+\frac{2\beta_{1}I_b}{\gamma_1+I_m})\nabla I_b+\frac{-\beta_{1} I^2_{b}}{(\gamma_1+I_m)^2}\nabla I_m\},$

$\textrm{div} \{\frac{-\beta_{2} I^2_{m}}{(\gamma_2+I_b)^2}\nabla I_b+ (d_2+\frac{2\beta_{2}I_m}{\gamma_{2}+I_b}+2\alpha_{2}I_m)\nabla I_m\}.$

The terms

$d_1+2\alpha_{1}I_b+\frac{2\beta_{1}I_b}{\gamma_1+I_m},\,\, d_2+\frac{2\beta_{2}I_m}{\gamma_2+I_b}+2\alpha_{2}I_m$

represent the self-diffusions and the terms

$\frac{-\beta_{1} I^2_{b}}{(\gamma_1+I_m)^2},\;\; \frac{-\beta_{2} I^2_{m}}{(\gamma_2+I_b)^2}$

represent the cross-diffusions. Here the term $\frac{-\beta_{1} I^2_{b}}{(\gamma_1+I_m)^2}\nabla I_m$ is due to the moving of the population of birds $I_b$ toward the location of increasing population of mosquitoes $I_m$ with rate $(-\frac{\beta_{1} I^2_{b}}{(\gamma_1+I_m)^2}<0)$. Similarly, the term $\frac{-\beta_{2} I^2_{m}}{(\gamma_2+I_b)^2}\nabla I_b$ is due to the moving of the population of mosquitoes $I_m$ toward the location of increasing population of birds $I_b$ with rate $(-\frac{\beta_{2} I^2_{m}}{(\gamma_2+I_b)^2}<0)$. The solution $(I_b(x),I_m(x))$ to problem $(4)$ is called a coexistence if $I_b(x)>0$ and $I_m(x)>0$ for every $x \in \Omega$.

Problem (1) is weakly-coupled parabolic system which only consider the random diffusion in a homogeneous environment. However, problem (4) implies that, in addition to the dispersive force, the diffusion also depends on population pressure from other population. This means that the population in (4) are not homogeneously distributed due to the consideration of self and cross diffusion terms. Moreover, the diffusive behavior in different populations also affect the distribution of resources. Thus, the consideration of diffusion and cross-diffusion effect is very reasonable and more close to reality, see for example, [29] for mixed-culture biofilm model, [17] for the tumor-growth model and [14,16,31,40] for the competition model. There are some valuable results about the roles of diffusion and cross-diffusion in the modeling of the dynamics of strongly coupled reaction-diffusion systems [5,11,13,16,19,22,25,31,30,38,40]. For instance, Shigesada et al. [31] proposed the strongly coupled elliptic system describing two species Lotka-Volterra competition model. Ko and Ryu studied a predator-prey system with cross-diffusion, representing the tendency of prey to keep away from its predators, under the homogeneous Dirichlet boundary conditions in [19]. Fu et al. investigated the global behavior of solutions for a Lotka-Volterra predator-prey system with prey-stage structure, under the homogeneous Neumann boundary conditions [11]. In 2014, Jia et al. [16] discussed a Lotka-Volterra competition reaction-diffusion system with nonlinear diffusion effects. In 2016, Braverman and Kamrujjaman [5] introduced a competitive-cooperative models with various diffusion strategies. More recently, Li et al. studied an effect of cross-diffusion on the stationary problem of a Leslie prey-predator model with a protection zone [22].

In recent years, researches on the existence and non-existence of the positive solutions for the dynamics of strongly-coupled elliptic systems have received comprehensive attention [16,19,20,28]. There are many standard approaches to derive the coexistence for the standard semi-linear parabolic system in mathematical models, such as construction of upper and lower solutions [12,16,18,28], bifurcation theory [6], fixed point theorem [19,42], ect. The upper and lower solutions method developed by Pao [27] is concise and effective to derive the coexistence. Based on the method, the coexistence for a general strongly-coupled system has been given in [28]. In [18], Kim and Lin studied the coexistence of three species in strongly coupled elliptic system. Gan and Lin in [12] considered the competitor-competitor-mutualist three species Lotka-Volterra model. Recently, Jia et al. [16] investigated the existence of the positive steady state solution of a Lotka-Volterra competition model with cross-diffusion.

Motivated by above problems, in this paper we are more interested in the nonnegative steady state solutions, that is, the coexistence of problem (4) describing a cross-diffusive WNv model in a heterogenous environment.

The plan of this paper is as follows: Section 2 is devoted to the basic reproduction number of problem (4) and its properties. The existence and non-existence of coexistence to (4) are discussed in Section 3. Finally, some numerical simulations and a brief discussion are given in Section 4.

---

2. Basic reproduction numbers

In this section, we first present the basic reproduction number for problem (4) and its properties for the corresponding system in $\Omega$. According to [10], the basic reproduction number is an expected number of secondary cases produced by a typical infected individual during its entire period of infectiousness in a completely susceptible population, and mathematically was defined as the dominant eigenvalue of a positive linear operator. Usually the basic reproduction numbers for the spatially homogenous models were calculated by the next generation matrix method [35], while for the spatially-dependent models, the numbers could be presented in the term of the principal eigenvalue of related eigenvalue problem [1] or the spectral radius of next infection operator [37,41].

Considering the linearized problem of (3), we have

$\left\{ \begin{array}{ll} \frac{\partial I_b}{\partial t}-d_1\Delta I_b = \alpha_{b}(x)\beta_{b}(x) I_m -\gamma_{b}(x)I_b, &x\in \Omega, t>0,\\ \frac{\partial I_m}{\partial t}-d_2\Delta I_m = \frac{A_m\alpha_{m}(x)\beta_b(x)}{N_b}I_b -d_{m}(x)I_m, &x\in \Omega, t>0,\\ I_b(x) = I_m(x) = 0,&x\in \partial\Omega. \end{array} \right. \label{a31b2}$

(5)

We now consider the following linear reaction-diffusion system

$\left\{ \begin{array}{ll} \frac{\partial u}{\partial t}-D\Delta u = F(x) u -V(x) u, &x\in \Omega, t>0,\\ u (x) = 0,&x\in \partial\Omega, \end{array} \right. \label{cc234641}$

(6)

where

$u = \begin{pmatrix}I_b\\ I_m\end{pmatrix}, \;\;\;D = \begin{pmatrix} d_1&0\\ 0&d_2 \end{pmatrix},$

$F(x) = \begin{pmatrix} 0&\alpha_{b}(x)\beta_{b}(x)\\ \frac{A_m\alpha_{m}(x)\beta_b(x)}{N_b}&0 \end{pmatrix},\;\;\; V(x) = \begin{pmatrix} \gamma_{b}(x)&0\\ 0&d_m(x) \end{pmatrix}.$

In addition, the interval evolution of individuals in the infectious compartments is governed by the following linear system

$\left\{ \begin{array}{ll} \frac{\partial u}{\partial t}-D\Delta u = -V(x) u, &x\in \Omega, t>0,\\ u (x) = 0,&x\in \partial\Omega. \end{array} \right. \label{cc3234641}$

(7)

Let $X_1: = \mathbb{C} (\overline{\Omega}, \mathbb{R}^2)$ and $X^+_1: = \mathbb{C}(\overline{\Omega}, \mathbb{R}^2_+)$. Set $T(t)$ be the solution semigroup on $X_1$ associated with system (7). We let $\Psi = (\phi,\psi)$ is the density distribution of $u$ at the spatial location $x\in\Omega$, we then see that $T(t)\Psi: = (T(t)\phi,T(t)\psi)$ represents the remaining distribution of infective birds and mosquitoes at time $t$. Therefore, the distribution of total new infective members is

$\int_{0}^{\infty} F(x)[T(t)\Psi](x)dt.$

Following the idea of [37,41], we define the linear operator

$L(\Psi)(x): = \int_{0}^{\infty} F(x)[T(t)\Psi](x)dt.$

It follows from the definition, we know that $L$ is a continuous and positive operator which maps the initial infection distribution $\Psi$ to the distribution of the total members produced during the infection period. Consequently, we define the spectral radius of $L$ as the basic reproduction number of system (5), that is,

$R_0^D = \rho(L).$

As in [23], to ensure the existence of the basic reproduction numbers we consider the following linear eigenvalue problem:

$\left\{ \begin{array}{ll} -d_1\Delta \phi = \frac{\alpha_{b}(x)\beta_{b}(x)}{R} \psi -\gamma_{b}(x)\phi+\mu \phi ,&x\in \Omega , \\ -d_2\Delta \psi = \frac{A_m\alpha_{m}(x)\beta_b(x)}{N_b R}\phi -d_{m}(x)\psi+\mu\psi,&x\in \Omega,\\ \phi(x) = \psi(x) = 0,&x\in \partial\Omega. \end{array} \right. \label{a31}$

(8)

For any $R>0$, the system is strongly cooperative, that is, $\alpha_{b}(x)\beta_{b}(x)>0$ and $\alpha_{m}(x)\beta_b(x) \frac{A_m}{N_b }>0$ for all $x\in \overline \Omega$. According to [2,9,33], for any $R>0$, there exists a unique value $\mu: = \mu_1(R)$, and called the principal eigenvalue, such that problem (8) admits a unique solution pair $(\phi_R,\psi_R)$ (subject to constant multiples) with $\phi_R>0$ and $\psi_R>0$ in $\Omega$. Moreover, $\mu_1(R)$ is algebraically simple and dominant, and the following properties hold.

Lemma 2.1. $\mu_1(R)$ is continuous and strictly increasing.

With the above definition, we have the following relation between the two principal eigenvalues.

Theorem 2.2. $sign(1-R_0^D)$ = $sign(\lambda_0)$, where $R_0^D = R_0^D(\Omega, \gamma_{b}(x), d_{m}(x))$ is the principal eigenvalue of the eigenvalue problem

$\left\{ \begin{array}{ll} -d_1\Delta \phi = \frac{\alpha_{b}(x)\beta_{b}(x)}{R^D_0} \psi -\gamma_{b}(x)\phi ,&x\in \Omega , \\ -d_2\Delta \psi = \frac{A_m\alpha_{m}(x)\beta_b(x)}{N_b R^D_0}\phi -d_{m}(x)\psi,&x\in \Omega,\\ \phi(x) = \psi(x) = 0,&x\in \partial\Omega \end{array} \right. \label{a311}$

(9)

and $\lambda_0$ is the principal eigenvalue of the eigenvalue problem

$\left\{ \begin{array}{ll} -d_1\Delta \phi = \alpha_{b}(x)\beta_{b}(x)\psi -\gamma_{b}(x)\phi+\lambda_0 \phi ,&x\in \Omega , \\ -d_2\Delta \psi = \frac{A_m\alpha_{m}(x)\beta_b(x)}{N_b }\phi -d_{m}(x)\psi+\lambda_0\psi,&x\in \Omega,\\ \phi(x) = \psi(x) = 0,&x\in \partial\Omega. \end{array} \right.$

(10)

Proof. In fact, $\lambda_0 = \mu_1(1)$. On the other hand, one can easily deduce from the monotonicity with respect to the coefficients in (8) that $\lim_{R\to 0^+}\mu_1(R)<0$ and $\lim_{R\to +\infty}\mu_1(R)>0$, therefore $R^D_0$ is the unique positive root of the equation $\mu_1(R) = 0$. The result follows from the monotonicity of $\mu_1(R)$ with respect to $R$.

Remark 1. Recalling that $\mu_1$ is monotonically increasing with respect to $\beta_b(x)$, in the sense that $\mu_1(\beta_{b,1}(x))<\mu_1(\beta_{b,2}(x))$ if $\beta_{b,1}(x)\geq \beta_{b,2}(x)$ and $\beta_{b,1}(x))\not \equiv \beta_{b,2}(x)$ in $\Omega$, we deduce from Lemma 2.1 that $R^D_0$ is monotonically increasing with respect to $\beta_b(x)$, and $R^D_0>1$ if $\beta_b(x)$ is sufficiently large.

If all coefficients are constant, we can provide an explicit formula for $R_0^D$, which is known as the basic reproduction number for the corresponding diffusive WNv model.

Theorem 2.3. If $\alpha_b(x) = \alpha^*_b$, $\alpha_m(x) = \alpha^*_m, \beta_b(x) = \beta^*_b$, $\gamma_b(x) = \gamma^*_b$ and $d_m(x) = d^*_m$, then the principal eigenvalue $R^D_0$ for (9), or the basic reproduction number for model (4), is expressed by

$R_0^D(\Omega) = \sqrt{\frac{ A_m\alpha^*_b(\beta^*_b)^2\alpha^*_m}{ N_b[d_1\lambda^*+\gamma^*_m][d_2\lambda^*+d^*_m]}}\ ,$

(11)

where $\lambda^*$ is the principal eigenvalue of $-\Delta$ in $\Omega$ with null Dirichlet boundary condition.

Proof. Let $\psi^*$ be the eigenfunction corresponding to the principal eigenvalue ($\lambda^*$) of $-\Delta$ in $\Omega$ with null Dirichlet boundary condition and

$P^* = \frac{ A_m\alpha^*_b(\beta^*_b)^2\alpha^*_m}{ N_b[d_1\lambda^*+\gamma^*_b][d_2\lambda^*+d^*_m]},$

$\phi^* = \frac{\alpha^*_b\beta^*_b}{\sqrt{R^*}[d_1\lambda^*+\gamma^*_b]}\psi^*.$

Then we know that $(\phi^*,\psi^*)$ is a positive solution of problem (9) with $R_0^D = \sqrt{P^*}$, and (11) follows directly from the uniqueness of the principal eigenvalue of (9).

---

3. Coexistence

In this section, inspired by [16,18,20,28], we first study the existence of a coexistence solution to problem (4) by constructing upper and lower solutions and then we establish the non-existence of the coexistence solution to problem $(4)$. For the convenience, we let

$\left. \begin{array}{ll} S = \{(I_{b},I_{m})\in C(\overline{\Omega})\times C(\overline{\Omega});\, \, (\hat{I_{b}},\hat{I_{m}})\leq (I_{b},I_{m})\leq (\tilde{I_{b}},\tilde{I_{m}}),x\in \overline{\Omega})\}, \end{array} \right.$

where $(\hat{I_{b}},\hat{I_{m}})$ and $(\tilde{I_{b}},\tilde{I_{m}})$ are given in the following definition.

Next we are going to give a sufficient condition for problem (4) to possess a positive solution by constructing upper and lower solutions as in [28]. To achieve this, we first give an equivalent form of problem (4):

$\left\{\begin{array}{ll} -\Delta [H_1(I_b,I_m)] = f_1(x,I_b,I_m) ,&x\in \Omega, \\ -\Delta [H_2(I_b,I_m)] = f_2(x,I_b,I_m) ,&x\in \Omega,\\ I_b(x) = I_m(x) = 0,&x\in \partial\Omega. \end{array} \right. \label{b2}$

(12)

where

$H_1(I_b,I_m) = (d_1+\alpha_{1}I_b+\frac{\beta_{1} I_b}{\gamma_1+I_m})I_b,$

$H_2(I_b,I_m) = (d_2+\frac{\beta_{2} I_m}{\gamma_2+I_b}+\alpha_{2}I_m)I_m.$

Taking

$u = H_1(I_b,I_m), \;\;v = H_2(I_b,I_m),$

then the Jacobian $J$ of the transformation $(I_b,I_m)\to (u,v)$ is given by

$\left.\begin{array}{rl} J = &\frac{\partial(u,v)}{\partial(I_b,I_m)} = \left|\begin{array}{cc} d_1+2\alpha_{1}I_b+\frac{2\beta_{1}I_b}{\gamma_1+I_m}&\frac{-\beta_{1} I^2_{b}}{(\gamma_1+I_m)^2} \\ \frac{-\beta_{2}I^2_{m}}{ (\gamma_2+I_b)^2}&d_2+\frac{2\beta_{2}I_m}{\gamma_2+I_b}+2\alpha_{2}I_m\\ \end{array}\right|\\\\ = &(d_1+2\alpha_{1}I_b+\frac{2\beta_{1}I_b}{\gamma_1+I_m})(d_2+\frac{2\beta_{2}I_m}{\gamma_2+I_b}+2\alpha_{2}I_m) -\frac{\beta_{1}\beta_{2}I^2_{b} I^2_{m}}{(\gamma_1+I_m)^2 (\gamma_2+I_b)^2}\\ \geq &d_1d_2+\frac{4\beta_{1}\beta_{2}I_b I_m}{(\gamma_1+I_m)(\gamma_2+I_b)}- \frac{\beta_{1}\beta_{2}I^2_{b} I^2_{m}}{(\gamma_1+I_m)^2 (\gamma_2+I_b)^2}\\ \geq &d_1d_2>0\;\textrm{for}\;\;(I_b,I_m)\geq(0,0). \end{array} \right.$

Therefore, the inverse $I_b = g_1(u,v)$, $I_m = g_2(u,v)$ exist whenever $(I_b,I_m)\geq(0,0)$. Hence, problem (4) reduces to the following equivalent form

$\left\{\begin{array}{ll} -\Delta u + k_1u = F_1(x, I_b,I_m),& x\in\Omega, \\ -\Delta v + k_2v = F_2(x, I_b,I_m),& x\in\Omega, \\ I_b = g_1(u,v),\, I_m = g_2(u,v),& x\in\Omega,\\ u(x) = v(x) = 0,\, \, &x\in\partial\Omega, \end{array} \right. \label{b3}$

(13)

where $F_i(x, I_b,I_m) = k_i H_i(I_b,I_m) + f_i(x, I_b,I_m)$ $(i = 1, 2)$ with $k_i>0$ $(i = 1, 2)$ chosen later.

In addition, from an elementary computations one can check that

$\frac{\partial I_b}{\partial u} = \frac{d_2+2\alpha_{2}I_m+\frac{2\beta_{2}I_m}{\gamma_2+I_b}}{J},\,\, \frac{\partial I_b}{\partial v} = \frac{\frac{\beta_{1}I^2_{b}}{(\gamma_1+I_m)^2}}{J},\,\,$

$\frac{\partial I_m}{\partial u} = \frac{\frac{\beta_{2}I^2_{m}}{(\gamma_2+I_b)^2}}{J},\,\, \frac{\partial I_m}{\partial v} = \frac{d_1+2\alpha_{1}I_b+\frac{2\beta_{1}I_b}{\gamma_1+I_m}}{J},$

which shows that $I_b = g_1(u,v)$ is nondecreasing in both $u$ and $v$, while $I_m = g_2(u,v)$ is also nondecreasing in both $u$ and $v$ for all $(I_b,I_m)\geq(0,0)$.

For the later analysis, we present the definition of upper and lower solutions to problem (13) as follows.

Definition 3.1. Assume that $F_1$ and $F_2$ are nondecreasing with respect to $I_b$ and $I_m$. A pair of $4$- nonnegative functions $(\tilde{I_{b}}, \tilde{I_{m}}, \tilde{u}, \tilde{v})$, $( \hat{I_{b}}, \hat{I_{m}}, \hat{u}, \hat{v})$ in $\mathcal{C}^2(\Omega)\cap \mathcal{C}(\overline{\Omega})$ are called ordered upper and lower solutions of $(13)$, if

$(0, 0)\leq( \hat{I}_b, \hat{I}_m) \leq (\tilde{I}_b,\, \tilde{I}_m)\leq (N_b, A_m), (\hat{u}, \hat{v}) \leq (\tilde{u},\, \tilde{v})$

and

$\left\{\begin{array}{lll} -\Delta \tilde{u}+ k_1\tilde{u} \geq F_1(x,\tilde{I_b},\tilde{I_m}) ,&x\in \Omega, \\ -\Delta \tilde{v}+ k_2\tilde{v} \geq F_2(x,\tilde{I_b},\tilde{I_m}) ,&x\in \Omega,\\ -\Delta \hat{u}+ k_1\hat{u} \leq F_1(x,\hat{I_b},\hat{I_m}) ,&x\in \Omega, \\ -\Delta \hat{v}+ k_2\hat{v} \leq F_2(x,\hat{I_b},\hat{I_m}) ,&x\in \Omega,\\ \tilde{I}_b\geq g_1(\tilde{u},\tilde{v}), \;\;\;\;\hat{I}_b\leq g_1(\hat{u}, \hat{v}),&x\in \Omega,\\ \tilde{I}_m\geq g_2(\tilde{u},\tilde{v}), \;\;\; \hat{I}_m\leq g_2(\hat{u}, \hat{v}),&x\in \Omega,\\ \tilde{u} (x)\geq 0\geq \hat{u} (x), \;\;\tilde{v} (x)\geq 0 \geq \hat{ v} (x),&x\in \partial\Omega. \end{array}\right. \label{b4}$

(14)

For definiteness, we select

$\tilde{I}_b = g_1(\tilde{u},\tilde{v}),~~\tilde{I}_m = g_2(\tilde{u},\tilde{v});$

$\hat{I}_b = g_1(\hat{u}, \hat{v}),~~\hat{I}_m = g_2(\hat{u}, \hat{v}),$

which is equivalent to

$\tilde{u} = H_1(\tilde{I}_b,\tilde{I}_m),~~\tilde{v} = H_2(\tilde{I}_b,\tilde{I}_m);$

$\hat{u} = H_1(\hat{I}_b,\hat{I}_m),~~\hat{v} = H_2(\hat{I}_b,\hat{I}_m).$

Then the requirements of $(\tilde{I}_b,\tilde{I}_m)$ and $(\hat{I}_b,\hat{I}_m)$ in (14) are satisfied and those of $(\tilde{u},\tilde{v}),(\hat{u}, \hat{v})$ are reduced to

$\left\{\begin{array}{lll} -\Delta [H_1(\tilde{I}_b,\tilde{I}_m)]+ k_1H_1(\tilde{I}_b,\tilde{I}_m) \geq F_1(x,\tilde{I_b},\tilde{I_m}) ,&x\in \Omega, \\ -\Delta [H_2(\tilde{I}_b,\tilde{I}_m)]+ k_2H_2(\tilde{I}_b,\tilde{I}_m) \geq F_2(x,\tilde{I_b},\tilde{I_m}) ,&x\in \Omega,\\ -\Delta [H_1(\hat{I}_b,\hat{I}_m)]+ k_1H_1(\hat{I}_b,\hat{I}_m) \leq F_1(x,\hat{I_b},\hat{I_m}) ,&x\in \Omega, \\ -\Delta [H_2(\hat{I}_b,\hat{I}_m)]+ k_2H_2(\hat{I}_b,\hat{I}_m) \leq F_2(x,\hat{I_b},\hat{I_m}) ,&x\in \Omega,\\ \tilde{I}_b (x)\geq 0\geq \hat{I}_b (x), \;\;\;\tilde{I}_m (x)\geq 0 \geq \hat{I}_m (x),&x\in \partial\Omega. \end{array}\right. \label{b5}$

(15)

Now we consider the monotonicity of $F_i$ $(i = 1,2)$. From direct computations it is easy to see that

$\frac{\partial F_1}{\partial I_b} = k_1(d_1+2\alpha_{1}I_b+\frac{2\beta_{1}I_b}{\gamma_1+I_m})-\alpha_{b}(x)\beta_{b}(x)\frac{I_m}{N_b}-\gamma_{b}(x),$

$\frac{\partial F_2}{\partial I_m} = k_2(d_2+\frac{2\beta_{2}I_m}{\gamma_2+I_b}+2\alpha_{2}I_m)-\alpha_{m}(x)\beta_b(x) \frac{I_b}{N_b}-d_{m}(x).$

If we choose

$k_1 = \max\limits_{x\in \overline{\Omega}}\{\frac{\alpha_b\beta_b A_m+N_b\gamma_b}{N_b d_1}\}(x), \;\;\; k_2 = \max\limits_{x\in \overline{\Omega}}\{\frac{\alpha_m\beta_b+d_m}{d_2}\}(x),$

then $F_1$ and $F_2$ are increasing with respect to $I_b$ and $I_m$, respectively, as long as $(0.0)\leq(I_b,\, I_m)\leq (N_b, A_m)$. On the other hand, the direct calculations show that

$\frac{\partial F_1}{\partial I_m} = -k_1\frac{\beta_{1}}{(\gamma_1+I_m)^2} I^2_{b}+\alpha_{b}(x)\beta_{b}(x)\frac{N_b-I_b}{N_b},$

$\frac{\partial F_2}{\partial I_b} = -k_2\frac{\beta_{2}}{(\gamma_2+I_b)^2} I^2_{m}-\alpha_{m}(x)\beta_b(x) \frac{A_m-I_m}{N_b}.$

Note $\frac{\partial F_1}{\partial I_m}(N_b,A_m)\leq0$ and $\frac{\partial F_2}{\partial I_b}(N_b,A_m)\leq0$ for any $\frac{\beta_{1}}{\gamma_1}$ and $\frac{\beta_{2}}{\gamma_2}$, respectively. To ensure that $\frac{\partial F_1}{\partial I_m}\geq 0$, $\frac{\partial F_2}{\partial I_b}\geq 0$, we have to modify the upper solution and we seek $((1-\delta_0)N_b, (1-\delta_0)A_m)$ as a new upper solution, where

$\delta_0\leq \min\limits_{x\in \overline{\Omega}}\{\frac{1}{2}, \frac{\gamma_{b}}{\alpha_{b}\beta_{b}A_m}, \frac{d_{m}}{\alpha_{b}\beta_{b}}\}$

from the first and second inequalities of (15). Let

$\beta^{*}_1 = \min\limits_{x\in \overline{\Omega}}\,\frac{\delta_0\alpha_b\beta_b\gamma^2_{1}}{N^2_{b} k_1}(x),\;\;\; \beta^{*}_2 = \min\limits_{x\in \overline{\Omega}}\,\frac{\delta_0\alpha_m\beta_b\gamma^2_{2}}{N_b A_m k_2}(x),$

then for $\beta_1\leq \beta^{*}_1$ and $\beta_2\leq \beta^{*}_2$, $F_i$ $(i = 1, 2)$ is monotone nondecreasing with respect to $I_b$ and $I_m$. Consequently, $(\tilde{I_{b}}, \tilde{I_{m}}, \tilde{u}, \tilde{v})$, $( \hat{I_{b}}, \hat{I_{m}}, \hat{u}, \hat{v})$ are a pair of ordered upper and lower solutions of problem (13).

To present the existence of a positive solution to (4), it suffices to find a pair of upper and lower solutions of (4). We seek such as in the form $(\tilde{I}_b,\,\tilde{I}_m) = (M_1,\, M_2)$, $(\hat{I}_b,\,\hat{I}_m) = (g_1(\delta d_1\phi,\,\delta d_2\psi), g_2(\delta d_1\phi,\,\delta d_2\psi))$ where $M_i$ $(i = 1,2)$ and $\delta$ are some positive constants with $\delta$ small enough, $(\phi, \psi)\equiv(\phi(x),\psi(x))$ is (normalized) positive eigenfunction corresponding to $\lambda_0$, and $\lambda_0$ is the principal eigenvalue of eigenvalue problem (10).

Indeed, $(M_1,\,M_2)$ and $(g_1(\delta d_1\phi,\,\delta d_2\psi), g_2(\delta d_1\phi,\,\delta d_2\psi))$ satisfy the inequalities in (15) if

${\small \left\{ \begin{array}{ll} -\Delta [(d_1+\alpha_{1}M_1+\frac{\beta_{1}M_1}{\gamma_1+M_2})M_1]\geq \alpha_{b}(x)\beta_{b}(x)\frac{(N_b-M_1)}{N_b} M_2 -\gamma_{b}(x) M_1 , \\ -\Delta [(d_2+\frac{\beta_{2}M_2}{\gamma_2+M_1}+\alpha_{2}M_2)M_2]\geq\alpha_{m}(x)\beta_b(x) \frac{(A_m-M_2)}{N_b} M_1 -d_{m}(x) M_2,\\ -\Delta[d_1\phi]\leq \alpha_{b}(x)\beta_{b}(x)\frac{(N_b-\hat{I}_b)}{N_b} d_2\psi/(d_2+\alpha_{2}\hat{I}_m+\frac{\beta_{2}\hat{I}_m}{\gamma_2+\hat{I}_b}) \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\;\;\;\;\;\;\; \;\;\;\;\; -\gamma_{b}(x)d_1\phi/(d_1+\alpha_{1}\hat{I}_b+\frac{\beta_{1}\hat{I}_b}{\gamma_1+\hat{I}_m}), \\ -\Delta[d_2\psi]\leq \alpha_{m}(x)\beta_b(x) \frac{(A_m-\hat{I}_m)}{N_b} d_1\phi/(d_1+\alpha_{1}\hat{I}_b+\frac{\beta_{1}\hat{I}_b}{\gamma_1+\hat{I}_m}) \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\;\;\;\;\; -d_{m}(x) d_2\psi/(d_2+\alpha_{2}\hat{I}_m+\frac{\beta_{2}\hat{I}_m}{\gamma_2+\hat{I}_b}). \end{array} \right. \label{b6} }$

(16)

The first two inequalities in (16) hold if we set

$\begin{array}{ll} (M_1,\,M_2) = ((1-\delta_0)N_b,(1-\delta_0)A_m). \end{array} \label{b7}$

(17)

Next, we notice that the relations

$\delta d_1\phi = (d_1+\alpha_{1}\hat{I}_b+\frac{\beta_{1}\hat{I}_b}{\gamma_1+\hat{I}_m})\hat{I}_b,\ \delta d_2\psi = (d_2+\alpha_{2}\hat{I}_m+\frac{\beta_{2}\hat{I}_m}{\gamma_2+\hat{I}_b})\hat{I}_m$

imply that $0 <\hat{I}_b\leq \delta\phi$ and $0 <\hat{I}_m\leq \delta\psi.$

If $R_0^D>1$, the principal eigenvalue of problem (10) is $\lambda_0 <0$, therefore we can choose $\delta$ sufficiently small such that the last two inequalities in (16) hold. Consequently, the pair $(\tilde{I}_b,\tilde{I}_m) = (M_1,M_2)$, $(\hat{I}_b,\hat{I}_m) = (g_1(\delta d_1\phi,\,\delta d_2\psi), g_2(\delta d_1\phi,\,\delta d_2\psi))$ are ordered upper and lower solutions of problem (4), respectively.

Using Theorem 2.1 of [28] leads to the following existence result :

Theorem 3.2. If $R_0^D>1$, problem $(4)$ admits at least one coexistence solution $(I_b(x),I_m(x))$ provided that $\beta_{1}$ and $\beta_{2}$ are sufficiently small.

To establish the non-existence of the coexistence solution to problem $(4)$, we have the following result.

Theorem 3.3. If $R_0^D(\Omega,\gamma_{b}(x) \frac{1}{1+(\frac{\alpha_{1}}{d_1}+\frac{\beta_{1}}{\gamma_1 d_1}) N_b}, d_{m}(x) \frac{1}{1+(\frac{\alpha_{2}}{d_2}+\frac{\beta_{2}}{\gamma_2 d_2}) A_m} )\leq1$, problem $(4)$ has no positive solution.

Proof. Suppose $(I_{b}^*(x), I_{m}^*(x))$ is a coexistence solution of problem $(4)$, that is, $(I_{b}^*(x), I_{m}^*(x))>(0,0)$ in $\Omega$ and satisfies

${\small \left\{\begin{array}{ll} -\Delta [(d_1+\alpha_{1}I_{b}^* +\frac{\beta_{1} I_{b}^*}{\gamma_1+I_{m}^*})I_{b}^*] = \alpha_{b}(x)\beta_{b}(x)\frac{(N_b-I_{b}^*)}{N_b} I_{m}^* -\gamma_{b}(x)I_{b}^*,& x\in \Omega , \\ -\Delta [(d_2+\frac{\beta_{2} I_{m}^*}{\gamma_2+I_{b}^*}+\alpha_{2}I_{m}^*)I_{m}^*] = \alpha_{m}(x)\beta_{b}(x)\frac{(A_m-I_{b}^*)}{N_b} I_{b}^* -d_m (x)I_{m}^*,&x\in \Omega,\\ I_{b}^*(x) = I_{m}^*(x) = 0,&x\in \partial\Omega. \end{array} \right. \label{a453} }$

(18)

First by the upper and lower solution method we know that $(I_{b}^*, I_{m}^*)\leq(N_b, A_m).$

Second, letting $w = (d_1+\alpha_{1}I_{b}^* +\frac{\beta_{1} I_{b}^*}{\gamma_1+I_{m}^*})I_{b}^*/d_1$ and $z = (d_2+\frac{\beta_{2} I_{m}^*}{\gamma_2+I_{b}^*}+\alpha_{2}I_{m}^*)I_{m}^*/d_2$, we have

$\left\{\begin{array}{ll} -d_1\Delta w = \alpha_{b}(x)\beta_{b}(x)\frac{(N_b-I_{b}^*)}{N_b} \frac{d_2 z}{d_2+\frac{\beta_{2} I_{m}^*}{\gamma_2+I_{b}^*}+\alpha_{2}I_{m}^*}&\\ \;\;\;\;\;\;\;\;\;\;\;\;-\gamma_{b}(x)\frac{d_1 w}{d_1+\alpha_{1}I_{b}^* +\frac{\beta_{1} I_{b}^*}{\gamma_1+I_{m}^*}},& x\in \Omega , \\ -d_2\Delta z = \alpha_{m}(x)\beta_{b}(x)\frac{(A_m-I_{b}^*)}{N_b} \frac{d_1 w}{d_1+\alpha_{1}I_{b}^* +\frac{\beta_{1} I_{b}^*}{\gamma_1+I_{m}^*}}&\\ \;\;\;\;\;\;\;\;\;\;\;\; -d_m (x)\frac{d_2 z}{d_2+\frac{\beta_{2} I_{m}^*}{\gamma_2+I_{b}^*}+\alpha_{2}I_{m}^*} ,&x\in \Omega,\\ w = I_{b}^*(x) = 0,&x\in \partial\Omega,\\ z = I_{m}^*(x) = 0,&x\in \partial\Omega, \end{array} \right. \label{a454}$

(19)

which means that

$\left\{\begin{array}{ll} -d_1\Delta w <\alpha_{b}(x)\beta_{b}(x)z -\frac{\gamma_{b}(x)}{1+(\frac{\alpha_{1}}{d_1}+\frac{\beta_{1}}{\gamma_1 d_1}) N_b} w,& x\in \Omega , \\ -d_2\Delta z<\alpha_{m}(x)\beta_{b}(x)\frac{A_m}{N_b} w -\frac{d_m (x)}{1+(\frac{\alpha_{2}}{d_2}+\frac{\beta_{2}}{\gamma_2 d_2}) A_m} z ,&x\in \Omega,\\ w = z = 0,&x\in \partial\Omega. \end{array} \right. \label{a455}$

(20)

On the other hand, the principal eigenvalue $\lambda_0$ in problem (10) meets

$\left\{ \begin{array}{ll} -d_1\Delta \phi = \alpha_{b}(x)\beta_{b}(x) \psi -\frac{\gamma_{b}(x)}{1+(\frac{\alpha_{1}}{d_1}+\frac{\beta_{1}}{\gamma_1 d_1}) N_b}\phi+\lambda_0\phi ,&x\in \Omega , \\ -d_2\Delta \psi = \alpha_{m}(x)\beta_b(x) \frac{A_m}{N_b} \phi -\frac{d_m (x)}{1+(\frac{\alpha_{2}}{d_2}+\frac{\beta_{2}}{\gamma_2 d_2}) A_m}\psi+\lambda_0\psi ,&x\in \Omega,\\ \phi(x) = \psi(x) = 0,&x\in \partial\Omega. \end{array} \right. \label{boo6}$

(21)

Comparing (20) with (21), we can easily deduce from the monotonicity with respect to the coefficients in (21) that $\lambda_0$ is monotone decreasing with respect to $\beta_{b}(x)$, which implies that $\lambda_0 <0$. Recalling Theorem 2.2 we can get that $R^D_0>1$, which is contrary to $R^D_0\leq1$.

Remark 2. Assume that all coefficients of $(4)$ are spatially-independent. $R^D_0$ is represented by (11). If $\alpha_b^*, \alpha^*_m$ or $\beta^*_b$ is big, then $R_0^D>1$ and problem $(4)$ admits at least one coexistence solution provided that $\beta_{1}$ and $\beta_{2}$ are sufficiently small. On the other hand, if $\alpha_b^*, \alpha^*_m$ or $\beta^*_b$ is small enough, then $R_0^D\leq1$ and problem $(4)$ has no positive solution.

Next we apply the monotone iterative schemes to construct the true solutions of (4). It follows from $R_0^D>1$, we know that $(M_1,M_2)$ and

$(g_1(\delta d_1\phi,\,\delta d_2\psi), g_2(\delta d_1\phi,\,\delta d_2\psi))$

are ordered upper and lower solution of problem (4), respectively. Using $(\bar{I}_b^{(0)},\bar{I}_m^{(0)}) = ((1-\delta_0)N_b,(1-\delta_0)A_m)$ and $(\underline{I}_b^{(0)}, \underline{I}_m^{(0)}) = (g_1(\delta d_1\phi,\,\delta d_2\psi), g_2(\delta d_1\phi,\,\delta d_2\psi))$ as two initial iterations, we can construct two sequences $\{(\bar{u}^{(n)},\bar{v}^{(n)})\}$ and $\{(\underline{u}^{(n)}, \underline{v}^{(n)})\}$ from the iteration process

$\left\{ \begin{array}{ll} -\Delta \bar{u}^{(n)}+k_1\bar{u}^{(n)} = F_1(x, \bar{I}_b^{(n-1)}, \bar{I}_m^{(n-1)}), &x\in \Omega,\\ -\Delta \bar{v}^{(n)}+k_2\bar{v}^{(n)} = F_2(x, \bar{I}_b^{(n-1)}, \bar{I}_m^{(n-1)}),&x\in \Omega,\\ -\Delta \underline{u}^{(n)}+k_1\underline{u}^{(n)} = F_1(x, \underline{I}_b^{(n-1)}, \underline{I}_m^{(n-1)}),&x\in \Omega,\\ -\Delta \underline{v}^{(n)}+k_2\underline{v}^{(n)} = F_2(x, \underline{I}_b^{(n-1)}, \underline{I}_m^{(n-1)}),&x\in \Omega,\\ \bar{I}_b^{(n)} = g_1(\bar{u}^{(n)},\bar{v}^{(n)}), \;\;\;\;\;\;\underline{I}_b^{(n)} = g_1(\underline{u}^{(n)}, \underline{v}^{(n)}),&x\in \Omega,\\ \bar{I}_m^{(n)} = g_2(\bar{u}^{(n)},\bar{v}^{(n)}), \;\;\;\;\;\underline{I}_m^{(n)} = g_2(\underline{u}^{(n)}, \underline{v}^{(n)}),&x\in \Omega,\\ \bar{u}^{(n)} (x) = \underline{u}^{(n)} (x) = 0,\;\;\; \bar{v}^{(n)} (x) = \underline{ v}^{(n)} (x) = 0,&x\in \partial\Omega, \end{array} \right. \label{b10}$

(22)

where $n = 1,2, \cdots$.

As in Lemma 3.1 of [28], the sequences $\{(\bar{u}^{(n)},\bar{v}^{(n)})\}$ and $\{(\underline{u}^{(n)}, \underline{v}^{(n)})\}$ governed by (22) are well-defined and possess the monotone property

$(\hat{u},\hat{v})\leq (\underline{u}^{(n-1)},\underline{v}^{(n-1)}) \leq (\underline{u}^{(n)}, \underline{v}^{(n)}) \leq (\bar{u}^{(n)},\bar{v}^{(n)})$

$\leq (\bar{u}^{(n-1)},\bar{v}^{(n-1)}) \leq (\tilde{u},\tilde{v})\ \textrm{for }\ n = 1,2, \cdots.$

Hence, the pointwise limits

$\lim\limits_{n \to \infty}(\bar{u}^{(n)},\bar{v}^{(n)})) = (\bar{u},\, \bar{v}),\, \;\;\;\lim\limits_{n \to \infty}(\underline{u}^{(n)},\, \underline{u}^{(n)}) = (\underline{u},\, \underline{v})$

exist and their limits possess the relation

$\begin{array}{ll} (\hat{u},\hat{v})\leq (\underline{u}^{(n)},\underline{v}^{(n)}) \leq (\underline{u},\underline{v}) \leq (\bar{u},\bar{v}) \leq (\bar{u}^{(n)},\bar{v}^{(n)}) \leq (\tilde{u},\tilde{v}) \end{array} \label{b01}$

(23)

for every $n = 1,2,\cdots.$

The last three equations of (22) give

$\bar{I}_b^{(n)} = g_1(\bar{u}^{(n)},\bar{v}^{(n)}), \;\;\; \underline{I}_b^{(n)} = g_1(\underline{u}^{(n)}, \underline{v}^{(n)}),$

$\bar{I}_m^{(n)} = g_2(\bar{u}^{(n)},\bar{v}^{(n)}), \;\;\; \underline{I}_m^{(n)} = g_2(\underline{u}^{(n)}, \underline{v}^{(n)}),$

which is equivalent to

$\begin{array}{ll} \bar{u}^{(n)} = H_1(\bar{I}_b^{(n)},\bar{I}_m^{(n)}), \;\;\;\underline{u}^{(n)} = H_1(\underline{I}_b^{(n)},\underline{I}_m^{(n)}),\\ \bar{v}^{(n)} = H_2(\bar{I}_b^{(n)},\bar{I}_m^{(n)}), \;\;\;\underline{v}^{(n)} = H_2(\underline{I}_b^{(n)},\underline{I}_m^{(n)}). \end{array} \label{b02}$

(24)

Now, by the above relation, letting $n\rightarrow \infty$ and applying the standard regularity argument for elliptic boundary problems, we derive that $(\bar{I}_b,\,\bar{I}_m)$ and $(\underline{I}_b,\,\underline{I}_m)$ satisfy

$\left\{ \begin{array}{ll} -\Delta [H_1(\bar{I}_b,\bar{I}_m)]+k_1H_1(\bar{I}_b,\bar{I}_m) = F_1(x, \bar{I}_b, \bar{I}_m), &x\in \Omega,\\ -\Delta [H_2(\bar{I}_b,\bar{I}_m)]+k_2H_2(\bar{I}_b,\bar{I}_m) = F_2(x, \bar{I}_b, \bar{I}_m),&x\in \Omega,\\ -\Delta [H_1(\underline{I}_b,\underline{I}_m)]+k_1H_1(\underline{I}_b,\underline{I}_m) = F_1(x, \underline{I}_b, \underline{I}_m),&x\in \Omega,\\ -\Delta [H_2(\underline{I}_b,\underline{I}_m)]+k_2H_2(\underline{I}_b,\underline{I}_m) = F_2(x, \underline{I}_b, \underline{I}_m),&x\in \Omega,\\ \bar{I}_b (x) = \underline{I}_b (x) = 0, \;\;\;\bar{I}_m (x) = \underline{ I}_m (x) = 0,&x\in \partial\Omega, \end{array} \right. \label{b003}$

(25)

which is equivalent to

$\left\{ \begin{array}{ll} -\Delta [H_1(\bar{I}_b,\bar{I}_m)] = f_1(x, \bar{I}_b, \bar{I}_m), &x\in \Omega,\\ -\Delta [H_2(\bar{I}_b,\bar{I}_m)] = f_2(x, \bar{I}_b, \bar{I}_m),&x\in \Omega,\\ -\Delta [H_1(\underline{I}_b,\underline{I}_m)] = f_1(x, \underline{I}_b, \underline{I}_m),&x\in \Omega,\\ -\Delta [H_2(\underline{I}_b,\underline{I}_m)] = f_2(x, \underline{I}_b, \underline{I}_m),&x\in \Omega,\\ \bar{I}_b (x) = \underline{I}_b (x) = 0, \;\;\;\bar{I}_m (x) = \underline{ I}_m (x) = 0,&x\in \partial\Omega. \end{array} \right. \label{b005}$

(26)

Therefore $(\bar{I}_b,\bar{I}_m)$ and $(\underline{I}_b,\underline{I}_m)$ are true solutions of (4). Moreover, $(\bar{I}_b,\bar{I}_m)$ and $(\underline{I}_b, \underline{I}_m)$ are maximal and minimal solutions in the sense that $(I_b,I_m)$ is any other solution of (4) in the sector

$< (\hat{I}_b,\hat{I}_m),(\tilde{I}_b, \tilde{I}_m)>: = \{(I_b, I_m)\in \mathcal{C}(\overline \Omega); (\hat{I}_b,\hat{I}_m)\leq (I_b, I_m)\leq (\tilde{I}_b, \tilde{I}_m)\, \textrm{on}\,\overline \Omega\},$

then $(\underline{I}_b, \underline{I}_m)\leq(I_b,I_m)\leq(\bar{I}_b,\bar{I}_m)$ on $\overline{\Omega}$. Furthermore, if $\bar{I}_b = \underline{I}_b$ or $\bar{I}_m = \underline{I}_m$, then $(\bar{I}_b,\bar{I}_m) = (\underline{I}_b,\underline{I}_m): = (I^*_b,I^*_m)$ and $(I^*_b,I^*_m)$ is the unique solution of $(4)$ in $\overline{\Omega}$. To achieve this, in fact, a subtraction of the third equation from the first equation in (26) yields that

$-\Delta [\frac{\beta_1 (I^*_b)^2(\bar{I}_m-\underline{I}_m)}{(\gamma_1+\bar{I}_m)(\gamma_1+\underline{I}_m)}] = \alpha_{b}(x)\beta_{b}(x)\frac{(N_b-I^*_b)}{N_b} (\bar{I}_m-\underline{I}_m)$

In light of $\beta_1>0$, $I^*_b>0$, $\alpha_{b}(x)\beta_{b}(x)>0$, $N_b-I^*_b>0$ and $(\bar{I}_m-\underline{I}_m) = 0$ on $\partial\Omega$, the above equation gives $\bar{I}_m\equiv\underline{I}_m$ in $\Omega$. Similarly as above one can show that $\bar{I}_b\equiv\underline{I}_b$ in $\Omega$. Therefore, $(\bar{I}_b,\bar{I}_m) = (\underline{I}_b,\underline{I}_m): = (I^*_b,I^*_m)$ which is unique solution of $(4)$ in $S$.

The above conclusions lead to the following theorem.

Theorem 3.4. Let $(\tilde{I}_b, \tilde{I}_m)$ and $(\hat{I}_b,\hat{I}_m)$ be a pair of ordered upper and lower solutions of $(4)$, respectively, then the sequences $\{(\bar{I}_b^{(n)},\bar{I}_m^{(n)})\}$ and $\{(\underline{I}_b^{(n)},\underline{I}_m^{(n)})\}$ provided from $(22)$ converge monotonically from above to a maximal solution $(\bar{I}_b,\bar{I}_m)$ and from below to a minimal solution $(\underline{I}_b,\underline{I}_m)$ in $S$, respectively, and satisfy the relation

$(\hat{I_{b}},\hat{I_{m}})\leq (\underline{I}_b^{(n)},\underline{I}_m^{(n)}) \leq (\underline{I}_b^{(n+1)},\underline{I}_m^{(n+1)})\leq(\underline{I}_b,\underline{I}_m)\leq(\bar{I}_b,\bar{I}_m) \leq (\bar{I}_b^{(n+1)},\bar{I}_m^{(n+1)})$

$\leq (\bar{I}_b^{(n)},\bar{I}_m^{(n)}) \leq (\tilde{I_{b}},\tilde{I_{m}})\ \textrm{for}\ n = 1,2, \cdots$

additionally, if $\bar{I}_b = \underline{I}_b$ or $\bar{I}_m = \underline{I}_m$, then $(\bar{I}_b,\bar{I}_m) = (\underline{I}_b,\underline{I}_m)( = (I^*_b,I^*_m))$ and $(I^*_b,I^*_m)$ is the unique solution of $(4)$ in $S$.

---

4. Numerical simulation and discussion

In this section, in order to illustrate our theoretical results, we simulate problem (4) with the following coefficients and parameters:

$d_1 = 0.2, d_2 = 0.4, \alpha_1 = 0.03, \alpha_2 = 0.04, \gamma_1 = 1, \gamma_2 = 1,$

$\alpha_b = 1+0.88\sin(\frac{\pi}{100} x), \;\;\; \alpha_m = 1+0.16\sin(\frac{\pi}{100} x),$

$\gamma_b = 1+0.6\sin(\frac{\pi}{100} x),\;\;\; d_m = 1+0.029\sin(\frac{\pi}{100} x),$

and $\beta_b = 1+0.09\sin(\frac{\pi}{100} x)$ and we also take the ratio $A_m/N_b = 20$ as in [21].

From Fig. 1, it is easy to see that there exist the upper solution sequence $\{(\overline {I}_b^{(n)}, \overline{I}_m^{(n)})\}$ which is monotone decreasing and the lower solution sequence $\{(\underline {I}_b^{(n)},\underline{I}_m^{(n)})\}$ which is monotone increasing, then one can see that there exists at least a coexistence solution of (4).

In this paper, to understand the impact of a cross-diffusion and environmental heterogeneity on the dynamics of WNv, we consider coexistence states of a cross-diffusive WNv model in heterogenous environments under Dirichlet boundary condition. This problem without spatially-dependent coefficients is similar to that has been studied in [16]. It is worth mentioning that problem (1) is weakly-coupled parabolic system which only involves the random diffusion in a homogeneous environment. However, in addition to the dispersive force, the diffusions of birds and mosquitoes are also interacted by each other and reaction depends on spatial heterogeneity of the environment. Therefore, we introduce the cross-diffusion terms $\Delta [(d_1+\alpha_{1}I_b+\frac{\beta_{1}I_b}{\gamma_1+I_m})I_b],$ $\Delta [(d_2+\frac{\beta_{2}I_m}{\gamma_2+I_b}+\alpha_{2}I_m)I_m]$ to model (1), which can better describe the interplay between birds and mosquitoes in diffusion.

The main result of this paper is twofold. Firstly, we introduce a definition of $R_0^D$, which is known as the basic reproduction number of problem (4) (Theorem 2.2). In the case that all coefficients are constants, we provide an explicit formula for $R_0^D$ (Theorem 2.3). Secondly, the coexistence of problem (4) is investigated by using method of upper and lower solutions and its associated monotone iterative schemes (Theorem 3.2 and Fig. 1) under condition $R_0^D>1$ provided that $\beta_{1}$ and $\beta_{2}$ are sufficiently small, whereas if $R_0^D\leq1$, problem (4) has no coexistence solution (Theorem 3.3). Our results show that no existence exists for small WNv transmission probabilities ($\alpha_m$ and $\alpha_b$), and small biting rate of mosquitoes on birds ($\beta_b$) (Remark 2). Moreover, the coexistence solution of problem (4) is between the maximal and minimal solution $(\bar{I}_b,\bar{I}_m)$ and $(\underline{I}_b,\underline{I}_m)$, respectively, and the true solution can be obtained by constructing the monotone iterative sequences (Theorem 3.4). However, the uniqueness of coexistence solution is still unclear.

We believe that the strongly-coupled problem (4) can produce much more complex dynamics of WNv than the weakly-coupled system (1). Such problems need further investigations. In fact, even for the corresponding parabolic problems with cross-diffusion, the existence of the solution is known only for some special cases, see [24,32] and references therein. To further investigate the effect of cross-diffusion in comparison to no cross-diffusion or small cross-diffusion, we come back to problem (3), Fig. 2 shows that the global solution of problem (3) exists and stabilizes to a positive steady-state for small cross-diffusion ($\beta_1 = 0.132$ and $\beta_2 = 0.11$), we can also see that the global solution of (3) exists for $\beta_1\leq 0.132$ and $\beta_2\leq 0.11$ by simulations. However, if we choose a little big cross-diffusion, for example, $\beta_1 = 0.133$ and $\beta_2 = 0.11$, we can see from Fig. 3 that the global solution of problem (3) does not exist. We leave it for future work.

---