疾病流行的模拟
作者：刘向明 林家瑞
单位：刘向明 林家瑞(湖北职工医学院 荆州434000)
关键词：流行病学；模拟；病毒
990301
摘 要 基于在流行病学数学模型中具有基本重要性的Reed-Frost模型，对疾病流行的模拟进行评价，首先介绍Reed-Frost模型及其实物模拟，然后论述模型的计算机模拟方法，最后将Reed-Frost模型的模拟方法推广到用于存在两种病毒的情形下疾病流行的模拟，使流行病学工作者及有关人员对疾病流行的模拟有一个简明扼要的了解，并可以借此利用计算机开始流行病学的模拟研究工作，也为一般生物医学数学模型研究人员在他们各自的领域里进行模拟研究提供参考与借鉴。
1 Reed-Frost模型与实物模拟
, 百拇医药
N.J.Reed和W.H.Frost早在1928年就在流行病学研究中引入了易感者与感染者交互作用模型。为此，首先选定适当的时间单位，这通常与疾病感染期的平均长度有关。然后，可令Sn表示第n个单位时间内易感者人数，Cn表示第n个单位时间内感染者人数，In表示第n个单位时间内免疫者人数，即到第n个单位时间已经历过感染阶段的人数，它们都是随机变量。在每一时刻个体可处在易感、感染和免疫这三种状态之一。易感者可随机地转变为感染者，但感染者则一定转变为免疫者，我们称这种转变为受到约束的转移。
假定人群大小不变，其数目为N，且是随机混匀的，即对于整个人群而言，任何易感者与任何感染者之间有效接触的概率相同。设P表示在第n个单位时间中易感者和感染者的有效接触率，在每个单位时间的终止时刻，考察人群并确定感染者、易感者和免疫者的新数目。在第n个单位时间中易感者和所有感染者均不接触的概率为(1-P)Cn，因此，
(1)
, 百拇医药
即Sn+1=k和Cn+1=m的概率为从Sn中选取k个个体的方式数乘上k个个体不与任一感染者接触的概率和m个个体均与感染者接触的概率。由于人群总数保持不变，In+1=N-K-m。
以上就是Reed-Feost模型的随机模式，因为
(2)
这启发我们把方程组
(3)
作为Reed-Frost模型的确定式模式，这里，S'n、C'n表示确定的人数。然而，不能期望有
E(Sn)=S'n
, 百拇医药
E(Cn)=C'n (4)
令a=-log(1-p),可将以上方程组重新写为
(5)
这就是所谓Kermack-Mckendrick模型。Reed-Frost模型确定了Sn和Cn为二项马尔可夫链。其中随时间变化的转移概率依赖于流行的状态，这个链难于分析，所以要借助于模拟进行研究。
Reed和Frost为了进行教学演示首先对模型作了实物模拟，他们用白、红、绿三种颜色的球表示所考察的人群。其中白球代表易感者，红球代表感染者，绿球代表免疫者。在各个单位时间里，把这些球与一定数量的黑球混匀后滚入一机械斜槽，将黑球想象为建筑物中房室的间隔或社区的边界，两个黑球间的部分为建筑物中房室和社区中的区块。任一易感者(白球)如和感染者(红球)在同一房室中相遇，则成为下一个单位时间的感染者；如白球所在房室里没有红球，则继续代表易感者；原来的红球在下一个单位时间里成为免疫者，用绿球代替；代表那些被感染了的易感者的白球用红球代替后加入到总体混匀。重复这个过程，直到不再出现白球或者红球。从上述可见，对于各个单位时间，可确定新的感染者和易感者数目，通过计算也能确定免疫者数目，并可列为表格形式，逐个区间的计数和累积结果。重复多次，就能确定流行规模的分布，观察流行的起伏。
, 百拇医药
为了解参数变化，譬如说降低有效接触率，所引起的流行过程的改变，只要在模拟时增加黑球的数目就可以了。其他一些因素的效应，例如感染者可恢复为易感者，也可类似的考察。这种模拟方法在不需要较多的数学与统计学的知识前提下成功地显示了简单模型所支配的疾病流行过程。但是，这种模拟方法是非常耗时的，缺乏实际应用价值。
2 计算机模拟
基于数学模型的模拟实际上是一个计算过程，在计算机技术不甚发达时模拟的实际可行性差。随着电子计算机技术的发展，数学模型方法在生物医学中获得了广泛而深刻的应用，基于数学模型的模拟可通过计算机技术实现，而不再依赖于实物。以下来说明基于Reed-Frost模型的疾病流行的模拟过程。
Reed-Frost模型中的Sn、Cn和In分别确定了第n个单位时间中的易感者、感染者和免疫者的数目，由于人群大小不变，故有N=Sn+Cn+In而Sn和Cn的条件密度函数由公式
(6)
, 百拇医药
递归地确定。这个模型可用三种方法分析研究：
①直接计算联合概率密度函数Pr[Sn=k,Cn=m];
②作模型的计算机模拟；
③采用各种近似方案。
虽然直接计算联合密度函数可以做到，但使人感到沉闷且结果受到限制，而计算机模拟却能发现引人注目的结果。
Reed-Frost链的计算机模拟是通过计算机生成一些可能的流行演化过程，称这些演化过程为样本路径。然后计算这些样本路径的多种统计量，用以描述疾病流行的动力学。具体方法如下：
首先确定初始数据，比方说S0=40,C0=1,I0=0
, 百拇医药
由Reed-Frost模型有
(7)
其中，Q=1-P,计算机按程序指令生成(0，1)区间上的均匀分布随机数X1,然后计算诸数值Q1,40,Q1,39,…,Q1,j1直到Q1,40+Q1,39+…+Q1,j1>X1,于是，取S1=j1,I1=41=j1;再生成另一随机数X2,按公式
(8)
计算诸数值Q2,j1,Q2,j1-1,…,Q2,j2,使Q2,j1+…+Q2,j2>X2,于是S2=j2,C2=j1-j2。将计算过程继续到得到某一代的易感者数目为止。
, 百拇医药
利用二项分布可得到生成样本路径的一个较快速的方法。在每个时刻n,生成jn-1个随机数ρ1,…,ρjn-1
令
(9)
则
重复以上过程多次，然后应用统计学方法分析处理这批数据，获得需要的信息。表1描述了易感者人群最终数量分布。
表1 易感者人群最终数量分布直方图数据表 j100
0
, 百拇医药 1
2
3…
34
35
36
37
38
39
40
相应的样本
路径数
4
2
, 百拇医药
1
0…
0
1
3
10
19
8
2
3 两种病毒情形下流行过程的模拟
在流行病研究中，一些问题涉及到相似肠道病毒之间的竞争干预，如polio病毒、coxsackie病毒和ECHO病毒。这些病毒争夺宿主体内相同的繁殖部位。于是，被一种病原体A感染的个体对另一种病毒B有暂时的干扰性。
, http://www.100md.com
现场观察提示但不能严格证实种痘或感染某种病毒可能抑制致病力更强的病原微生物所引起的流行的出现。当然不可能直接以人群作实验来对这种效应进行观察研究，也不存在适合的非人类模型，但是可进一步利用模拟研究这个问题。为了反映由此引发的各种公共卫生问题，有必要增加模型的复杂性。然而可利用对Reed-Frpst模型稍作推广得到较简单的模型，来获得观察干预的可能性的一些信息。
假定有两种病毒，例如coxsackie和polio，分别用A与B来表示。对于一种病毒引起的疾病，存在四种可能的状态，即：易感状态S，感染状态C，免疫状态I，以及受到一种病毒的感染后的对另一种病毒的非易感状态T。由于存在两种病毒，这里用两个字母表示一个别的状态，例如SS表示双重易感者，要注意的是不是任何字母都表示可能的状态。例如T仅仅只在TC和CT中出现才有实际意义，其他可能的状态是SI、IS、CI、IC和II。这样，在存在两种病毒的前提下，原来Reed=Frost模型中三个状态增加到八个状态，原来只有一个从感染状态到免疫状态的有约束的转移，现在有了四个有约束的转移，随机转移也从一个增加到三个，见表2。
, 百拇医药
表2 两种病毒情形下的状态转移表 有约束的转移
随机转移
不可能向其他状
态转移的状态
TC→SI
SI→(SI或CI)
CT→IS
IS→(IS或IC)
II
CI→II
SS→(SS或CI或IC)
IC→II
, http://www.100md.com
表2中前两种随机转移与Reed-Frost模型中随机转移情形类似，在模拟中可以相同的方式处理。第三种随机转移涉及到要作出两种决策。首先，要决定个体是否能避免与所有两种病毒的感染者接触。用QA表示两种病毒的易感者不与任一A病毒感染者接触的概率，QB有类似的解释。那么，这个易感者不和这两种病毒的任一感染者接触的概率为：QT=QA.QB
现在可使用(0，1)中的均匀随机数X作出决策，如果X≤QT不成立，则就要进一步决定个体是成为A病毒感染者还是B病毒感染者。如果成功接触的概率PA和PB相同，那么，受A病毒感染的分数由公式
fA=nc,i,A/(nc,i,A+nc,i,B) (10)
, 百拇医药
给出是合理的。其中nc,i,A和nc,i,B分别是第i个单位时间中A病毒与B病毒的感染者。然而，如PA和PB不同，则应使用其他的竞争风险模型，如有人选择的模型导致
fA=nc,i,A.log(QA)/[nc,i,Alog(QA)+nc,i,Blog(QB)] (11)
这里QA=1-PA,QB=1-PB。无论选择什么样的模型，都要求
fA+fB=1 QT+fA(1-QT)
成立，则出现转移SS→TC。
Reed-Frost模型及其简单推广的计算模拟方法是对更复杂模型进行计算机模拟的基础，人们可籍此编写计算机程序，实施模拟。在熟悉了基本模型的模拟方法后，进而可以设计处理更加复杂的实际问题的研究方案，进行模拟。
参考文献
1 方积乾.微积分初步与生物医学应用.北京：北京医科大学与中国协和医科大学出版社，1990.
收稿日期：1998-10-02, 百拇医药
单位：刘向明 林家瑞(湖北职工医学院 荆州434000)
关键词：流行病学；模拟；病毒
990301
摘 要 基于在流行病学数学模型中具有基本重要性的Reed-Frost模型，对疾病流行的模拟进行评价，首先介绍Reed-Frost模型及其实物模拟，然后论述模型的计算机模拟方法，最后将Reed-Frost模型的模拟方法推广到用于存在两种病毒的情形下疾病流行的模拟，使流行病学工作者及有关人员对疾病流行的模拟有一个简明扼要的了解，并可以借此利用计算机开始流行病学的模拟研究工作，也为一般生物医学数学模型研究人员在他们各自的领域里进行模拟研究提供参考与借鉴。
1 Reed-Frost模型与实物模拟
, 百拇医药
N.J.Reed和W.H.Frost早在1928年就在流行病学研究中引入了易感者与感染者交互作用模型。为此，首先选定适当的时间单位，这通常与疾病感染期的平均长度有关。然后，可令Sn表示第n个单位时间内易感者人数，Cn表示第n个单位时间内感染者人数，In表示第n个单位时间内免疫者人数，即到第n个单位时间已经历过感染阶段的人数，它们都是随机变量。在每一时刻个体可处在易感、感染和免疫这三种状态之一。易感者可随机地转变为感染者，但感染者则一定转变为免疫者，我们称这种转变为受到约束的转移。
假定人群大小不变，其数目为N，且是随机混匀的，即对于整个人群而言，任何易感者与任何感染者之间有效接触的概率相同。设P表示在第n个单位时间中易感者和感染者的有效接触率，在每个单位时间的终止时刻，考察人群并确定感染者、易感者和免疫者的新数目。在第n个单位时间中易感者和所有感染者均不接触的概率为(1-P)Cn，因此，
(1), 百拇医药
即Sn+1=k和Cn+1=m的概率为从Sn中选取k个个体的方式数乘上k个个体不与任一感染者接触的概率和m个个体均与感染者接触的概率。由于人群总数保持不变，In+1=N-K-m。
以上就是Reed-Feost模型的随机模式，因为
(2)这启发我们把方程组
(3)作为Reed-Frost模型的确定式模式，这里，S'n、C'n表示确定的人数。然而，不能期望有
E(Sn)=S'n
, 百拇医药
E(Cn)=C'n (4)
令a=-log(1-p),可将以上方程组重新写为
(5)这就是所谓Kermack-Mckendrick模型。Reed-Frost模型确定了Sn和Cn为二项马尔可夫链。其中随时间变化的转移概率依赖于流行的状态，这个链难于分析，所以要借助于模拟进行研究。
Reed和Frost为了进行教学演示首先对模型作了实物模拟，他们用白、红、绿三种颜色的球表示所考察的人群。其中白球代表易感者，红球代表感染者，绿球代表免疫者。在各个单位时间里，把这些球与一定数量的黑球混匀后滚入一机械斜槽，将黑球想象为建筑物中房室的间隔或社区的边界，两个黑球间的部分为建筑物中房室和社区中的区块。任一易感者(白球)如和感染者(红球)在同一房室中相遇，则成为下一个单位时间的感染者；如白球所在房室里没有红球，则继续代表易感者；原来的红球在下一个单位时间里成为免疫者，用绿球代替；代表那些被感染了的易感者的白球用红球代替后加入到总体混匀。重复这个过程，直到不再出现白球或者红球。从上述可见，对于各个单位时间，可确定新的感染者和易感者数目，通过计算也能确定免疫者数目，并可列为表格形式，逐个区间的计数和累积结果。重复多次，就能确定流行规模的分布，观察流行的起伏。
, 百拇医药
为了解参数变化，譬如说降低有效接触率，所引起的流行过程的改变，只要在模拟时增加黑球的数目就可以了。其他一些因素的效应，例如感染者可恢复为易感者，也可类似的考察。这种模拟方法在不需要较多的数学与统计学的知识前提下成功地显示了简单模型所支配的疾病流行过程。但是，这种模拟方法是非常耗时的，缺乏实际应用价值。
2 计算机模拟
基于数学模型的模拟实际上是一个计算过程，在计算机技术不甚发达时模拟的实际可行性差。随着电子计算机技术的发展，数学模型方法在生物医学中获得了广泛而深刻的应用，基于数学模型的模拟可通过计算机技术实现，而不再依赖于实物。以下来说明基于Reed-Frost模型的疾病流行的模拟过程。
Reed-Frost模型中的Sn、Cn和In分别确定了第n个单位时间中的易感者、感染者和免疫者的数目，由于人群大小不变，故有N=Sn+Cn+In而Sn和Cn的条件密度函数由公式
(6), 百拇医药
递归地确定。这个模型可用三种方法分析研究：
①直接计算联合概率密度函数Pr[Sn=k,Cn=m];
②作模型的计算机模拟；
③采用各种近似方案。
虽然直接计算联合密度函数可以做到，但使人感到沉闷且结果受到限制，而计算机模拟却能发现引人注目的结果。
Reed-Frost链的计算机模拟是通过计算机生成一些可能的流行演化过程，称这些演化过程为样本路径。然后计算这些样本路径的多种统计量，用以描述疾病流行的动力学。具体方法如下：
首先确定初始数据，比方说S0=40,C0=1,I0=0
, 百拇医药
由Reed-Frost模型有
(7)其中，Q=1-P,计算机按程序指令生成(0，1)区间上的均匀分布随机数X1,然后计算诸数值Q1,40,Q1,39,…,Q1,j1直到Q1,40+Q1,39+…+Q1,j1>X1,于是，取S1=j1,I1=41=j1;再生成另一随机数X2,按公式
(8)计算诸数值Q2,j1,Q2,j1-1,…,Q2,j2,使Q2,j1+…+Q2,j2>X2,于是S2=j2,C2=j1-j2。将计算过程继续到得到某一代的易感者数目为止。
, 百拇医药
利用二项分布可得到生成样本路径的一个较快速的方法。在每个时刻n,生成jn-1个随机数ρ1,…,ρjn-1
令
(9)则
重复以上过程多次，然后应用统计学方法分析处理这批数据，获得需要的信息。表1描述了易感者人群最终数量分布。
表1 易感者人群最终数量分布直方图数据表 j100
0
, 百拇医药 1
2
3…
34
35
36
37
38
39
40
相应的样本
路径数
4
2
, 百拇医药
1
0…
0
1
3
10
19
8
2
3 两种病毒情形下流行过程的模拟
在流行病研究中，一些问题涉及到相似肠道病毒之间的竞争干预，如polio病毒、coxsackie病毒和ECHO病毒。这些病毒争夺宿主体内相同的繁殖部位。于是，被一种病原体A感染的个体对另一种病毒B有暂时的干扰性。
, http://www.100md.com
现场观察提示但不能严格证实种痘或感染某种病毒可能抑制致病力更强的病原微生物所引起的流行的出现。当然不可能直接以人群作实验来对这种效应进行观察研究，也不存在适合的非人类模型，但是可进一步利用模拟研究这个问题。为了反映由此引发的各种公共卫生问题，有必要增加模型的复杂性。然而可利用对Reed-Frpst模型稍作推广得到较简单的模型，来获得观察干预的可能性的一些信息。
假定有两种病毒，例如coxsackie和polio，分别用A与B来表示。对于一种病毒引起的疾病，存在四种可能的状态，即：易感状态S，感染状态C，免疫状态I，以及受到一种病毒的感染后的对另一种病毒的非易感状态T。由于存在两种病毒，这里用两个字母表示一个别的状态，例如SS表示双重易感者，要注意的是不是任何字母都表示可能的状态。例如T仅仅只在TC和CT中出现才有实际意义，其他可能的状态是SI、IS、CI、IC和II。这样，在存在两种病毒的前提下，原来Reed=Frost模型中三个状态增加到八个状态，原来只有一个从感染状态到免疫状态的有约束的转移，现在有了四个有约束的转移，随机转移也从一个增加到三个，见表2。
, 百拇医药
表2 两种病毒情形下的状态转移表 有约束的转移
随机转移
不可能向其他状
态转移的状态
TC→SI
SI→(SI或CI)
CT→IS
IS→(IS或IC)
II
CI→II
SS→(SS或CI或IC)
IC→II
, http://www.100md.com
表2中前两种随机转移与Reed-Frost模型中随机转移情形类似，在模拟中可以相同的方式处理。第三种随机转移涉及到要作出两种决策。首先，要决定个体是否能避免与所有两种病毒的感染者接触。用QA表示两种病毒的易感者不与任一A病毒感染者接触的概率，QB有类似的解释。那么，这个易感者不和这两种病毒的任一感染者接触的概率为：QT=QA.QB
现在可使用(0，1)中的均匀随机数X作出决策，如果X≤QT不成立，则就要进一步决定个体是成为A病毒感染者还是B病毒感染者。如果成功接触的概率PA和PB相同，那么，受A病毒感染的分数由公式
fA=nc,i,A/(nc,i,A+nc,i,B) (10)
, 百拇医药
给出是合理的。其中nc,i,A和nc,i,B分别是第i个单位时间中A病毒与B病毒的感染者。然而，如PA和PB不同，则应使用其他的竞争风险模型，如有人选择的模型导致
fA=nc,i,A.log(QA)/[nc,i,Alog(QA)+nc,i,Blog(QB)] (11)
这里QA=1-PA,QB=1-PB。无论选择什么样的模型，都要求
fA+fB=1 QT+fA(1-QT)
成立，则出现转移SS→TC。
Reed-Frost模型及其简单推广的计算模拟方法是对更复杂模型进行计算机模拟的基础，人们可籍此编写计算机程序，实施模拟。在熟悉了基本模型的模拟方法后，进而可以设计处理更加复杂的实际问题的研究方案，进行模拟。
参考文献
1 方积乾.微积分初步与生物医学应用.北京：北京医科大学与中国协和医科大学出版社，1990.
收稿日期：1998-10-02, 百拇医药