具有扩散项的互惠–寄生耦合模型图灵斑图研究

Analysis of Spatial Pattern in Mutualistic-Parasitic Coupled System with Diffusion

摘要: 本文将扩散因子引入互惠–寄生耦合系统，建立互惠–寄生耦合空间模型，研究其斑图动力学性态。利用斑图理论，给出耦合系统产生Hopf分支和Turing分支的临界表达式，并由此确定出图灵斑图出现的参数区域，运用有限差分法模拟模型的斑图结构。结果表明，当扩散存在时互惠–寄生耦合系统具有黑眼(点状)斑图结构，寄生物在空间呈现孤立的高密度分布。

Abstract: In this paper, we establish mutualistic-parasitic coupled spatial model incorporating diffusion, and study the dynamical behaviors of this spatial model. By spatial pattern theory, we obtain the conditions for Hopf bifurcation and Turing bifurcation, and give the exact parameter space for Turing domain. By utilizing difference approximation method, we find that the spatial model ex-hibits black eye spotted pattern, which shows that the diffusion can result in an isolated high den-sity of parasite in the space.

1. 引言

寄生物与宿主协同进化的研究，不仅在种群生态学领域受到广泛关注，而且对农作物培育、畜牧养殖、致病细胞的过度复制等领域的预测研究有重要启示 [1] - [7] 。目前研究寄生物与宿主相互作用的经典模型都是基于寄生物将导致宿主种群减少的前提下建立的。然而，寄生物在很多情况下也会导致宿主种群的增长，如大肠杆菌和乳酸杆菌与其宿主的关系 [8] [9] 。若系统中寄生物与宿主之间既存在互惠关系也存在寄生关系，则可称该系统为互惠–寄生耦合系统。文献 [10] 在不考虑扩散因素的情况下，建立了互惠–寄生耦合模型，并研究了其动力学性态。事实上，扩散广泛存在于生物系统中，自然界各种生物通过扩散方式扩大栖居范围，以便获得更多的资源。本文拟在互惠–耦合模型的基础上，引入扩散因子，建立具有扩散项的互惠–寄生耦合模型空间模型，研究其斑图动力学性态，进而得出扩散对宿主与寄生物密度空间分布的影响。

2. 主要结果

2.1. 建立空间模型

首先给出互惠–寄生耦合模型 [10]

$\{\begin{array}{l}\frac{\text{d}H}{\text{d}t}=H\left[r_1-{\alpha }_{1}H+\left(\mu -\frac{P}{H}\right)\beta P\right]\equiv F\left(H,P\right)\\ \frac{\text{d}P}{\text{d}t}=P\left(r_2-{\alpha }_2\frac{P}{H}\right)\equiv G\left(H,P\right)\end{array}$ (1)

其中 $H\left(t\right),P\left(t\right)$ 分别表示宿主与寄生物的种群密度， $\mu$ 表示寄生物对宿主的固有促进水平，当系统内寄生物与宿主的数量之比不超过 $\mu$ 时，寄生物对宿主具有促进作用；当该数量之比超过 $\mu$ 时，寄生物对宿主具有抑制作用；其它参数的含义见表1。

注：这里的寄生物的固有增长率 $r_2$ 是指其不受密度制约情况下的增长率，以上各个系数均为正实数。

由Fick法则我们假定个体的扩散是随机的，这导致如下的表达：

$\frac{\partial H}{\partial t}=d_1{\nabla }^{2}H$ ， $\frac{\partial P}{\partial t}=d_2{\nabla }^{2}P$ ，

其中 ${\nabla }^2=\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}$ 是二维空间的Laplacian算子， $d_1$ ， $d_2$ 分别表示扩散系数，结合模型(1)，可得如下空间模型：

Table 1. Definitions and the symbols used in model (1)

表1. 模型1中的符号及含义

$\{\begin{array}{l}\frac{\partial H}{\partial t}=F\left(H,P\right)+d_1{\nabla }^{2}H\\ \frac{\partial P}{\partial t}=G\left(H,P\right)+d_2{\nabla }^{2}P\end{array}$ (2)

其中 $F\left(H,P\right)=H\left[r_1-{\alpha }_{1}H+\left(\mu -\frac{P}{H}\right)\beta P\right]$ ， $G\left(H,P\right)=P\left(r_2-{\alpha }_2\frac{P}{H}\right)$ 。从生物意义上考虑，假设这里的所有系数都是正的。

2.2. 空间模型的分支分析

设模型(1)存在唯一的均匀定态解 $\left(H^{*},P^{*}\right)$ ，即 $H^{*}$ 与 $P^{*}$ 满足

$\{\begin{array}{l}F\left(H^{*},P^{*}\right)=0\\ G\left(H^{*},P^{*}\right)=0\end{array}$ .

对均匀态解 $\left(H^{*},P^{*}\right)$ 作一个微扰，令 $H=H^{*}+h$ ， $P=P^{*}+p$ ，代入方程(2)，经泰勒级数展开并去掉高阶项，可得线性微扰方程

$\{\begin{array}{l}\frac{\partial h}{\partial t}=a_{11}h+a_{12}p+d_1{\nabla }^{2}h\\ \frac{\partial p}{\partial t}=a_{21}h+a_{22}p+d_2{\nabla }^{2}p\end{array}$ (3)

其中

$a_{11}={\frac{\partial F}{\partial H}|}_{\left(H^{*},P^{*}\right)}=\frac{\left(\beta {\left(r_2/{\alpha }_2\right)}^2-{\alpha }_1\right)\cdot r_1}{{\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta \cdot r_2/{\alpha }_2}$ ,

$a_{12}={\frac{\partial F}{\partial P}|}_{\left(H^{*},P^{*}\right)}=\frac{\left(\mu -\text{2}{r}_2/{\alpha }_2\right)\cdot \beta \cdot r_1}{{\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta \cdot r_2/{\alpha }_2}$ ,

$a_{21}={\frac{\partial G}{\partial H}|}_{\left(H^{*},P^{*}\right)}=r_2^2/{\alpha }_2$ , $a_{22}={\frac{\partial G}{\partial P}|}_{\left(H^{*},P^{*}\right)}=-r_2$ .

将微扰变量在傅里叶空间展开，令

$\left(\begin{array}{l}h\\ p\end{array}\right)={\displaystyle \underset{\kappa }{\sum }\left(\begin{array}{l}{c}_{\kappa }^1\\ c_{\kappa }^2\end{array}\right){\text{e}}^{{\lambda }_{\kappa }t+i\kappa h}}$ (4)

代入微扰方程(3)得特征方程

${\lambda }_{\kappa }\left(\begin{array}{l}{c}_{\kappa }^1\\ c_{\kappa }^2\end{array}\right)=\left(\begin{array}{cc}{a}_{11}-{\kappa }^2{d}_1& a_{12}\\ a_{21}& a_{22}-{\kappa }^2{d}_2\end{array}\right)\left(\begin{array}{l}{c}_{\kappa }^1\\ c_{\kappa }^2\end{array}\right)$ (5)

解特征方程(5)得如下色散关系：

${\lambda }_{\kappa }^2-t{r}_{\kappa }{\lambda }_{\kappa }+{\Delta }_{\kappa }=0$

其中

$t{r}_{\kappa }=a_{11}+a_{22}-{\kappa }^2\left(d_1+d_2\right)=t{r}_0-{\kappa }^2\left(d_1+d_2\right)$

$\begin{array}{c}{\Delta }_{\kappa }=a_{11}{a}_{22}-a_{21}{a}_{12}-{\kappa }^2\left(a_{11}{d}_2+a_{22}{d}_1\right)+{\kappa }^4{d}_1{d}_2\\ ={\Delta }_0-{\kappa }^2\left(a_{11}{d}_2+a_{22}{d}_1\right)+{\kappa }^4{d}_1{d}_2\end{array}$ ,

解色散关系得

${\lambda }_{\kappa }=\frac{t{r}_{\kappa }\pm \sqrt{t{r}_{\kappa }^2-4{\Delta }_{\kappa }}}{2}$ .

由斑图动力学理论及(4)知，当所有的特征值 ${\lambda }_{\kappa }$ 的实部都小于0时，微扰随时间衰变为0，因而系统是稳定的。而当系统有一个特征值 ${\lambda }_{\kappa }>0$ 时，微扰变量将随时间不断增加，系统的均匀定态解 $\left(H^{*},P^{*}\right)$ 失稳。如果该反应扩散系统是有界的，由于系统的唯一均匀定态解失稳，它一定会向某种时空非均匀态转化，而这种转化必然伴随着系数的某种时空对称性破缺 [11] 。设 ${\lambda }_{\kappa }$ 是系统控制参量 $\mu$ 的函数： ${\lambda }_{\kappa }={\lambda }_{\kappa }\left(\mu \right)$ ，在临界点 ${\mu }_c$ ，有 ${\lambda }_{\kappa }\left({\mu }_c\right)=0$ ， $\frac{\partial {\lambda }_{\kappa }}{\partial \mu }\ne 0$ ，则 $\mu ={\mu }_c$ 是系统的一个动力学分岔点 [11] 。最常见的分岔有：跨临界分岔，Hopf分岔和图灵分岔。

下面，类似文献 [12] 的分析，给出三种分支临界线的表达：

1) 在跨临界分支中，非扩散模型的一个特征值消失。因此，解 ${\lambda }_1\cdot {\lambda }_2={\Delta }_{\kappa }=0$ ， $\kappa =0$ ，可得跨临界分支参数 $\mu$ 的临界值表达式：

$\mu ={\alpha }_1{\alpha }_2/\left(r_2\beta \right)+r_2/{\alpha }_2$ .

由于参数 $\mu$ 的取值范围为 $0\le \mu <{\alpha }_1{\alpha }_2/\left(r_2\beta \right)+r_2/{\alpha }_2$ ，故系统(3)不会出现跨临界分支。

2) Hopf分支不稳定的出现对应于下面情况：对于非扩散模型来说，当根从负到正跨越实轴时产生了一对纯虚根。从数学角度看，当满足下列条件：

$Im\left({\lambda }_{\kappa }\right)\ne 0$ , $Re\left({\lambda }_{\kappa }\right)=0$ .

在 $\kappa =0$ 系统产生了Hopf分支，即

$t{r}_0=a_{11}+a_{22}=0$ ,

$\left(\beta {\left(r_2/{\alpha }_2\right)}^2-{\alpha }_1\right)\cdot r_1/\left({\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta \cdot r_2/{\alpha }_2\right)-r_2=0$ .

方程两边同时乘以 ${\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta \cdot r_2/{\alpha }_2$ 得

$r_1\left(\beta {\left(r_2/{\alpha }_2\right)}^2-{\alpha }_1\right)-r_2\left({\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta \cdot r_2/{\alpha }_2\right)=0$ ,

$r_1\left(\beta {\left(r_2/{\alpha }_2\right)}^2-{\alpha }_1\right)-r_2\left({\alpha }_1+\beta \cdot {\left(r_2/{\alpha }_2\right)}^2\right)+(r_2\cdot \beta \cdot r_2/{\alpha }_2)\mu =0$ ,

即可得Holf分支参数 $\mu$ 的临界值表达式：

$\mu =r_1/r_2\cdot {\alpha }_1{\alpha }_2/\left(r_2\beta \right)+{\alpha }_1{\alpha }_2/\left(r_2\beta \right)+r_2/{\alpha }_2-r_1/{\alpha }_2$ .

3) 令

$B\left({\kappa }^2\right)=\left(\begin{array}{cc}{a}_{11}-{\kappa }^2{d}_1& a_{12}\\ a_{21}& a_{22}-{\kappa }^2{d}_2\end{array}\right)$ ,

则 $\det \left(B\left({\kappa }^2\right)\right)$ 是关于 ${\kappa }^2$ 的二次多项式，它在某些 ${\kappa }^2$ 值处取得极小值。首先假设非扩散模型在平衡点处是稳定的(特别的 $\det \left(B\left(0\right)\right)>0$ )，且图灵不稳定的发生需满足下面两个条件：(a) ${\kappa }^2>0$ (如果这个条件不成立，对于所有的 ${\kappa }^2=0$ ，由 $\det \left(B\left(0\right)\right)>0$ 可推出 $\det \left(B\left({\kappa }^2\right)\right)>0$ )是图灵不稳定的必要条件；(b) 图灵不稳定的充分条件是 $\det \left(B\left({\kappa }^2\right)\right)<0$ ，线性分析给出了图灵斑图产生的条件如下：

$t{r}_0=a_{11}+a_{22}<0$ , ${\Delta }_0=a_{11}{a}_{22}-a_{21}{a}_{12}>0$ ，

$a_{11}{d}_2+a_{22}{d}_1>0$ ， ${\left(a_{11}{d}_2+a_{22}{d}_1\right)}^2>4d_2{d}_1\left(a_{11}{a}_{22}-a_{21}{a}_{12}\right)$ ；

图灵分支产生的条件如下：

$Im\left({\lambda }_{\kappa }\right)=0$ , $Re\left({\lambda }_{\kappa }\right)=0$ , $\kappa ={\kappa }_T\ne 0$ ,

且波数 ${\kappa }_T$ 满足

${\kappa }_T^2=\sqrt{\frac{a_{11}{a}_{22}-a_{21}{a}_{12}}{d_1{d}_2}}$ .

由 ${\lambda }_{{\kappa }_T}=\frac{t{r}_{{\kappa }_T}\pm \sqrt{t{r}_{{\kappa }_T}^2-4{\Delta }_{{\kappa }_T}}}{2}$ 得

$\{\begin{array}{l}t{r}_{{\kappa }_T}=0\\ t{r}_{{\kappa }_T}^2-4{\Delta }_{{\kappa }_T}=0\end{array}$

其中，

$t{r}_{\kappa }=a_{11}+a_{22}-{\kappa }_T^2\left(d_1+d_2\right)$ ,

${\Delta }_{\kappa }=a_{11}{a}_{22}-a_{21}{a}_{12}-{\kappa }_T^2\left(a_{11}{d}_2+a_{22}{d}_1\right)+{\kappa }_T^4{d}_1{d}_2$ ,

即

$\{\begin{array}{l}{a}_{11}+a_{22}-{\kappa }_T^2\left(d_1+d_2\right)=0\\ a_{11}{a}_{22}-a_{21}{a}_{12}-{\kappa }_T^2\left(a_{11}{d}_2+a_{22}{d}_1\right)+{\kappa }_T^4{d}_1{d}_2=0\end{array}$ ,

整理得

$\{\begin{array}{l}{a}_{11}+a_{22}-{\kappa }_T^2\left(d_1+d_2\right)=0\left(6\right)\\ 2\left(a_{11}{a}_{22}-a_{21}{a}_{12}\right)-{\kappa }_T^2\left(a_{11}{d}_2+a_{22}{d}_1\right)=0\left(7\right)\end{array}$

由 $\left(7\right)\times \left(d_1+d_2\right)-\left(6\right)\times \left(a_{11}{d}_2+a_{22}{d}_1\right)$ 得

$2\left(a_{11}{a}_{22}-a_{21}{a}_{12}\right)\left(d_1+d_2\right)-\left(a_{11}+a_{22}\right)\left(a_{11}{d}_2+a_{22}{d}_1\right)=0$ ,

整理

$\left(a_{11}-a_{22}\right)\left(a_{22}{d}_1-a_{11}{d}_2\right)-2a_{12}{a}_{21}\left(d_1+d_2\right)=0$ ,

即

$\begin{array}{l}\left(\frac{\left(\beta {\left(r_2/{\alpha }_2\right)}^2-{\alpha }_1\right)\cdot r_1}{{\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta \cdot r_2/{\alpha }_2}+r_2\right)\left(-r_2\cdot d_1-\frac{\left(\beta {\left(r_2/{\alpha }_2\right)}^2-{\alpha }_1\right)\cdot r_1}{{\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta \cdot r_2/{\alpha }_2}\cdot d_2\right)\\ -2\frac{\left(\mu -\text{2}{r}_2/{\alpha }_2\right)\cdot \beta \cdot r_1}{{\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta \cdot r_2/{\alpha }_2}\cdot r_2^2/{\alpha }_2\cdot \left(d_1+d_2\right)=0\end{array}$ ,

整理得

$\begin{array}{l}\left(\beta {\left(r_2/{\alpha }_2\right)}^2-{\alpha }_1\right)\cdot r_2{d}_1+{\left(\beta {\left(r_2/{\alpha }_2\right)}^2-{\alpha }_1\right)}^2\cdot r_1/\left({\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta \cdot r_2/{\alpha }_2\right)\cdot d_2\\ +r_2^2\cdot d_1\left({\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta \cdot r_2/{\alpha }_2\right)/r_1+\left(\beta {\left(r_2/{\alpha }_2\right)}^2-{\alpha }_1\right)\cdot r_2{d}_2\\ +2\beta \left(\mu -2r_2/{\alpha }_2\right)\cdot r_2^2/{\alpha }_2\left(d_1+d_2\right)=0\end{array}$ ,

两边同时乘以 $\left({\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta \cdot r_2/{\alpha }_2\right)/\left(r_2{d}_2\right)$ 得

$\begin{array}{l}\left(\beta {\left(r_2/{\alpha }_2\right)}^2-{\alpha }_1\right)\cdot \left({\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta \cdot r_2/{\alpha }_2\right)\cdot d_1/d_2\\ +{\left(\beta {\left(r_2/{\alpha }_2\right)}^2-{\alpha }_1\right)}^2\cdot r_1/r_2+r_2/r_1\cdot d_1/d_2\cdot {\left({\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta \cdot r_2/{\alpha }_2\right)}^2\\ +\left(\beta {\left(r_2/{\alpha }_2\right)}^2-{\alpha }_1\right)\cdot \left({\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta \cdot r_2/{\alpha }_2\right)\\ +2\beta \left(\mu -2r_2/{\alpha }_2\right)\cdot r_2/{\alpha }_2\cdot \left(1+d_1/d_2\right)\left({\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta \cdot r_2/{\alpha }_2\right)=0\end{array}$

令

$A=\beta \cdot {\left(r_2/{\alpha }_2\right)}^2-{\alpha }_1$ , $D=d_1/d_2$ , $r=r_1/r_2$ ,

上式可变为

$\begin{array}{l}\left(1+D\right){\alpha }_1-\mu \beta r_2/{\alpha }_2\cdot \left(1+D\right)+\beta {\left(r_2/{\alpha }_2\right)}^2\left(1+D\right)+Ar\\ +D/r\cdot \left({\alpha }_1^2-2{\alpha }_1\left(\mu -r_2/{\alpha }_2\right)\beta r_2/{\alpha }_2+{\left(\mu -r_2/{\alpha }_2\right)}^2{\beta }^2{\left(r_2/{\alpha }_2\right)}^2\right)\\ +2\mu \beta r_2/{\alpha }_2\cdot \left(\left(1+D\right)/A\right)\left({\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta r_2/{\alpha }_2\right)\\ -4\beta {\left(r_2/{\alpha }_2\right)}^2\left(\left(1+D\right)/A\right)\left({\alpha }_1-\left(\mu -r_2/{\alpha }_2\right)\beta r_2/{\alpha }_2\right)=0\end{array}$ .

进一步，令

$a={\beta }^2\cdot {\left(r_2/{\alpha }_2\right)}^2\cdot \left(D/r-2\left(1+D\right)/A\right)$ ,

$b=-\beta \left(r_2/{\alpha }_2\right)\cdot \left(1+D\right)+2{\alpha }_1\beta \cdot r_2/{\alpha }_2\cdot \left(\left(1+D\right)/A-D/r\right)+2{\beta }^2{\left(r_2/{\alpha }_2\right)}^3\cdot \left(3\left(1+D\right)/A-D/r\right)$ ， $c=\left(1+D\right)\cdot \left({\alpha }_1+\beta {\left(r_2/{\alpha }_2\right)}^2\right)+Ar+D/r\cdot {\left(\left({\alpha }_1+\beta {\left(r_2/{\alpha }_2\right)}^2\right)\right)}^2+4\beta {\left(r_2/{\alpha }_2\right)}^2\cdot \left(1+D\right)$ .

整理得

$a{\mu }^2+b\mu +c=0$ ,

由此得图灵分支参数 $\mu$ 的临界表达式为

$\mu =\frac{-b\pm \sqrt{b^2-4ac}}{2a}$ .

根据 $\mu$ 的取值范围，上述两个值里面取正，且要保证 $b^2-4ac\ge 0$ 。

下面我们将利用以上Hopf分支和Turing分支的临界表达式，确定出图灵斑图出现的区域，即Turing区域。由区域图可以清晰的看出系统的稳定与失稳区域。

在图1中，红线和黑线分别表Hopf和Turing分支线，参数空间中区域I表示稳定区域；区域II表示Hopf没有失稳，Turing失稳区域，即Turing区域；区域III表示Hopf和Turing同时失稳区域。参数取为 $r_1=1$ ， $r_2=0.5$ ， ${\alpha }_1=0.9$ ， ${\alpha }_2=0.4$ ， $d_1=0.1$ ， $d_2=1.2$ 。

3. 数值模拟

一个具体的系统将会出现哪种斑图结构是斑图选择的问题，它基于比较难的非线性分析理论，故使用分析法或一般方法对空间模型的动力学性态进行研究是比较困难的。而借助于计算机进行数值模拟是解决该问题的一个有效方法。用计算机求解一个复杂的方程组，事实上就是将连续问题转化为离散问题。二维空间中定义的连续反应扩散系统的求解问题可以转变为由 $m\times n$ 格子组成的离散区域问题来进行求解。每两个格子间的空间距离定义为步长 $\Delta h$ 。在离散系统中，描述扩散的拉普拉斯算子 ${\nabla }^2={\partial }^2/\partial x^2+{\partial }^2/\partial y^2$ 用有限差分法进行离散，即导数是 $\Delta h$ 上的差分逼近，当 $\Delta h\to 0$ 时，差分与求导近似相等。同样时间演化也离散化，时间步长为 $\Delta t$ ，可以用Euler法解决，下一个时间步长的状态依赖于前一个时间步长。

下面我们首先给出空间模型离散化后的差分格式，并进行数值模拟。模拟过程中采用零流边界条件(Neumann边界条件)，离散区域为 $120\times 120$ 个格子。此外，对于空间模型的数值模拟还依赖于下面两个初值条件：一种是种群密度随机分布，另一种是种群在水平空间呈不同的密度分布，而在垂直空间呈同一密度分布。在这里，我们选择第一种情况，种群密度随机分布。首先利用有限差分法给出模型(3)的差分格式：

$\{\begin{array}{l}{H}^{k+1}\left(i,j\right)=H^k\left(i,j\right)+\Delta t\cdot F_1\left(H,P\right)+\frac{\Delta t\cdot d_1}{{\left(\Delta h\right)}^2}\left[H^k\left(i+1,j\right)+H^k\left(i-1,j\right)+H^k\left(i,j+1\right)+H^k\left(i,j-1\right)-4H^k\left(i,j\right)\right]\\ P^{k+1}\left(i,j\right)=P^k\left(i,j\right)+\Delta t\cdot F_2\left(H,P\right)+\frac{\Delta t\cdot d_2}{{\left(\Delta h\right)}^2}\left[P^k\left(i+1,j\right)+P^k\left(i-1,j\right)+P^k\left(i,j+1\right)+P^k\left(i,j-1\right)-4P^k\left(i,j\right)\right]\end{array}$

Figure 1. The parameter space of bifurcation

图1. 系统(3)的分支参数图

Figure 2. The time evolution of spatial pattern for parasitism

图2. 寄生物在空间分布的时间演化图

其中

$F_1\left(H,P\right)=H^k\left(i,j\right)\cdot \left[r_1-{\alpha }_1{H}^k\left(i,j\right)+\left[\mu -\frac{P^k\left(i,j\right)}{H^k\left(i,j\right)}\right]\beta P^k\left(i,j\right)\right]$ . $F_2\left(H,P\right)=P^k\left(i,j\right)\cdot \left[r_2-{\alpha }_2\frac{P^k\left(i,j\right)}{H^k\left(i,j\right)}\right]$

零流边界条件离散为：

$H_n\left(x_0,y_j\right)=H_n\left(x_1,y_j\right),H_n\left(x_{N+1},y_j\right)=H_n\left(x_N,y_j\right)$ ,

$H_n\left(x_i,y_0\right)=H_n\left(x_i,y_1\right),H_n\left(x_i,y_{N+1}\right)=H_n\left(x_i,y_N\right)$ ,

$P_n\left(x_0,y_j\right)=P_n\left(x_1,y_j\right),P_n\left(x_{N+1},y_j\right)=P_n\left(x_N,y_j\right)$ ,

$P_n\left(x_i,y_0\right)=P_n\left(x_i,y_1\right),P_n\left(x_i,y_{N+1}\right)=P_n\left(x_i,y_N\right)$ .

利用上述差分格式进行数值模拟，选取空间步长 $\Delta h=1$ ，时间步长 $\Delta t=0.01$ 。由于寄生物与宿主的空间分布很类似，在这里我们只给出寄生物的空间斑图。各图表示寄生物种群在空间状态下的分布，图中各颜色表示寄生物种群的密度(数量)，不同的空间位置种群的密度不同，即图中呈现出不同的颜色。

图2给出了在迭代0，1000，10,000，100,000次时寄生物的斑图的空间演化，从图中可以看出初始不规则的斑图最终演化为规则的黑眼点状斑图，且动力学性态不再发生任何改变。模拟过程初始条件为正平衡点的小扰动，参数选取为 $r_1=1$ ； $r_2=0.5$ ； ${\alpha }_1=0.9$ ； ${\alpha }_2=0.4$ ； $\beta =1$ ； $\mu =0.91$ 。

基金项目

国家自然科学基金(31600299, 31560134)；陕西省自然科学基础研究计划项目(2017JQ3020)；陕西省高校科协青年人才托举基金(20160234)；宝鸡文理学院重点项目(ZK2017021)。