具交叉扩散和弱Allee效应的反应扩散方程图灵斑图研究

Turing Pattern Study of Reaction Diffusion Equations with Cross Diffusion and Weak Allee Effect

摘要: 本文提出了一类具有交叉扩散的反应扩散方程，用于研究一类具有弱Allee效应的捕食食饵反应扩散模型在交叉扩散驱动下的图灵斑图动力学现象。研究发现，交叉扩散是该模型产生图灵斑图必不可少的条件，如果没有交叉扩散，则模型不会发生图灵失稳。同时，通过选择不同参数，我们发现模型的斑图现象十分丰富，可以出现点斑图、点–线混合斑图、线斑图三种类型。

Abstract: This article proposes a class of reaction diffusion equations with cross diffusion to study the Turing pattern dynamics of a predator-prey reaction diffusion model with weak Allee effect driven by cross diffusion. Research has found that cross diffusion is an essential condition for the model to generate Turing patterns. Without cross diffusion, the model will not experience Turing instability. At the same time, by selecting different parameters, we found that the pattern phenomenon of the model is very rich, which can appear in three types: spots pattern, spots-stripes mixed pattern and stripes pattern.

1. 引言

在自然界中，许多生物总是混杂在一起，它们之间的关系非常复杂，包括捕食、寄生、共存和竞争等，并将形成一个庞大而复杂的食物网。其中，捕食关系最为常见。如果一种生物以另一种为食，我们称前者为“掠食者”，后者为“猎物”。比如，蛇–鹰，斑马–狮子，虫–鸟等就是典型的捕食者与猎物关系。1926年，Lotka [1] 和Volterra [2] 提出了著名的Lotka-Volterra捕食–被捕食模型，该模型是微分方程理论在生物科学中应用的重要成果。尽管该模型本身比较简单，不适用于现实生活中的许多问题，但它促进了生物数学的发展。如今，许多学者 [3] [4] [5] [6] [7] 研究捕食和被捕食模型的相关问题，这些模型最初是从Lotka-Volterra模型慢慢演化而来的。正如数学家和生物学家所探索的，他们提出了两种功能性反应：猎物依赖和捕食者依赖，这可以用来更准确地描述捕食者和猎物之间的关系。在食饵依赖型中，Holling-II功能型反应 $\frac{au}{1+mu}$ 已被广泛研究。其中，张 [8] 研究了如下具有的第二类功能反应函数的捕食食饵模型：

$\{\begin{array}{l}{u}^{\prime }=ru\left(1-\frac{u}{K}\right)-\frac{au}{1+mu}-\alpha uv\\ v^{\prime }=\beta uv-dv\end{array}$ (1)

其中， $u,v$ 分别代表食饵和捕食者的密度， $r,K,a,\alpha ,\beta ,d$ 都是正参数，分别代表食饵的内禀增长率、环境最大容纳量、捕食者的攻击速率、消耗率、转化率、死亡率。张 [8] 研究该模型并将它分为四个参数区域：强Allee区、弱Allee区、无Allee区和灭绝区。其中，弱Allee区的参数条件满足：

$\mathcal{D}=\left\{\vartheta =\left(a,r,m,K\right)\in R_{+}^4:armK>0,\frac{1}{mK}<\frac{a}{r}<1\right\}$ (2)

然而，上面描述的ODE模型只反映了种群随时间的演化，但在现实世界中，有一个绝对不能忽视的因素——空间。例如，近年来，新型冠状病毒肆虐。由于病毒很容易通过空气传播，如果数学家和生物学家在研究病毒时只考虑时间因素，没有考虑病毒在空间传播造成的空间异位，那么研究所得到的结果就不太真实。同样，在动物世界中，总是存在着空间位置的变化，为了生存、为了捕食、为了逃避捕食、为了躲避自然灾害，等等，他们必须有空间上的扩散。因此，在研究问题时，最好考虑空间因素，在捕食–被捕食ODE模型中引入所谓的“扩散项”，使研究结果更加准确。这样所建立的方程称为反应扩散方程，近年来受到国内外学者的广泛关注。例如，Subrata Dey等 [9] 研究了一个具有广义捕食者和自由扩散的捕食者–食饵系统：Arancibia-Ibarra，Claudioet等 [10] 研究了一个扩散Holling-Tanner捕食者–食饵模型，为捕食者提供了一个替代的食物来源：Hu等 [11] 研究了一个强耦合的扩散捕食者–食饵系统：Tang和Song [12] 研究了一个具有群集行为和交叉扩散的捕食者–食饵模型：

$\{\begin{array}{ll}{\partial }_{t}u=d_{11}\Delta u+d_{12}\Delta v+u\left(1-u\right)-\sqrt{u}v\hfill & x\in \Omega ,t>0\hfill \\ {\partial }_{t}v=d_{21}\Delta u+d_{22}\Delta v+\gamma v\left(-\beta +\sqrt{u}\right)\hfill & x\in \Omega ,t>0\hfill \\ {\partial }_{\nu }u={\partial }_{\nu }v=0\hfill & x\in \partial \Omega \hfill \\ u\left(x,0\right)=u_0\left(x\right)\ge 0,v\left(x,0\right)=v_0\left(x\right)\ge 0\hfill & x\in \Omega \hfill \end{array}$ (3)

这里， $d_{11},d_{12},d_{21},d_{22}>0$ 。 $d_{11},d_{22}$ 表示自由扩散系数，它表示物种因自身原因的扩散速率， $d_{12},d_{21}$ 代表交叉扩散系数，它代表物种因另一物种而交叉形成的扩散速率。因此，受到张 [8] 和Tang [12] 的启发，我们提出了如下的模型：

$\{\begin{array}{ll}\frac{\partial u}{\partial t}=d_{11}\Delta u+d_{12}\Delta v+ru\left(1-\frac{u}{K}\right)-\frac{au}{1+mu}-\alpha uv\hfill & x\in \Omega ,t>0\hfill \\ \frac{\partial v}{\partial t}=d_{21}\Delta u+d_{22}\Delta v+\beta uv-dv\hfill & x\in \Omega ,t>0\hfill \\ \frac{\partial u}{\partial \nu }=\frac{\partial v}{\partial \nu }=0\hfill & x\in \partial \Omega \hfill \\ u\left(x,0\right)=u_0\left(x\right)\ge 0,v\left(x,0\right)=v_0\left(x\right)\ge 0\hfill & x\in \Omega \hfill \end{array}$ (4)

这里，我们将区域Ω限制至二维空间 $\Omega =\left[0,L\right]\times \left[0,L\right]$ ： $\Delta =\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}$ 是拉普拉斯算子： $\frac{\partial u}{\partial \nu }=\frac{\partial v}{\partial \nu }=0$ 指Neuman边界条件，代表着物种不会跑出这个边界： $u\left(x,0\right)=u_0\left(x\right)\ge 0$ ， $v\left(x,0\right)=v_0\left(x\right)\ge 0$ 则代表着物种正的初值，符合现实意义。

本文中，我们将在第二部分简述ODE模型(1)的一些结论，并分析PDE模型(4)在弱Allee条件(2)下发生图灵失稳的条件，第三部分则是在这个图灵失稳条件下进行数值模拟，展现丰富的图灵斑图现象，第四部分对模型提出总结与展望。

2. 理论分析

2.1. ODE模型稳定性分析

为了方便讨论问题，我们定义如下的非线性函数

$H\left(u\right)=r\left(1-\frac{u}{K}\right)-\frac{a}{1+mu}=\frac{h\left(u\right)}{K\left(1+mu\right)}$ (5)

这里，

$h\left(u\right)=-rm{u}^2+r\left(mK-1\right)u+\left(r-a\right)K$ (6)

容易知道，当 $\vartheta =\left(a,r,m,K\right)\in \mathcal{D}$ 时，我们有 $mK-1>0$ ， $\left(r-a\right)K>0$ ，因而，我们知道 $h\left(u\right)$ 是一个经典的二次函数，在第一象限是先单调递增再单调递减， $h\left(0\right)=\left(r-a\right)K>0$ ，并且有且仅有一个正根p满足 $h\left(p\right)=0$ ，因此， $H\left(u\right)$ 在第一象限也仅有一个正根p。

由(5)，我们可以将方程(1)简写成：

$\{\begin{array}{l}{u}^{\prime }=uH\left(u\right)-\alpha uv\\ v^{\prime }=\beta uv-dv\end{array}$ (7)

为了研究(7)的平衡点稳定性，我们必须深入研究非线性函数 $H\left(u\right)$ 的性态。对 $H\left(u\right)$ 求导得：

$H^{\prime }\left(u\right)=\frac{am}{{\left(1+mu\right)}^2}-\frac{r}{K}$ (8)

令 $H^{\prime }\left(u\right)=0$ ，则有 $\left|1+mu\right|=\sqrt{\frac{amK}{r}}>1$ ，因而， $H^{\prime }\left(u\right)=0$ 有且仅有唯一的正根，记为

$b=\frac{\sqrt{\frac{amK}{r}}-1}{m}$ (9)

综上所述，我们知道函数 $H\left(u\right)$ 在 $\left(0,p\right)$ 之间大于0，在 $\left(p,+\infty \right)$ 之间小于0，并且 $H\left(u\right)$ 在 $\left(0,b\right)$ 之间单调递增，在 $\left(b,+\infty \right)$ 之间单调递减。诸如，当我们取 $a=1$ ， $r=2$ ， $m=2$ ， $K=2$ ，则 $b=0.2071$ ， $p=1.7808$ ，那么 $H\left(u\right)$ 的函数图像如图1所示。

由张 [8] 对于(7)的弱Allee部分的平衡点的分析，我们知道，模型(7)永远有两个边界平衡点： $\left(0,0\right)$ ， $\left(p,0\right)$ 。而(7)的正平衡点 $E^{\ast }\left(u^{\ast },v^{\ast }\right)$ ， $u^{\ast }=\frac{d}{\beta }$ ， $v^{\ast }=\frac{H\left(\frac{d}{\beta }\right)}{\alpha }$ 只有在当 $0<\frac{d}{\beta }<p$ 时才会存在。下面的定理揭示了正平衡点 $E^{\ast }\left(u^{\ast },v^{\ast }\right)$ 的稳定性。

Figure 1. Example of $H\left(u\right)$ image

图1. $H\left(u\right)$ 图像示例

定理1.

1) $E^{\ast }\left(u^{\ast },v^{\ast }\right)$ 在 $0<u^{\ast }<b$ 时是不稳定的，在 $b<u^{\ast }<p$ 时是稳定的。

2) 当 $u^{\ast }=b$ 时，模型(7)发生Hopf分岔。

证明：

1) 容易知道，模型(7)在 $E^{\ast }\left(u^{\ast },v^{\ast }\right)$ 的雅可比矩阵为：

$J=\left(\begin{array}{cc}{j}_{11}& j_{12}\\ j_{21}& j_{22}\end{array}\right)=\left(\begin{array}{cc}{u}^{\ast }{H}^{\prime }\left(u^{\ast }\right)& -\alpha u^{\ast }\\ \beta v^{\ast }& 0\end{array}\right)$ (10)

那么，J的特征方程为：

${\lambda }^2-T\lambda +D=0$ (11)

其中，T，D分别是J的迹和行列式，表达式为：

$T=u^{\ast }{H}^{\prime }\left(u^{\ast }\right),D=\alpha \beta u^{\ast }{v}^{\ast }>0$ (12)

由上面对于 $H^{\prime }\left(u^{\ast }\right)$ 的分析，我们知道当 $0<u^{\ast }<b$ ， $H^{\prime }\left(u^{\ast }\right)>0$ ， $T\left(J\right)>0$ ；当 $b<u^{\ast }<p$ ， $H^{\prime }\left(u^{\ast }\right)<0$ ， $T\left(J\right)<0$ 可以立即得到结论。

2) 当 $u^{*}=b$ ，我们知道 $T\left(J\right)=u^{\ast }{H}^{\prime }\left(u^{\ast }\right)=0$ ，并且横截条件为

${\frac{\partial \left(u^{\ast }{H}^{\prime }\left(u^{\ast }\right)\right)}{\partial u^{\ast }}|}_{u^{\ast }=b}={\left(u^{\ast }{H}^{″}\left(u^{\ast }\right)+H^{\prime }\left(u^{\ast }\right)\right)|}_{u^{\text{*}}=b}<0$ (13)

由著名的Poincare-Andronov-Hopf分支定理，结论显然成立。

2.2. PDE模型稳定性分析

为了得到在 $E^{\ast }\left(u^{\ast },v^{\ast }\right)$ 发生图灵失稳的参数条件，我们总是假定 $E^{\ast }\left(u^{\ast },v^{\ast }\right)$ 在ODE模型(7)中总是存在并且稳定的，即下列条件总是成立的：

$b<u^{*}<p$ (14)

现在假设

$\left(\begin{array}{c}u\\ v\end{array}\right)=\left(\begin{array}{c}{u}^{\ast }\\ v^{\ast }\end{array}\right)+\left(\begin{array}{c}{u}_k\\ v_k\end{array}\right){\text{e}}^{\lambda t+i\left(k\cdot r\right)}$ (15)

是模型(4)在 $E^{\ast }\left(u^{\ast },v^{\ast }\right)$ 的线性化系统的非常数解。这里， $k=\left(k_x,k_y\right)$ 且 $k=\sqrt{k_x^2+k_y^2}$ 也叫做波数， $r=\left(x,y\right)$ 是二维空间行向量， $x,y$ 是空间变量。则模型(4)的线性化系统的特征方程为：

${\lambda }^2-T_k\lambda +D_k=0$ (16)

这里，

$T_k=T-k^2\left(d_{11}+d_{22}\right)$ (17)

$D_k=\mathcal{M}{k}^4+\mathcal{N}{k}^2+D$ (18)

这里：

$\mathcal{M}=d_{11}{d}_{22}-d_{12}{d}_{21},\mathcal{N}=d_{12}{j}_{21}+d_{21}{j}_{12}-d_{22}{j}_{11}-d_{11}{j}_{22}$ (19)

容易知道，当 $b<u^{*}<p$ 时， $T<0$ ，故而 $T_k<0,\forall k$。因此，能够使得发生图灵失稳现象的只能是当对于适当的k，有 $D_k<0$ 发生的情况。现记 $\mu =k^2>0$ 并定义如下函数：

ℊ $\left(\mu \right)=\mathcal{M}{\mu }^2+\mathcal{N}\mu +D$ (20)

为了得到图灵失稳的条件，我们需要寻找恰当的 $\mathcal{M},\mathcal{N}$ 使得ℊ $\left(\mu \right)$ 在某些 $\mu >0$ 上小于0，由于 $D>0$ (见(12))，通过对二元函数ℊ $\left(\mu \right)$ 结合图像分析，我们很容易知道当满足以下任一条件时，那么“ℊ $\left(\mu \right)$ 在某些 $\mu >0$ 上小于0”：

(i) $\mathcal{M}<0$ ；或(ii) $\mathcal{M}\ge 0$ ， $\mathcal{N}<0$ 且 ${\mathcal{N}}^2-4\mathcal{M}D>0$ (21)

综上所述，模型(4)发生图灵失稳的参数条件为：(14) & (21)。

现在假设交叉扩散 $d_{12}=d_{21}=0$ ，那么 $\mathcal{M}>0$ ， $\mathcal{N}=-d_{22}{j}_{11}-d_{11}{j}_{22}=-d_{22}{u}^{\ast }{H}^{\prime }\left(u^{\ast }\right)>0$ ，依据判别条件(21)，我们知道，如果没有交叉扩散系数，则有 $T_k<0,D_k>0$ 恒成立，换句话说，模型(4)将不会发生图灵失稳现象。因此，在该模型中，交叉扩散系数至关重要。

3. 数值模拟

现在，我们对模型(4)在条件(14) & (21)下进行数值模拟。通过如下的五点差分格式

$\Delta u_{ij}^k=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}=\frac{u_{i,j+1}^k+u_{i,j-1}^k+u_{i+1,j}^k+u_{i-1,j}^k-4u_{ij}^k}{h^2}$ (22)

可以实现对拉普拉斯算子做离散化，这里h指的是空间步长。再者，通过设置在正平衡点附近的随机小扰动作为初值。鉴于捕食者与食饵所形成的图灵斑图是一致的，因而我们仅展示食饵的斑图情况。以下三个例子中，我们用Matlab进行数值模拟，并将网格固定为50 × 50，空间步长为1，时间步长为0.01，演变时间为50,000。

例1：如果我们选择参数 $a=1$ ， $r=2$ ， $m=2$ ， $K=2$ ， $\alpha =1$ ， $\beta =1$ ，并且 $d=1\in \left(b,p\right)=\left(0.2071,1.7808\right)$ 满足条件(14)。至于扩散系数，我们选取 $d_{11}=d_{22}=d_{12}=1$ 。此时，若选取 $d_{21}=2$ ，则 $\mathcal{M}=-1<0$ 满足条件(21)的(i)，此时，可以得到如下图2线斑图。我们可以看到，食饵在此时以线的形式分布着，并且是交叉纵横的。在线上，食饵的密度高。但在线外，食饵的密度低。

Figure 2. Stripes pattern

图2. 线斑图

例2.如果我们选择参数 $a=1$ ， $r=2$ ， $m=2$ ， $K=2$ ， $\alpha =1$ ， $\beta =1$ ，并且 $d=1\in \left(b,p\right)=\left(0.2071,1.7808\right)$ 满足条件(14)。至于扩散系数，我们选取 $d_{11}=3$ ， $d_{12}=1$ ， $d_{21}=6$ ， $d_{22}=3$ ，此时， $\mathcal{M}=3>0$ ， $\mathcal{N}=-3<0$ ， ${\mathcal{N}}^2-4\mathcal{M}D=5>0$ ，符合条件(21)的(ii)。此时，图灵斑图如图3所示。可以看到，当时间较大时，不同于例1的是，食饵不再以线形分布，而是聚集成了一个个的斑点，并且越靠近斑心中心，种群密度就越高，越远离斑心中心则密度越低。这种斑图称为“点斑图”。

例3.如果我们选择参数 $a=1$ ， $r=2$ ， $m=2$ ， $K=2$ ， $\alpha =1$ ， $\beta =1$ ，并且 $d=1\in \left(b,p\right)=\left(0.2071,1.7808\right)$ 满足条件(14)。至于扩散系数，我们选取 $d_{11}=3$ ， $d_{12}=1$ ， $d_{21}=7$ ， $d_{22}=3$ ，此时， $\mathcal{M}=2>0$ ， $\mathcal{N}=-4<0$ ， ${\mathcal{N}}^2-4\mathcal{M}D=13.333>0$ ，符合条件(21)的(ii)。此时，图灵斑图如图4所示。可以看到，当时间较大时，不同于例1，2的是，食饵不仅仅是单纯以圆点状或线型分布，而是圆点与线互相竞争，最终点、线共同存在。这叫作“点–线”混合斑图。

Figure 3. Spots pattern

图3. 点斑图

Figure 4. Mixed pattern

图4. 混合斑图

4. 总结展望

本文通过对一个具有交叉扩散和弱Allee效应的反应扩散方程捕食食饵模型进行研究，我们发现交叉扩散系数是产生图灵失稳必不可少的条件。同时。通过选取不同参数，可以生成丰富多彩的图灵斑图(Turing pattern)，包括点斑图、线斑图和点–线混合斑图，这反映了在不同的参数条件下种群在空间上的分布情况，对种群生态学的认识、研究与防控是具有一定的意义的。

值得展望的是，本文只讨论了弱Allee效应情况下的斑图情况，但对于强Allee效应、无Allee效应的情况又有什么情况发生？有什么异同点？这是值得后续进行研究的。此外，如果参数选择得当，还可以出现螺旋波斑图等更加美丽的斑图现象，这同样也值得关注。