复杂网络上流行病传播动力学的爆发阈值解析综述

李睿琪, 王　伟,舒盼盼, 杨　慧,潘黎明,崔爱香,唐　明

(1.电子科技大学互联网科学中心,成都611731; 2.北京师范大学系统科学学院,北京 100875)



复杂网络上流行病传播动力学的爆发阈值解析综述

李睿琪1, 2, 王伟1,舒盼盼1, 杨慧1,潘黎明1,崔爱香1,唐明1

(1.电子科技大学互联网科学中心,成都611731; 2.北京师范大学系统科学学院,北京 100875)

摘要：对流行病传播爆发阈值的理论解析方法进行总结，主要介绍平均场、点对近似、主方程、边渗流、空穴理论、边划分以及谱分析这7种常用的动力学解析方法的前提假设、具体思路、步骤及其应用局限，并且梳理总结了SIS与SIR模型爆发阈值的异同。

关键词：复杂网络；流行病传播；爆发阈值；理论解析

0引言

流行病的预警与干预作为公共安全的重中之重，已受到社会各界的广泛关注[1-2]。然而想要进行流行病传播规律的实体研究是不可接受、不可行的，数理模型自然成为流行病传播研究方向的希望与重点，通过对数理模型的分析及与实证的对比，让我们对于流行病传播的机制与防控的效率能够有更为深入而客观的认识。早在1760年，Daniel Bernoulli就提出了第一个流行病传播模型[3]；1927年，Kermack和McKendrick确立了流行病传播的现代数理模型框架[4]，激励了无数科学家的后续工作；1992年，Anderson和May的重要著作[1]则是对过去流行病传播研究工作的精到总结。自从1998年复杂网络理论诞生以来[5-7]，网络思维逐步深入人心。网络化的思考方式，可以将流行病传播过程中的重要因素(例如个体交互的结构、移动性、交互模式)进行很好的建模[7-11]，复杂网络上的传播动力学研究成为了流行病传播研究新的方向和希望。十几年间，复杂网络上的流行病传播研究受到了极大关注，取得了飞跃性的发展[12-13]。早在年，Pastor-Satorras等率先研究了因特网上的病毒传播问题[14-15]，通过理论解析发现在热力学极限情况下病毒爆发的阈值接近于零，这暗示无论多小的传染概率都将导致病毒长期存在于网络中。2004年德国马普所的Hufnagel等在全球性流行病传播方面做出了开创性工作[16]，他们利用全球航空网络及相关航班数据重现了2003年SARS的时空演化斑图。美国印第安纳大学的Vespignani等主持了一项名为全球性流行病传播建模(Global Epidemic and Mobility Modeler)的项目。这个项目基于复杂网络理论，根据病毒传播的特点和世界交通的数据，计算病毒可能的传播情况，从而使我们能够对将来可能发生的疫情进行预警[17]。令人振奋的是，他们较为准确地预测出了2009年甲型H1N1流感的时空演化斑图[18]。2013年，Brockman等利用航空网络的流量数据，定义了地点之间的有效距离，将原先观察到的较为复杂的流行病传播演化过程重新映射回反应扩散过程，并且可以通过这一方法有效地反推出流行病的源头[19]。

在流行病研究中，上述这类模型通常被称作复合种群模型，它们的基本假设是将地点看作节点，将个体看作节点中的粒子。目前很多研究已经足够精细、足够定量化，已然能够为政府决策提供一定建议与支持。然而更为基础的流行病研究则是将个体看作网络中的节点，分析其上传播动力学特征。真实的网络其结构通常具有时效性、动态演化性，然而静态网络是这些更为复杂情况的基础，所以如何能够将静态网络传播动力学有效地分析出来，是理论研究的第一步，也是本文关注的重点。

为了能够解释并预测流行病传播中的一些重要问题，许多研究者利用计算机模拟进行了许多相应的尝试；但科学家们对这样的结果往往并不满意，想要找寻出其内在的规律，用更为准确的方式来更加精确地进行预测。另一方面，即便对于静态网络上传播动力学的解析，想要得出一般的解析解也是极其困难甚至不可能的，网络层级性、聚集性等复杂性因素也迫使研究者们发展相适应的数理解析方法、计算模型。而且值得一提的是，虽然不同的动力学具体机制不尽相同，流行病传播动力学模型及其解析方法还被应用、借鉴到了更为广泛、更为一般的动力学行为研究上，例如信息传播、观点传播、文化规范和社会行为、电网级联失效[11]等等。

在流行病传播动力学的研究中，流行病的爆发阈值无疑是非常重要的，它对于流行病的评估、预警、干预策略的选择都有极为重大的影响。研究者发现网络上的流行病爆发这类宏观涌现行为的研究方法路线与统计物理中的非平衡相变非常接近[20]。受此启发，在过去十年，诞生了众多流行病传播解析方法：从经典的平均场方法到更为严格的定量化数值解析方法[21-22]。然而目前复杂网络研究中，阈值解析的方法繁多、鱼龙混杂，部分文章往往只是对某些细节进行了考虑、可以适用于某些略为特殊的情况，但却未必会有较好的拓展性与普适性。我们进行这样的解析方法综述，一个重要的目的就是对目前各类繁多的理论方法进行一个梳理，探寻不同方法的解析思路、适用条件、准确性等重要问题，让我们对于这些相应的方法能有一个更为直观而全面的认识和理解。

本文将以流行病传播的解析方法为主要研究对象，力求对传播动力学爆发阈值的解析方法进行较为全面的分析。按照各类方法使用的数学工具，将其分为3大类：微分方程法，生成函数法以及谱分析法。具体介绍平均场(也会在这部分简介SIS与SIR模型以及流行病传播对应的更为广义的相变)、点对近似、主方程、边渗流、空穴理论、边划分以及谱分析这7种常用的动力学解析方法，力求较为全面地覆盖传播动力学研究中重要的解析方法，并探寻各类解析方法内在的机制原理、前提假设、解析思路、物理含义及其局限性,以期为初学者提供一个较为清晰的学习脉络，为广大复杂网络研究人员提供一定的参考。

1平均场方法

1.1平均场理论概述

在复杂网络研究中，进行解析时最为简洁的方法当数平均场近似。平均场方法计算过程简单且原理易于理解，基本上只要熟悉常微分方程的人都具备运用平均场方法的能力，这也大大降低了理论解析的门槛。平均场方法被广泛地应用于流行病传播动力学解析及预测[15,23-29]、投票模型[30]以及同步现象[31]的解析当中，其解析结果可较好地反应实际情况。概言之，平均场是把环境对对象的作用进行集体处理，以平均作用效果替代单个作用效果的加和的方法。虽然其精确适用的限制条件非常严格，但是由于其解析思想简洁、过程简明，平均场近似被广泛地应用。然而不能完全满足平均场应用的假设条件时，其理论解析结果的准确性和如何进行相应的改进等问题就受到了众人的关注，也得出了许多重要的结论[24-29]。本节兼顾最新结论对于原始平均场方法的指导与更正，但并不主要讨论这些最新的发现，只在后面部分大略涉及一些关于平均场最新进展的讨论。我们力争对平均场方法进行一个完整而清晰的描述并给出具体的解析案例，使初学者对于平均场有一个整体的把握。我们将主要关注平均场方法在全接触SIR模型阈值解析中的应用，并以小世界网络与无标度网络为例给出具体的解析过程及相应分析结论。

1.2平均场方法提出的背景

平均场理论最初是由朗道为解释二级相变而提出[32-35]。对于所有的连续相变，都可以定义一个序参量，使它在临界点的一边为零，表示相对的无序状态，而在临界点的另一边不为零，表示相对的有序状态[32-35]。与一级相变(即不连续相变，如物质气、液、固三态间的转变)不同，二级相变(连续相变，如铁磁相变、超导相变、超流相变、以及合金的有序无序相变等)不具有两相共存的临界变化阶段和清晰的界面(例如在冰融化这一一级相变过程中，存在冰与水两种相位共存的情形)，相变在物质各处同时地突然发生。在二级相变过程中，各处同时出现新相与旧相互相嵌套(而非清晰的相位变化界面，此时很难判定它究竟处于哪一个相位)，形成所谓“花斑”[32]。这说明二级相变中，物质分子与其近邻的其它分子如何互相作用、结合成什么相，主要取决于跨越一切尺度的所有分子之间的相互作用的总体效果，亦即“在连续相变的临界点，关联长度趋于无穷”；相比之下，对于一级相变至关重要的分子的周围环境的“局部信息”(如晶格常数、晶格结构等)对二级相变的影响并不重要。因此，在二级相变中，可以把跨越一切尺度的所有分子之间的相互作用的总体效果等价于一个“平均场”，不去计算局部的、处处不同的相互作用情况——而这正是平均场方法的核心思想。例如对于二维Ising模型的相变求解，在低温时，自旋关联很大，热强度涨落较大(亦即偏离平均值较多)，平均场方法很难正确描述相应现象[36]。

虽然连续相变是一个理想的平衡过程(即不存在能量与物质的流动，或者说流动极其缓慢、在局部可认为无变化的物理过程)，然而，以广义的序参量随驱动量变化的幂律为标志的临界现象涉及更为广泛的领域。对于流行病的传播过程，可以将其等价为一种广义的相变。在流行病停止传播的状态，节点的感染密度ρ为零时，系统相对无序；而在流行病传播的状态，ρ大于零，系统相对有序。因此可以把ρ定义为广义的序参量。对于流行病的传播，大量的实验和理论结果都证明对于有效感染率λ存在一个阈值λc，只有在λ-λc>0时，流行病才能全局传播，因此可以把λ-λc定义为热力学驱动量。本节将介绍如何运用平均场近似方法的思想，写出流行病传播的动力学方程，并在一定限制条件下进行求解分析。

1.3SIR模型的假设及其微分方程表达形式

经典的SIR(Susceptible-Infected-Recovered/Removed)与SIS(Susceptible-Infected-Susceptible)模型是Reed和Frost在1920年的一篇未发表论文中首先提出的，这也是首次开创性地利用微分方程来进行动力学的描述[1,37]。他为个体定义了3种状态：易感态(Susceptible)、感染态(Infected)、恢复态(Recovered/Removed)，任意个体必处于这3种状态中的某一种。SIR模型假设流行病传播的周期远远小于个体的寿命，所以不考虑节点的死亡，整个流行病传播的过程中，网络节点总数恒定，亦即网络为静态网络。一个易感态节点(S)以某一概率λ被任一个感染态节点(I)感染, 一个感染态节点(I)以概率μ康复且变为终生免疫的恢复态(R)。SIR模型中，节点状态演化过程为S→I→R，或一直保持S状态；当传播过程结束后，只会有S与R两种状态存在。

对于随机规则网络(节点度相同)和均匀网络(度分布服从正态或者泊松分布)，设网络中节点的平均度为〈k〉(它实则是对任意节点在某一时刻接触能力的描述，对于单点接触(Contact Process, 又译作接触过程, 但笔者认为单点接触这一说法更加形象直观)或者有限接触能力的传播情形，只须将此处修改为相应数值即可，稍后在1.6节中还会详细介绍)，流行病的感染率为λ，节点的恢复概率为μ，用S(t)表示t时刻处于S态节点的密度；ρ(t)、R(t)依次对应I态与R态的节点密度。则SIR模型的微分方程组如式(1)～(4)[15]。

(1)

(2)

(3)

1=S(t)+ρ(t)+R(t)

(4)

第一个方程其含义是t时刻S态节点被感染为I态的概率，亦即S态节点密度减少的变化率。为了得到这样一个变化率，就必须知道一个S态节点会连接到多少个处于感染态的节点，以确定其被感染的概率。平均地说来，一个节点会有〈k〉个邻居，其邻居在t时刻处于感染态的概率为ρ(t)，进而可知其处于感染态的邻居平均说来有〈k〉ρ(t)个，所以第一个方程的精确的形式应当是dS(t)/dt=-S(t)(1-(1-λ)〈k〉ρ(t))；但是这样一个方程较难处理，将S(t)(1-(1-λ)〈k〉ρ(t))使用泰勒展开后，舍掉λ的二阶及以上的高次项而得到的线性近似即可得到S(t)λ〈k〉ρ(t)这一线性表达，而这样的线性近似就自然而然地会引入误差，只有在其近似适用范围内(λ〈k〉ρ(t)≪1)，才能认为较为精确。

第二个方程是表征t时刻处于I态节点密度变化的概率，这当中有原先处于S态被感染后新增的个体，也有由I态以概率μ恢复为R态减少的个体，故由两部分组成。而R态t时刻的变化率则只可能是由I态节点转化而来，故第三个方程仅有μρ(t)一项。

分析SIR模型对应的微分方程，可以看出它隐含了许多近似假设以及平均场的思想。〈k〉的引入是平均场思想的重要体现，也是平均场方法最大的假设所在：亦即认为波动可以忽略，每一个节点的度都可以用平均度〈k〉来代替。上面提到的线性近似在处理阈值时适用的原因，是由于在低于传播阈值λc的情况下，当网络中节点感染密度远远小于1时，对S(t)(1-(1-λ)〈k〉ρ(t))进行泰勒展开，得到S(t)λ〈k〉ρ(t)是较为精确的；然而当感染密度较大(大略为ρ(t)>1/〈k〉)时，上述的近似结果从理论上讲将不再十分精确(偏差将增大,而且通常解析值会高于实际值)，但已有的工作大多都表明它仍能较好地描述其趋势、定性解释许多现象，甚至定量地来看也可以得到不错的结果[12,13,28,38-41]。

1.4平均场理论应用的约束条件及准确性分析

平均场方法从平均的、概率的角度去考量问题，免去了对于细节的计算，故而易于使用；但原则上来讲，其适用条件是非常严格的，历史上著名的科学案例(如二维Ising模型的相变问题)明确警示了滥用平均场方法的不良后果。数学上，平均场方法的适用范围只能是完全图，或者说系统结构是充分混合的：在这种情况下，系统中的任何一个个体能够以等概率接触其他个体。只有满足这一条件，平均场方法得到的解析结果才是完全精确的。

在复杂网络研究中，与其说平均场理论是一种技术手段，把它看作一种思想方法或许更为恰当，即从平均的、概率的角度去考量问题，进行宏观的把握。其技术手段无外乎常微分方程、泰勒级数等等。从平均场方法自身发展的历史来看，它也是不断地缩小近似的范围、进行更为合理的近似，著名的异质平均场方法(Heterogeneous Mean-Field Method)[15,23]就是将度相同的节点进行平均近似，而不再把所有节点看作相似的。平均场的这种思想方法也渗透到了点对近似(Pairwise Appropriation)、主方程(Master equation)[28]等更为精确的方法中，这些方法近似范围更小；当然这些方法超越平均场近似之处在于，考虑了点对情况以捕捉动力学相关性[28]、关注度的演化[42]等等，使得相应微分方程更为精确。

在复杂网络研究中，标准的平均场方法一般隐含了3个前提假设[28](我们会在1.5中结合具体的平均场方程来依次具体阐述这些假设是如何被体现的)：

1)无局部聚类：对于任意度为k的节点A，其邻居的状态互相独立：亦即A某一邻居B1的状态只受A的影响，而不受A其它邻居(B2-Bk)的影响。故而当邻居之间存在连边时，这一假设并不完全满足；然而应当注意的是，实证发现在感染率小于或等于阈值时，流行病的传播过程基本上是树形分支过程，产生的回路极少[27]，而这也能部分地说明为何平均场方法在预测阈值上较为准确。

2)无模体：所有度相同的节点，其动力学行为可用其平均情况代替。网络中如果存在大量模体，即便两个节点度相同，但由于可能处于不同模体的不同位置中，它们对于传播所产生的影响将非常不同。例如三元模体中，度相同的节点在不同的模体中可能处于不同的位置、发挥不同的作用：有些可能是接收者，而有些则是作用者[43]；有些是纯然的传播者，有些则是完全的被感染者。这类区分在有向网络中更为明显，例如图1a中被双环选中的度为2的两个节点，在此种情形下，虽然二者度相同但其动力学行为必将大不相同。姑且不论有向网络，即便是在无向网络中，图1b中4个被选中的节点虽然度均为3，然后它们在传播过程中，第二个时间步所导致的影响将极不相同。在许多文献中通常会看到平均场方法的一个基本假设是波动可忽略[23,41,44]，它实则是指所有的情形与平均情况相近，其波动不大：具体说来即网络中所有节点度及动力学行为相近；或者起码在度相同的情形下，其动力学行为相近。

3)无动力学相关性：更新某一节点A的状态时，A的状态与其邻居(B1-Bk)的状态相互独立，即A状态的更新变化不会作用于Bi的状态。这一假设实则是忽略了邻居节点间的动力学相关性，有理论研究发现，假设 3)正是平均场理论误差产生的主要来源[28](点对近似恰恰是捕捉到了这种动力学相关性；对于这一假设的包含与否，恰是平均场方法与点对近似最大的区别之所在)。

但是当聚类及模体作用较强时，点对近似同样无法精确描述相应动力学行为[28]，退火网络邻接矩阵中存储的为其连接概率)由于每个时间步都会以相应概率进行边的重连，故而它取消了动力学相关性[12,45](但是应当注意的是，退火网络仍然有拓扑相关性，对于任意节点其连接情况还是明确的，只是从0-1二元值变为了概率值，所以同配网络和异配网络的退火邻接矩阵的整体数值分布情况还是很不同的)；另外，对于动态网络，当节点移动的速度较大时，也能够较好地满足平均场的假设[46]。故而一般认为网络中节点平均度较大时解析结果相对准确：因为当平均度较大时，网络更靠近于全连通图，所以这种情况下动力学相关性反而较弱，不再起到决定性作用；当平均度较小时，聚类、模体、动力学相关性的作用往往更为明显而无法被忽略。

1.5SIR模型阈值解析

1.5.1无网络结构时的模型求解

首先对于方程(1)～(4)进行初步分析。由于流行病的有效传播率等于λ/μ，这一分式可化为x/1的形式，μ值的变化只会影响感染演化的时间尺度[1]，故而假设μ=1可不失一般性。λ/μ也是单点接触中的基本再生数R0，基本再生数R0表示一个感染个体在充分混合的S态群体中一次接触平均能感染的个体数；当有网络结构时，SIS过程的平均基本再生数R0=〈k〉λ/μ。通过这一概念可以初步定义流行病爆发阈值：只有当R0>1时流行病才可能传播开来，R0<1时流行病会在有限时间步内消亡。

(5)

要考查流行病最终是否在全局爆发(或者说流行病能否长时间的存在于网络中)，可直接分析传播过程结束达到稳态时的情况

R∞=1-S∞=1-e-λ〈k〉R∞

(6)

λc=1/〈k〉

(7)

但是稍后会发现它并不是SIR模型阈值的精确解，而是其上界。

除阈值这一问题外，我们往往还很关心最终的感染密度。由函数导数性质可知，f′(0)的值越大，其函数曲线在此处斜率越大，最终与x轴交点的值也可能更大，亦即R∞值较大，且此时随λ〈k〉的增大f(1)=-e-λ〈k〉→0。故而从平均场角度来看，λ〈k〉数值越大，则最终感染的比例R∞的数值也越大。

在λ=λc时，R∞的值较小趋近于0，|-λ〈k〉R∞|<1，故进行泰勒展开，可以得到

R∞=1-e-λ〈k〉R∞

(8)

(9)

又因为λc=1/〈k〉，可得

(10)

然而当|-λ〈k〉R∞|>1时，式中的泰勒展开将不再适用；虽然方程R∞=1-e-λ〈k〉R∞形式非常简单，然而不幸的是对于R∞并没有一个简单的闭合形式解即可由有限个初等函数与算术运算、指数、对数、根式表示，一定不包含无穷级数、连分数、积分、微分或极限)。但是仍然可以得到相对精确的数值解：将方程看作y=R∞与y=1-e-λ〈k〉R∞联立，在[0,1]区间内求解即可，这时通过程序模拟便可以较为容易地获得较为精确的结果，而且这一方法具有较好的普适性，应用范围也更广。

1.5.2复杂网络上的求解

方程组(1)～(4)从平均的角度出发，认为每个节点度均为〈k〉；然而当网络度分布取值较广时，这样的假设将无法刻画真实情况。在无标度网络中，节点度分布呈现幂律，且往往存在胖尾效应，〈k〉的近似与实际情况大相径庭，无法用平均度较好地描述节点的度值，而在ER(Erdös-Rényi)、WS(Watts-Strogatz)网络中，度分布往往呈现指数有界特性，概率密度在均值〈k〉两侧迅速减小。

为了更准确地刻画异质网络上的流行病传播过程，Pastor-Satorras和Vespignani提出了异质平均场方法，假定度为k的所有节点所处的作用环境相近、遵从相同的动力学规律(这正是在1.4节中提到的假设1)[14-15]，这样就不能在宏观的层面而须要在介观层面单独考虑度不同的节点相应的行为，此时引入变量Xk(t)表示t时刻度为k的节点中处于X(可为S、I、R)状态节点的密度。相应的SIR模型微分方程如下[23]：

(11)

(12)

(13)

1=Sk(t)+ρk(t)+Rk(t)

(14)

对于方程组(11)～(14)含义的解释与方程组(11)～(14)相同，不同之处在于对网络中节点按度分类后，考虑不同的度对于传播的影响时，ρ(t)不再能有效描述任意度为k的给定节点的一条边连接到一个感染态节点的概率；所以引入一个新的变量Θ(t)来表征这一概率。根据相应物理意义可得其表达式：

(15)

当网络无度关联性时，即度为k的节点并不以特定的概率而是会随机地连向度为k′的节点。

(16)

对于同配与异配网络条件概率P(k′|k)表达式将有很大不同。文献[24]详细讨论了网络存在度关联时的情况。

同时，式(16)还包含了平均场近似的假设3：亦即认为节点A的k个邻居(B1-Bk)之间状态互相独立，A与其任何两个邻居Bi、Bj不会构成三角形，故而B1-Bk感染与否只与A有关，是相互独立事件，所以式(16)中才能够对所有ρk′进行加和。

由式(11)～(16)可以解得未感染节点的密度函数：

(17)

其中辅助函数

(18)

Φ(t)的物理意义是到t时刻、网络中任意一条边的另一端节点处于R态的概率。当网络上的动力学在到达稳态时，流行病究竟有没有在全局爆发，就可以通过Φ(t)得到清晰的认识。为了得到感染节点密度的表达式，我们发现观察Φ(t)的均值随时间的变化率会更为方便一些，结合式(13)可以得到：

(19)

(20)

(21)

1.5.3同质网络的求解

当考虑具体的网络结构之后，可以得到更为具体的解析情况：

先来看随机网，当其规模N足够大(N≫k)时，将其度分布近似为泊松分布是较为严格、精确的：由于ER随机网络中任意两边都会以p的概率进行重连，所以某一节点其度为k的概率

(22)

通过这样一个简单的推导，在热力学极限情况下，ER随机网络的度分布就是服从均值为〈k〉的泊松分布。

将P(k)带入到等式(21)中，进而可知

(23)

网络最终感染密度

R∞=1-S∞=1-e-λ〈k〉Φ∞≈λ〈k〉Φ∞

(24)

方程中省略了二次及以上的高阶项。为了得到较高的精确度，保留最相关的项，在dΦ(t)/dt中将泰勒级数展开到二阶：

(25)

(26)

其中A为形如ec的常数。将(10)代入(8)即得

(27)

通过计算机模拟也验证了这一幂律关系(见图2)。在临界点附近λc存在幂律行为，表明这一过程是连续相变过程。

图1　模体示意图[43]

图2　阈值点附近的感染密度变化情况[23]

1.5.4异质网络的求解

(28)

其中kmax～N1/(γ-1)[51]，在热力学极限下(亦即当N→∞)，kmax→∞，分母发散，进而有λc→0。在这里也可以明显看出，在热力学极限情况下，其阈值会趋近于0。虽然现在有文章指出在SIS模型下SF网络的阈值为0，对于γ>2.5的情形阈值消失的原因主要是由于度最大的中心节点在任何传染概率之下均会处于活跃状态，而当γ<2.5时则主要是由于进行k-core分解后，k-core值最高的团簇处于活跃态所致[25,52]；但原有的解释在SIR模型中的解释仍是准确的[25]。

相应地

(29)

Φ∞在稳态时被看作一个常数，所以可将其提到积分号外部，而在下面Φ(t)的求解中，则无法采用如此假设，因为它是与k直接相关的函数：

(30)

(31)

其中,γE是相应的Γ函数。对其进行积分后可得：

(32)

其中A为一积分常数。当r→∞时，

(33)

因而结合式(29)和(33)，可以得到

(34)

对于非零的λ值，最终的R态节点密度都不为0；而这一分析也与λc的解析相一致(见图3)。

图3　感染密度与1/λ的函数关系[23]

1.6平均场方法的应用

可以发现，应用平均场方法的求解过程充满了这样或那样的近似与假设，但这也大大简化了解析的难度，让我们得以抓住最主要的影响因素，得到相对准确的结论。在上述的解析过程中，大多都是进行线性假设，舍掉了相应的高阶项；在对于阈值的解析中，由于感染率一般很小，故而舍弃高阶项产生的误差会相对较小；然而在最终感染密度的解析中，这样的线性近似会使其解析结果偏大：改进过的非微扰平均场(NonPerturbative Heterogeneous Mean-field method)保留了这些高阶项，使得解析的准确度大大提高[38]，但由于它还是假设度相同的节点动力学行为相同，进行了相应的平均代换，所以其数值结果仍与真实模拟值略有差异，而且通常比实际值略高。

普遍认为平均场方法主要有本文介绍的3个假设，但在文中也可以看到标准的异质平均场方法也没有明确地去考虑度关联性和节点的接触能力(即感染节点会接触自己邻居的多大比例，是全部抑或只是一小部分)；然而这两个因素并不是其应用的前提，因为这两种因素都可以被平均场方法捕捉到，改进后较为完整的微分方程组如下[53]：

(35)

(36)

(37)

1=Sk(t)+ρk(t)+Rk(t)

(38)

方程组(35)～(38)中表征节点的接触能力，亦即某一节点只能接触到目标节点比例的边；当A=1时即为单点接触(Contact Process：每次只选择某一邻居进行感染，若邻居已然处于感染态，则跳过感染过程)。这样一来，就将度关联性与接触能力这两个因素纳入了考虑范畴，届时只须根据具体情况给出P(k′|k)以及相应的A值即可。对于有权网络，其上单点接触往往存在依赖权重进行偏好选择，对于这一情况，只须将接触概率与感染概率分开考虑并将其相乘即可[38]。

另外对于SIS模型，由于I态节点不会变为R态，而会以一定概率μ恢复为S态，其微分方程如下[15]：

(39)

(40)

1=Sk(t)+ρk(t)

(41)

此外，当SIR传播模型中感染节点不以某一概率μ而是经过固定时间步τ后恢复时，其传播过程是可以严格对应于渗流过程的：流行病的传播过程相当于以概率p(p=1-(1-λ)τ)从源节点不断进行加边[27,54]；对于时间步为变量(或者以概率恢复，二者效果基本相同)的情况，由于引入了更强的随机性和时序关联性，而无法被渗流方法所准确描述[55]。渗流方法会在后面的章节详细阐述。

1.7平均场方法小节

纵观平均场方法，它更多的是一种思想方法而非技术手段：亦即从平均的、概率的角度去考虑问题，进行宏观的把握，把环境对物体的作用进行集体处理，以平均作用效果替代局域作用效果的加和，其具体的数学技术手段无外乎率方程、常微分方程、泰勒级数等等。

其特点正是只考虑最为主要的影响因素，考查其速率变化，略去其它相比之下影响较小的因素：如不同度之间的转换[42]、点对选取后的情况[29]、子节点对父节点的重感染概率[27]等相应因素。平均场方法正是在一系列假设与近似基础上，牺牲了一定的准确性来换取其极大的简洁性。然而是否存在整合上述多种影响因素的精确方法仍然是一个极大的挑战；是否一定是考虑的因素越多其结果就越精确，同样是一个重要的问题；探究不同因素之间的互相影响以及不同因素对于传播动力学作用的影响程度，将是非常有意义与挑战性的工作。

平均场理论被广泛地应用到流行病传播阈值及最终感染密度的解析中[15,23]，它出现之后也成为其它理论方法争相参照的对象,将平均场方法看作复杂网络动力学解析中最基础的方法应该也会被大多数人所接受。通过对于平均场方法更深入的研究和应用，也得出了许多重要结论，诸如：对于SIS模型，无标度网络在二阶矩发散的情况下，其阈值均趋近于0[24]。而对于SIR模型，则无论标度因子的大小，都是k-core值最高的团簇在影响着流行病的爆发阈值和传播过程[52]。另外高集群系数和负度关联会保护相应的无标度网络(抑制流行病的传播)[26]。

平均场方法在稳态时感染密度的解析中也占有主导地位[15]；同时在谣言传播解析中也被广泛应用[44,56]，文献[30]中在信息传播中引入了信息之间的互相竞争，其数学解析也是通过平均场方法完成的。在真实流行病的预测应用中，许田[57]应用平均场方法下的SIR模型有效地模拟了北京2003年4月18日到6月16日每天新增SARS病例统计数据的研究结果，作者得到的结果清晰地显示出数据模拟结果相当于实际数据的一种“平滑化”，意即其结果相当于全局的、平均的描述。而且在新兴的多层网络研究中，平均场方法也为进一步理解其上传播阈值减小这一动力学行为提供了有力解释[58]。

2点对近似

2.1点对近似方法概述

在上一节中，我们已经指出平均场方法误差产生的主要原因即是忽略了节点间的动力学相关性；为了更准确地描述流行病传播动力学过程，许多研究者不再假设节点间动力学无关(即节点状态与邻居节点状态无关，只考虑处于不同态节点数目的变化)，开始考虑节点状态之间的关联性[59-61]，也就是将处于不同状态节点的连边数量或比例作为一个动态参量，考虑不同类型连边数目的变化。想要得出不同类型节点对数目变化的微分方程，需要依赖于网络中不同类型的三元组数目(三元模体只是三元组的一种特例：三元组只要保证AB、BC分别相连即可，无论AC相连与否均可；三元模体则必须保证三者均两两相连)。为使相应的微分方程闭合，需要用节点对的数目来近似网络中三元组的数目，这种方法被称做点对近似(Pairwise Approximation)[7,61]。之所以用节点对数目来进行相应的近似，原因在于在非线性模型中，描述第n阶矩的微分方程会依赖于第n+1阶矩，这样得出的一系列方程是不可解的；为了使方程可解，必须使方程在某一阶矩上闭合。假设要使方程在第n阶闭合，必须用n阶及n阶以下的矩来估计第n+1阶矩，这个过程被称为矩闭合近似。

点对近似方法同样是用一组常微分方程来描述简单空间模型的动力学行为，例如简单的宿主—寄生虫模型[60,72-73]、疾病传播模型[7,42,61,64]、投票模型[65]、空间博弈[74-76]等等。它不仅描述网络中节点的状态变化，而且还同时考虑“存在连接的节点对”(下文中为简便起见,将简称“节点对”)中两个节点状态之间的相关性。这样就解决了困扰平均场方法的动力学相关性，其假设条件相对平均场方法更弱，进而可以更为准确地描述网络中的动力学过程[42,64-65]。

点对近似方法不但提高了解析的精度，而且它还具有更广的适用范围和拓展性：它既适用于静态网络，包括同质网络[42,62]和异质网络[42,64-65]；又适用于动态网络[66-68]。当网络结构并非固定不变时，例如节点或边存在增加或删除的情形时，节点与其邻居的连边会不断发生变化，此时网络结构的变化对传播过程的影响很大，所以必须考虑网络中不同状态的连边数量随时间的变化，因此点对近似方法很适于描述动态网络上的动力学过程。

下面将分别讨论在不同类型的网络上，如何使用点对近似方法来解析不同传播模型的动力学过程。

2.2静态网络上的求解

2.2.1点对近似方法预备知识

在真正给出描述动力学过程相应的微分方程之前，需要一些有别于平均场方法的预备知识。

与平均场方法相似，考虑处于不同态节点数量的变化率时，仍然使用微分方程进行描述。定义抽象函数f用于描述网络中节点状态(例如用f表示处于某种状态的节点数目(如[A])或某种类型的连边数目(如[AB]))，当网络趋近于无穷大时，可将f近似看作是连续的，那么其变化率可用如下微分方程表示为

(42)

接下来，问题就集中在如何求解[ABC]和[ABA]上。

当A≠C时，

(43)

其中,

(44)

当A=C时，

(45)

其中,

(46)

接下来分别求解Γσ(A|B|A)和Γσ(A|B|C)[62]：

=EBj=1{[Qj(A|B)-EBj=1(Qj(A|B))]2}=varBj=1[Qj(A|B)]

(47)

其中EBj=1[]表示状态为B的节点上相应状态节点数目的平均值，varBj=1[]表示方差。也就是说，Γσ(A|B|A)是网络中所有B态节点所连接的A态邻居节点的个数的方差。

下面分别来分析Qj(A|B)服从不同分布时适用的情况：

(48)

2)Qj(A|B)服从二项分布，这种情况主要适用于规则网络(每个节点有m个固定的邻居)：

(49)

与求解Γσ(A|B|A)相类似，

(50)

由此可知Γσ(A|B|C)是Qj(A|B)和Qj(C|B)的协方差。下面来考虑在B态节点其邻居中C态A和态节点的分布情况：

[ABC]=[AB][BC]/[B]

(51)

(52)

(53)

到这里，就求解出了网络中的三元组相应的数目，然而上述这些计算仍然没有涉及到网络上的动力学。下面针对不同的传播模型(SIS/SIR)来给出具体的解析。假设流行病的感染概率为β，恢复概率为γ。

2.2.1.1SIS模型上的求解

首先来看节点数目变化率的微分方程组,依然应用类似于平均场方法的线性近似：假设一个S态节点有n个I态邻居，则其被感染概率的准确值是1-(1-β)n，进而在满足约束条件(nβ≪1)时近似为nβ,进而有：

(54)

(55)

[S]+[I]=N

(56)

接下来再来看点对近似方法所特有的关于不同连边类型的变化率方程：

(57)

(58)

(59)

[SS]+2[SI]+[II]=M

(60)

其实从上述的方程组中，仍然可以发现这里面仍然有许多必要的假设以及近似：首先仍然是把每个感染或者恢复看作独立事件，这也正是每个式子中都有求和项的原因；而且由于这一假设，真实情况下可能存在的II→SS这一类的转化就没有被考虑其中。

[S]+[I]=N

(56)

[SS]+2[SI]+[II]=M

(60)

(61)

(62)

(63)

接下来，如若求解流行病爆发阈值，只需设定好相应的初始值，然后不断地迭代方程，当方程的最终解由0变为非0时，就可以数值地确定出相应的阈值。

对于SI模型，只需要将SIS模型中的恢复概率置为0即可，此处不作赘述。

2.2.1.2SIR模型上的求解

SIR模型解析思路与SIS模型情形相似，只是多了一个R态，下面只给出相应的微分方程组，不再展开后续的求解：

(64)

(65)

(66)

(67)

(68)

(69)

(70)

(71)

(72)

2.2.2异质网络上的求解

这里静态的异质网络[42,64-65]是指网络度分布较广、大多数节点的度值较明显偏离于网络平均度的网络，例如无标度网络；在此种情况下，用平均度来代表所有节点的度往往无法得出正确的估计，仍须在介观的层次上来对不同的度分别进行考虑，写出各自相应的微分方程，在这一解析思路上它与异质平均场方法是相似的。

2.2.2.1SIS模型上的求解

对于感染个体的动力学，包括感染个体以概率γ恢复为易感个体和易感个体被感染态邻居感染；这就须要把邻居节点的状态信息包含到节点的状态方程中：

(73)

对于易感个体的动力学变化，相类似地可构建出方程(76)。

(74)

为了能够求解方程(73)～(74)，接下来需要构建相应的方程来模拟节点对的动力学：

(75)

式(75)中的各项表示[SnIm]的增加量(由[SnSm]中Sm的被感染、[InIm]中In的恢复所导致)和减少量(由[SnIm]中的Sn被Im或其他与此Sn相连的I态节点感染、[SnIm]中Im的恢复所导致)。

类似地可得其他类型点对的动力学方程：

(76)

(77)

理论上来说，方程中的三元组可以用更完整的四元组来表示，更高阶的处理可依此类推。但是，这样会导致方程变得非常复杂甚至不可解，所以为了简化，类似于均匀网络的分析，进行点对近似，也就是用不同类型边(二元组)的数量来近似估计三元组的数目[61,69]：

(78)

进而依照同质网络中相应的求法，即可解出网络上流行病爆发的阈值。

2.2.2.2SIR模型求解

其相应的微分方程为

(79)

(80)

(81)

(82)

(83)

(84)

(85)

(86)

(87)

其求解过程与SIS模型相似，此处不再赘述。

2.3动态同质网络上的求解

相比于静态网络，动态网络最大的特点在于网络的拓扑结构不再一成不变，而是随着时间的演化、周围环境的变化，可能会发生相应的改变。Gross等[67]使用点对近似方法来模拟动态网络上的传播，其模型为：随机图的规模为N，平均度为〈k〉，每个节点处于S态或I态之一，在每个时间步对于每条SI类型的边，S节点以概率β被感染为I态节点，感染节点以概率γ恢复为易感态，此外对每条SI边，S节点以概率w(w∈[0,1]自适应地断开与I态节点的连边，并将断边随机重连到一个S态节点上。在这个模型中，零阶矩是感染态和易感态节点的数目，即[I]和[S]；一阶矩是平均到每个节点上的SS、SI和II边数目，即[SS]、[SI]和[II]；二阶矩是三元组的数目，即[ABC](A,B,C∈{I,S})，这里的边和三元组都是无向的(仍可将其看作双向边，这里只会影响到[M]的数值以及相应方程前面的系数“2”)。[I]及[SS]和[II]的动力学微分方程加上两个守恒条件[I]+[S]=N和[SS]+[SI]+[II]=N〈k〉/2就可以完整地刻画零阶矩和一阶矩的动力学。

由此得出的方程组(88)～(90)。

(88)

每个S态节点的邻居中I态节点的平均数量为[SI]/[S]，所以每个S态节点被感染的概率为β[SI]/[S]，S态节点密度为[S]，因此由感染引起的I态节点增加的速率为β[SI]；因为每个I态节点都会以概率γ恢复为健康态，所以由于感染态节点恢复引起的I态节点减少的速率为γ[I]。

(89)

在一次感染事件中，流行病通过一条SI边来传播，会将SI边变为II边，因此一次成功感染事件会产生至少一条II边。但是如果除了感染它的那个I态节点之外，这个新被感染节点的邻居中还有其他I态节点，那么这次感染事件还会产生其他的II边；因此，一次感染事件总共产生的II边的数目是1+[ISI]/[SI]，其中“1”表示感染发生的那条SI边，第二项表示由这条SI边形成的ISI三元组的数目。由此可得，II边产生的速率为β[SI](1+[ISI/[SI]=β([SI]+[ISI])。如果II边上一点以概率β恢复，那么II边减少；一个I态节点形成的II边的平均数量为2[II]/[I]，其中的“2”是由于一条II边连接了两个I态节点、在统计边数时须要将一条II边重复计数一次,因此II边减少的速率为2γ[II]。

(90)

SI边中的I态节点恢复会使SS边增加，相应的速率为γ[SI]，同时，每个重连事件都会增加一条SS边，所以由重连引起的SS边增加的速率为w[SI]；因为每条SS边中的S节点可能同其他I态节点相连，所以若这条SS边中的S节点被感染就会造成SS边减少，每条SS边形成的SSI类型的三元组数目为[SSI]/[SS]，由此可知，SS边减少的速率为β[SSI]。

此时，这个方程组仍不是一个闭合的模型，它们还依赖于两个未知二阶矩[SSI]和[ISI]。而方程中的一阶矩已经能够很好地反映重连的效果，因此使用点对近似方法来使系统闭合：用零阶和一阶矩来近似地估计二阶矩。

ISI三元组的一半是一条SI边，其密度为[SI]。为了近似一条给定SI边形成的ISI三元组数量，必须计算出连向此S节点的其余I态节点的平均数量。假设除去给定SI边之外，S节点平均还有〈q〉条边，其中〈q〉表示平均剩余度，即通过一条随机选择的边到达的节点它除去此边之外拥有的平均边数。又因为这〈q〉条边中每条边是一条SI边的概率为[SI]/〈k〉[S]。已知SI边的密度、SI边中S节点有其他SI边的概率，所以[ISI]=κ[SI]2/[S]，其中κ=〈q〉/〈k〉。初始网络是随机网络，剩余度分布q(k)=(k+1)P(k+1)/〈k〉，平均剩余度

(91)

点状图对应数值模拟结果，实线为方程

又因为P(k)服从均值为λ的泊松分布(P(k)=e-λλk/k!)，因此相应的均值E(k)和方差D(k)=E(k2)=[E(k)]2=〈k2〉-〈k〉2的关系为E(k)=D(k)=λ=〈k〉，由以上两式可知，〈k2〉=〈k〉2=〈k〉，因此〈q〉=〈k〉；在此情况下κ=1，兼由式(48)可得[ISI]=[SI]2/[S]。类似的，由式(51)可得，[SSI]=2[SS][SI]/[S]，其中的“2”表示一条SS边中的两个S节点(在有向图中将没有这个系数“2”)。由此可得到一个闭合系统的微分方程组：

(88)

(92)

(93)

其中圆圈表示实际模拟的结果，方框表示用点对近似解出的结果，此处初始感染比例i(0)=0.3，感染概率p=0.008，重连概率w=0.3

2.4点对近似方法适用条件及准确性概析

点对近似方法适于对没有较短回路的网络结构进行解析，它对拓扑多为局域树形结构的网络动力学过程近似非常准确。当网络中没有强连通模体时，若i和j直接相连，同时节点i与j存在另一条非直接相连的路径，即i和j之间存在的回路中较长的那条路径，这条路径越长，由它导致的节点i和j的状态之间的关联性就越小，越可以忽略，因此点对近似方法的近似较准确。而且点对近似方法扩展性极强，可以处理几乎各种类型的网络。但是对于动态网络中，当在演化过程中会出现较强的关联性时，点对近似就不能比较准确地刻画其演化情况，如图5对比了自适应模型[67]在双稳态参数下感染i(t)比例随时间演化的模拟结果和点对近似分析的结果，由图可知点对近似此时解析不准确，因为网络的重连导致网络中出现了很强的社区性，关联性也随之变强。

然而许多实际社会网络中都存在很多回路，存在大量模体时[70-71]，点对近似方法同样无法精确刻画，其误差仍然存在；对于模体这一因素，平均场方法与点对近似都无能为力。点对近似对于平均场方法最大的改进就在于其加入了对动力学相关性的考虑，使得分析结果更为准确；同时，点对近似方法具有很好的拓展性，适用于许多模型(例如简单的宿主-寄生虫模型[60,72-73]、流行病传播模型[7,42,61,68]、投票模型[65]、空间博弈[74,76]等等)。

3主方程方法

3.1主方程方法概述

到目前为止，应用微分方程的解析方法中，最为精确的当数主方程方法，这是马尔可夫过程概率分布的演化方程，研究概率分布随时间的演化。主方程是统计物理里最重要的方法之一，它几乎是普遍适用的，广泛地并应用于化学、生物学、金融、人口动力学、复杂网络动力学等诸多学科领域，以及布朗运动、流体、半导体等具体问题的解析当中。

在前面介绍的平均场和点对近似方法都假设度相同的节点动力学行为相似，而点对近似却更近一步考虑了节点对间的动力学相关性。为了更精确地描述传播动力学过程，Gleeson更深入一步，不仅考虑度为k的节点的状态变化，而且还考虑其邻居的状态变化，提出了相应的主方程方法来对复杂网络上SIS模型传播动力学进行解析[42]；该方法基于Marceau等[77]提出的思路，并进行了相应改进。实验结果表明主方程方法比平均场和点对近似方法得到的解析结果更准确；而且可以通过适当的近似，退化为点对近似方法和平均场方法的演化方程；另外此方法还可用于对自旋模型等其他只包含两个状态转换的动力学过程进行解析。

3.2对SIS模型的求解

给定网络的度分布，利用配置模型[78]可以生成相应的静态无向网络。在SIS模型中，每个节点要么处于易感状态S，要么处于感染状态I。

Gleeson提出的主方程方法假定节点的感染概率和康复概率依赖于节点的度k和感染邻居数m，令Fk，mdt表示一个度为k且在t时刻具有m个感染邻居的易感节点，在t+dt时刻变为感染态的概率，令Rk，mdt表示一个度为k且具有m个感染邻居的感染节点，在dt时间间隔内恢复为易感态的概率。对于SIS模型，有Fk，m=λm和Rk，m=μ[1]，其中λ和μ分别表示感染率和康复率。

对于每个m值(0≤m≤k)，sk,m(t)和ik,m(t)演化的主方程为

(94)

(95)

其中，每个方程右边的前两项分别表示度为k的易感节点被感染和感染节点康复为易感态；其余4项考虑其邻居的状态变化(S态邻居节点被感染或者I态邻居节点康复)：其中βs(k-m)sk，m表示度为k的易感节点其k-m个易感邻居节点变为感染态的个数，βs表示在dt时间间隔内易感节点的易感邻居变为感染态的比例；相应地，γs表示dt时间间隔内易感节点的感染态邻居康复为易感态的比例；βi表示dt时间间隔内感染节点的易感邻居变为感染态的比例；γi表示dt时间间隔内感染节点的感染态邻居康复为易感态的比例。其余的4项旨在考虑度相同但邻居状况不同的节点，它们产生哪些变化后会对sk，m的数目产生影响。在解析相应邻居节点变化率的过程中，研究者实际上应用了点对近似和平均场的思路：式中βs表示在dt时间间隔内整个网络中所有sk，m中度为k的易感节点被感染后，所导致的原有S-S边变为S-I边占原有S-S边的比例，亦即对于网络中任意一条S-S边其在dt时间间隔后变为S-I边的概率；相应地，γs表示dt时间间隔内整个网络中所有ik，m中度为k的感染态节点恢复后，所导致的原有S-I边变为S-S边占原有S-I边的比例；βi表示dt时间间隔内整个网络中所有sk,m中度为k的易感节点被感染后，所导致的原有I-S边变为I-I边占原有I-S边的比例；γi表示dt时间间隔内整个网络中所有ik,m中度为k的感染态节点恢复后，所导致的原有I-I边变为I-S边占原有I-I边的比例。

我们发现，βs虽然是关注sk,m中那个度为k的s态节点被感染，但是由于是对于全网中所有S态节点进行了考虑，所以可以用这样一个概率来描述任意一条S-S边中某一端S节点被感染的概率。

而且通过这样一个式子，也可以发现为什么系统(1)～(2)中在考虑sk,m的变化率时只考虑了sk,m-1,sk,m+1，而没有考虑sk,m-2,sk,m+2及之后的项，因为后面的项它们本身的转换概率阶数更高，将其忽略所产生的影响将非常小。该方法对节点的分类较点对方法更为细致，而且其所进行的近似也非常合理，这也是它高精确度的原因所在。

上述这4个参数的值可以通过跟踪网络中每条边的两个节点的状态变化计算得到：例如对于βs，首先计算t时刻两个节点均是易感节点的边数(即S-S边数)，其次计算在dt时间内由S-S边变为S-I边的数目，并且对于所有不同的度都进行计算后求取平均，后者与前者的比值即是βs；其它3个参数γs、βi、γi通过相同的计算方法得到。βs、γs、βi、γi这4个参数的计算公式为

(96)

(97)

(98)

(99)

给出一定的初始条件，数值迭代上述方程即可得出其感染密度，当感染密度从0变为非0时即可求得其爆发阈值。

3.3主方程方法近似条件

主方程方法与点对近似、平均场理论都使用微分方程描述相应的动力学变化，而且它也应用了前两者的一些解析思想，通过适当的近似，就可以得到标准点对近似以及平均场的方程。

我们仍然可以通过考虑点对的情况，亦即pk(t)(对于任意度为的易感节点,随机选择的一个邻居在t时刻处于感染态的概率，亦即某一点度为k的S-I边的比例)和qk(t)(对于度为k的感染节点,随机选择的一个邻居在t时刻处于感染态的概率，亦即某一点度为k的I-I边的比例)，来对系统(94)，(95)进行近似推导出点对近似方法。

(100)

(101)

(102)

更进一步，可以从平均场的角度，利用先前1.5小节中所介绍过的概率Φ(在SIS模型中的含义为：在任意时刻,随机选择一条边、其另一端节点处于感染态的概率)对系统(94)，(95)进行近似，推导出平均场近似方法；当网络不存在度关联性时，Φ=P(k)kρk/k=〈ρkk/〈k〉〉。利用概率Φ代替系统(100)～(102)中的pk(t)和qk(t)，亦即将sk，m≈(1-ρk)Bk，m(Φ)和ik，m≈ρkBk，m(Φ)代入系统(94)，(95)即可推导出平均场近似方法[79]：

(103)

系统(103)的初始条件为ρk(0)=ρ(0)。当k的变化范围是[0,kmax]时，系统(103)至多包括kmax+1个方程。进行进一步展开，即可得到

(104)

这与1.5节中式(12)是一致的(SIS与SIR模型中ρk变化率的方程是相同的)。

3.4主方程方法小节

Gleeson提出的对SIS传播模型进行解析的主方程方法，是目前基于微分方程组解法中最为精确的方法，它不仅考虑了节点状态的变化，而且还考虑节点了的邻居节点状态的变化，将动力学演化过程考虑得更为细致。主方程方法比点对近似方法和平均场方法的解析结果更为准确(如图6所示)，它解析出的感染密度随时间的演化情况和模拟结果非常接近。而且通过适当的近似，可由主方程推导出点对近似方法和平均场方法。另外该方法不仅可以解析SIS传播动力学，它还适合于其它只包括两个状态的动力学，例如自旋模型、投票模型[79]等。

主方程方法的不足是系统中演化方程的个数随最大度的平方规模增长，当最大度的值较大时，此方法的计算复杂度较高。

4边渗流方法

4.1边渗流方法概述

从图论的角度来看，渗流理论描述的是在一个随机图中的连通簇的行为。在物理、化学和材料科学中，渗流理论被用来描述有孔材料中液体的流动，例如岩石中的石油等。同时渗流理论还可以用来研究路由器网络故障、森林火灾、网络上流行病的免疫、流行病传播。在复杂网络动力学研究中，渗流方法通常分作边渗流与点渗流，点渗流即为删除网络中的节点(同时删除连接到这个节点的边)的过程，例如网络中路由器的故障、在流行病传播过程中个体的免疫，都可以用点渗流来描述。而边渗流即为删除节点之间连边的过程。本节将主要介绍边渗流的方法在求解SIR及SIS模型阈值上的应用，边渗流方法从网络中边的角度去考虑网络中流行病的传播。

对于SIR模型，系统的长期行为(例如网络中最终的感染密度等不包含时间演化过程的性质)可以近似等价地映射为网络中的边渗流过程。用渗流来研究流行病传播的方法最早由Mollison[80]和Grassberger[81]提出，并由Newman等[82-84]发展完整。本节会简要回顾Newman在任意度序列的随机图上得到的一些基本结论[7]。另外，先前的理论认为只有SIR模型传播动力学才可以映射为渗流过程，然而有研究者发现只需要引入边的重感染概率[27]就可以解析SIS模型的传播动力学。

4.2对SIR模型的求解

在SIR模型中，感染节点以λ的传播率传播流行病，并且假设节点的染病周期为τ，即节点经历确定的τ时间步后从I状态恢复为R态(这与节点以1/τ的概率恢复有较大的差异，以概率恢复会导致更大的随机性，一定程度上破坏了平均假设，会使得边的渗流概率呈现出一定的异质性)。假设动力学演化过程中时间步是连续的(对于离散情况引入一个无穷小量Δt→0即可)，那么对于任意一条边，流行病在τ时间内没有通过其成功传播的概率为

(105)

那么，流行病通过任意一条边成功传播的概率为

φ=1-e-λτ

(106)

如果令网络中的边分别以概率φ和1-φ被占据和不占据，那么流行病传播过程就相当于一个边渗流概率为φ的过程：被占据即意味着这条边足以传播流行病，不被占据则反之。当一条边被占据，意味着当流行病达到这条边任意某一端的时候，这条边足以传播流行病，而非一定会传播流行病。因此，被占据的边意味着潜在的传播可能性。当流行病从初始感染节点开始传播，可以看到流行病最终可以传播到的节点恰好就是那些与初始节点之间有被占据的边组成通路的那些节点。因此，传播的最终状态就是初始节点所在的边渗流簇中的节点处于恢复态(R)。由此，流行病爆发的阈值与相对应的边渗流的阈值是相同的。

以任意度序列的随机网络为例。设网络的度分布为pk，u为一个节点不通过某条边连接到最大簇的平均概率。这里假定每条边的占有概率是相等的，那么，计算u时有两种情形须要考虑：一是这条边没有被占据；二是这条边被占据，但是这条边另一端的节点同样不属于最大连通簇。而第二种情况当且仅当另一端的节点不通过任意一条其他的边连接到最大连通簇才成立。如果另一端的节点有k条剩余边(亦即除了连接到节点的这条边之外其它边的数目)，那么发生的概率为uk。因此不通过某条边连接到最大簇的总概率为1-φ+φuk。因此，针对k取平均，就可以得到一个自洽方程。

(107)

qk=(k+1)pk+1/〈k〉

(108)

根据式(106),(108)及式(107)导函数可以图像地解出传播相变的阈值。式(107)的解为y=u与y=1-φ+φg1(u)这两个函数的交点。显然u=1恒成立为一个平凡解。只有当存在u=1之外的解时，网络中才存在一个最大连通簇使得流行病爆发。因为qk为概率必有qk大于0，所以对于u≥0，g1(u)及其各阶导数必然非负。因此，1-φ+φg1(u)与u在u=1处恰相切(此即为相变点)，对式(107)等号两边进行求导可得

(109)

再将式(106)带入到式(109)中，即可得到

(110)

因此当λτ高于这个阈值ln(〈k2〉-〈k〉)/(〈k2〉-2〈k〉)的时候，即很有可能发生流行病的爆发，因为初始感染者有可能没有落在最大连通簇里。而当λτ小于这个值的时候，流行病则一定不会爆发。

当然，在上面介绍的经典渗流方法中，将SIR模型映射为边渗流只能反映系统长期的行为，而无法包含随时间实时演化的信息，这也是上述方法最大的局限所在[83,85]。

4.3对SIS模型的求解

先前鲜有理论将SIS传播动力学映射到渗流过程，原因在于节点被感染后会以一定概率恢复，这样流行病的传播就不再能精确地等效于一个简单的加边过程。但是Parshani和Havlin等[27]创造性地运用子节点对父节点的重感染概率，有效地克服了这一问题，对于SIS模型传播阈值的预测也比一般的平均场方法更加准确。具体做法为

首先计算每一条边被占据的概率

p=1-(1-λ)τ

(111)

然后从一个相对微观的角度来考虑流行病传播的过程：当网络不存在度关联性时，对于任意一个节点j，它被度为k的感染节点i(节点i会有1条入边，k-1条出边；此处入边与出边是针对动力学过程而言的，网络仍然是无向网络)触及到的概率是kP(k)/〈k〉。只要网络是树形结构(或者环极少进而可以被忽略时)，在SIR模型中，平均说来被节点i所感染的节点数目[85]

(112)

为了简便起见，在后面的公式中沿用这一代换κ-1=〈k2〉/〈k〉-1。

然而对于SIS模型，上述的方程不再准确。在SIS传播模型中，流行病传播过程不再是一个有向无环的树形分支过程：网络中的连边可能会多次地传播流行病，更为重要的是先前渗流的边也会恢复成非占有态。因此，直观地说来这应当会导致最终稳态下某条边渗流的难度增大，使SIS模型的阈值在同等参数条件下相较SIR模型要小些，也就是说SIS模型中的流行病更容易大范围地传播开。对应到具体的数学解析上，研究者引入子节点对于父节点重感染概率π来克服上述的问题，这样相应的渗流方程为

〈ni〉SIR=p(κ-1)+π

(113)

这里对于重感染概率的求解是关键点。大道至简，复杂的问题往往能够溯源到一个最简单的模型：

以图7中的i节点为例，它在感染j1节点之后，会在之后一定时间步之后恢复，具体值取决于i感染j1发生在i被感染后的第几个时间步，设这一变量为s，那么在接下来的(τ-s)步时间内，i仍处于感染态，而i恢复之后j1尚处于感染态，因此i、j1在感染态上存在一个时间差τ-(τ-s)=s。在这s时间步内，j1节点在每一时间步都会以概率λ去感染i节点。据此写出子节点j1感染其父节点i的概率的方程

图6　感染密度的时序演化[42]

图7　感染过程示意图[27]

(114)

从式(114)可以看出其过程，由于i节点已将j1节点感染所以整个式子首先是一个条件概率P(j→i|i→j)，节点i感染j1节点的概率为1-(1-λ)τ，故而在条件概率中这一项会出现在分母上；再来看分子，既然i在第s步感染了节点j1，那么在前s-1步，都未感染成功：(1-λ)s-1λ；再接下来，就是j1在余下的s步时间内感染节点i：1-(1-λ)s。式子展开后是一个典型的等比数列，化简后即得右式。

然而这一过程仍然略为理想，因为除去j1节点，i节点的其它子节点j2—jk以及其父节点a也会以πt的概率去感染i节点，所以在i的邻居节点之间存在一个竞争

(115)

其中πt表示了j1感染其父节点i的概率，求和符号中每一项表示除了节点j1外，i感染了它的k′个邻居并被它们重感染的概率。对于这k′+1个重感染节点，仅有第一个感染i的节点会产生作用，其他k′个邻居便不会再有任何影响。既然这些节点都处于等同地位，那么将相应值除以k’+1就可以得到先前被i感染的子节点重复感染i的概率。

当〈ni〉SIS=1时，由式(111),(113)～(115)即可解得阈值λc。并且在图8中依次比较了不考虑重感染概率(p对应的曲线)、考虑重感染但不考虑竞争(πt对应的曲线)以及考虑竞争效应的重感染情况(π对应的曲线)。另外Parshani和Havlin等[27]还研究了不同网络结构、以及SIS与SIR模型不同的动力学对于传播阈值的影响(见图9)。

图8　3种近似方法对应的理论观测值(相应曲线)与模拟结果(点状图)[27]

图9　不同情况下流行病爆发阈值与的关系[27]

4.4边渗流方法的适用性分析及最新进展

在原有的理论框架下，当节点感染后不以一定概率而是经过特定时间步后恢复为健康态，则SIR模型下的流行病传播可以映射为边渗流[55]；当节点以一定概率恢复时，其情况较为复杂；在4.3中介绍的方法认为重感染率π是此种概率恢复情形中的上界，故而仍可较准确地描述此类的传播问题[27]。

总体来说边渗流方法相较于传统的平均场方法在解析SIR模型下传播阈值时更为精确。但是Newman在文献[82]中的经典方法也有一定的局限性：例如对于“网络大小有限、网络中存在簇会导致较强的动力学相关性)、流行病传播过程随时间演化、感染和恢复概率异质不恒定”等问题都不能求解。然而Newman在文献[82]中的方法中的这些局限性在研究者们的不断努力下，一定程度上都得到了克服：Marder[85]将边渗流的方法做了适当的改变，使得能观察到流行病随时间的演化过程；Pierre-André No⊇l等在文献[86]的基础上研究了有限大小网络上流行病的演化过程；Miller[88]研究了对于异质的感染和恢复概率的问题，求出了感染比例的上界和下界。Newman[89-90]给出了网络中存在聚类时的解析方法，文献[83],[91]研究了两种流行病在单个网络上传播的情形。Allard等[92]则研究了多重属性网络上的流行病传播，即使网络中每条边的占有概率都不同，也能精确求解其阈值。Funk等[93]将此方法推广到重叠网络上，研究了两种流行病在网络上的传播过程。Dickison, Wang Y等[94-95]则将渗流的方法应用到理解双层互连接网络上的流行病动力学。可以说渗流方法自身的扩展性较强，能够克服自身先前的各种局限，在传播动力学的解析方面也发挥着越来越重要的作用。

5空穴理论

空穴理论最早由Perl[96]提出，尔后由Mézard等[97-98]将其应用到统计力学中。物理学中的空穴是指电子挣脱价键的束缚成为自由电子后，其价键中所留下的空位。空穴理论在物理、化学、生物、计算机科学中有着广泛的应用，仅在网络科学领域中，就有用于解析流行病传播动力学[99]、计算有向网络中环的个数[100]、求解图的二分问题[101]、计算最大连通子图[102]、分析网络中节点的作用[103-104]等等。

利用空穴理论的思想求解SIR流行病传播动力学最早是由Karrer等[99]提出的。他们在求解过程中，会将原先的无向网络转换成有向网络，其中边的方向表示流行病的传播方向。由于处于S态节点不能感染它的邻居节点、只能被它的I态邻居节点感染，也就是说S态节点只有入边而没有出边。借用空穴理论的思想该S态节点处于空穴态，对应到物理过程中处于空穴态的节点只能作为受体接触电子。因此，求解一个节点t在时刻处于S态的概率，也就是计算它在t时刻处于空穴态的概率，再根据SIR模型中3种状态的转换关系就可以求解出节点处于每种状态的概率以及流行病传播的爆发阈值。

此前研究SIR模型的传播动力学时，为了简化解析过程，通常会假设恢复时间服从指数分布，恢复过程为泊松过程[23,105]，然而这一假设与现实生活很不相符。举例来说，季节性流行病恢复时间服从泊松分布而并非指数分布，其均值通常在一周左右[99]。这使得人们对流行病传播过程的描述无论从定性理解还是定量分析上都有一定的偏差，Karrer等[99]利用空穴理论的思想解决了这一难题。

图10　空穴理论网络示意图

5.1SIR流行病模型求解过程

本小节介绍Karrer等[99]如何利用空穴理论求解无向网络上SIR模型流行病动力学。如图10a所示，流行病在传播过程中，节点既可能被邻居节点感染，也可能感染邻居节点。为了描述流行病的传播方向，他们将网络中的无向边变为两条有向边(如图10b)，边的方向表示流行病的传播方向。由于处于S态节点不能感染邻居节点，只能被邻居节点感染。因此，删除从它出发的所有边，则该节点处于空穴态(如图10c中节点i为空穴态)；此时，原来的无向网络已变为有向网络。求解任意节点在t时刻处于S态的概率时，只需要考虑它在初始时刻处于S态(空穴态)的概率且在t时间步内没有被邻居节点感染的概率即可。接下来将考察两种网络上的流行病传播动力学：树形网络和非树形网络。

5.1.1树形网络上的求解

所谓树形网络，是指网络中不存在有限环即环的边数为有限大小)。例如配置网络中节点之间随机连边，且网络规模趋于无穷时，就可以近似地看作树形结构或者局部树形结构[7,82]。相对于非树形网络而言，树形网络上的流行病传播研究过程要简单一些。如果节点i在t时刻处于S态(空穴态)，那么只可能是因为它在初始时刻处于S态并且在之后t时间步内没有被邻居感染，因而相应概率为

(116)

其中，Si(t)表示节点i在t时刻处于S态，P(Si(t))表示相应概率，简记为P(Si)；z表示节点在初始时刻处于S态的概率；N(i)表示节点i的邻居构成的集合；Hi→j(t)表示节点在j时间t内没有成功感染节点i的概率。等式(116)给出了节点i在t时刻处于S态的概率，但是式子含有未知量Hi→j(t)。

根据动力学机制可知，只有I态节与它的S态邻居接触才可能传播流行病。因此，一个节点在被感染后在t到dt内，感染它的一个邻居的概率为

(117)

对式(117)左右两边对t积分，可得一个I态节点在恢复之前将流行病传给它的S态邻居的概率，即Newman等[82]所说的有效传播率。其中，s(t)dt表示被节点感染后，t到t+dt内感染一个S态邻居的概率；r(t)dt表示节点被感染后，t到t+dt内恢复的概率。

(118)

(119)

对等式(119)进行数值求解，就可以得到节点i在t时刻处于I态的概率。由于节点在任意时刻只能处于3种状态中的某一种，即P(Si)+P(Ii)+P(Ri)=1成立，再联立等式(116)和等式(119)就可求出节点i在t时刻处于各种状态的概率。

5.1.2非树形网络

在实际网络中，由于边的连接并非随机，树形网络几乎不存在[7,99]。即使在配置网络中，当网络规模有限时也会导致网络中出现有限环[108]。空穴理论只能精确求解树形网络上SIR模型，对于非树形网络只能给出流行病传播范围的上界。下面将介绍利用空穴理论求解非树形网络上SIR模型的传播动力学。

类似于树形网络上流行病传播的求解过程，首先将原来的无向网络变为有向网络，边的方向表示流行病传播方向，S态节点在传播过程中任意时刻都为空穴态。此时，还需要对节点和边都赋予一定的权值。其中，节点的权值τi表示它被感染后所需要的恢复时间为τi；有向边的权值wij表示j节点被感染后，感染它所指向的节点i所需要耗费的时间为wij。根据恢复时间分布函数r(τ)对每个节点赋予权值，同理根据分布函数s(τ)对每条边赋予权值。在对边赋予权值时需要注意，若wij<τj表明节点i被节点j感染时，节点j还没有恢复为R态，这表明节点j可以成功感染节点i，此时为有向边j→i赋值wij；否则，表示节点j在其处于感染态的时间步内无法将节点i感染，此时令wij=∞。

假设网络中任意两点l，k存在一条从l到k的有向路径，其路径长度Lkl表示节点l被感染后节点被k感染需要等待的时间为Lkl。因此，求解节点i在t时刻处于S态的概率就变得相当容易，即只需要满足：它在初始时刻没被感染，相应概率为z；不存在初始时刻被感染且它到节点i的有向路径长度小于t的节点j。

不妨以节点i为中心，将所有到节点i有向路径小于t的节点构成集合记为Ω(ni表示集合Ω中节点数目)。因此，集合Ω中任意一个节点k到节点i的有向路径都小于t。若集合Ω中所有节点在初始时刻都没有被感染，那么节点i在t时刻也不会被Ω中的节点感染。此时，节点i在t时刻处于S态的概率为zni+1。对所有的wij和τj进行平均，则节点i在t时刻处于S态的概率为

P(Si)=z〈zni〉

(120)

其中，〈…〉表示对所有wij和τj进行平均。

(121)

图11　网络上的流行病传播示意图

由切比雪夫不等式[99]，式(121)可化为

(122)

其中Hi←j(t)=〈znij〉，表示时间t内节点j没有将流行病传给节点i的概率。

在计算P[Sj(t′)|i处于空穴态]时，可以利用不等式(122)。但是，此时节点i和节点j都处于空穴态，需要做适当的处理，即添加从节点i出发到除节点j之外的所有邻居的向边(图10d)，相当于原网络中只删除了从节点i到节点j的有向边。在做此处理后，势必增加了到达节点j的路径数，从而增加了节点j被感染的概率。于是有

P[Sj(t′)|i处于空穴态]≥P[Sj(t′)|删除从i到j的有向边]

(123)

再联立不等式(122)，有

(124)

成立。不等式(124)和等式(118)形式很相像，只是这里是不等式。因此，无法精确计算出网络中S态节点的比例。

(125)

经过一次迭代后，令不等式(125)右边为

(126)

它对于任意的τ都有f(τ)≥0成立，根据不等式(126)有表达式

(127)

成立。再根据不等式(126)有表达式

(128)

成立。经过m次迭代，得到关系式

(129)

(130)

利用等式(122)，有

(131)

表达式(131)给出了在非树形网络中P(Si)满足的关系，进而可以求解出P(Si)最小值和P((Ii)+P(Ri)的最大值，即网络中被感染节点的最大比例。

5.1.3配置网络上的流行病传播

配置网络的度分布P(k)的生成函数为

(132)

剩余度分布q(k)的生成函数为

(133)

因此，在任意节点在时刻t，节点没有被它的某个邻居感染的平均概率

(134)

进一步可知，t时刻网络中处于S态节点的比例为

(135)

(136)

即为没有被邻居感染的概率。根据表达式(131)t时刻节点i处于S态的概率为

(137)

当t→∞时，网络中没有节点处于I态，只有S态和R态节点，即P(Ri)=1-P(Si)成立。对于配置网络，当网络规模很大时可以近似看成树形网络，此时Fi←j和等式(116)的Hi←j含义相同，而非树形网络根据表达式(130)都有Fi←j≤Hi←j成立，即Fi←j是Hi←j的下界。因此对于树形网络，可精确地求出节点在稳态时处于S态的概率；然而对于非树形网络而言，只能近似求出节点i在稳态时处于S态的最小概率，即处于R态的最大概率。

当t→∞时，等式(134)变为

H1=1-p+pzG1(H1)

(138)

等式(135)变为

P(S)=zG0(H1)

(139)

其中，H1=H1(∞)表示一条边在整个传播过程中没有被传播流行病的概率。根据等式(138)和等式(139)就可得稳态时处于S态和R态的比例。根据等式(138)，可用图解法求出pc=〈k〉/(〈k2〉-〈k〉)，这和Newman等[82]利用渗流理论得到的阈值相同。当p>pc时，网络中就会有一定比例的节点被感染；反之，则几乎没有更多的节点被感染。

最后，实验模拟结果也非常精准地验证了之前理论分析结果的正确性与准确性。图12给出了在ER网络上感染密度随时间变化的理论值和实验结果。其中γ和β分别表示恢复率和感染率，因此r(t)=γe-γ t，s(t)=e-β t[7]。根据式(117)有f(t)=βe-(β+γ)t，将t′=t-τ带入等式(134)可得

(140)

在等式(140)两边同时对t微分有

(141)

网络节点数量N=105，平均度为3，模拟重复100次以上。红色实线对应根据等式(135)在每个t时刻求得的理论值，点状图对应模拟值

初始条件H1(0)=1，因此等式(141)的解

所以数值地求解出H1的值，将其带入等式(138)即可求解出网络稳态时S态节点的比例，进而也能求出I态和R态节点的比例。

5.2空穴理论总结及展望

对于流行病动力学的研究，最为重要的问题便是考察在什么条件下疾病会成为流行病，以及流行病的影响范围。此节介绍的Karrer等[99]利用空穴理论求解SIR模型的方法，可以求解感染时间和恢复时间服从任意分布的情况，且能精确求解树形网络上的流行病传播范围以及动力学演化过程。空穴理论无疑成为人们认识流行病动力学的一个新工具。当然也应该清楚地认识到，对于非树形网络空穴理论只能求解出最大感染比例，而无法给出准确的结果。

张昊等[107]在Karrer等[99]和Shiraki等[102]工作基础之上，对节点度之间存在关联性时的流行病控制问题做了研究。他们利用空穴理论得到的解析结果和仿真结果非常吻合，在一定程度上加深了我们对于流行病传播的认识。Altarelli等[108]最近利用空穴理论解决了传播源优化问题。当然，空穴理论对于SI、SIR、SEIR流行病传播模型也都适用[99]。

6边划分方法

在研究复杂网络上的流行病传播时，人们往往会关注动力学过程、稳态时感染密度和传播阈值。异质平均场、马尔可夫过程、点对近似等建模方法已成为研究流行病传播的主流工具。然而，这些方法在描述动力学过程时，微分方程组维数很高，难以求解和找到对应的解析表达式。Volz[109]利用边划分的方法很好地解决了这一难题。他用低维非线性方程组描述了复杂网络上SIR流行病传播动力学模型，可求解在任意时刻流行病感染传播范围。

6.1配置网络上SIR流行病传播

在研究流行病传播模型时，非常重要的问题之一就是随机选择一节点u，它处于S、I和R态的概率分别是多少。由于节点u是随机选取的，当网络规模很大时它处于S态的概率就为网络中S节点所占的比例S(t)，处于I态和R态的概率也是类似的。如果知道S(t)，就能很容易求出I(t)和R(t)。

节点u处于S态只有它的所有邻居节点到此时刻都没有成功感染它，其概率大小取决于它的邻居数量、I态邻居的变化率和邻居处于I态的概率。由于节点u的二层邻居数通常大于它的平均度，因此即使知道网络中I态节点比例也不知道节点u的邻居处于I态的概率。因此，求解节点u的邻居节点处于I态的概率也就成为重点之一。一旦求出此概率，就很容易得到节点u处于S态的概率。

假设节点u的每个邻居相互独立。然而，节点u的一个邻居被感染，它可能感染节点u，而节点u又可能会感染邻居节点w，使得假设不成立。这里只需要假设节点u不传播流行病即可保证它的邻居之间相互独立。虽然做了如此假设这并不影响节点处于各种状态的概率。因为，节点u处于S态时，它的邻居之间没有关联性，因此不影响它处于S态的概率。当它处于I态时，邻居节点不会相互独立，但是这也不影响它处于I态的时间和恢复概率，因此不影响它处于I态和R态的概率。因此，不允许节点u传播流行病，并假设它的邻居相互独立并不会影响它处于各种状态的概率。

在介绍边划分求解方法时，还有必要说明只有当初始时刻有足够多的I态节点使得流行病在传播的过程中不会受到随机扰动而受影响时才会很准确。当然，这里足够多的I态节点相对于整个网络而言也是非常少的。

(142)

(143)

再联立

(144)

为在配置网络上SIR流行病的演化方程。根据等式(143)和等式(144)就可以求出在任意时刻节点处于各种状态的概率。在这里只需要求解一个方程，其求解时间相对于异质平均场、点对近似和马尔可夫等方法更少。

(145)

(146)

当R0>1时，等式(146)有两个解，其中有一个解为1，另一个解小于1，小于才是此时寻求的解。再根据等式(144)可求出稳态时流行病传播范围R(∞)=1-G(θ(∞))。

初始时刻选取1%的节点作为传播源。理论值为虚线，模拟值为实线

6.2最新研究进展

在Volz 提出边划分方法后，Miller[110-111]对此方法做了进一步简化，得到方程和。最近，Miller、Slim和Volz在综述文献[112],[114]中详细介绍了利用边划分所解决的流行病传播问题。他们在文献[112]中主要介绍了利用边划分如何求解不同度分布、不同接触模式下的SIR流行病传播。文献[113]中主要介绍对于不同的情况如何选取不同边划分模型，而文献[114]重点在于介绍在其他情况下的应用。比如，有向网络上、异质的感染率和恢复率[88]、节点类型多样性网络[115]、动态网络[116]上的流行病传播。Volz等[117]还利用边划分方法研究了异质网络和时空聚集模式，以及初始感染比例对动力学的影响[118]。不仅如此，Valdez等还利用此方法对自适应流行病传播模型[119]、动态渗流[120-121]做了相应的研究。然而，这里介绍的边划分方法并不适用于SIS模型，杨紫陌与周涛[48]利用另外的边划分方法求解权重网络上SIS单点接触流行病传播模型。

7马尔可夫链方法

7.1马尔可夫链概述

马尔可夫链方法因安德烈·马尔可夫(Markov A. A.)得名，在数学中是指具有马尔可夫性质的离散时间随机过程：亦即在这一随机过程中，预测未来的状态只依赖于先前的部分状态信息，而与更早之前的历史状态无关。通常可分为一步条件转移概率与N步(或称多步)条件转移概率，前者下一时刻状态只与当前状态和前一时间步有关；后者则与之前N个时间步状态都有关,例如物理学中矩磁的相变过程。对于服从一步条件转移概率的指定个体，已知当前状态信息，根据马尔可夫过程便能预测下一时刻的状态信息。另外特别值得一提的是，马尔可夫方法可以处理离散时间，不需要像其它方法假设时间连续。

马尔可夫链是描述对象的一种状态序列，其每个时刻的状态值取决于前面的状态，是具有马尔可夫性质的随机变量{X1，X2，…，Xt}的一个数列。这些变量的取值范围(即它们所有可能取值的集合)，被称为“状态空间”(或“相空间”)；Xn的值表示相应变量在t时刻所处的状态，如果Xt+1对于过去状态的条件概率分布仅是Xt的函数：亦即P(Xt+1=x|X1=x1，X2=x2，…，Xt=xt)=P(Xt+1=x|Xt=xt)，其中x的不同取值对应过程中的不同状态；那么上面这个恒等式可以被称作具有马尔可夫性质，更为具体的是指“一步条件转移概率”。

马尔可夫链通常用于排队理论和统计学中的建模，还可作为信号模型用于熵编码技术(例如算术编码)。马尔可夫链还有众多生物学方面的应用，模拟生物人口过程是最为常见的。隐蔽马尔可夫模型还被用于生物信息学，用以编码区域或进行基因预测。此外，马尔可夫链因其特有的马尔可夫性质，被广泛应用在病毒传播的预测研究上，例如模拟植物病毒传播[135]和解析流行病的爆发[136-137]。

P(Xu=xu|{Xv},v≠u}=P(Xu=xu|Xv,v∈Vu)

(147)

其中Vu表示节点u的邻居节点所构成的集合。在马尔可夫随机网络中，只有节点对中两节点间存在关联性，因而可以用度分布P(k)和条件概率分布P(k′|k)这两个参量来刻画其上的动力学。这类网络的马尔可夫特征意味着，所有高阶相关性都可以表示为P(k)和P(k′|k)的函数。

7.2基于马尔可夫链的SIS传播动力学建模

考虑无向连通网络G=(V,E)，感染概率为β，感染节点的恢复概率为u。定义pi，t是节点i在t时刻处于感染态的概率,ξi，t是节点i在t时刻没有被其邻居感染的概率(在SIS模型中，将所有节点的状态看成是一个状态组合，由于每个节点状态可能为易感态或感染态，所以N个节点的状态组合共有2N种。换言之，可以将传播过程看成是包含2N种不同状态的马尔可夫链：节点在t时刻的状态只取决于t-1时刻的状态)，

(148)

其中Vi表示节点i的邻居节点集合，节点i在t时刻没有收到来自邻居节点的感染包括两种可能：1)邻居节点j处于感染态，但是在t时刻没能感染节点i，即式中pj，t-1(1-β)；2)邻居j在t-1时刻处于易感态，即式中(1-pj，t-1)。

进一步，若满足以下3种情况之一，可推知节点i在t时刻处于健康态：

1)i在t时刻之前是健康态，并且在t时刻没有被其邻居感染；

2)i在t时刻之前已被感染，并且在t时刻恢复同时没有再次被其邻居感染；

3)i在t时刻之前已被感染，在t时刻接触到了感染态邻居，在这些无效的接触(因发生接触时处于感染态,因而不会被其邻居感染,称这些接触是无效的)之后恢复为健康态。

因此节点i在t时刻处于健康状态的概率为

1-pi，t=(1-pi，t-1)ξi，t+μpi，t-1ξi，t+μpi，t-1(1-ξi，t)/2，(i=1，2，…，N)

(149)

本节将主要关注基于上述马尔可夫链模型的传播动力学(我们前面各个章节所涉及的传播动力学模型均可用马尔可夫链来进行描述)解析方法：特征值分析法和微观尺度上的离散时间马尔可夫方法。两种方法都基于节中的数学模型，其差别在于后续的解析过程中：前者主要通过非线性系统解的稳定性进行相应分析；而后者则是一种更为数值化的非线性系统迭代解析过程。

7.2.1特征值分析

7.2.1.1特征值分析方法

定义为网络的邻接矩阵，AT是矩阵A的转置，At是矩阵A的t次幂，λi，A表示矩阵A的第i大特征值。下面不妨以一种更为数学化、更为严谨的形式来给出对于相应结论的证明。

定理: 若一种流行病消亡，则有β/μ<τ=1/λi，A，其中τ表示流行病传播阈值。

证明：由公式(149)可知，

(150)

基于近似公式，当时a≪1，b≪1时，省略二阶项(1-a)(1-b)≈1-a-b,上述证明中进行了两处近似：

由公式(150)可得到

(151)

将式(151)转换为矩阵表示法

Pt=((1-μ)I+βA)Pt-1=SPt-1=StP0

(152)

其中t时刻的各节点状态空间Pt=(p1，t，…，pN，t)T，S=(1-μ)I+βA称为系统矩阵。

引理：矩阵S的第i大特征值λi，S=1-μ+βλi，A，而且S的特征项量与A的完全相同。

引理证明：vi,A为矩阵A的特征值λi，A对应的特征向量，由定义可知Avi,A=λi，Avi,A进而有

Svi,A=(1-μ)vi,A+βAvi,A=(1-μ)vi,A+β λi，Avi,A=(1-μ+β λi，A)vi,A

7.2.1.2稳定性分析

在均匀网络中，节点的度都比较集中在平均度附近(即k≈〈k〉，∀k)，ρ(t)表示感染个体的平均密度。其平均场速率方程应用平均场理论可写作

dρ(t)/dt=-ρ(t)+λ〈k〉ρ(t)(1-ρ(t))

(153)

此近似方程在ρ(t)→0时(即处于传播的初始阶段)成立，原因已在平均场部分进行了解释。当t→∞时，由∂ρ(t)=0即可解得传播阈值λc=〈k〉-1[14]。当时λ<λc时，ρ=0；若λ≥λc，则有ρ～(λ-λc)。

对于异质网络，节点度值的波动性很大，无法再使用全局的平均近似；此时定义度为k的感染个体密度为ρk(t)，相应的反应速率方程为[14-15]

dρk(t)/dt=-ρk(t)+λk(1-ρk(t))Θk(t)dt

(154)

然而对于节点之间有关联性的网络，上面的推导是不正确的，因为在Θk(t)中没有考虑边的连出节点的度k。对于马尔可夫网络，网络的关联性是完全由条件概率P(k′|k)来定义的，其中

(155)

kP(k′|k)P(k)=k′P(k|k′)P(k′)≡〈k〉/P(k,k′)

(156)

很明显ρk=0是一个平凡解；在阈值附近，只有少许节点被感染(即ρk→0)，省略掉ρkρk′这类高阶项，我们可以将式(149)线性化为

(157)

如果雅可比矩阵L至少有一个正特征值，那么ρk=0这个解就不稳定[7]。依据相应的定义，连接矩阵C中元素Ckk′=kP(k′|k)，由式(156)可得，若vk是C某一特征值对应的特征向量，那么P(k)vk是C的转置矩阵CT里对应相同特征值λi，C的特征向量，由此可知C的所有特征值都为实数。假设λ1，A是C的最大特征值，那么只要雅可比矩阵L的最大特征值-1+λλ1，C>0，即λ>1/λ1,C时，ρk=0这个解就不稳定，就会有另一个非零解成为实际的稳态，即网络处于染病态，由此可知传播阈值λc=1/λ1，C。

可以发现，无论使用网络的何种矩阵表示方式，都得到了定性相同且定量相同的结论。然而真实网络往往都存在一定程度的度关联性(例如同配性或异配性导致的度关联)，而且通常很难用一个闭合形式的条件概率P(k′|k)来刻画。因此Wang等[123]提出了一个普适的病毒传播分析模型，能够刻画底层网络拓扑对传播的影响，并且适合于任何网络结构。

对于均匀网络或随机网络，其邻接矩阵的最大特征值是平均度，因此传播阈值为λc=1/〈k〉，这与平均场得出的结果一致；对于服从幂律分布且大小有限的网络，文献[7]通过对实际网络和模型网络中的阈值模拟来对比邻接矩阵特征值方法和平均场方法确定阈值的准确性，模拟结果明显表明λc=1/λ1，A更准确。Chung等[124]更从理论上计算出了服从幂律度分布且大小有限的网络的邻接矩阵最大特征值

(158)

(159)

具体计算中需要注意以下3点：1)如果初始向量x(0)刚好与最大特征值对应的特征向量正交，那么此方法无效。由佩龙—弗洛比尼斯定理[7]可知，一个非负实矩阵的最大特征值非负，且最大特征值对应的特征向量的元素都为非负，因此若初始向量x(0)的元素都为正数，则x(0)不会和最大特征值对应的特征向量正交，所以选择初始向量时应该使其元素都是正数；2)向量的元素每次迭代之后都会增长，因此最终这些元素在计算机中会因为过大而溢出，因此需要周期性地对向量进行归一化；3)判断x(t)是否收敛为最大特征值对应的特征向量时可以针对不同的初始向量计算两次，由此在允许的误差范围内判断其何时收敛。

图14　流行病爆发阈值随传播率的演化情况[14-15]

模拟结果均是在俄勒冈州自治系统网络上进行的，其度分布服从幂律N=10 900,

在得到最大特征值对应的特征向量后，再次乘以邻接矩阵，由矩阵的特征值和特征向量的关系可知，两个向量的任意元素的比值即为最大特征值，通常取几个较大的元素的比值的平均值，这样可以减小误差、提高准确性。

7.2.1.3特征值方法的准确性及适用性分析

通过观察对于阈值和感染密度的具体理论预测与数值模拟情况(详见图14与图15)，通过求网络邻接矩阵的最大特征值来确定传播阈值的方法，其结果比平均场的解析更为准确，并且此方法既适用于实际网络也适用于模型网络。

通过以上分析，网络邻接矩阵的最大特征值的大小决定了阈值，此外最大特征值对于诸如信息传播[127]、同步[128-129]、渗流[130-132]等其他动力学解析的影响也非常大。通过求解网络邻接矩阵的最大特征值来解析传播阈值虽然普遍适用于各种网络结构，但当网络较大时很难直接快速求出其最大特征值，因此只能用近似的方法，例如上面提到的幂乘法，这会使结果的准确性有所降低。

7.2.2离散时间的非线性系统数值迭代解析法

这一方法同样基于7.2中的数学模型，但与7.3.1中分析解的稳定性完全不同。

对于二状态SIS传播模型，目前的理论研究多集中于全接触又译作反应过程,每个时间步,感染节点会接触其全部邻居节点)和单点接触(每个时间步,感染节点只随机接触它的一个邻居节点进行试图的感染,若邻居已处于感染态则跳过感染过程)两种特殊的接触过程；但实际上，节点之间的接触过程则与具体问题本身密切相关，具有不确定性。文献[134]采用基于接触过程的概率离散马尔可夫链模型，并将其分别应用到含权和无权网络，研究不同的接触过程对流行病动态传播的影响。

pi(t+1)=(1-qi(t))(1-pi(t))+(1-μ)pi(t)+μ(1-qi(t))pi(t)

(160)

对于由N个式(151)组成的非线性系统，应用不动点迭代(或牛顿迭代)就可以得到对应的数值解(若为便捷,可使用Matlab中的fsolve函数进行相应的求解)，其与蒙特卡洛模拟的对比结果如图16和图17所示。图16和图17表明，不论是在微观尺度(插图中的具体节点概率)上，还是在不同接触情形下，基于离散马尔可夫链方法得到的数值解都与模拟结果吻合得很好，充分验证了该方法的准确性。

其中N=10 000，μ=1，λ→∞，

图17　接触个数取不同值时,航空网络上

7.2.2.1上述方法相应的拓展变形

为了能够更准确地描述流行病传播动力学以及自适应网络的相应情况，文献[29]将退火邻接矩阵(Annealed Adjacency Matrix, AAM,该矩阵中的每个元素Aij都是0-1之间的实数,表示节点i-j相连的概率)引入到离散马尔可夫链方法[126]中,得到节点i在t+1刻为健康态的退火方程：

si(t+1)=si(t)+μ(1-si(t))-

si(t)(1-qi(t))

(161)

将邻接矩阵换成是与传播过程一起随时间演化的退火矩阵，文献[29]进一步考虑了自适应网络上的马尔可夫链SIS传播模型，如图18所示，结果再一次验证了退火近似的准确性，而且对于度较小的节点预测性也很高。

如果应用平均场假设，认为度相同的节点具有相同的感染概率，那么∀i，当ki=k时，有si=ρk。从而可对方程(161)进行一个简单代换，得出异质平均场方程的一个变形(姑且用NewHMF与先前的方法进行一个区分)：

(162)

其中Nk=NP(k)。

通过对比发现，相较于一般的异质平均场方程(HMF)，退火方程(Annealed,由式(161)迭代解析可得)和上述异质平均场方法(New HMF,由式(162)迭代解析可得)的数值解与蒙特卡罗模拟结果(MC)有着很好的吻合(见图19)。

一般来讲，P(k，k′)=kP(k′|k)/(NP(k′))。从而，上述异质平均场方程可以写成如下形式：

(163)

当对于任意的k，当λ(1-ρk)≪1时，将式(154)线性化后可得出与经典HMF方法相同的解析式：

(164)

换句话说，HMF方程只在流行病爆发阈值附近才是精确的，在λ较大时它所使用的近似都会使得误差产生。

相应地，如果式(162)中使用的邻接矩阵，对应得到的方法被称作非微扰异质平均场方法(Non-Perturbative HMF)[38]。

7.2.3小结

回顾前面介绍的7种方法，可以发现大多数方法是在从网络层面、在较宏观的尺度分析问题，或者也只是精确到度层阶(以度对网络中的节点进行分类)的介观尺度，只有点对近似和主方程精确到了点对的状态这一较为微观的视角；那么有没有方法可以从微观的角度去研究宏观尺度上的一些性质呢？马尔可夫链方法带给我们一个新的视角。它与常见的微分方程(以平均场方法为代表)有较大差异。从群体角度来看，经典的平均场方法的确能够较好地描述、预测传播动力学的临界现象，但往往不能给出指定个体的感染信息，然而在许多实际应用和理论分析中, 特定个体的信息往往是人们所需要和关注的。马尔可夫方法恰好能解决这一问题。

N=104，μ=1

N=5 000，〈k〉=4

7.3.2小节中介绍的数值迭代方法准确性非常之高，其优点在于几乎没有进行近似，保留了网络结构的几乎所有信息，这也正是它明显优于其它方法的特点。然而，也正是由于这样的特性，使得计算的时间复杂度较高，比如马尔可夫链方法的状态空间随着节点个数呈指数增长,很难用于处理实际的大型网络。这或许也是它最大的问题所在。另外，它的另一显著优点是可以处理离散时间，无须连续时间假设；这也克服了之前许多方法的一些局限，使得结果更为准确。而基于退火矩阵的马尔可夫链方法具有更好的拓展性，可以对基于节点状态的动态演化网络实行动力学分析，且对于大多数节点预测十分准确；虽然所做的改动并不大，只是将邻接矩阵变为了退火矩阵，但其收效却是颇有些“意想不到”的：基于退火邻接矩阵(矩阵的每个元素值代表节点对之间的连接概率)的马尔可夫链方法还能够很好地用于动态网络中的传播动力学分析，如预测自适应网络中的感染密度分析。

8总结与讨论

本文前述各节分析了如何具体求解具体网络在特定动力学规则下的爆发阈值，细心的读者或许已经发现SIS和SIR的临界感染分布是否存在本质区别是一个聚讼已久的话题。接下来，我们将尽量清晰地梳理一下这一问题的发展脉络并给出确定阈值的方法。

8.1SIS与SIR模型阈值求解的演进

(165)

(166)

很接近，而QMF阈值只对γ<5/2有效。局域化分析的相关工作[143-144]则证明度分布指数γ>5/2时SIS真实阈值有限，与异质平均场结论近似。文献[145]指出任何度分布衰减小于指数的小世界网络上都不存在流行病。

较早的SIR阈值结论等同于SIS[1]，即无关联的均匀网络上的SIR传播阈值为λc=1/〈k〉。Moreno等[147]考虑节点度的异质性，推得复杂网络上SIR传播阈值与先前SIS相同λc=〈k〉/〈k2〉。如此得出，二阶矩〈k2〉对传播阈值起着决定性作用。在无标度网络结构中，当N趋于无穷大时二阶矩发散，传播阈值为0。Boguna等[47]对此进行了进一步补充，发现考虑节点度的异质性时，在无关联网络上的SIR阈值为

(167)

关联网络上的SIR阈值为最大特征值的倒数

(168)

上述研究都是针对全接触模型，有研究者单独研究了单点接触情况下的阈值行为。Zanette等[151]考虑基于单点接触的SIR，发现当小世界网络的随机性p比较大时，SIR感染范围分布形式为幂律分布尾部截断。幂律指数约为-1.5，分布的尾部呈现二次凸起，且直观认为最大凸起对应的感染范围为Rmax～0.25N，随着网络规模N线性增长。Ben-Naim等[152-153]基于均匀混合假设考虑SIR阈值点的爆发规模。得出阈值点的感染范围服从幂律分布P(n)～n-3/2。对于规模为N的网络，SIR的最大感染范围n～N2/3，平均感染范围〈n〉～N1/3，流行病最长生存时间t*～N1/3，平均生存时间〈t〉～lnN。基于随机游走的思想，文献[154]再次验证了上述结论。Khaleque[155]考虑空间网络(节点距离异质分布)上的SIR传播，模拟得出在阈值之下感染范围呈指数衰减；阈值点为幂律分布，阈值之上呈双峰分布。Petter Holme[156]发现，SIR的流行病最长生存时间对应的感染概率大于且独立于传播阈值。

8.2对于SIS和SIR模型传播阈值的模拟判定方法

鉴于理论传播阈值的重要性，对理论模拟阈值进行准确检测便显得尤为重要。然而如何确定一个传播模型的模拟阈值，这一问题研究者们一直争论不休，提出了许多不同的指标。我们同样对其进行一个梳理。

鉴于易感性对估计SIS传播阈值的有效性，我们尝试将其进一步应用于SIR传播模型。基于随机规则网络，我们对SIS和SIR传播过程进行了大量模拟，并根据易感性

χ=N(〈ρ2〉-〈ρ〉2)/〈ρ〉[141]

(169)

的峰值估计传播阈值。基于易感性能够准确估计SIS的传播阈值[12]，但对估计SIR传播阈值失效。

在磁化系统中，通常采用序参量的比值来确定平衡相变的临界点，其一般定义[159]为

(170)

显然，当q=s=1时即对应可变性Δ：

(171)

另外我们的研究小组发现，实际网络的理论阈值与模拟阈值存在明显区别。绝大多数网络中，由可变性和易感性两种方法确定的SIS模拟阈值差不多(具有明显负关联性的网络除外)，但二者确定的SIR模拟阈值存在很大差别。特别是，由易感性确定的SIR阈值总是偏离理论阈值很多，由可变性确定的模拟阈值则更偏向于异质平均场(整体看来，实际网络的特征值阈值要小于异质平均场阈值)。当特征值阈值与模拟阈值差很多时，对应的异质平均场和特征值阈值也差很多。通过计算IPR参数[143]可以发现，模拟阈值和特征值阈值较吻合时IPR值较小，二者不一致时IPR值较大[161]。也就是说，特征值阈值和模拟阈值不一致应该是特征向量局域化所致。另外，当网络呈负关联(如ego-facebook,AS2000)时，平均场和特征值阈值都与模拟阈值不一致，这与无标度网络上的结果分析基本一致。

8.3讨论

网络动力学作为复杂网络研究中的重要课题，具有极其重要的意义。我们在文章中介绍的方法大多都可以应用于其它动力学过程，如投票者模型、网络同步等等，而不只局限于流行病的传播。

在这篇综述中，纵观我们介绍的这七种方法，其所应用的数理方法大致可分为3类：微分方程、生成函数、基于马尔可夫链方法。这些看似不同的方法，却都能比较准确地描述现实网络中动力学的行为特征；然而也应当看到不同的方法所对应的解析难度以及分析角度着实有较大不同。

基于微分方程的3种方法(平均场近似、点对近似和主方程)里，主方程方法考虑的动力学过程、网络拓扑等情况最为全面，也最为精确；然而它的计算机复杂性也是最高的，其维数很高，无法像平均场方法那样在人力可解的情况下给出一个还算不错的预测。这3种方法都可以描述动力学的时间演化过程，相比之下，经典的渗流理论和矩阵分析在这一方面略逊一筹。然而这3种方法它们都假设时间是连续的，而马尔可夫方法则很好地克服了这一问题，而且通过它的框架也可以衍生出时间离散的平均场方法。

另外，可以发现影响传播过程的网络拓扑因素有很多：例如度相关性、动力学相关性、局域效应(如聚类和模体等)、有限尺度效应、网络异质性、传播率的异质性、感染的非对称性、拓扑随时间的演化、邻居节点度转化等等；另外网络上个体的自适应行为、人类动力学与网络拓扑的共演化[162]、网络多层耦合性与连边的多关系特性也是极其重要、近年研究较热的问题。那么是否能够有一种考虑所有重要因素的综合方法是一个有趣的问题。想要达成这样一个远大的目标，有很多问题依然摆在我们面前：在传播动力学中如此繁复的影响因素里，它们各自对动力学产生的影响究竟有多大？哪些对于传播动力学来说是至关重要不可或缺、哪些又是相对无足轻重的？甚至更深一步，基于网络的流行病传播动力学解析方法其根本局限在哪里？这些都是我们未来最终须要回答的问题。

我们运用和拓展本文介绍的7种常用理论方法来描述网络传播动力学，取得了一系列研究成果。对于异质平均场，我们将其拓展后用于描述信息-疾病非对称耦合传播动力学，发现信息扩散有利于抑制流行病传播[165-166]。杨慧等[167]拓展点对近似方法来分析个体异质性对自适应传播所带来的影响。Xu等[164]利用主方程方法准确地刻画了多关系社会网络中的信息传播。王伟等[168-170]将边划分理论拓展用于描述非马尔可夫社会传播动力学，发现网络结构和动力学参量都对相变类型有很大的影响。舒盼盼等[160]利用可变性方法对常用理论方法得到的阈值有效性进行了分析。

我们去评判一种方法的优劣与否，不只要看它在某些特定情景下的解析是否准确，也须要考虑这一方法是否具有较好的可用性与扩展性，能够加入对于更多重要因素的考虑。一个理论方法所基于的假设越少，其适用性往往越广；同时，还要看到，理论解析不同于数值模拟，解析一定是建立在一定的假设上的，它对于数据的需求性不应当过强，最好只需要基于较少的初始数据(如度分布、度关联性)，其更多的应当强调理论解析性，应较少的加入过多的数值模拟；有时从某种意义上讲，数值求解未必会比蒙特卡罗模拟更具优势。可解释性以及易懂性，是理论解析方法所应具有的；如果一种解析方法自身就很复杂、很难为我们提供一些直观的定性理解，那么即使它十分准确也未必能称作好的方法。理论解析事关对于机制的理解，哪些假设和近似是合理的关系到是否能够描述出动力学过程中的关键因素，这是非常本质而且极其重要的方面。

参考文献:

[1]AndersonRM,MayRM.InfectiousDiseasesofHumans:dynamicsandcontrol[M].Oxford:OxfordUniversityPress, 1992.

[2]DaileyDJ,GaniJ.EpidemicModeling:AnIntroduction[M].Cambridge:CambridgeUniversityPress, 2001.

[3]BernoulliD.Essaid'unenouvelleanalysedelamortalitecauseeparlapetiteveroleetdesavantagesdel'inoculationpourlaprevenir[J].MemMathPhysAcadRoySci, 1760, 1-45.

[4]KermackWO,McKendrickAG.Acontributiontothemathematicaltheoryofepidemics[J].ProcRSocLondA, 1927, 115, 772: 700-721.

[5]AlbertR,BarabsiAL.Statisticalmechanicsofcomplexnetworks[J].RevModPhys, 2002, 74:47.

[6]BoccalettiS,LatoraV,MorenoY,etal.Complexnetworks:structureanddynamics[J].PhysRep, 2006,424: 175.

[7]NewmanMEJ.Networks:AnIntroduction[M].Oxford:OxfordUniversityPress, 2010.

[8]ButtsCT.Revisitingthefoundationsofnetworkanalysis[J].Science, 2009, 325: 414.

[9]JacksonM.SocialandEconomicNetworks[M].Princeton:PrincetonUniversityPress, 2010.

[10] Vespignani A. Predicting the behavior of techno-Social systems[J]. Science, 2009, 325 (5939): 425.

[11] Vespignani A. Modeling dynamical processes in complex socio-technical systems[J]. Nature Physics, 2012, 8: 32-39.

[12] Dorogovtsev S N,Goltsev A V. Critical phenomena in complex networks[J]. Rev Mod Phys, 2008, 80: 1275.

[13] Barrat A, Barthelmy M, Vespignani A.Dynamical Processes on Complex Networks[M]. New York: Cambridge University Press, 2008.

[14] Pastor-Satorras R, VespignaniA. Epidemic spreading in scale-free networks[J]. Phys Rev Lett, 2001, 86: 32001.

[15] Pastor-Satorras R, Vespignani A. Epidemic dynamics and endemic states in complex networks[J]. Phys Rev E, 2001, 63(6): 066117.

[16] Hufnagel L, Brockmann D, Geisel T. Forecast and control of epidemics in aglobalized world[J]. Proc Natl Acad Sci, 2004, 101: 15124.

[17] Colizza V, Barrat A, Barthélemy M, et al. The role of the airline transportation network in the prediction and predictability of global epidemics[J]. Proc Natl Acad Sci, 2006, 103: 2015.

[18] Balcan D, Hu H, Goncalves B,et al. Seasonal transmission potential and activity peaks of the new influenza A(H1N1): a Monte Carlo likelihood analysis based on human mobility[J]. BMC Medicine, 2009, 7: 45.

[19] Brockmann D, Helbing D.The hidden geometry of complex, network-driven contagion phenomena[J]. Science, 2013, 342 (6164): 1337-1342.

[20] Henkel M, Hinrichsen H, Lubeck S. Non-equilibrium phase transition: Absorbing Phase Transitions[M]. Netherlands: Springer Verlag, 2008.

[21] Danon L,Ford A P,House T,et al. Networks and the epidemiology of infectious disease[J]. Interdisciplinary Perspectives on Infectious Diseases, 2011:284909.

[22] Keeling M, Rohani P. Modeling Infectious Diseases in Humans and Animals[M]. Princeton: Princeton University Press, 2010.

[23] Moreno Y, Pastor-Satorras R, Vespignani A. Epidemic outbreaks in complex heterogeneous networks[J]. Eur Phys J B, 2002, 26: 521-529.

[24] Boguna M, Pastor-Satorras R, Vespignani A.Absence of epidemic threshold in scale-free networks with degree correlation[J]. Phys Rev Lett, 2003, 90: 2.

[25] Castellano C, Pastor-Satorras R. Thresholds for epidemic spreading in networks[J]. Phys Rev Lett, 2010, 105(21): 218701.

[26] Eguiluz V M, Klemm K. Epidemic threshold in structured scale-free networks[J]. Phys Rev Lett, 2002, 89: 10.

[27] Parshani R, Carmi S, Havlin S. Epidemic threshold for the susceptible-infectious-susceptible model on random networks[J]. Phys Rev Lett, 2010, 104(25): 258701.

[28] Gleeson J P, Melnik S, Ward J A,et al. Accuracy of mean-field theory for dynamics on real-world networks[J]. Phys Rev E, 2012, 85(2): 026106.

[29] Guerra B, Gómez-Gardees J. Annealed and mean-field formulations of disease dynamics on static and adaptive networks[J]. Phys Rev E, 2010, 82(3): 035101.

[30] Sood V,Redner S. Voter model on heterogeneous graphs[J]. Phys Rev Lett, 2005, 94(17): 178701.

[31] Pikovsky A S, Rosenblum M G, Kurths J. Synchronization in a population of globally coupled chaotic oscillators[J]. EPL, 1996, 34: 165.

[32] 于渌, 郝柏林. 相变和临界现象[M]. 北京: 科学出版社, 1984.

[33] 林宗涵. 热力学与统计物理学[M]. 北京: 北京大学出版社, 2007.

[34] Landau L D, Lifshitz. Statistical Physics[M]. 3rd Edition London: Pergamon Press, 1986.

[35] Stanley H E. Introduction to Phase transition and Critical Phenomena[M]. New York: Oxford Univ Press, 1971.

[36] 苏汝铿. 统计物理学[M].第2版,上海:高等教育出版社, 2004.

[37] Bailey N T J. The Mathematical Theory of Infectious Diseases and Its Applications[M]. New York: Hafner Press, 1975.

[38] Gomez S, Gómez-Gardenes J, Moreno Y,et al. Nonperturbative heterogeneous mean-field approach to epidemic spreading in complex networks[J]. Phys Rev E, 2011, 84(3): 036105.

[39] Durrett R. Some features of the spread of epidemics and information on a random graph[J]. Proc Natl Acad Sci, 2010, 107: 4491.

[40] Baronchelli A, Loreto V. Ring structures and mean first passage time in networks[J]. Phys Rev E, 2006, 73(2): 026103.

[41] Baronchelli A, Pastor-Satorras R. Mean-field diffusive dynamics on weighted networks[J]. Phys Rev E, 2010, 82(1): 011111.

[42] Gleeson J P. High-accuracy approximation of binary-state dynamics on networks[J]. Phys Rev Lett, 2011, 107(6): 068701.

[43] Milo R, Shen-Orr S,Itzkovitz S,et al. Network motifs: simple building blocks of complex networks[J]. Science, 2002, 298: 824.

[44] Zanette D H, Negro R, Argentina. Dynamics of rumor propagation on small-world networks[J]. Phys Rev E, 2001, 65(4): 041908.

[45] Baronchelli A, Castellano C, Pastor-Satorras R. Voter models on weighted networks[J]. Phys Rev E, 2011, 83(6): 066117.

[46] Liu Z, Wang X, Wang M. Inhomogeneity of epidemic spreading[J]. Chaos, 2010, 20(2): 023128.

[48] Yang Z, Zhou T. Epidemic spreading in weighted networks: an edge-based mean-field solution[J]. Phys Rev E, 2012, 85(5), 056106.

[49] Watts D J, Strogatz S H. Collective dynamics of ‘small-world’ networks[J]. Nature, 1998, 393: 4.

[50] Barabási A L, Albert R. Emergence of scaling in random networks[J]. Science, 1999, 286: 509-512.

[51] Newman M E J. Power laws, Pareto distributions and Zipf's law[J]. Contemporary Physics, 2005, 46: 5.

[52] Castellano C, Pastor-Satorras R. Competing activation mechanisms in epidemics on networks[J]. Sci Rep, 2012, 2: 371.

[53] Yang R, Wang B H, Ren J,et al. Epidemic spreading on heterogeneous networks with identical infectivity[J]. Phys Lett A, 2007, 364: 189-193.

[54] Alexandrowicz Z. Critically branched chains and percolation clusters[J]. Phys Lett A, 1980, 80: 4.

[55] Kenah E, Robins J M. Second look at the epidemic spread[J]. Phys Rev E, 2007, 76(3): 036113.

[56] Karp R,Schindelhauer C,Shenker S,et al. Randomized rumor spreading[Z]. Foundations of Computer Science. Proceedings of 41st Annual Symposium, 2000:565-574.

[57] 许田, 张培培, 姜玉梅, 等. 流行病传播模型与SARS[J]. 自然杂志, 2004, 26(1): 20-25.

Xu Tian, Zhang Peipei, Jiang Yumei, et al. Epidemic spreading models and SARS[J]. Ziran Zazhi, 2004, 26(1): 20-25.

[58] Saumell-Mendiola A, Serrano M A, Boguna M. Epidemic spreading on interconnected networks[J]. Phys Rev E, 2012, 86(2): 026106.

[59] Rand D A. Correlation equations and pair approximations for spatial ecologies[J]. Advanced Ecological Theory: Principles and Applications, 2009, 12 (34):329-368.

[60] Keeling M J, Rand D A. Spatial correlations and local fluctuations in host parasite systems[J]. From Finite to Infinite Dimensional Dynamical Systems, 2001: 5-57.

[61] Keeling M J.The effects of local spatial structure on epidemiological invasions[J].Proc R SocLond B, 1999, 266 (1421): 859-867.

[62] Morris A J. Representing spatial interactions in simpleecological models[D]. UK: University of Warwick, 1997.

[63] Benoit J, Nunes A, Teloda G M. Pair approximation models for disease spread[J]. Eur Phys J B, 2006, 50 (1/2): 177-181.

[64] Eames K T D, Keeling M J. Modeling dynamic and network heterogeneities in the spread of sexually transmitted diseases[J]. Proc Natl Acad Sci, 2002, 99: 13330-13335.

[65] Pugliese E,Castellano C. Heterogeneous pair approximation for voter models on networks[J].EPL, 2009, 88 (5): 58004.

[66] Shaw L B,SchwartzI B. Fluctuating epidemics on adaptive networks[J]. Phys Rev E, 2008, 77 (6): 066101.

[67] Gross T, D’Lima C J D, Blasius B. Epidemic dynamics on an adaptive network[J].Phys Rev Lett, 2006, 96 (20):208701 .

[68] Shaw L B, Schwartz I B. Enhanced vaccine control of epidemics in adaptive networks[J].Phys Rev E,2010, 81 (4): 046120.

[69] Keeling M J, Rand D A, Morris A J. Correlation models for childhood epidemics[J].Proc R Soc B, 1997, 264 (1385): 1149-1156.

[70] CaldarelliG, Pastor-Satorras R, Vespignani A. Structure of cycles and local ordering in complex networks[J].Eur Phys J B,2004, 38 (2): 183-186.

[71] GleissP M, Stadler P F, Wagner A, et al.Relevant cycles in chemical reaction networks[J]. Adv Complex Syst, 2001, 4 (5): 207-226.

[72] Sato K, Matsuda H, Sasaki A. Pathogen invasion and host extinction in lattice structured populations[J]. J Math Biol, 1994, 32 (3): 251-268.

[73] Keeling M J. The ecology and evolution of spatial host-parasite systems[D]. UK: University of Warwick, 1995.

[74] Nakamaru M, Matsuda H, IwasaY. The evolution of cooperation in a lattice-structured population[J].J Theor Biol, 1997, 184 (1): 65-81.

[75] Szabó G,Hauert C. Evolutionary prisoner’s dilemma games with voluntary participation[J].Phys Rev E,2002, 66 (6): 062903.

[76] Perc M,Marhl M. Evolutionary and dynamical coherence resonances in the pair approximated prisoner'sdilemma game[J]. New J Phys, 2006, 8: 142.

[77] Marceau V, No⊇l P A, Hébert-Dufresne L,et al. Adaptive networks: coevolution of disease and topology [J]. Phys Rev E, 2010, 82(3): 036116.

[79] Durrett R, Gleeson J P, Lloyd A L,et al. Graph fission in an evolving voter model[J]. Proc Natl Acad Sci, 2012, 109: 3682-3687.

[80] Mollison D. Spatial contact models for ecological and epidemic spread[J]. J Roy Stat Soc B, 1977, 39: 283-326.

[81] Grassberger P. On the critical behavior of the general epidemic process and dynamical percolation[J]. Math Biosci, 1982, 63: 157-172.

[82] Newman M E J. Spread of epidemic disease on networks[J]. Phys Rev E, 2002, 66(1): 016128.

[83] Newman M E J. Threshold effects for two pathogens spreading on a network[J]. Phys Rev Lett, 2005： 95(10):108701.

[84] Meyers L A, Newman M E J, Martin M,et al. Applying network theory to epidemics: Control measures for mycoplasma pneumoniae outbreaks[J]. EID, 2003, 92: 204.

[85] Madar N, Kalisky T, Cohen R, et al. Immunization and epidemic dynamics in complex networks[J]. Eur Phys J B, 2004, 38: 269.

[86] No⊇l P A, Davoudi B, Brunham R C, et al. Time evolution of epidemic disease on finite and infinite networks[J]. Phys Rev E,2009, 79(2): 026101.

[87] Marder M. Dynamics of epidemics on random networks[J]. Phys Rev E, 2007, 75: 066103.

[88] Miller J C. Epidemic size and probability in populations with heterogeneous infectivity and susceptibility[J]. Phys Rev E, 2007, 76(1): 010101.

[90] Newman M E J. Random graphs with clustering[J]. Phys Rev Lett, 2009, 103(5): 058701.

[91] Newman M E J. Competing epidemics on complex networks[J]. Phys Rev E, 2011, 84(3): 036106.

[92] Allard A, No⊇l P A, DubéL J. Heterogeneous bond percolation on multitype networks with an application to epidemic dynamics[J]. Phys Rev E, 2009, 79(3): 036113.

[93] Funk S, Jansen V A. Interacting epidemics on overlay networks[J]. Phys Rev E, 2010, 81(3): 036118.

[94] Dickison M, Havlin S, Stanley H E. Epidemics on interconnected networks[J]. Phys Rev E, 2012, 85(6), 066109.

[95] Wang Y, Xiao G. Epidemics spreading in interconnected complex networks[J]. Phys Lett A, 2012, 376: 42-43.

[96] Pearl J. Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference[M]. San Francisco:Morgan Kaufmann,1988.

[97] Mézard M, Parisi M, Virasoro M. Spin Glass Theory and Beyond World Scientific[M]. Singapore:Pergamon Press, 1987.

[98] Mézard M, Parisi G. The Bethe lattice spin glass revisited[J], Eur Phys J B,2001, 20(2):217-233.

[99] Karrer B, Newman M E J. Message passing approach for general epidemic models[J]. Phys Rev E, 2010,82(1):016101.

[100]Bianconi G, Gulbahce N. Algorithm for counting large directed loops[J]. J Phys A Math Theor, 2008,41(22):224008.

[101]Sulc P, Zdeborov L. Belief propagation for graph partitioning[J]. J Phys A MathTheor, 2010,43(28):285003.

[102]Shiraki Y, Kabashima Y. Cavity analysis on the robustness of random networks against targeted attacks: influences of degree-degree correlations[J]. Phys Rev E, 2010, 82(3): 036101.

[103]Reichardt J, Alamino J, Saad D. The interplay of microscopic and mesoscopic structure in complex networks[DB/OL].[2014-06-30]. arXiv.org/abs/1012.4524.

[104]Yoon S, Goltsev A V, Dorogovtsev S N, et al. Belief-propagation algorithm and the ising model on networks with arbitrary distributions of motifs[J]. Phys Rev E,2011,84(4):041144.

[105]Pastor-Satorras R, Vespignani A. Epidemic dynamics in ?nite size scale-free networks[J]. Phys Rev E,2002,65(3): 035108.

[106]Shao J, Buldyrev S V, Braunstein L A, et al. Structure of shells in complex networks[J]. Phys Rev E, 2009, 80(3):036105.

[107]张昊,陈超,王长春. 基于空穴理论的复杂网络流行病传播控制[J]. 电子科技大学学报,2011, 40:4.

Zhang Hao, Chen Chao, Wang Changchun. Epidemic spread and control on complex networks using cavity theory[J]. JUESTC, 2011, 40:4.

[108]Altarelli F, Braunstein A, Dall'Asta L, et al. The spread optimization problem[DB/OL].[2014-06-30]. arXiv/org/abs/1203.1426v1.

[109]Volz E. SIR dynamics in random networks with heterogeneous connectivity[J]. Math Biol, 2008, 56:293-310.

[110]Miller J C. A note on the derivation of epidemic final sizes[J]. Bull Math Biol, 2012, 74:2125-2141.

[111]Miller J C. A note on a paper by Erik Volz: SIR dynamics in random networks[J]. Math Biol, 2011, 62:349-358.

[112]Miller J C, Slim A C, Volz M E. Edge-based compartmental modeling for infectious disease spread part I: an overview[J]. Journal of the Royal Society Interface,2011,9(20):890-906.

[113]Miller J C, Volz M E. Edge-based compartmental modeling for epidemic spread Part II: model selection and hierarchies[J]. Journal of Mathematical Biology,2013,67(4):869-899.

[114]Miller J C, Volz M E. Edge-based compartmental modeling for infectious disease spread Part III: disease and population structure[J]. Eprint Arxiv,2011,9(70):890-906.

[115]Volz E, Frost S D W, Rothenberg R, et al. Epidemiological bridging by injection druguse drives an early HIV epidemic[J]. Epidemics,2010,2(3):155-164.

[116]Volz E, Meyers L A. Epidemic thresholds in dynamic contact networks[J]. J R Soc Interface, 2009, 6: 233-241.

[117]Volz E M, Miller J C, Galvani A, et al. Effects of heterogeneous and clustered contact patterns on infectious disease dynamics[J]. PLoS Computational Biology, 2011, 7(6): e1002042.

[118]Miller J C. Epidemics on networks with large initial conditions or changing structure[J]. Plos One,2014,9(7):e101421.

[119]Valdez L D, Macri P A, Braunstein L A. Intermittent social distancing strategy for epidemic control[J].Phys Rev E, 2012, 85(3): 036108.

[120]Valdez L D, Macri P A, Braunstein L A. Temporal percolation of the susceptible network in an epidemic spreading[J].PLoS ONE,2012,7(9): e44188.

[121]Valdez L D, Macri P A, Braunstein L A. Temporal percolation of a susceptible adaptive network[J]. Physica A, 2013, 392(18), 4172-4180.

[123]Wang Y, Chakrabarti D, Wang C, et al. Epidemic spreading in real networks: an eigenvalue viewpoint[C]//Proceedings of the 22nd Symposium on Reliable Distributed Systems, Florece, Italy,2003: 25-34.

[124]Chung F, Lu L, Vu V. Spectra of random graphs with given expected degrees[J]. Proc Natl Acad Sci. 2003, 100 (11): 6313-6318.

[126]Kephart J O, White S R. Directed-graph epidemiological models of computer viruses[C]//Proceedings of the 1991 IEEE Computer Science Symposium on Research in Security and Privacy, New York, USA, 1991: 343-359.

[127]Chakrabarti D, Leskovec J, Faloutsos C, et al. Information survival threshold in sensor and P2P networks[C]//26th IEEE International Conference on Computer Communications, Anchorage, USA, 2007:1316-1324.

[128]Restrepo J G, Ott E, Hunt B R. Emergence of synchronization in complex networks of interacting dynamical systems[J], Physica D, 2006, 224 (1/2): 114-122.

[129]Feng J, Jirsa V K, Ding M. Synchronization in networks with random interactions: theory and applications[J], Chaos, 2006, 16 (1): 015109.

[130]Cohen R, Erez K, Ben-Avraham D, et al. Resilience of the internet to random breakdowns[J]. Phys Rev Lett, 2000, 85 (21): 4626-4628

[131]Vazquez A, Moreno Y. Resilience to damage of graphs with degree correlations[J]. Phys Rev E, 2003, 67 (1): 015101(R).

[133]Molloy M, Reed B. The size of the giant component of a random graph with a given degree sequence[J]. Combinatorics, Prob Comput, 1998, 7(3):295.

[134]Gómez S, Borge-Holthoefer J, Meloni S,et al. Discrete-time Markov chain approach to contact-based disease spreading in complex networks [J]. EPL, 2010, 89 (3): 38009.

[135]Gibson G J. Markov chain Monte Carlo methods for fitting spatiotemporal stochastic in plant epidemiology[J]. Appl Stat, 1997, 46(2):215-233.

[136]Ganesh A, Massoulié L, Towsley D. The effect of network topology on the spread of epidemics[C]∥ Proceeding of IEEE INFOCOM 2005, 2005: 1455-1466.

[137]O’Neill P D. A tutorial introduction to Bayesian inference for stochastic epidemic models using Markov chain Monte Carlo methods [J]. Math Biosci, 2002, 180(1/2):103-114.

[138]Dorogovtsev S N, Mendes J F F. Evolution of networks[J]. Adv Phys,2002, 51, 1079.

[139]Chakrabarti D,Wang Y, Wang C, et al. Epidemic thresholds in real network[J]. ACM Trans Inf Syst Secur,2008, 10: 1.

[140]Mieghem P V, Omic J, Kooij R. Virus spread in network[C]. IEEE ACMTrans Netw,2009, 17： 1.

[141]Ferreira S C,Castellano C,Pastor-Satorras R. Epidemic thresholds of the susceptible-infected-susceptible model on networks:a comparison of numerical and theoretical results[J]. Phys Rev E,2012, 86(4): 041125.

[142]Cator E, Mieghem P V. Second-order mean-field susceptible-infected-susceptible epidemic threshold[J]. Phys Rev E,2012, 85(5): 056111.

[143]Goltsev A V,Dorogovtsev S N,Oliveira J G, et al. Localization and spreading of diseases in complex networks[J]. Phys Rev Lett,2012, 109(12): 128702.

[144]Lee H K,Shim P S,Noh J D. Epidemic threshold of the susceptible-infected-susceptible model on complex networks[J]. Phys Rev E,2013, 87(6): 062812.

[145]Bogun M, Castellano C, Pastor-Satorras R. Nature of the epidemic threshold for the susceptible-infected-susceptible dynamicsin networks[J]. Phys Rev Lett,2013, 111(6): 068701.

[146]Givan O,Schwartz N,Cygelberg A, et al. Predicting epidemic threshold on complex networks: Limitations of mean-field approaches[J]. J Theor Biol,2011, 288, 21.

[147]Moreno Y,Pastor-Satorras R,Vespignani A. Epidemics outbreaks in complex heterogeneous networks[J]. Eur Phys J B, 2002, 26: 521-529.

[148]Callaway D S, Newman M E J,Strogatz S H,et al.Network robustness and fragility: percolation on random graphs[J]. Phys Rev Lett, 2000, 85, 5468-5471.

[149]Youssef M,Scoglio C. An individual-based approach to SIR epidemics in contact networks[J]. J Theor Biol, 2011, 283.136-144.

[150]Prakash B A, Chakrabarti D, Faloutsos M, et al. Get the flu (or mumps) Check the eigenvalue![DB/OL].[2014-06-30].arxiv.org/abs/10040060.

[151]Zanette D H. Critical behavior of propagation on small-world networks[J].Phys Rev E, 2001, 64(5)： 050901(R).

[152]Ben-Naim E,Krapivsky P L. Size of outbreaks near the epidemic threshold[J]. Phys Rev E, 2004, 69(5): 050901(R).

[153]Ben-Naim E,Krapivsky P L. Scaling behavior of the threshold epidemics[J]. Eur Phys J B, 2012, 85: 1.

[154]Kessler D A, Shnerb N M. Solution of an infection model near threshold[J]. Phys Rev E, 2007, 76(1): 010901(R).

[155]Khalleque A, Sen P. The susceptible-infected-recovered model on a Euclidean network[J]. J Phys A: Math Theor,2013, 46(9): 095007.

[156]Holme P. Extinction times of epidemic outbreaks in networks[J]. PloS ONE,2013, 8: e84429.

[157]Hong H,Ha M,Park H. Finite-size scaling in complex networks[J].Phys Rev Lett, 2007, 98(25): 258701.

[158]Binder K, Heermann D W. Monte carlo simulation in statistical physics,[M]. 5th ed Berlin: Springer-Verlag, 2010.

[159]Ferreira S C,Ferreira R S, Castellano C, et al. Quasistationary simulations of the contact process on quenched networks[J]. Phys Rev E, 2011, 84(6): 066102.

[160]Shu P, Wang W, Tang M, et al, Numerical identification of epidemic thresholds for susceptible-infected-recovered model on finite-size networks [J]. Chaos, 2015, 25 (6), 063104.

[161]Wang W, Liu Q H, Zhong L F, et al, Predicting the epidemic threshold of the susceptible-infected-recovered model [DB/OL]. [2014-06-30]. www.nature.com/articles/srep24676.

[162]Li R, Wang W, Di Z. Effects of human dynamics on epidemic spreading in Cote d'Ivoire[DB/OL]. [2014-06-30]. arxiv.org/abs/160500899v1.

[163]Li R, Tang M, Hui P M. Epidemic spreading on multi-relational networks[J]. Acta Phys Sin, 2013, 62(16): 168903.

[164]Xu E H W, Wang W, Xu C, et al, Suppressed epidemics in multirelational networks [J]. Phys Rev E， 2015, 92 (2): 022812.

[165]Wang W, Tang M, Yang H, et al, Asymmetrically interacting spreading dynamics on complex layered networks [J]. Sci Rep, 2014, 4: 5097.

[166]Liu Q H, Wang W, Tang M, et al, Impacts of complex behavioral responses on asymmetric interacting spreading dynamics in multiplex networks [J]. Food & Nutrition Sciences, 2015,6(6):555-561.

[167]H Yang, M Tang, T Gross, Large epidemic thresholds emerge in heterogeneous networks of heterogeneous nodes [J]. Sci Rep. 2015, 5, 13122.

[168]Wang W, Tang M, Zhang H F, et al, Dynamics of social contagions with memory of nonredundant information [J]. Phys Rev E, 2015, 92(1): 012820.

[169]Wang W, Tang M, Shu P, et al, Dynamics of social contagions with heterogeneous adoption thresholds: crossover phenomena in phase transition [J]. New J Phys, 2016, 18(1): 013029.

[170]Wang W, Shu P, Zhu Y X, et al, Dynamics of social contagions with limited contact capacity [J], Chaos, 2015, 25(10): 103102.

(责任编辑耿金花)

Review of Threshold Theoretical Analysis About Epidemic Spreading Dynamics on Complex Networks

LI Ruiqi1,2, WANG Wei1, SHU Panpan1,YANG Hui1,PAN Liming1, CUI Aixiang1, TANG Ming1

(1.Web Science Center, University of Electronic Science and Technology of China, Chengdu 611731,China;2.School of Systems Science, Beijing Normal University, Beijing 100875, China)

Abstract:In this Review, we mainly focus on solving the threshold theoretically, and introduce seven common methods, including Mean-Field theory, Pair-wise Approximation, Master equation, Generating Function, Edge Percolation, Cavity method, Edge Classification and Spectral analysis. We also summarize the difference of thresholds between the SIS and SIR model. We are aiming to provide a clear picture for beginners and a good reference for researchers.

Key words:complex networks; epidemic spreading;outbreak threshold; theoretical analysis

文章编号：16723813(2016)01000139;

DOI:10.13306/j.1672-3813.2016.01.001

收稿日期：2014-12-24

基金项目：国家自然科学基金(11105025,11575041)

作者简介：李睿琪(1992-)，男，河北张家口人，博士研究生，主要研究方向为复杂网络传播动力学、城市演化动力学、交通分析、人类动力学、风险投资网络分析、大数据分析与建模。

中图分类号：N93;N94

文献标识码：A