CN102945310A - 基于自治计算的流行病传播网络建模与推断 - Google Patents

Abstract

基于自治计算的流行病传播网络建模与推断，包括下列步骤：采用多自治体建模流行病传播网络和流行病传播过程，采用蒙特卡罗模拟和负反馈机制从流行病监控数据中推断出传播网络结构及与流行病有关的生物学参数。本发明为流行病监控与防治提供了一种新方法，与现有方法相比，本发明具有如下主要优点：1)提供了一种适应范围更广的流行病传播网络建模方法；2)提供了一种有效的流行病传播网络推断方法，可从流行病监控数据中推断出流行病传播网络结构以及与流行病有关的生物学参数；3)本发明所提供的方法可用于流行病风险评估和流行病防治策略的有效性验证。

Description

基于自治计算的流行病传播网络建模与推断

技术领域

本发明属于信息技术和公共卫生的交叉领域，尤其涉及流行病控制、复杂网络和数据挖掘等领域。

背景技术

流行病的每次爆发都会给人类社会带来巨大损失。建立理论模型理解和模拟流行病的传播与爆发过程，建立监控系统收集流行病的传播数据，进而基于模型与监控数据进行准确的风险评估，帮助决策者制定相应的预防和应急方案，是有效控制流行病传播与爆发、最大限度减少生命财产损失的新途径。

流行病的传播过程由流行病模型和传播网络共同决定。医学领域对流行病的研究有较长历史，提出了多种数学模型，其中最具代表性的是基于动力系统的SIR模型和SIS模型，以及基于它们的多种改进模型。对这类模型的早期研究大都假设流行病的传播媒介均匀分布，没有考虑人口在空间上分布不均对流行病传播过程的影响。人口等流行病传播媒介在空间上的分布可采用网络结构进行抽象与建模。1983年，Grassberger将流行病传播问题建模为键逾渗(bond percolation)问题，最先在网络框架下研究了SIR模型的基本性质。1998年小世界网络模型提出后，Newman最早研究了小世界网络结构对SIR模型的影响。在此基础上，后续工作深入研究了小世界网络模型的各种参数(如随机长连概率)对SIR模型的影响，小世界网络传播动力方程的分形、不动点、混沌和周期分叉等多种数学性质。一些工作还将SIS和SEIR等模型推广到网络上，讨论这些流行病模型在小世界网络上的传播特性。1999年无标度网络模型提出后，人们发现包括社会网络在内的很多实际网络都是非均匀网络，其度分布满足幂率而非泊淞分布。此后，各种传播模型在无标度网络上的临界值分析成为新的研究热点。

以上研究在分析网络结构对流行病传播产生的影响时，大都事先假定传播网络结构已知(如某种规则网络)或者服从某种特定的网络生成模型(如小世界网络和无标度网络模型)，进而分析流行病模型在不同网络结构上的临界值及各种宏观动力学性质。然而实际问题所涉及的传播网络具有更高的复杂性，往往不能满足这些假设。实际传播网络与假设网络结构上的差异会显著影响分析的准确性与可靠性。随着信息技术在医疗领域的广泛应用，记录各种流行病传播过程的监控数据(surveillance data)与日俱增，这些数据记录了与流行病传播和爆发相关的时空信息，包括入院、确诊、出院、死亡等医疗过程的统计数据，地理位置信息以及当地人口的统计数据。如果能从这些监控数据中推断出导致传播现象发生的传播网络结构，就有望解决真实传播网络难以建造的问题。

基于观测到的传播数据推断传播网络结构是研究流行病传播规律的新思路，在流行病等医学领域还未见到相关研究，类似问题最先受到计算机领域研究者的关注。该方面的早期工作主要针对信息传播现象，提出能够根据传播现象推断出隐含信息传播网络的学习算法。2004年，HP公司的Adar等人较早开展了这方面的尝试，基于博客(blog)被话题(topic)感染时间的先后次序、博文的文本信息以及博客主页间的超链接结构推断出博客间的隐含影响关系。有监督的学习方法需要大量正确标注的样本，而在实际问题中这些样本往往难以获取。IBM公司Gruhl等人提出了模拟博客阅读和撰写博文行为的主题传播模型，进而提出了类似EM算法的非监督参数估计算法，根据博文发表时间的先后次序估计出博客间的影响概率，建立隐含的主题传播网络。2010年，斯坦福计算机系的研究小组发表了从多组时间序列中推断出隐含传播网络结构的论文，该工作推广了信息传播网络推断问题，使之能描述流行病传播，并针对独立级联模型给出一种基于子模块(submodular)贪心优化的近似学习算法，可推断出近似最优的隐含传播网络。在此基础上，该小组又相继提出了缺失观测数据的处理方法和基于凸规划的传播网络的推断算法。

以上工作从算法和实验的角度验证了从传播数据中学习传播网络的可行性，可为研究流行病传播网络推断问题提供参考。然而现有的传播网络推断方法不能有效应用于流行病监控数据。其主要原因如下：①现有方法所能处理的传播数据与可获得的流行病监控数据形式不同。现有方法处理的传播数据都是基于个体(individual based)的时间序列，以“某个体在某时间被感染”的方式依次记录信息在网络上的传播过程。然而如此细粒度的流行病监控数据通常很难获得，各级疾病控制中心(CDC)和各地区医院提供的流行病监控数据通常是基于空间位置(location based)的时间序列，以“某区域在某个时间段内感染、出院或死亡多少人”的方式记录流行病的传播和爆发过程。②现有传播网络推断方法假设被观察到的传播数据足够多，需要多次传播过程对应的时间序列数据才能较为准确的推断传播网络结构。然而在实际应用中，可获得的流行病监控数据通常十分有限。③流行病监控数据具有空间和时间上的多尺度(multi-scale)特性，这些特性都为现有的传播网络推断方法带来困难。④现有传播网络推断方法主要针对信息传播过程，大多采用级联和阈值模型建模信息的传播机制，难于扩展到SIR、SIS或SEIR等模型，而这些模型非常适合于刻画流行病的传播和爆发过程。

流行病的传播过程是复杂系统上的动力学过程。自治计算(AOC，Autonomy Oriented Computing)为建模和分析复杂系统及发生在其上的动力学行为提供一种新的理论和方法，采用“自底向上”的“白盒”方法分析和发现隐藏在复杂现象背后的基本机制和规律，为此我们采用自治计算方法建模和推断流行病传播网络。

发明内容

现有的传播网络结构推断方法大都面向信息传播过程，所能处理的数据与可获得的流行病监控数据形式和特性均不相同，不适合处理具有粗粒度、时空多尺度和数据缺失等特性的流行病监控数据。针对以上问题，提出了基于自治计算的流行病传播网络建模方法和网络结构推断方法。该方法采用多自治体建模传播网络结构和流行病传播过程，采用蒙特卡罗模拟和负反馈机制调节系统参数，以缩小模拟系统产生的数据与实际监控数据间差异为目标，改变自治体的行为，促使模拟系统向真实系统逐步演化，以此方式推断出传播网络结构及与流行病相关的生物学参数。

附图说明

图1所示的流程图给出本发明提供的基于自治计算的流行病传播网络建模与推断方法的基本步骤。

图2所示自治体的状态转换规则。

图3所示的流程图给出了模拟系统参数推断的具体步骤。

图4所示2009年香港地区H1N1猪流感传播的真实数据与模拟系统产生的模拟数据。

图5(a)所示采用模拟系统模拟流行病传播过程；图5(b)所示模拟系统的参数估计过程。

图6所示缺失数据对爆发趋势估计的影响。

图7所示时间窗大小对爆发趋势估计的影响。

图8所示基于真实数据的流行病传播网络推断过程。图8(a)所示香港地图；图8(b)所示香港岛的区域划分和真实数据；图8(c)所示香港岛的流行病传播网络。

图9所示基于传播网络推断的流行病感染风险评估。图9(a)所示第40天时累积感染病例的空间分布图；图9(b)所示第100天时感染风险的空间分布图；图9(c)所示第40天至第100天期间的感染风险的时空分布图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细说明。实施例仅用于说明本发明，但不用来限制本发明的应用范围。

图1，一种基于自治计算的流行病传播网络建模与推断方法的流程开始于步骤101。

步骤102构造模拟流行病传播过程的自治体系统(命名为D-AOC系统)，具体如下：

首先，建模流行病传播网络。

流行病传播网络包含N个节点。节点可表示具有不同大小的空间区域(如户、村、乡、镇、县、市)，也可表示复合群体(如不同年龄段的人口)。节点间的链接表示流行病传播途径，链接权重表示传播强度。流行病传播网络表示为一个N×N矩阵A＝[a_ij]_N×N，a_ij表示节点i与节点j的链接权重，满足0≤a_ij≤1和

其次，建模自治体。

D-AOC系统中的自治体定义如下：

流行病传播网络中的每个节点可包含多个自治体，节点i中自治体定义为如下的一个5元组：

(i，s，β_i，γ_i，τ_i)

其中，i表示自治体所在网络节点的标识符；s表示自治体当前所处的状态，自治体在每个时刻处于S(易感)、E(潜伏)、I(感染)和R(恢复)四种状态之一；β_i，1/γ_i，1/τ_i分别表示节点i中自治体对流行病的抵抗力、流行病的潜伏期和患病之后的康复期，满足0≤β_i，γ_i，τ_i≤1。令B＝(β₁，…，β_N)，H＝(γ₁，…，γ_N)，F＝(τ₁，…，τ_N)，向量B，H和F表示与流行病相关的生物学参数。

图2给出自治体的状态转换规则，图2中虚线圈和实线圈分别对应不可观测和可观测状态。自治体的状态转换规则描述如下：

规则1：当位于节点i中的自治体处于S状态时，将以

的概率进入E状态。

称为节点i在t时刻的感染强度，由自治体对流行病的抵抗力和流行病传播网络共同决定，计算如下：

${\mathrm{\λ}}_i^{\left(t\right)}={\mathrm{\β}}_i{\mathrm{\Σ}}_{j=1}^N\frac{I_j^{\left(t\right)}}{P_j^{\left(t\right)}}{a}_{\mathrm{ij}}$

其中β_i表示节点i中自治体对流行病的抵抗力，

表示t时刻节点j中处于I状态的自治体数目，

表示t时刻节点j中自治体的总数，a_ij表示节点i和节点j之间的网络链接权重。

规则2：自治体处于E状态已经持续1/γ_i时间后进入I状态；

规则3：自治体处于I状态已经持续1/τ_i时间后进入R状态；

规则4：自治体进入R状态后则永久停留在该状态，意味着它获得了对流行病的终身免疫。

步骤103给出了采用蒙特卡罗模拟和负反馈机制，从流行病监控数据中推断出D-AOC系统参数的方法，具体如下：

令D表示实际观测到的流行病监控数据，定义为一个N×T矩阵D＝[d_it]_N×T，d_it表示t时刻节点i中新出现的流行病患者数目，T表示监控时间长度。在实际应用中，流行病监控数据D通常具有如下几个特性：

①时空特性。流行病监控数据D由N个时间序列构成，D矩阵的第i行对应在节点i上观察到的时间序列，记录了节点i所表示的空间区域或复合群体新增病例随时间的变化情况。

②时空多尺度特性。流行病传播网络中的节点可具有不同的空间大小，一些节点对应“村”和“乡”等较小的空间区域，另一些节点对应“县”和“市”等较大的空间区域。不同空间尺度的区域对应的时间序列可能具有不同的时间粒度。如“村”和“乡”等小尺度区域以“天”或者“周”为时间单位普查和统计流行病数据，而“县”或“市”等大尺度区域以“月”或者“季度”为时间单位普查和统计新增病例。

③数据缺失。实际的监控数据常常存在数据缺失的情况，大体可分为两种情况。第一，为了降低获取监控数据的成本，通常不对所有节点而仅对部分节点进行监控和普查，因而会造成某些节点对应的时间序列数据整组缺失。第二，如②中讨论，空间范围大的节点仅提供时间粒度较大的数据，造成细粒度数据缺失。如某些节点仅能提供以周或月为单位的统计数据，而不是每天的统计数据。

基于D，图3模拟系统参数推断的过程开始于步骤301。

步骤302：初始化系统参数。随机生成2m组系统参数，令P_i表示第i组系统参数，V_i表示与P_i对应的一组系统参数，其中1≤i≤m；令迭代次数Q等于0；令LP和GP分别表示局部最优和全局最优的系统参数；初始时，令LP和GP均为P₁；初始时，令E_max为一个充分大的整数。

步骤303：运行D-AOC系统的蒙特卡洛模拟，获取流行病模拟数据

对i从1到m，基于P_i，运行D-AOC系统的蒙特卡洛模拟，获取模拟过程产生的流行病监控数据

其中

表示t时刻节点j中新出现的处于I状态的自治体数目，T表示模拟过程持续的时间。

蒙特卡洛模拟，按如下方程进行：

$\left\{\begin{array}{c}{\mathrm{\ΔS}}_i^{(t+1)}=-{\mathrm{\λ}}_i^{\left(t\right)}{S}_i^{\left(t\right)}\\ {\mathrm{\ΔE}}_i^{(t+1)}={\mathrm{\λ}}_i^{\left(t\right)}{S}_i^{\left(t\right)}-{\mathrm{\γ}}_i{E}_i^{\left(t\right)}\\ {\mathrm{\ΔI}}_i^{(t+1)}={\mathrm{\γ}}_i{E}_i^{\left(t\right)}-{\mathrm{\τ}}_i{I}_i^{\left(t\right)}\\ {\mathrm{\ΔR}}_i^{(t+1)}={\mathrm{\τ}}_i{I}_i^{\left(t\right)}\end{array}\right.$

其中，γ_i和τ_i含义如上所述；

分别表示t时刻节点i中处于S状态、E状态和I状态的自治体数目；

和

分别表示t+1时刻节点i中处于S状态、E状态、I状态和R状态的自治体数目的变化量。特别的，

表示节点i中新出现的处于I状态的自治体数目。

步骤304：计算

偏离D的程度

其中

表示矩阵X＝[x_it]_N×T的二阶范式。

步骤305：如果

足够逼近真实监控数据D，则停止算法并返回GP作为从D推断出的系统参数。

步骤306：更新LP和GP。根据以上计算出的m个偏差E_i(1≤i≤m)，更新LP和GP，具体如下：

首先，令E_u是m个偏差中的最小者；将LP更新为P_u；

其次，如果E_u小于E_max，则将GP更新为P_u，同时将E_max更新为E_u。

步骤307：更新P_i和V_i。对i，从1到m执行如下操作：

根据LP和GP采用负反馈机制更新P_i和V_i：

首先，将V_i更新为wV_i+c₁·rand·(LP-P_i)+c₂·rand·(GP-P_i)；

其次，将P_i更新为P_i+V_i；

其中，rand表示随机数发生器，返回0到1范围内的随机数；w，c₁和c₂是3个预先设置的参数。

步骤308：令迭代次数Q增1。

步骤309：如果Q小于规定的迭代次数则返回步骤303执行，否则停止算法，返回GP作为从D推断出的系统参数。

例1：基于真实数据的模型验证

检验上述流行病传播网络模型的合理性。实验数据采用2009年H1N1猪流感在香港爆发的真实监控数据，该数据记录了2009年5月至2010年3月香港地区累计病例的变化情况。图4红色实线所示，5月1日至6月12日期间已出现H1N1患者，但并未开始大规模传播，6月12日之后感染人数呈快速增加。在H1N1爆发的中后期，香港政府进行了患者隔离和疫苗接种等措施，对H1N1传播进行了控制，新增病例逐渐减少。由于本发明提供的流行病传播模型不考虑人为控制措施对流行病传播的影响，因此选用H1N1爆发过程中人为干预较少的前3个月数据来验证模型的合理性。

根据年龄，香港受H1N1影响最大的人口可划分为5个年龄组：5-14岁(A1)，15-24岁(A2)，25-44岁(A3)，45-64岁(A4)，65岁及以上(A5)。相应的，流行病传播网络节点表示为由对应年龄组构成的复合群体，因此该传播网络包含5个网络节点。网络链接权重表示复合群体间的交互概率。D-AOC系统的参数如表1所示，其中包括不同年龄段的感染率(β_i)、潜伏期(γ_i)和恢复期(τ_i)以及不同年龄组之间的交互频率。基于这些参数运行D-AOC系统模拟H1N1在香港的传播过程，得到累计新增病例的近似期望值，图4蓝色实线所示。从中可见，D-AOC系统能够较好的模拟H1N1前100天感染人数的增长趋势，但由于数据噪声和人为干预的影响，模拟结果没有与真实数据完全一致。

表1

例2爆发趋势和模型参数估计

检验D-AOC系统的参数估计方法，具体步骤如下：①给定一组系统参数θ＝{A，θ₁，θ₂，θ₃，θ₄，θ₅}，其中θ_i＝<β_i，γ_i，τ_i>，1≤i≤5；②基于θ运行D-AOC系统模拟流行病的传播过程，记录此过程中感染病例数目随时间的变化情况，作为观测到的流行病监控数据D；③基于D，采用图3中的方法估计出参数

和流行病爆发趋势④对比实际值与估计值，分析参数和爆发趋势的估计情况。图5(a)所示，其中的条形图表示通过上述步骤②得到的流行病监控数据D(新增病例的时间序列)，红色实线表示通过上述步骤③估计出的爆发趋势(新增病例估计值的时间序列)。为度量估计出的爆发趋势

与监控数据D是否吻合，采用公式 $E(D,\hat{D})=1-\frac{D\&CenterDot;{\hat{D}}^T}{{\left|\right|D\left|\right|}_2\&CenterDot;{\left|\right|\hat{D}\left|\right|}_2},$ 其中 ${\left|\right|X\left|\right|}_2=\sqrt{{\mathrm{\&Sigma;}}_{i=1}^N{\mathrm{\&Sigma;}}_{t=1}^T{x}_{\mathrm{it}}^2}$ 表示矩阵X＝[x_it]_N×T的二阶范式，并且有

实验采用表1中参数作为待估计的参数θ。参照图5(b)给出了系统参数的推断过程，随着迭代步数的增加，当前参数下对应的爆发趋势逐步逼近观察值D，直至两者的偏差小于规定阈值。从图中可看出，通过负反馈机制推断出的D-AOC系统参数能够较好的模拟监控数据所描述的流行病爆发趋势(偏差为10^-5级)。表2给出了参数推断过程结束时得到的参数估计值

比较表1和表2可看出，图5中方法能够较为准确的估计出流行病传播网络结构及与流行病有关的生物学参数。

表2

例3缺失数据对爆发趋势估计的影响

实际可得的监控数据可能会有不同程度的缺失，缺失数据虽不可见，对流行病传播的影响却实际存在。基于缺失监控数据估计出流行病的爆发趋势后，可从估计数据中猜测出缺失的传播数据。本实验研究不同情况下的数据缺失对爆发趋势估计产生的影响，首先采用表1中参数运行D-AOC系统，模拟流行病的传播过程，采用不同时间尺度从中采样得到观测到的监控数据；基于采样所得的监控数据估计参数和爆发趋势；根据公式计算估计所得的爆发趋势与实际监控数据的偏差，进而分析不同情况的数据缺失对爆发趋势估计产生的影响。

实验中讨论的采样策略分别是：①对每个网络节点每天采样；②对每个网络节点以2天为间隔采样(即获取2天数据之和)；③对每个网络节点以7天为间隔采样(即获取7天数据之和)；④对5个网络节点分别以1天、2天、7天、15天和30天为间隔采样。第一种采样策略对应无数据缺失的情况，作为基准对比数据缺失情况。图6给出了不同采样策略对应的参数估计过程，从中可看出：①相对于无数据缺失的情况(对应图中红色曲线)，基于缺失数据估计出的爆发趋势与实际爆发趋势具有一定偏差；②采样的时间尺度越大，估计出的爆发趋势与实际爆发趋势的偏差越大，但都控制在10^-2以下，意味着即使监控数据具有不同程度缺失，D-AOC的参数推断过程也可较为准确的估计出流行病爆发趋势。

例4时间窗大小对爆发趋势估计的影响

流行病传播具有实时性，因而希望能快速且准确的估计出流行病的爆发趋势，从而尽早制定相应的防范措施。[1，T]表示获取流行病监控数据的时间窗，如果能以较小的时间窗估计出较为准确的爆发趋势，则意味着能够在爆发高峰之前准确的估计出流行病的爆发趋势。首先采用表1中参数运行D-AOC系统，模拟流行病传播过程，进而采用不同时间窗[1，T]从中采样得到监控数据；基于采样所得的监控数据估计[T，T+Δt]时期内的爆发趋势；根据公式

计算估计所得的爆发趋势与监控数据的偏差，进而分析不同T和Δt的取值对爆发趋势估计产生的影响。

图7中黑色实线表示流行病的爆发趋势，图中其它三条曲线分别表示Δt取值60天、90天和120天对应的估计偏差。从中可看出：①可获得的监控数据越多(时间窗越大)，估计出的传播趋势越准确；②预测的时间越短(Δt取值越小)，估计出的传播趋势越准确；③在传播初期(T＜114)，感染人口比例少，爆发趋势不明显，基于当前监控数据估计出的传播趋势具有较大偏差(大于0.3)；④而当传播过程刚刚进入开始明显增长的阶段(T＞134)，基于当前监控数据能够较好的估计出今后2到3个月的爆发趋势(偏差小于0.2)，这意味着本算法能够在爆发高峰到来之前较为准确的估计出今后一段时间内的流行病爆发趋势。⑤当传播过程离开高峰期后(T＞214)，监控数据累计的足够充分，无论Δt取值大小，都能够较好的估计出后期的爆发趋势。

例5流行病感染风险评估

流行病传播网络建模与推断的一个重要应用是流行病感染的风险评估，即基于当前的监控数据估计出不同地区的人口在后续时间内被感染的概率。图8和9以香港地区为例说明了基于传播网络推断评估感染风险的基本原理。①首先利用经纬线对全香港地区进行分区，形成层次网络的第一层，如图8(a)所示。②由于该分区粒度大，不能满足“同质”假设，为此需要对第一层中的每个分区再次划分得到层次网络的第二层。选择香港岛，将其进一步细分为12×12个小区域，使得每个小区域近似满足“同质”假设，如图8(b)所示。图中每个区域中心的实心点表示布设的监控点(医院或诊所)，提供实时监控数据(图中表示为相应的时间序列)。③融合各监控点的监控数据推断出传播网络结构(如图8(c)所示)和模型参数，进而执行系统模拟一定时间内(如2个月内)流行病在香港岛的传播过程。④利用模拟数据和公式

计算出各地区人口在不同时刻被感染的风险

并在地图的相应区域用不同深浅的颜色表示感染风险的大小，如参照图9(c)所示。图8(b)中，对香港岛细分为144个小区域中包含85个有效区域，因此图8(c)所示的传播网络包含85个节点，网络链接权重表示自治体从一个区域移动到其它区域的概率。图9(a)给出了该示例中流行病爆发至40天时累计病例的空间分布，其中点表示病例。基于前40天的监控数据，采用上述过程计算出未来60天各区域中个体被感染的风险，绘制出感染风险的时空分布图。图9(b)给出了未来第100天时香港岛各地区的感染风险，左下角的颜色条表示不同灰度对应的感染风险数值。图9(c)给出了未来60天内香港岛各地区的感染风险随时间演化的趋势。

Claims (7)

1.一种基于自治计算的流行病传播网络建模与推断方法，其特征在于，包括如下步骤：

101：构造模拟流行病传播过程的自治体系统(命名为D-AOC系统)；

102：采用蒙特卡罗模拟和负反馈机制，从流行病监控数据中推断出D-AOC系统的参数，包括：流行病传播网络和与流行病有关的生物学参数。

2.根据权利要求书1所述的D-AOC系统，其特征在于，按照如下步骤构造：

201：建模流行病传播网络；

202：建模自治体；

3.根据权利要求书2所述的流行病传播网络，其特征在于，按照如下方式建模：

流行病传播网络包含N个节点。节点可表示具有不同大小的空间区域(如户、村、乡、镇、县、市)，也可表示复合群体(如不同年龄段的人口)。节点间的链接表示流行病传播途径，链接权重表示传播强度。流行病传播网络表示为一个N×N矩阵A＝[a_ij]_N×N，a_ij表示节点i与节点j的链接权重，满足0≤a_ij≤1和

4.根据权利要求书2所述的自治体，其特征在于按照如下方法定义自治体和自治体的状态转换规则：

401：定义自治体

流行病传播网络中的每个节点可包含多个自治体，节点i中自治体定义为如下的一个5元组：

(i，s，β_i，γ_i，τ_i)

其中，i表示自治体所在网络节点的标识符；s表示自治体当前所处的状态，自治体在每个时刻处于S(易感)、E(潜伏)、I(感染)和R(恢复)四种状态之一；β_i，1/γ_i，1/τ_i分别表示节点i中自治体对流行病的抵抗力、流行病的潜伏期和患病之后的康复期，满足0≤β_i，γ_i，τ_i≤1。令B＝(β₁，…，β_N)，H＝(γ₁，…，γ_N)，F＝(τ₁，…，τ_N)，向量B，H和F表示与流行病相关的生物学参数。

402：定义自治体的状态转换规则

规则1：当位于节点i中的自治体处于S状态时，将以

的概率进入E状态。

称为节点i在t时刻的感染强度，由自治体对流行病的抵抗力和流行病传播网络共同决定，计算如下：

${\mathrm{\&lambda;}}_i^{\left(t\right)}={\mathrm{\&beta;}}_i{\mathrm{\&Sigma;}}_{j=1}^N\frac{I_j^{\left(t\right)}}{P_j^{\left(t\right)}}{a}_{\mathrm{ij}}$

其中β_i表示节点i中自治体对流行病的抵抗力，

表示t时刻节点j中处于I状态的自治体数目，

表示t时刻节点j中自治体的总数，a_ij表示节点i和节点j之间的网络链接权重。

规则2：自治体处于E状态已经持续1/γ_i时间后进入I状态；

规则3：自治体处于I状态已经持续1/τ_i时间后进入R状态；

规则4：自治体进入R状态后则永久停留在该状态，意味着它获得了对流行病的终身免疫。

5.根据权利要求书1所述的系统参数推断，其特征在于，按如下步骤进行：

令向量<A，B，H，F>表示待推断的系统参数，各分量含义如上所述。

令D表示实际观测到的流行病监控数据，定义为一个N×T矩阵D＝[d_it]_N×T，d_it表示t时刻节点i中新出现的流行病患者数目，T表示监控时间长度。基于D，按照如下方式推断出D-AOC的系统参数：

501：随机生成2m组系统参数，令P_i表示第i组系统参数，V_i表示与P_i对应的一组系统参数，其中1≤i≤m；

502：令迭代次数Q等于0；

503：令LP和GP分别表示局部最优和全局最优的系统参数；

初始时，令LP和GP均为P₁；初始时，令E_max为一个充分大的整数。

504：对i从1到m执行如下操作：

基于P_i，运行D-AOC系统的蒙特卡洛模拟，获取模拟过程产生的流行病监控数据

其中

表示t时刻节点j中新出现的处于I状态的自治体数目，T表示模拟过程持续的时间；

计算

偏离D的程度 $E_i={\left|\right|D-{\hat{D}}_i\left|\right|}_2,$ 其中 ${\left|\right|X\left|\right|}_2=\sqrt{{\mathrm{\&Sigma;}}_{i=1}^N{\mathrm{\&Sigma;}}_{t=1}^T{x}_{\mathrm{it}}^2}$ 表示矩阵X＝[x_it]_N×T的二阶范式；

505：根据以上计算出的m个偏差E_i(1≤i≤m)，更新LP和GP，具体步骤如下：

令E_u是m个偏差中的最小者；

将LP更新为P_u；

如果E_u小于E_max，则将GP更新为P_u，同时将E_max更新为E_u；

506：对i从1到m执行如下操作：

根据LP和GP采用负反馈机制更新P_i和V_i；

507：令迭代次数Q增1。如果Q小于规定的迭代次数则返回步骤504执行；

508：返回GP作为从D推断出的系统参数。

6.根据权利要求书5所述的D-AOC系统的蒙特卡洛模拟，其特征在于，按如下方程进行：

$\left\{\begin{array}{c}{\mathrm{\&Delta;S}}_i^{(t+1)}=-{\mathrm{\&lambda;}}_i^{\left(t\right)}{S}_i^{\left(t\right)}\\ {\mathrm{\&Delta;E}}_i^{(t+1)}={\mathrm{\&lambda;}}_i^{\left(t\right)}{S}_i^{\left(t\right)}-{\mathrm{\&gamma;}}_i{E}_i^{\left(t\right)}\\ {\mathrm{\&Delta;I}}_i^{(t+1)}={\mathrm{\&gamma;}}_i{E}_i^{\left(t\right)}-{\mathrm{\&tau;}}_i{I}_i^{\left(t\right)}\\ {\mathrm{\&Delta;R}}_i^{(t+1)}={\mathrm{\&tau;}}_i{I}_i^{\left(t\right)}\end{array}\right.$

其中，

γ_i和τ_i含义如上所述；

分别表示t时刻节点i中处于S状态、E状态和I状态的自治体数目；

和

分别表示t+1时刻节点i中处于S状态、E状态、I状态和R状态的自治体数目的变化量。特别的，

表示节点i中新出现的处于I状态的自治体数目。

7.根据权利要求书5所述的根据LP和GP采用负反馈机制更新P_i和V_i，其特征在于，按如下步骤进行：

701：将V_i更新为wV_i+c₁·rand·(LP-P_i)+c₂·rand·(GP-P_i)；

702：将P_i更新为P_i+V_i；

其中，rand表示随机数发生器，返回0到1范围内的随机数；w，c₁和c₂是3个预先设置的参数。

Images