交叉扩散驱动的SI 模型空间斑图*

陆源杉 肖敏† 万佑红 丁洁 蒋海军

1) (南京邮电大学自动化学院,人工智能学院,南京 210023)

2) (新疆大学数学与系统科学学院,乌鲁木齐 830047)

目前国内外关于SI 模型空间格局的研究大多数局限在自扩散以及系统参数对斑图模式的影响,而关于交叉扩散对空间格局的演化机理研究成果较少.本文建立了一个具有自扩散和交叉扩散的空间流行病模型,研究了在有无自扩散驱动系统不稳定的情况下,交叉扩散对SI 模型的稳定性、稳定速度以及斑图结构的影响.研究发现,在无自扩散驱动系统不稳定的情况下,引入交叉扩散能够激发Turing 斑图的产生;在自扩散驱动系统不稳定的情况下,交叉扩散可以实现斑图结构的改变;对于SI 模型的稳定速度,不论有无自扩散驱动系统不稳定,交叉扩散都影响了其到达稳定所需时间,且在不同的交叉扩散系数下,所需时间也不同.因此,交叉扩散对于SI 模型的稳定性、稳定速度、斑图结构都有重要的影响.

1 引言

1952 年,Turing[1]发现在反应扩散系统中,在没有扩散的情况下对小时间扰动稳定的齐次稳态在有扩散项的情况下变得不稳定,这种现象被称作Turing 不稳定,这种不稳定在空间域中产生了一种称为Turing 斑图的模式.在这一发现之后,许多研究者受到启发,将Turing 斑图应用于物理、化学、生物和生态过程等各个领域[2-8].

反应扩散系统是空间运动建模的重要工具之一,已广泛应用于化学反应[9]、谣言传播[10]、图像处理[11]、流行病传播[12]等许多空间分布系统.反应扩散过程提供了由于局部反应和扩散运动而发生的多种物质的浓度变化的最佳模型.以疾病传播为例,将不同位置的个体划分为不同的种类,表示个体相对于模型疾病的易感状态、感染状态.局部反应过程代表同一地点的个体相互接触并根据感染动态改变状态.扩散过程代表空间上个体的自由扩散,这可能是疾病从一个区域扩散到另一个区域的原因.

在流行病学领域,研究人员提出了反应扩散模型,Wang 等[13]研究了包含人口与流行病过程的反应扩散模型以捕捉环境的空间异质性和个体的移动对疾病动态的影响.阮中远[14]介绍了在流行病学领域基于复杂网络的反应扩散模型,以此加深人们对流行病传播的理解.Wang 等[15]研究了基于网络化亚种群的空间流行病的建模,促进了人们对突发疾病大规模传播动力学的认识.空间格局是疾病传播的特征,可以预测空间中的流行病传播动态,为研究和管理流行病的空间动态提供了重要的理论依据,进而指导政策制定[16].

大多数反应扩散流行病模型考虑了传染病病原体在空间中的自发运动或扩散过程,也就是说在这些模型中只考虑到了自扩散项存在的情况[17-21].比如Wang 等[22]研究了具有饱和感染力的反应扩散流行病模型自扩散系数对Turing 斑图的影响.Sun 等[23]提出了具有非线性发病率的流行病模型的感染率β 不同对空间格局的影响.

交叉扩散是指考虑不同物种或类型之间相互影响的扩散过程,即一个隔间的浓度梯度引起另一个隔间的流通,而不是仅存在隔间本身的流通.它首先由Kerner[24]提出.以流行病学领域为例,在现实生活中易感者有能力识别感染者并且选择远离或接近他们[25,26],随后研究人员对流行病学模型的空间格局进行了一系列的研究.Aly 等[27]研究了一种具有自扩散和交叉扩散项SIS 模型的稳定性问题,并通过数值模拟研究了Turing 模式.Triska等[28]提出了一种具有负值的感染者交叉扩散SI模型,并研究了负值感染者交叉扩散对斑图种类的影响.

然而,当前针对交叉扩散对空间流行病模型的研究大多都局限在Turing 不稳定以及对斑图数值仿真的初步结果,但对于交叉扩散对空间流行病模型的稳定性、稳定速度以及斑图结构的影响机制研究相对较少.为了更好地了解这些机制,本文研究空间流行病SI 模型的空间格局的形成等问题.

本文的结构组织如下: 第2 节提出了一个具有自扩散和交叉扩散项的空间流行病模型;第3 节分析了该模型发生Turing 不稳定的条件;第4 节通过一系列的斑图仿真模拟,阐明了交叉扩散对系统的稳定性、稳定速度以及斑图结构的影响;第5 节对本文工作进行总结并展望未来的工作.

2 模型介绍

常见的空间流行病模型有SIS 模型,SIR 模型和SI 模型,这些模型被广泛应用于描述和分析传染病在人群中的传播过程,对理解和控制传染病的传播具有重要意义.SIS 模型将人群分为两类: 易感人群(S)与感染人群(I),该模型假设易感人群(S)在被感染为感染人群(I)后能够恢复成易感人群(S),通常用于研究一些康复后不会获得长期免疫的疾病,如普通感冒、疟疾等[29].SIR 模型包含3 种人群: 易感人群(S)、感染人群(I)、康复人群(R),与SIS 模型不同的是,SIR 模型假设感染人群能够恢复成康复人群,即从感染人群恢复并具有永久免疫的人群,该模型适用于研究感染后能够获得长期免疫的疾病,如白日咳、麻疹[30].

本文考虑空间流行病模型中的一种基本模型:SI 模型.一般来说,SI 模型的病原体活跃的人群包括两个亚群体: 一类是易感人群(S),另一类是感染人群(I).其中S 和I 都是关于时间与空间的函数,它们都可以将病毒传播给健康个体.与SIR 或SIS 模型不同,SI 模型假设易感人群一旦被感染,他们就会持续处于感染状态,不会恢复成易感状态,该模型一般用于研究一旦感染后就永久带病原体的疾病,许多性传播疾病,如艾滋病,都遵循SI控制的流行病学[31].

文献[23]考虑了以下具有非线性发病率的空间流行病SI 模型:

其中,d为被感染者的疾病相关死亡率;A是人口的吸纳率;μ为疑似感染者和被感染者的自然死亡率;βSpIq是非线性的发病率,相比于双线性发病率βSI具有更丰富的动力学行为[32],其中β 是传播率,本文假设p=1,q=2;∇2是二维空间中的拉普拉斯算子,x,y表示空间;DS与DI分别表示易感个体和感染个体的自扩散系数.

本文在此研究基础上引入交叉扩散.模型如下:

其中,D1为易感人群交叉扩散系数,D1∇2I代表感染人群对易感人群扩散的影响;D2为感染人群的交叉扩散系数,D2∇2S代表易感人群对感染人群扩散的影响.交叉扩散系数可以为正、零或负[17].易感人群交叉扩散系数D1为正值,表示易感人群向感染人群密度较低的方向移动(易感人群移动后仍保持易感状态),即易感人群倾向于远离感染人群.感染者的交叉扩散系数D2为正值,表示感染人群向易感人群密度低的方向移动(感染人群移动后仍保持感染状态),即感染人群倾向于远离易感人群.模型(2)需要用初始种群进行分析:

3 稳定性分析

为了找出反应扩散系统(2)发生Turing 不稳定的条件,首先研究无扩散系统平衡点的稳定性.当无扩散时,即D1=D2=DS=DI=0,系统(2)退化为

通过计算可得系统(3)有3 个平衡点:

E3对应于流行病灭绝的情况,并且E2是不稳定的[23],所以本文只研究平衡点E1处的稳定性,记E1=(S1,I1).系统(3)在平衡点E1处的雅克比矩阵如下:

其中,fq与gq(q=S,I) 分别是f(S,I),g(S,I) 在平衡点E1处关于q的偏导数.系统(3)的特征方程为

其中,m0=-(fS+gI),h0=fSgI-fIgS.

做如下假设:

定理1若(H1) 成立,则系统(3)在平衡点E1处局部渐进稳定.

证明若(H1) 成立,根据赫尔维茨判据可得特征方程(4)所有特征根实部均小于0,即系统(3)在平衡点E1处渐近稳定.

注1Turing 不稳定本质是由扩散引起,因此下文对于含扩散系统的Turing 不稳定分析均建立在无扩散系统(3)在平衡点E1处稳定的前提下,即(H1) 成立.

考虑仅含自扩散存在的情况,即D1=D2=0,DS0,DI0.系统(2)退化为

若平衡点E1处的微扰形式为eλtcos(kxx)×cos(kyy),则在平衡点E1处得到系统(5)的雅克比矩阵如下:

其中,k为波数,并且满足k2=+由雅克比矩阵J1得到系统(5)的特征方程如下:

定理2对系统(5)有如下结论:

1)若(H1),(H2) 均成立,则系统(5)在平衡点E1处局部渐进稳定;

2)若(H1),(H3) 均成立,则系统(5)在平衡点E1处发生Turing 不稳定.

证明

1)若(H1) 成立,则有m1＞0.当(H1),(H2)均成立时,对于所有的k2都有h1＞0 .此时方程(6)无正实部特征根,系统(5)在平衡点E1处局部渐进稳定.

2)若(H3) 成立,存在k2使得h1＜0 .此时方程(6)存在正实部特征根,则系统(5)在平衡点E1处发生Turing 不稳定.

从概念的描述看，第一学段的6个版本是结合图形给出概念名称；第二学段是结合图形给出概念的文字定义；第三学段不仅给出概念的文字定义，还结合图形说明概念的符号表示方法．

注2当系统处于局部渐进稳定时,易感人群(S)与感染人群(I)种群密度在空间的各个位置趋于平衡态.当系统处于Turing 不稳定时,易感人群(S)与感染人群(I)种群在空间中出现非均匀的密度分布.

考虑自扩散与交叉扩散均存在的情况,在平衡点E1处得到系统(2)的雅克比矩阵如下:

其中,k为波数,并且满足k2=+.由雅克比矩阵J2得到的色散关系如下:

做如下假设:

定理3若(H1),(H4) 均成立,则系统(2)在平衡点E1处发生Turing 不稳定.

证明若(H1),(H4) 均成立,总有合适的波数k2使得α ＜0,β ＜0,此时特征方程(7)存在正实部的特征根,即系统(2)在平衡点E1处发生Turing不稳定.

4 数值模拟

在二维空间中由反应扩散系统定义的连续问题在具有M=N=200 个网格位置的离散域中求解.Δh为网格点之间的间距,时间步长为Δt.在离散系统中,用有限差分计算描述扩散的拉普拉斯函数.本文设置Δh=1,Δt=0.01 .离散域初始条件选取为稳态周围0.01 量级的均匀分布随机扰动,其形式为

其中ξ1表示在200 × 200 的网格中添加0—1 之间的随机值.取参数A=1,µ=1.8,d=1,β=35时,(H1) 成立,由定理1 可知,系统(2)在无扩散时是局部渐进稳定的.分别取自扩散系数DS=6,DI=6(无自扩散驱动系统不稳定);DS=6,DI=1(自扩散驱动系统不稳定)两种情况进行对比分析.

注3本文主要参考了文献[23],故仿真中也基于该文献中的系统参数以及自扩散系数取值,并且所取值的系统参数及自扩散系数均满足上文的稳定性条件.在稳定性分析部分,我们发现改变交叉扩散系数会影响系统在平衡点处特征方程的特征值,即改变交叉扩散系数可能会对系统动力学产生一定影响,故在本文的仿真中选取不同的交叉扩散系数以此研究交叉扩散对系统动力学的影响.

注4当DS=6,DI=6时,(H1),(H2)均成立,由定理2 可知,系统(5)在平衡点E1处局部渐进稳定,即无自扩散驱动系统不稳定.当DS=6,DI=1时,(H1),(H3) 均成立,系统(5)在平衡点E1处发生Turing 不稳定,即自扩散驱动系统不稳定.

4.1 交叉扩散诱导Turing 斑图的形成

本节讨论当DS=6,DI=6 时,交叉扩散系数D1,D2对系统(2)空间格局的影响机理.

绘制系统(2)的色散关系曲线,如图1 所示.当D1=D2=0 时(蓝线)时,系统(2)退化为仅含自扩散的系统(5),观察到方程(7)所有特征根实部均为负数,此时系统(2)是局部渐近稳定的.当交叉扩散系数D2=6(黄线),D2=8 (紫线)时,观察到引入交叉扩散使得方程(7)在一定波数k2范围内存在正实部特征根,则系统(2)发生Turing不稳定.

图1 色散关系曲线Fig.1.Dispersion relationship curve.

验证色散关系所得出的结论,对系统(2)进行斑图仿真模拟.图2为DS与DI比值为1(无自扩散驱动系统不稳定),在未引入交叉扩散时不同迭代步数下的空间格局.在迭代200 步时,由于初值的随机扰动,此时系统(2)的空间格局处于一个不稳定的瞬态模式;当迭代到第10000 步,此时空间格局已是稳定的纯色均匀态.

图2 不同迭代步数下的空间格局Fig.2.Spatial pattern under different iteration steps.

由图1 色散关系曲线易知,引入交叉扩散能够改变系统(2)的局部稳定性,即分布均匀的纯色空间格局会发生改变.图3(a)—(d)为DS与DI比值为1 时,引入交叉扩散系统(2)出现的空间格局.图3(a)和图3(b)分别为交叉扩散系数D2=6,8时出现的孔洞-条纹斑图与条状斑图.图3(c)为交叉扩散系数D1=-6,D2=6 时出现的条状斑图.图3(d)为交叉扩散系数D1=-8,D2=8 时出现的点-条状斑图.

图3 不同交叉扩散取值的空间格局(a)D2=6 ;(b) D2=8;(c) D1=-6,D2=6;(d)D1=-8,D2=8Fig.3.The spatial pattern of of different cross diffusion values: (a)D2=6; (b) D2=8; (c) D1=-6,D2=6 ;(d) D1=-8,D2=8 .

注5上述模拟结果表明,在无自扩散驱动系统不稳定的情况下,引入交叉扩散能够使原本局部稳定的系统转变为不稳定的状态,即交叉扩散能够诱导Turing 斑图的形成.

4.2 交叉扩散改变系统的空间格局

首先讨论当DS=6,DI=1,交叉扩散系数D2=0,易感染者交叉扩散系数D1对系统(2)空间格局的影响机理.

图4(a)为不同交叉扩散系数D1下的色散关系曲线.当DS与DI的比值为6,且D1=0 (蓝线)时,方程(7)已经存在正实部特征根,系统(2)发生Turing 不稳定.交叉扩散系数D1为负值,即D1=-0.2,-1(色散关系图对应红、黄曲线)时,能够看出随着D1的减小,方程(7)特征根的正实部越来越大,此时系统(2)仍然会发生Turing 不稳定.当交叉扩散系数D1为正值,即D1=1,3.5 (色散关系图对应紫、绿曲线)时,能够看出随着D1的增大,方程(7)特征根的实部变小直至出现负值,即引入交叉扩散使原本在平衡点E1处不稳定的系统(2)转变为稳定的状态.

图4 不同取值的交叉扩散系数 D1 的色散关系曲线及空间格局(a)色散关系曲线;(b) D1 =0;(c) D1=-0.2;(d) D1=-1 ;(e) D1=1;(f)D1=3.5Fig.4.Dispersion relationship curves and spatial patterns of cross diffusion coefficients D1 with different values: (a)Dispersion relationship curve;(b) D1=0;(c) D1=-0.2;(d) D1=-1;(e) D1=1;(f) D1=3.5 .

图4(b)为当D1=0 时的空间格局为点-条状斑图;图4(c)和图4(d)分别为负交叉扩散系数D1=-0.2,-1的空间格局;当D1=-0.2,虽然此时的空间格局仍为点-条状斑图,但可以明显看出相较于D1=0 的空间格局点状的斑图数量更多,而当D1=-1 时,此时的空间格局已经全部变为点状斑图.图4(e)和图4(f)分别为正交叉扩散系数D1=1,3.5的空间格局.当D1=1,此时空间格局为迷宫状斑图;当D1=3.5,此时的空间格局变为了分布均匀的纯色图,即系统(2)转变为局部稳定状态.

注6图4 所示结果表明,在自扩散驱动系统不稳定的情况下引入易感染者交叉扩散系数D1改变了系统(2)的Turing 斑图结构,并且使自扩散驱动不稳定的系统转变为局部稳定状态.

下面考虑当交叉扩散系数D1=0 时,感染者交叉扩散系数D2对系统(2)的空间格局影响机理.同样取DS=6,DI=1 .

图5(a)为不同交叉扩散系数D2下的色散关系曲线.当正交叉扩散系数D2=0.5,2 时(色散关系图对应橙、黄曲线),能够看出随着D2的增大,方程(7)特征根的实部越来越大,此时系统(2)发生Turing 不稳定.当负交叉扩散系数D2=-1,-2时,方程(7)特征根的实部越来越小直至出现负值,即交叉扩散引入使得系统(2)在平衡点E1处转变为局部稳定的状态.

图5(b)为D2=0 时的空间格局为点-条状斑图.图5(c),(d)分别为正交交叉扩散系数D2=0.5,2 时出现的点-条状斑图,能够看出随着D2的增大点状斑图数量逐渐增加,且条状斑图的形状越来越不规则.图5(e)和图5(f)分别为负交叉扩散系数D2=-1,-2时的空间格局.当D2=-1,此时的空间格局为洞孔状斑图;当D2=-2 时,此时的空间格局变为了分布均匀的纯色图.

注7图5 所示结果表明,在自扩散驱动系统不稳定的情况下引入感染者交叉扩散系数D2改变了系统的Turing 斑图结构,并且使自扩散驱动不稳定的系统转变为局部稳定状态.

4.3 交叉扩散改变系统的稳定速度

本节研究在有无自扩散驱动系统不稳定的情况下,讨论交叉扩散对系统(2)到达稳定状态所需时间的影响机理.

首先考虑在DS=6,DI=6 的情况下,当交叉扩散系数D2=0时,易感者交叉扩散系数D1对系统(2)到达稳定状态所需时间的影响机理.

由图6 可知,虽然引入交叉扩散没有改变系统(2)局部稳定状态,但是能观察到交叉扩散对系统(2)到达稳定状态所需要的时间具有一定影响.在迭代200 步时,系统(2)的空间格局由于初值的随机扰动,当D1=-9,-9.8,-10.2 时都处于一个不规则的瞬态模式.当迭代到第450 步时,三组参数的斑点轮廓都逐渐减弱并且能够看出D1越小,斑点的轮廓越明显,即系统(2)到达稳定状态所需要的时间越多.当D1=-9,迭代到第800 步时,系统(2)的空间格局进入了均匀的纯色分布状态;当D1=-9.8,-10.2 时,分别到第1500,3000 步时系统(2)的空间格局才到达均匀的纯色状态.

图6 不同 交叉扩散系数对应的系统演化图(a)D1=-9;(b) D1=-9.8;(c)D1=-10.2Fig.6.System evolution diagrams corresponding to different cross diffusion coefficients: (a)D1=-9;(b) D1=-9.8;(c) D1=-10.2 .

考虑不同迭代步数时空间离散域中感染者I到达平衡点浓度的网格个数.当传播率β=35 时,此时平衡点处的浓度I1=0.2362.选择M=N=200,共有40000 个网格.图7 展示了在不同的交叉扩散系数下,不同迭代步数稳定到平衡点处浓度的网格数.在相同的迭代步数下D1的值越大,稳定到平衡点浓度的网格个数就越多,即在无自扩散驱动系统不稳定的情况下,负交叉扩散系数D1的值越大,稳定的速度越快,系统(2)到达稳定状态的时间越少.

下面考虑在DS=6,DI=1 的情况下,当交叉扩散系数D2=0时,易感染者交叉扩散系数D1对系统(2)到达稳定状态所需时间的影响机理.

由图8 可知,引入交叉扩散不仅改变系统(2)原本不稳定的状态,而且对系统(2)到达稳定状态所需要的时间也具有一定影响.当迭代到第700 步时,对于D1=4,斑点的轮廓已经较模糊,而当D1=3.5,3.2时,仍有较清晰的斑点轮廓存在,且当D1=3.2时斑点的轮廓更为清晰.当D1=4 迭代到1800 步时,系统(2)的空间格局进入了均匀的纯色分布状态,而当D1=3.5,3.2 时,分别到2600 步和5300 步才能够进入均匀的纯色状态.

图8 不同交叉扩散系数对应的系统演化图(a)D1=4;(b) D1=3.5;(c)D1=3.2Fig.8.System evolution diagrams corresponding to different cross diffusion coefficients: (a)D1=4;(b) D1=3.5;(c) D1=3.2 .

当传播率β=35时,平衡点处的浓度I1=0.2362,考虑空间离散域中感染者I到达平衡点浓度的网格个数.图9 展示了在不同交叉扩散系数下,不同迭代步数稳定到平衡点浓度的网格数.在相同迭代步数下D1的值越大,稳定到平衡点的网格个数就越多,即在自扩散驱动系统不稳定的情况下,正交叉扩散系数D1越大,稳定速度越快,系统(2)到达稳定状态的时间越少.

图9 当 D1=4,3.5,3.2 时,不同迭代步数稳定到平衡点的网格个数Fig.9.The number of grids that stabilize to the equilibrium point with different iteration steps when D1=4,3.5,3.2.

5 结 论

本文研究了交叉扩散对具有非线性发病率SI模型的空间格局影响机理等问题,如斑图的形成与斑图结构改变、系统稳定速度等.分析了无扩散系统平衡点的稳定性,并且阐明了存在扩散项时发生Turing 不稳定的条件.通过数值仿真得到了一些结论: 在无自扩散驱动系统不稳定的情况下,引入交叉扩散能够改变系统的局部稳定性,诱导斑图的形成.在自扩散驱动系统不稳定的情况下,引入交叉扩散能够改变斑图的结构.对于易感染者交叉扩散系数D1,当D1的取值为负时,斑图的结构从点-条状斑图转变为点状斑图,而当D1的取值为正时,斑图结构从点-条状斑图转变为迷宫状斑图,最后转变为分布均匀的纯色图.对于感染者交叉扩散系数D2,当D2取值为正时,斑图的转变与易感者交叉扩散系数D1取值为负时相同,逐渐转变为点状斑图,而当D2取值为负时,斑图的结构却出现孔洞状结构,最后转变为分布均匀的纯色图.对于系统稳定速度,在无自扩散驱动系统不稳定的情况下,交叉扩散会影响系统的稳定速度.易感者交叉扩散系数D1为负值且D1越大,系统稳定的速度越快.当存在自扩散驱动系统不稳定时,交叉扩散使系统从不稳定的状态转变为局部稳定状态且对其稳定速度也具有一定影响.易感者交叉扩散系数D1为正值且D1越小,系统(2)稳定的速度越慢.

基于本文的研究结果,未来将考虑将二维的空间格局研究推广至三维的空间格局,并且将整数阶模型推广至分数阶模型.扩散项在非线性动力学系统中有着重要的作用,未来也将考虑将扩散项转换为复杂网络,研究复杂网络所诱导的Turing 斑图.