基于无标度网络的埃博拉病毒传播模型研究

毛海舟（浙江工贸职业技术学院基础部，浙江 温州 325000）

赵利彬（闽江学院数学系，福建 福州 350108）

埃博拉病毒的爆发与传播给非洲乃至全球的经济社会带来了巨大的影响，细致地对病毒传染病模型进行研究和分析，不仅有利于消除民众不必要的恐惧或盲目乐观的情绪，也可以为政府的应急决策提供参考。

传统的数学模型往往用一组确定性的微分方程来描述传染病传播的动力学行为［1，2］，这使得研究者很容易分析模型中各个参数对系统演化行为的影响，但这些参数所对应的实际意义却难以精确刻画。因此，一种控制行为到底会对哪些参数产生影响，产生多大的影响都是很难确定的，这也是微分方程模型解释能力强于控制能力的原因。近年来，复杂系统理论提出了一种崭新的超越还原论思想的建模方法，即为系统设定一些规则后，让系统在一定的环境下自发的演化，然后考察系统演化过程中表现出来的若干性质。复杂网络上的传播模型能够很好地体现复杂系统建模的思想，并且由于考虑了传播网络拓扑结构本身对传播行为的影响，因此可以得到很多传统模型不能得到的结果［3］。下面，笔者基于复杂网络理论建立了埃博拉病毒传播的随机演化模型，并将该模型外推得到仿真模型，通过试验分析来了解埃博拉病毒在一个封闭系统（如一个国家或一个地区等）中的传播模式。

1 复杂网络（Barabasi－Albert模型）理论

自然界中的许多系统都可以用复杂网络来描述，如Internet是由路由器和传输介质组成的复杂网络，人脑是由神经元通过互连形成的复杂网络，而社会则是人与人通过各种各样的关系联系起来的。复杂网络的广泛存在，使得众多科学家致力于复杂网络自身性质特别是动力学性质的研究［3～5］，如社会人际网络怎样促成了传染病的传播，如何在Internet上最有效的防止病毒的攻击，什么样的拓扑结构在冲击下最为稳定。目前，有关复杂网络系统的研究得益于数据采集能力的提高和高性能计算工具的出现，在此基础上，各种实际网络复杂的拓扑结构被揭示出来。

1998年，Watts等［6］提出了小世界网络的概念，称之为WS模型，该模型具有很小的平均距离和较大的簇系数，即通常意义上的小世界效应［3，4］。WS模型尽管具有小世界效应，但不足于用来刻画真实系统中的复杂网络，其与早期的一些随机网络一样，其顶点度的分布符合Possion分布，但针对真实系统的复杂网络统计分析表明，复杂网络的顶点度满足Power－law分布［5，7］。因此，用于模拟现实的传播网络至少应当具备小平均距离、大簇系数和Power－law度分布3个特征。1999年，Barabasi［7］提出了构造复杂网演化算法，可以同时实现复杂网络的3大特征，因而笔者进行试验分析时均采用Barabasi－Albert模型。

大量的实证研究表明，真实网络几乎都具有小世界效应［8～10］，但同时研究者也发现许多网络的度分布都服从幂率分布［9］，也就是说具有某个特定度的节点数目与这个特定的度之间的关系能够用一个幂函数近似表示。幂函数是一条下降速度相对缓慢的曲线，这就使得度很大的节点可以在真实网络中存在。由于幂函数具有这种标度不变性，因此把节点度服从幂率分布的网络叫做无标度网络。

Barabasi［7］考察了网络的生成机制，发现增长和择优连接是实际网络演化过程中的2个基本特性，同时还提出了能够产生无标度特性的第1个网络模型，即Barabasi－Albert模型。该模型的构造算法如下：①增长。从一个包含m0个节点的初始网络开始，每个时间步引入一个新的节点，并且连接到m个已经存在于网络中的节点上，这里有m≤m0。②择优连接。新节点与一个已经存在的节点i相连接的概率Πi与节点i的度数ki满足如下关系：

经过t时间步后，该算法产生一个具有m0＋t个节点和mt条边的网络。

BA网络模型具有以下特性：

1）度分布。Barabasi等［11］用平均场的方法得到BA网络的度分布为：

这表明BA网络的度分布可以由幂指数为－3的幂律函数来近似描述。

2）平均路径长度。BA网络的平均路径长度为：

由式（3）可见，BA网络具有小世界特性。

3）聚集系数。BA网络的聚集系数为：

由式（4）可见，当网络规模充分大时，BA网络的聚集特性已不再明显［7］。

4）BA网络中存在中枢节点（Hub节点）和富者愈富现象，幂律度分布的胖尾特性导致无标度网络中存在少数含有大量连接边的中枢节点，择优连接特性必然引起富者愈富现象。

5）BA网络既具有鲁棒性又具有脆弱性。面对节点的随机失效，网络具有鲁棒性；面对蓄意攻击时，因为中枢节点的存在，网络变得极其脆弱。

2 埃博拉病毒传播模型的建立与分析

2.1 基本假设

埃博拉病毒是一种传染性强的呼吸系统疾病，其传播的动力学特征除了受传播网络拓扑结构的影响外，还与病毒活性、易感人群抵抗能力、隔离措施等诸多因素有关。为便于研究，对基于无标度网络的埃博拉病毒的传播模型作如下假设。

1）每个病人的传染能力相同，即不考虑超级传播事件的影响。

2）病毒的活性在演化过程中是不变的，也不考虑病毒的变异行为，即忽略影响埃博拉病毒活性的通风条件、温度等因素。

3）群体中每个个体的免疫能力相同，并不受年龄结构的影响，即忽略不同个体的不同免疫力，这是因为模型所关心的是整体表现出来的统计性质而非特定个体的情况。

4）每个病人的潜伏期相同，且在潜伏期中不具有传染性。

5）在没有隔离控制措施下每个病人从发病开始到具有传染性的时间是一样的。

6）病人从发病到传染性消失的这段时间内传染能力强弱不随时间变化。

2.2 符号说明

用图G（V，E）表示埃博拉病毒的传播网络，其中V是顶点集，E是边集，每个顶点代表系统中的1个人，2个人之间有1条边相连，表示这2个人有可能发生接触，从而可能导致病毒扩散，这里记为相邻个体。相邻个体在一个单位时间（d）内的接触概率用r表示，每一次事实接触都属于以下3种情况之一：病人与病人的接触、健康者之间的接触、病人与健康者之间的接触。显然，其中病人与健康者之间的接触值得关注，因为这种接触是否造成健康者染病，与病人的传染能力和健康人的抵抗能力相关。假设每个病人的传染能力是一样的，用参数η表示，而易感人群的免疫能力也是一样的，用参数ζ表示，λ＝η／ζ表示该次传播的有效致病因子，健康者每次和病患的接触都会累积有效致病因子，若其超过致病阈值χ，则认为该健康者已经染病。健康者染病后会经过一段潜伏期，潜伏期长度用tq表示，感染者在潜伏期里不具有传染性。病人在病发后的传染期内具有传染性，传染期的长度记为ty。初始化时，系统随机指定一个个体处于疾病初发时期，然后系统根据既定规则自发地以单位时间（d）迭代演化。

埃博拉病毒的传播大致可以分为3个阶段，即自由传播阶段、过渡阶段和控制阶段。自由传播阶段时人们完全没有认识到埃博拉病毒的危害，一切行为如常，政府也没有积极控制的行为；控制阶段时政府已出台严格的隔离、诊治和预防措施，民众对埃博拉病毒的危害性有充分的认识，开始广泛采取喷洒消毒液、戴口罩、服用中药、较少与人接触等预防措施，控制阶段与自由传播阶段的不同表现为：ty锐减（由于隔离政策），r下降和ζ上升（由于民众防范意识上升）；过渡阶段是指由于埃博拉病毒染病人数积累到一定数量后，民间已经有了一些消息，但政府还没有正式的疫情解释和行之有效的防范措施，这时候ty会有少许的下降，主要是因为部分患者较早意识到问题的严重性，从而主动到医院就医，同时r也会有一些下降并伴随着ζ的少许上升，但远不如控制阶段明显，因此叫做过渡阶段。

2.3 SIRSLS模型的构建

真实流行病学中，低易感个体如果和感染个体接触不会马上被感染，只是失去这部分免疫力变成易感个体，针对这种感染机制，医学界提出了一种新的符合实际的传染病传播模型——具有低易感的传染病模型（即SIRSLS模型）。

对于无标度网络，由于度分布呈幂律形式，即网络中的度很不均匀，因而对于不同度值的节点需要作不同的处理。用sk（t）、ik（t）、rk（t）、lsk（t）分别表示网络中度为k的易感节点、感染节点、康复节点和低易感节点在t时刻占总节点的比例，且满足如下关系式：

而在整个BA无标度网络中t时刻感染节点的比例为：

式中，ik（t）表示度为k的染病节点在t时刻的感染密度；P（k）表示节点的度分布情况的分布函数。

根据平均场理论，SIRSLS模型在BA无标度网络中的传播过程可以用下面的微分方程组表示：

式中，Θk（t）在t时刻度为k的节点有一条边指向染病节点的概率；ik（t）表示度为k的染病节点在t时刻的感染密度；λ表示通过一条边与染病节点相连时的感染概率；δ表示免疫丧失率；μ表示免疫保留率。

则在t时刻度为k的节点有一条边指向染病节点的概率为：

式中，〈k〉表示网络的平均度，即网络上所有节点的平均值。

方程组（7）的解析解难以求解，当传播过程达到稳定状态（或传播结束）时，传染病的传播机理有特殊意义，因此求其特解。

当t→ ∞ 时，ik（t）、rk（t）、lsk（t）不再发生变化，此时，由此可得方程组的稳态解。

其传播阈值为：

与小世界网络一样，在无标度网络中SIRSLS模型的传播阈值也是和网络的拓扑结构、免疫丧失率和免疫保留率有关。但是在无限大的无标度网络中，〈k〉→∞，此时传播阈值趋于零。当μ＝0时，其传播阈值为

2.4 系统演化

试验中取tq＝13，ty＝20，χ＝0.99，η＝ζ＝1.0，网络中有500人。利用建立的无标度网络中SIRSLS模型进行系统演化500d后累计染病人数的时间序列（在不同接触概率r下）如图1所示。

若任由疾病自由传播，即使接触概率较低，也会波及系统中的大多数人，但并非所有的人都会被传染，这一点是传统的传染病模型较难解释的。研究发现，染病人数上限存在的原因，不是源于某些拥有特强免疫能力的个体，而是网络自身从高扩散能力结构逐步演化为低扩散能力结构，在不同τ（网络扩散能力因子）值下，其网络扩散能力演化过程如图2所示。从图2可以看出，不同的网络扩散能力因子都存在感染人数上限，且随着τ值的下降感染人数上限也随之下降，因而通过网络结构自身的这种免疫机制可以阻止埃博拉病毒的传播。

因为只有处于传染期的病人与健康人之间的接触，才可能导致埃博拉病毒的传播，所以用连接正处于传染期的病人与健康人的边的数目来衡量该网络目前的扩散能力。可以发现，在传播过程中，τ值经历了从小到大、又从大到小的过程。初始时患病者很少，所以τ值也小，然后随着患者的增多，τ值也随之增加；在传播了一段时间后，由于一些患病者失去了传染能力，而且健康者人数下降，因此τ随之下降。当τ值下降到某个很小的值后，新的传染事件成为一个小概率事件，如果在一段时间内该事件不发生（所有对τ值有贡献的病人经过了传染期），τ值将变成0值，从而累计患病人数达到上限。图1中2条曲线分别表示接触概率为0.020和0.015的2次试验中τ值的演化，从中可以看到类似单峰的形态。

图1 自由传播时染病人数的时间曲线图

图2 网络扩散能力演化示意图

3 几内亚埃博拉病毒传播模型的仿真试验

根据世界卫生组织发布数据显示，截止2014年8月，几内亚共有495例感染，死亡363例；利比里亚感染516例，死亡282例；塞拉利昂感染691例，死亡286例。可以看出，几内亚死亡人数居首位。下面，以几内亚作为研究对象，分析埃博拉病毒在该区域的传播情况。

虽然基于复杂网络的埃博拉病毒传播模型讨论了埃博拉病毒传播的动力学性质，但该模型处理的只是一个规模为500人的封闭系统，而几内亚的人口数目超过千万，超出处理能力数个量级。因此，该模型并不能解释真实数据。下面，笔者对该模型做外推处理，使其能够在一定程度上解释真实数据。

在模拟真实的仿真模型中，整个系统被分为医院外和医院内2个部分。在医院外选择了7000个节点，每个节点代表总人数为2000的一个相对封闭的子系统。任意2个子系统都有可能以一个较小的概率产生“接触”，在单位时间里，2个子系统的接触使得分别位于2个子系统内部的一对节点状态互换。一个健康的子系统是指系统内部所有个体都是健康的，其他情况的系统都被认为是不健康的，埃博拉病毒从一个不健康系统向健康系统的传播源于该系统中的染病个体在与健康系统接触中被选中并进行了状态互换。子系统内部的埃博拉病毒传播情况可以用随机演化模型处理，由于人数都是2000，实际上可以统一处理。系统初始时随机选择某个系统中的某个个体为初发疾病的患者，然后按照既定的规则自发演化。病人在发病一段时间后将被移入医院，这段时间的长短与政府隔离政策的执行力度有关。医院系统共有10000个医生，病人在进入医院后立刻与其中5名医生形成接触关系，并且每天与这些医生接触1次，当然，由于严格的消毒措施，医生的免疫防护能力与一般易感人群时不一样的。需要说明的是，在这个模型的初步讨论中，依然不考虑超级传播者事件。同时假设治愈出院的病人数量不对系统的整体行为产生影响。

3.1 数据预处理

在仿真试验之前，首先需要进行一些数据挖掘和预处理的工作。以几内亚的数据为例，可以得到的数据是累积到某天的确诊人数，这仅是模型中已被医院收入的染病者的一部分，另一部分应从真实数据中累积到某天的疑似病例中产生，这里就需要估计疑似病例中有多少是埃博拉病毒感染者。为此，笔者给出了如下处理办法。

定义y（t）表示t时刻后发现的疑似病例总数，z（t）表示t时刻后疑似病例转为确诊病例的总人数。设原确诊人数的时间序列为｛x（t）｝，修正后的时间序列为｛x＊（t）｝，时间序列的起点为t0，t时刻从人群中发现的确诊病例为l（t），新增疑似病例为h（t）。这里l（t）和h（t）可以看作初始时刻确诊病例总数和疑似病例总数。y、z、x、h、l是已知的，x＊是待求的。下面给出一种简单算法：

几内亚埃博拉病毒病死亡人数随时间变化修正前后死亡人数对比图如图3所示。需要说明的是，仿真试验时使用修正后的数据。

由于无法获知最初感染者的发病时间，因此最初感染者发病时间距离公布数据的t0时刻（几内亚埃博拉首例报告期）的长度也是有待确定的参数之一。为此选择不同的参数进行了大量试验，尽可能使得因埃博拉病毒死亡人数的时间序列与真实数据相吻合，再利用最小二乘法从中选择最优拟合参数来计算几内亚在未来一年的感染人数。

3.2 仿真结果分析

在无疫苗条件下几内亚埃博拉病毒感染者每日死亡人数预测曲线图如图4所示。从图4可以看出，爆发初期，死亡人数是每天7人（自由传播阶段），此时政府采取了相关措施，如隔离、灭菌和保护等。由于埃博拉病毒通常在死亡后的感染者体内存活达数个小时，传播源往往是死亡的动物或病人尸体。直接或间接物理接触死亡病人的体液、血液、污染物（衣服、被单）等都有可能导致感染，因而必须对病人尸体进行火化处理才能有效防止病毒感染。在经历一段时间的隔离措施后，随着生活环境的改善，病毒的传染速度将会得到控制，因而每天的死亡人数下降到2人左右（过渡阶段）。但是，由于人群对防疫措施的不够重视，比如允许病人外出等原因，导致死亡人数又会很快上升。因为政府重新出台防疫条例并加强改善居住环境，使得病毒传染重新得到控制，因而每日死亡人数逐渐稳定在2～3人（控制阶段）。当然，后期如果要彻底消灭埃博拉病毒，还必须加速研制疫苗和药物。

图3 几内亚埃博拉病毒感染者随时间变化修正前后死亡人数对比图

图4 在无疫苗条件下几内亚埃博拉病毒感染者每日死亡人数预测曲线图

4 结语

城市化和全球化使得人类交往更加密集、快速和频繁，这增加了传染病扩散的可能性，加之传染病传播途径增多以及部分病原体变异毒性增强等因素，更加大了传染病防控的难度。复杂网络理论为研究传染病的防控提供了新的方法，为控制真实网络中的病毒传播提供了重要的理论依据。以几内亚为例，利用区域埃博拉传播模型对埃博拉病毒在该区域的传播情况进行分析，认为埃博拉病毒传播率在小规模范围内有所波动，但是最终会下降并趋于稳定状态。针对上述情况，可以采取相关措施，即通过改进免疫策略来阻止传染病传播、加速疫苗与病毒治疗药物的研究和加大疫情防控知识的宣传，由此减少并最终阻断埃博拉病毒的传播。

［1］Baile N T J.The Mathematical Theory of Infectious Diseases［M］.London：Griffin，1975.

［2］Murray J.Mathematical Biology［M］.Berlin：Springer－Verlag，1993.

［3］Strogatz J D.Exploring complex networks［J］.Nature，2001，410：268～276.

［4］Albert R，Barabasi A L.Statistical mechanics of comples networks［J］.Review of Modern Physics，2002，74：47～91.

［5］Wang X F.Complex networks：Topology，dynamics and synchronazation［J］.International Journal of Bifurcation ＆ Chaos，2002，12：885～916.

［6］Watts D J，Strogatz S H.Collective dynamics of“small world”networks［J］.Nature，1998，393：440～442.

［7］Barabasi A L，Albert R.Emergence of scaling in random networks［J］.Science，1999，286：509～512.

［8］Albert R，Barabási A L.Statistical mechanics of complex networks［J］.Reviews of Modern Physics，2002，74：47～49.

［9］Ebel H，Mielsch L I.Borbholdt S.Scale－free topology of email networks［J］.Physical Review E，2002，66：1162～1167.

［10］Sen P，Dasgupta P，Chatterjee A，et al.Small－world properties of the Indian railway networks［J］.Physical Review，2003，67：210～215.

［11］Barabasi A L，Albert R，Hong H.Mean－field theory for scale－free random networks［J］.Physical，1999，272：173～187.