融入迁移限制与双层集合种群网络的传染病模型

顾秋阳琚春华张双竹

(1.浙江工业大学管理学院,浙江杭州 310023; 2.浙江工业大学中国中小企业研究院,浙江杭州 310023;3.浙江工商大学管理工程与电子商务学院,浙江杭州 310018)

1 引言

传染病会在人群中迅速产生大范围传播,这是由于此类疾病容易通过复杂的人口网络进行传播.2020年初,新型冠状病毒肺炎(corona virus disease-19,COVID-19)疫情爆发,其病原体被国际病毒分类委员会命名为严重急性呼吸综合征冠状病2(severe acute respiratory syndrome coronavirus 2,SARS-CoV-2).新冠病毒以高传播效率、严重感染后果以及捉摸不定的流行时间对人类健康构成了持续的威胁[1].此次疫情发生以来,党中央高度重视,始终把人民群众生命安全和身体健康放在第一位,中央政治局常委会两次召开会议进行专题研究.由于全球疫情的全面爆发,如何制定流行病动态控制策略和模型备受学界关注.到目前为止,已有学者提出了各种用于分析流行病传播模型[2-4],并对人群中的疫苗接种和流行病传播状况进行了研究.学界通常将传染病传播表示为基于节点(即agent)的模型,链接(路径或边)表示个体间的交互作用;另一种方法是使用反应扩散等微分方程描述种群模型.已有学者在生态建模中建立了集合种群模型[5-6],集合种群是由某物种在空间上分离的子种群(斑块)组成的.个体可以在斑块间进行迁移[7-8].随机游走或迁移是与许多传播过程相关的基本机制,故研究迁移对流行病传播的影响具有重要意义[9-12].而异质连接和移动模式对传染病传播有重要影响.Colizza等[13]通过复杂网络中的数值模拟,探讨了传染病在集合种群模型中的传播行为.结论表明,由于网络结构的复杂性,在多种群动态模型中求解微分方程并不容易.Nagatani等[8]研究了完全图、环形图和星图上的集合种群模型中随机游走节点的传染病传播模型.本文使用不同的解析方法推导出了上述种类图的传染病阈值.且由于网络异质性对传染病传播与阈值产生了重要影响,传染病动力学与网络拓扑的耦合问题已引起了学界的广泛关注[14-17].Onaga等[18]认为可通过控制被感染者的活动范围,以减缓传染病的蔓延.且其认为被感染者的活动范围(即迁移范围)与易感染者的活动范围(即迁徙范围)有很大不同,被感染会大大降低人口的流动性.

已有学者基于经典易感-感染-易感(susceptible infected susceptible model,SIS)模型,在静态和时序图中构建了集合种群动态模型进行研究[19],在集合种群模型中,所有种群(或个体)都沿同一迁移路径移动,此路径由单一网络图表示.然而,当某种群的作用范围明显不同于与其他种群时,各种群的迁移路径都不同.从一个种群到另一种群,用不同的图描述种群间的迁移路径是必要的.由于受感染个体的迁移路径与易感染个体不同,故本文使用双层网络方法对受感染个体的限制迁移对传染病传播的影响进行研究.

本文在双层网络上建立了传染病传播的集合种群模型.其中,当个体受感染时,与易感染个体相比其迁移能力会受到限制.本研究用双层网络代替传统的单层网络来表示受感染个体的限制性迁移.在双层网络上,本文研究了传染病传播是如何随着迁移限制进行变化的,并在双层网络的集合种群模型中,建立了易感-感染-易感(SIS)模型.本研究将双层网络上的SIS模型与随机游走模型相结合,并在双层网络上推导出了集合种群动态模型中的反应扩散方程,并对其进行数值求解.实验结果表明,在集合种群动态模型中,传染病的传播很大程度上取决于受感染个体的迁移限制.

2 预备知识

在SIS传染病模型中,可将流动个体分为易感染节点(susceptible,记为S)或感染节点(infected,记为I).并将节点的总密度设置为ρ0,将子种群的密度设置为ρI,故可得如式(1)所示:

故可将易感染节点(S)和感染节点(I)间的交互过程表示如式(2)-(3)所示:

并将其转化过程如式(2)所示的传染率设置为β,而如式(3)所示的转化过程的恢复率可设置为γ.如系统中仅包含单个斑块,则均值场理论(mean field theory,MFT)成立,具体如式(4)所示:

其中ρI(t)表示t时刻感染个体的密度.达到稳态时的解可表示为如式(5)所示:

本研究还构建了一个由N个斑块(节点)(N≥2)组成的网络,其中每个斑块的平均场动态保持不变,而个体的动态由反应扩散方程表示.集合种群的分布结构由网络图进行表示,子种群i在图中表示为节点i,且其中许多个体通过链接在节点间随机游走,在子种群i中的随机游走者总密度可由如式(7)所示计算:

其中:ki为节点i的度,Ni表示节点i的最近邻.而式(7)右侧从左至右第1项表示最近邻的进入个体数量,第2项表示子种群中的外流个体的数量.参考文献[20]可知,即使系数取小于1的值,相变点也不会改变(由于要设置动态的外流参数较难,故在此只考虑固定系数),故在此假设系数为1.故在集合种群模型中,所有子种群的感染率和恢复率相同.且在集合种群模型中,子种群i中受感染节点密度可表示为如式(8)所示:

其中:式(8)右侧从左至右的第1项和第2项分别表示扩散引起感染密度的增加和减少,第3项表示如式(2)所示的感染密度(I)增加,而第4项表示如式(3)所示的感染密度(I)减少.在集合种群模型中,感染率和恢复率不依赖于子种群,但其数值相同.

当个体数量趋于无穷时,集合种群模型才为有效;如感染个体的密度是负数,则将密度设为零.具有3个节点的集合种群模型是同时具有同质图和异质图的最简单的模型,且在同质和异质图中传染病动态不一致.如图1所示,图1(a)报告了完全图/环形图,节点1-3的度都为2;图1(b)报告了星图,其中节点1的度为2,节点2和节点3的度为1.其中,完全图/环形图为同质图,星图为异质图.

图1 不同3节点网络结构示意图Fig.1 Diagram of different three-node network structures

图2为不同三节点网络结构下的平衡密度与恢复率关系图.

图2 不同三节点网络结构下的平衡密度与恢复率关系图Fig.2 Relations diagram between equilibrium density and recovery rate under different three-node network structures

图2(a)报告了在恢复率γ下的易感染节点与感染节点的平衡密度的分布图,其对应图1(a)中的环形图;且参考文献[21]的做法,将初始条件设置为ρS,1(0)=0.5,ρS,2(0)=0.2,ρS,3(0)=0.2,ρI,1(0)=0.05,ρI,2(0)=0.01,ρI,3(0)=0.04,β=1.初始状态下,节点1-3处的感染密度和易感染密度相同.随着恢复率γ从0增加到时,被感染密度呈线性下降.如式(5)所示平均场由在图中由MET表示.随着恢复率γ从0增长到,易感密度呈线性上升.当γ=时,易感密度饱和;当γ >时,易感密度恒定为.图2(b)报告了在恢复率γ下的易感染节点与感染节点的平衡密度的分布图,其对应图1(b)中的星图;且与图2(a)具有一致的初始值.节点2和节点3处的感染密度一致.而中心1的感染密度为节点2和节点3 的两倍.节点1-3的感染密度与如式(5)所示平均场理论存在较大偏差,阈值γC=0.383.当恢复率γ从0增长到0.4时,虽然易感染节点与感染节点的变动趋势与图2(a)中一致,但其变化速度存在区别.

3 融入集合种群的传染病模型

3.1 双层网络集合种群传染病模型

本文提出了一个上融入双层集合种群网络和迁移限制的传染病模型.由于感染者的活动范围受到限制,感染者和易感染者的迁移路径明显存在不同.在此情况下,传染病模型应从单层图上如式(7)-(8)所示转化为在双层图上描述其动态.子种群(即斑块)由双层网络上的节点表示,第1张图上的链接表示易感个体斑块间的迁移路径,第2张图上的链接表示受感染个体斑块间的迁移路径.由于受感染者的迁移受到限制,第2张图与第1张图有所不同.易感个体以随机游走的方式移经过第1张图上的链接,而受感染个体以随机游走的方式移动通过第2张图上的链接.在具有N个节点的双层集合种群传染病模型中,节点i处的易感染个体和受感染个体密度可由式(9)-(10)推导得到

3.2 三斑块集合种群传染病模型

本文首先考虑具有3个斑块的集合种群传染病模型.受感染者与易感染者的迁移途径不同,故在双层图中推导SIS传染病模型动态,并展示了传染病传播是如何受感染个体迁移限制而发生变化的.图3中的第1张图的易感染个体迁移路径用黑线表示,第2张图的受感染个体迁移路径用灰线表示.图3(a)中的易感染个体迁移路径用循环图表示,受感染个体的迁移路径用星型图表示.由于感染个体迁移受到限制,感染个体从不直接在子种群2和3间进行迁移.易感染个体在亚群1,2和3之间随机游走.

图3 当N=3时的双层网络示意图Fig.3 Diagram of the two-level network when N=3

本研究将如式(9)-(10)应用到图3(a)的双重图中.对节点1-3子种群中易感染和受感染个体的密度可由如式(11)-(16)反应扩散方程推导得到

在图3(b)中,第1张图上的易感染个体在循环图上迁移,与图3(a)一致,第2张图上的子种群2和子种群3间只有一个链接(即子种群2和子种群3中的受感染个体往返于子种群3和子种群2,子种群1上的受感染个体不能迁移到其他子种群中);而易感染个体可以在任意节点中随机游走.本研究将式(9)-(10)运用到图3(b)的双层网络中,子种群1-3中易感染和感染个体的密度可由如式(17)-(22)反应扩散方程推导得到

在图3(c)中,易感染个体的迁移路径由星图表示.受感染个体的迁移路径由节点1和节点2间的链接表示;由于迁移限制,受感染个体仅在子种群1和节点2间进行迁移,个体随机游走在双层网络之上.本研究将式(9)-(10)运用于图3(c)的双层网络中.子种群1-3中易感染和被感染个体密度由如式(23)-(28)所示反应扩散方程推导得到

在图3(d)中,第1张图上的易感染个体在星图上迁移,第2张图上的子种群2和3间只有一个链接,子种群2和子种群3中的受感染个体往返于子种群3和子种群2,子种群1中的受感染个体不会迁移到其他子种群种.本研究将式(9)-(10)运用到图3(d)的双层网络中.子种群1-3中易感染和感染个体的密度可由如式(29)-(34)所示反应扩散方程推导得到

图4 不同双层网络三节点拓扑结构下平衡密度对恢复率的影响Fig.4 Effect of equilibrium density on recovery rate under three-node topology of different two-level networks

3.3 四斑块集合种群传染病模型

对于存在4个节点的图,共有6种不同配置模式.本研究在6种不同拓扑结构的单层图上对传染病传播动态进行了研究,并在双层网络上研究了受感染个体迁移受限的传染病传播动态[8].本研究在图5设置了两种具有4节点的双层网络传染病动态如图5(a)-(b)所示.第1张图上易感染个体的迁移路径用黑线表示,第2张图上的受感染个体的迁移路径用灰线表示.在图5(a)中,易感染个体的迁移路径由星图表示,受感染个体的迁移路径也由同一星图表示;在图5(b)中,易感染个体的迁移路径由星图表示,受感染者的迁移路径由路径图表示.

图5 当N=4时的双层网络示意图Fig.5 Diagram of the two-level network when N=4

本研究将等式(9)-(10)应用于图5(a)中的星图.种群1-4中易感染和感染个体的密度由式(35)-(42)所示反应扩散方程推导得到

易感染个体在式(35)-(38)中的扩散项与式(39)-(42)中被感染个体相符.这是由于被感染个体与易感染个体的移动路径相同.本研究将式(9)-(10)运用到如图5(b)所示的双层网络中.子种群1-4中易感染和受感染个体的密度可由如式(43)-(50)所示反应扩散方程推导得到

其中,式(47)中没有扩散项,这是由于节点1处的受感染个体不会迁移到其他节点.式(47)-(50)中受感染个体的扩散项明显不同于式(43)-(46)中的易感个体的扩散项.

图6(a)报告了在初始值为

且β=1时,图5(a)上双层网络中恢复率与平衡稳态易感和受感染密度的关系.由于节点2-4的拓扑等价,节点3和节点4处的受感染密度与节点2处的保持一致.受感染密度随恢复率的上升呈线性下降.节点2-4的易感密度相同.随着恢复率从0(γ=0)增加到左右时(γ=0.34),节点1处的易感密度呈线性上升;当γ=0.34时呈饱和状态,当γ>0.34时,节点1处的易感密度为0.5.随着恢复率从0(γ=0)增加到时(γ=0.34),节点2-4处的平衡稳态易感密度呈线性上升,当γ=0.34时呈饱和状态,当γ >0.34时,节点2-4处的易感染密度保持恒定值

图6(b)报告了与图6(a)具有相同初始值,且β=1时,图5(b)上双层网络中恢复率γ与平衡稳态易感和受感染密度的关系.由于节点2和节点4的拓扑等价,故节点4处的受感染密度与节点2处相同.随着恢复率从1(γ=0)转变为0(γ=0.5),节点1(中心)处的受感染密度呈线性下降趋势.节点2-4处受感染密度为0,受感染个体都集中在节点1(中心).节点2-4处易感密度相同,随着恢复率从0(γ=0)上升到(γ=0.5),节点1的易感染密度呈线性上升,当γ=0.5时呈饱和状态,当γ >0.5时节点1处易感密度保持恒定值0.5.随着恢复率从0(γ=0)上升到(γ=0.5),节点2-4处平衡稳态易感密度呈线性上升,当γ=0.5时呈饱和状态,当γ >0.5时,节点2-4处易感染密度保持恒定值.故双层网络的拓扑结构不同,节点1-4处的易感个体和受感染个体的平衡密度有很大差异.

图6 不同双层网络四节点拓扑结构下平衡密度对恢复率的影响Fig.6 Effect of equilibrium density on recovery rate under four-node topology of different two-level networks

3.4 环形星图中的传染病模型

最后,本文对传染病在环形星图上的传播进行研究.易感染个体在第1张星图上移动,受感染个体在第2张环形图上移动,且第1张图与第2张图的链接不重合.图7(a)报告了当N=6时星图的示意图,节点编号为1-6.第1张星图用黑线表示,第2张环形图用灰线表示.且节点1为第1张图的中心,易感染个体通过中心节点1在星图上游走,中心节点1的受感染个体不会迁移到节点2-6,但节点2-6处的感染者可以在环形图中的节点2-6间迁移.受感染个体的迁移路径不会与易感染个体重叠.图7(b)报告了N=6的星图,其中易感染个体在第1张图上移动,感染个体只能在灰色路径上移动.

图7 当N=6时的双层环状星型网络示意图Fig.7 Diagram of the two-level ring star network when N=6

在图7(b)报告了的双层网络上的数值计算出的平衡密度结果与图8(a)一致,在此不再赘述.且本研究发现切断第2张循环图中的任何链接,平衡密度都不会发生改变.本研究通过计算分析得出了环形星图上的平衡密度.其在平衡状态下满足如式(57)-(58)所示:

图8 单层/双层星图拓扑结构下平衡密度对恢复率的影响Fig.8 Effect of equilibrium density on recovery rate under two-level and one-level star-graph topology

在拓扑等效节点上,处于平衡状态的种群数量相同,具体可表示为如式(59)所示:

除了中心节点1 外,节点1上没有受感染的个体,受感染的个体只能在中心节点1上存活.平衡密度的解析式(64)-(67)与图3(d)和图5(b)中具有3个和4个节点的环形星图的数值解一致.且对于具环形星图而言,其解与图8(a)中的数值解一致.受感染个体的平衡密度与节点数N无关.传染病阈值科表示为γc=通过分析可知阈值与环形星图的节点数N无关,故所得结论具有普适性.且为比较在单层网络传染病模型与本文所提双层网络模型中平衡密度对恢复率的影响的区别,参考文献[23]等的做法在星图中设置单层网络模型参数,并得到如图8(b)所示在单层星图拓扑结构下的仿真结果.比较可知,本文所提双层网络结构可有效降低感染节点密度,将感染节点限制在中心节点(度值最高的子种群)中,较单层网络结构更优.

为进一步体现本文所提模型在真实传染病环境中的可用性和科学性,本研究参考文献[24]计算得到的新冠疫情中的真实动态参数值(传染率、恢复率和基本再生数等),将本文所提模型在星图中进行了进一步仿真实验,得到如图9所示实验结果.比较图8-9结果可知,本文所提融入迁移限制与双层集合种群网络的传染病模型能够较好的模拟真实情况中的传染病传播,并能对传染病传播达到一定的控制效果.

图9 真实参数环境下双层星图拓扑结构中平衡密度对恢复率的影响Fig.9 Effect of equilibrium density on recovery rate in the topology of bilayer star map under real parameter environment

4 结论

文献[25-26]等在社会和空间层网络上提出了具有异质混合和传播行为的种群流行模型,研究了异质混合和传播行为对传染病传播的影响.而其所使用的双层网络中的集合种群模型类似于本文使用的双层网络模型.然而,本文所提基于经典SIS模型提出的融入迁移限制和集合种群概念的传染病模型,用于表示易感染个体和受感染个体的两次迁移.由于受感染个体的迁移路径与易感染个体的迁移路径有明显区别,故本文使用双层网络以描述易感染和受感染个体的迁移路径.在双层网络的种群模型中研究了受到迁移限制的个体对传染病传播情况的影响.本研究得出了具有3,4和6个节点的双层集合种群网络传染病模型的数值解.实验结果表明,受感染个体的迁移限制对模型的平衡密度具有重要影响,特别是在星图(或环形星图)上,受感染个体仅在中心节点生存,而不会在其他节点出现.故可知集合种群图(网络)结构在控制传染病的传播中起着重要的作用.且如果受感染个体的迁移仅限于低度节点间的链接,受感染个体仅在中心节点生存,在其他节点消失;如感染个体在低度节点的迁移受到限制,则感染个体会聚集在高度节点上.星图是典型的异质图,在现实生活中具有异质性的斑块环境类似于星图拓扑,故其能够较好的模拟现实生活中的传染病传播环境.如受感染个体的迁移被限制在较低程度的区域内,那么在环形星图上推导出的结论则可以应用于现实网络中的传染病传播研究.

尽管本研究已提出了上述具有重要意义的发现,但还是具有一些局限性,其中一些可能会为未来的进一步研究指明方向:首先,本文所提模型是在节点充分混合的假设下构建了融入集合种群传染病模型.但当斑块内部存在强烈的异质性时假设不成立,而对于具有强异质性的系统可尝试使用Agent-Based模型进行建模.其次,本文基于SIS模型进行建模,没有考虑到免疫者的要素,故本文认为将来可以尝试使用SIR模型、SEIR模型等更复杂的模型进行建模.并可参考文献[27]中所提多层集合种群网络结构与个体周期性迁移,对传染病模型进行改进,进一步提升模型的科学性.再次,可以先考虑外流参数动态化的情景下的双层集合种群网络的传染病模型.最后,可以基于机器学习、深度学习等算法进行模型改进和参数的动态调整,进一步提升模型的科学性.