跳转至

第 13 章:流行病的地理传播与控制(Geographic Spread and Control of Epidemics)

本章论述流行病的地理空间传播问题。作者开门见山地指出,与时间维度上对疾病/流行病的发生发展与控制已被充分研究相比,地理空间维度的传播还远未被理解和研究透。然而,对传染病、药物滥用、谣言或错误信息等时空传播过程建立逼真模型是非常必要的。核心问题是如何在模型中纳入并量化空间效应。本章按"简单空间流行病模型 → 历史案例(1347–1350 年欧洲黑死病)→ 当代现实案例(欧洲狐狸中的狂犬病)"的顺序展开,最后讨论控制策略,并强调这些模型类型并不局限于单一疾病。

13.1 流行病空间传播的简单模型(Simple Model for the Spatial Spread of an Epidemic)

本节给出卷 I 第 10.2 节非空间流行病模型的简化空间版。问题是在均匀易感人群 S₀ 中引入一些感染者,确定疫情的时空扩散。人口只分两类:位置 x、时刻 t 的感染者 I(x, t) 与易感者 S(x, t)。两类都用相同扩散系数 D 进行简单扩散(即 S、I 都游走),从易感到感染的转化率正比于 rSI(r 为常数参数,反映疾病的传播效率),感染者因病死亡率为 aI,1/a 是感染者寿命。模型 (13.1) 是 ∂S/∂t = −rIS + D∇²S、∂I/∂t = rIS − aI + D∇²I。这正是卷 I 第 10 章 (10.1)(10.2) 加扩散项。

作者集中于一维问题;二维研究在 13.8 节给出。引入 I = I/S₀、S = S/S₀、x = √(rS₀/D) x、t = rS₀ t、λ = a/(rS₀)(S₀ 为代表密度),把方程化为 (13.3):∂S/∂t = −IS + ∂²S/∂x²、∂I/∂t = I(S − λ) + ∂²I/∂x²。三个维度参数 r、a、D 缩减为唯一无量纲组合 λ。基本传染率 R₀ = 1/λ(参卷 I 10.2 节)有几种等价含义:(i) 一名原发感染者在完全易感人群中产生的二代感染数;(ii) 一个时间尺度的度量——传染期时间尺度 1/(rS₀) 与感染者寿命 1/a 的相对值。

具体研究的是均匀易感人群中的传染病行波传播:寻找存在行波解的条件和波速。设 I(x, t) = I(z)、S(x, t) = S(z)、z = x − ct 代入 (13.3) 得到常微分系统 (13.5):I″ + cI′ + I(S − λ) = 0、S″ + cS′ − IS = 0,c 为待定波速。边界条件 I(−∞) = I(∞) = 0、0 ≤ S(−∞) < S(∞) = 1 描述一个感染者脉冲向未感染人群传播。图 13.1 给出由 (13.3)(初值紧支集、与 (13.6) 相容)计算出的 λ = 0.75 行波;图 13.5 是另一例(仅感染者扩散)。(13.5) 是 4 阶相空间系统。

S(z) 必须是 z 的单调递增函数(因为若在某 z 处 S 有局部极大,则 S′ = 0、(13.5) 第二式给出 S″ = IS > 0,对应极小,矛盾)。在波前 S → 1、I → 0 附近对 (13.5) 第一个方程线性化得到 I″ + cI′ + (1 − λ)I ≈ 0(参卷 I 第 13 章 13.2 节中对 Fisher–Kolmogoroff 方程行波解的处理),解为 I(z) ∝ exp[((−c ± √{c² − 4(1 − λ)}) z/2)]。要求 I(z) > 0 且 → 0 故不能振荡,由此得波速 c 与参数 λ 必须满足 c ≥ 2√(1 − λ)、λ < 1,即 (13.9)。λ > 1 时无波解,故 λ = 1 是疫情波传播的阈值条件。维度形式 a/(rS₀) < 1((13.10)),恰好与卷 I 10.2 节在无空间情况下得到的存在流行病的阈值条件相同。

仿照 Fisher–Kolmogoroff 方程的经验,作者预期由完整非线性系统算出的行波几乎总是演化为最小波速 c = 2√(1 − λ) 的波(除少数特殊情形)。维度形式波速 V = 2√(rS₀D)·√(1 − a/(rS₀))((13.11))。S(z) → 1 的衰减也是指数级:把 (13.5) 第二式在 z → ∞ 处 S = 1 − s、s 小时线性化得 s″ + cs′ − I = 0;与 (13.8) 一起得 S(z) ∼ 1 − O(exp[−(c ∓ √{c² − 4(1 − λ)})z/2]),即指数趋近 1。

阈值 (13.10) 蕴含几条重要控制推论:(a) 存在最小临界易感密度 Sc = a/r,低于此则不能形成传染病波;(b) 给定 S₀ 与 a,存在临界传播系数 rc = a/S₀,超过它才能传播;(c) 给定 r 与 S₀,存在临界死亡率 ac = rS₀,超过它则抑制传播——疾病越快致死,流行越难扩散。这些对控制策略都有启示。易感者种群可通过疫苗接种或扑杀减少(13.6–13.9 节将详细讨论免疫与扑杀的效果)。给定死亡率与种群密度 S₀,若能通过隔离、医疗干预等降低传播系数 r,可能破坏 (13.10) 从而阻止疫情传播。对阈值 a/(rS₀) < 1 还应注意:易感人口的突然涌入把 S₀ 抬过 Sc 也可触发新疫情。

本节末尾指出,二物种模型可自然推广为三物种 SIR 模型。13.5–13.9 节将把这一推广用于当前欧洲狂犬病疫情的具体研究。

13.2 1347–1350 年黑死病在欧洲的传播(Spread of the Black Death in Europe 1347–1350)

14 世纪中叶的黑死病是鼠疫大流行,作者引用 Langer (1964) 的统计与 McEvedy (1988) 的史评。鼠疫(bubonic plague)由鼠疫杆菌(Bacillus pestis)引起,主要通过黑鼠身上的跳蚤传染给人;感染后 2–3 天内约 80% 死亡,总共死掉欧洲四分之一到三分之一的人口。1347 年 12 月鼠疫由东方船只带入意大利,随后几年内以约 200–400 英里/年的速度向北穿过欧洲,图 13.2 给出其时空传播的波前。

1350 年主疫情之后,1356 年德国爆发第二次鼠疫;之后周期性的小型疫情持续出现,但规模再不及 1347 年那一次。13.4 节会给出"部分免疫"扩展,使主波后出现周期性小幅复发,对应图 13.6、13.7。13.5 节的三物种狂犬病模型(图 13.9、13.10)则展示更显著的周期性流行波——那里的复发周期可以解析地估计。

中世纪欧洲对鼠疫的反应各式各样:成群结队的"鞭笞者"半裸苦行、宣讲末世论;富人的精致象牙柄皮鞭有遗存。15 世纪末威斯特法伦(Westfalen-Lippe)发现一份"治愈"方子:将快孵蛋切去尖端、放走胚胎、混生藏红花、煎至深褐,再加白芥末、莳萝、鹤嘴与 theriak(当时的江湖万灵药),吞后 7 小时禁食——疗效"未见记录"。

鼠疫有腺型、肺型、败血型三种。败血型为杆菌在血液中极快速繁殖,即使现代也几乎必死,肺型为空气传播极易传染。肺型死者有的"坐下即倒地身亡"(败血型)。伦敦 1664–1666 年大瘟疫期间流行英国童谣 "Ring-a ring o'roses..."(玫瑰花环/口袋里花束/喷嚏/我们全倒地")被认为创作于该时期;当时用洋葱、大蒜塞鼻或"鸟嘴面具"阻挡被认为致病的"恶臭"。

作者比较了 1665 年伦敦大瘟疫(数据比黑死病多得多)。Samuel Pepys 描述牛津是"放荡不羁行为",Daniel Defoe 1722 年日记生动写照当时情景。鼠疫并未因伦敦大瘟疫而绝迹:Gregg (1985) 详述了 1850 年始于中国云南、1959 年才被 WHO 正式宣告结束的最近一次鼠疫大流行,导致 1300 多万人死亡;1959 年后报告病例仍持续。越战期间(1965–1975)成千上万人死亡。旧金山 1906 年地震后的三年疫情中约 200 人死亡,造成美国西部(特别是新墨西哥)至今是世界上两个最大残余鼠疫灶区之一(另一个在俄罗斯)。鼠疫杆菌自西海岸稳步东移,1984 年已在中西部发现,传播速率约 35 英里/年。鼠疫存在于近 30 种哺乳动物中(松鼠、花栗鼠、郊狼、土拨鼠、小鼠、田鼠、家养宠物、蝙蝠等),并非只有鼠。纽约估计每人一鼠的鼠密度让东海岸爆发风险巨大。鼠疫症状常被误诊或迟诊,肺型极易感染。

回到建模,作者引用 Noble (1974) 的工作估计参数。1347 年欧洲约 8500 万人,平均密度 S₀ ≈ 50/英里²。设瘟疫与谣言都按扩散传播,扩散距离 L 的时间 O(L²/D);按当时通讯条件新闻约 100 英里/年得 D ≈ 10⁴ 英里²/年。传播系数 r 由跳蚤-鼠-人接触决定,Noble 估 r ≈ 0.4 英里²/年;平均感染期约两周(可能偏长),死亡率 a ≈ 15/年。代入得 λ = a/(rS₀) ≈ 0.75;由 (13.11) 得波速 V ≈ 140 英里/年,与 Langer (1964) 引用的 200–400 英里/年相比偏低,但与粗略的参数估计范围相符。

模型没有纳入人口密度不均、随机性等因素,但能把握地理传播的某些全局特征。Raggett (1982) 对英国 Eyam 村 1665–66 年疫情的随机模型给出的对照并不优于确定性模型。但在涉及小数目的空间模型里,随机性更关键。Keeling 与 Gilligan (2000) 提出一个有意思的鼠疫时空模型,鼠疫是人畜共患病,他们把鼠种群与人种群同时纳入(含随机性),用元胞自动机处理空间随机性;模型显示疾病可在鼠亚群中持续多年,并给出鼠种群中疾病向人群传播的判据,数据取自北美的啮齿动物种群与现行参数估计。

13.3 狂犬病简史:事实与神话(Brief History of Rabies: Facts and Myths)

狂犬病或为最可怕的疾病:患者在漫长而恐怖的病程中经历最可怕的噩梦式体验最终死去。希波的圣奥古斯丁把它列入灾祸清单(含疯狂、监禁、破产)。尽管已有有效疫苗(接触后立即接种几乎百分百有效),但若已发病(出现临床症状)则无药可治。Patricia Morrison 博士(私人通讯,1992)讲述 8 世纪列日主教圣于贝尔(St. Hubert)的故事:于贝尔死后百年尸体被掘出(据说未腐),运至阿登山区一座贫困修道院;修道院急需圣物吸引朝圣者,于是于贝尔成为"圣物",引来许多"奇迹治愈";圣于贝尔镇(St. Hubert des Ardennes)因此在修道院周围兴起。于贝尔本是年轻贵族,打猎时见一白鹿角间有十字架,遂皈依。此情节直接抄袭圣尤斯塔斯传。于贝尔是猎人主保圣人;今天镇上各种仪式几乎从不提他与狂犬病的故事。

11 世纪某僧侣记载,被狂犬(或疑似狂犬)动物(多为狗或狼)咬伤的人会被带到于贝尔神龛。仪式由神父在朝圣者前额开一口、塞入取自已贝尔主教披带的一根线;修士们说披带由圣母编织、亲手由天使送达。"线"仪式(la taille)之所以被联系到狂犬病,可能是因当时认为狂犬由狗肛门或舌下的虫引起,圣披带的线由此被视为"接种"。中世纪那件据说永不缩小的披带,至今可在教堂的圣髑盒中看到。18 世纪初(吸血鬼传说盛期)记录显示有 1956 人接受过 la taille;1920 年代仍有信众;据 Morrison 1988 年的访问今天或许仍有。教堂近圣于贝尔祭坛的墙上有一大铁环,被绑在上面的是些"狂吠、嚎叫、抽搐"的可怜人;九日待满后若仍"治愈"则为圣徒显灵,若死了就是"未蒙垂怜"。圣徒与修道院稳赢。狂犬病歇斯底里(rabies hysteria,指以为自己得狂犬病而出现大部分"狂躁"症状)在文献中记录良好;铁环九日后仍未死的"可怜人"死了,那就说圣于贝尔已决定不出手。

吸血鬼在 17 世纪最后 25 年首次被广泛提及,认为是复活的尸体爬出坟墓、吸睡人血。18 世纪在巴尔干尤其被恐惧,Voltaire 称 1730–35 年间为"唯一话题"。Gómez-Alonso (1998) 的有趣论文提出一种可能解释,认为狂犬病可能是源头,本文的主要元素取自其工作;原文含有新旧参考文献的全面清单。传说起于 17 世纪末一份报告:尸体中充满液态血。村民们说看见狗形的鬼或可怕的人攻击人、抓人喉咙。报告提供其他血淋淋细节。故事在 1731–1732 年塞尔维亚 Medvedja 村达到顶点:农民死亡被归因于吸血鬼。17 具被发掘的尸体被木桩刺穿、斩首、火化——这正是 Voltaire 提到的"唯一话题"和众多启蒙运动知名人物的来源。文学与电影中著名的 Dracula 实际在 19 世纪末才出现。

吸血鬼据说有种种习性:离墓去性交或猎杀无辜者取血。被吸血鬼袭击、吃过吸血鬼杀死的动物、死于狂犬病或鼠疫、甚至"伟大的情人"皆可变成吸血鬼;预防方式也多种多样。两具"吸血鬼"尸体的可见迹象是:突出的生殖器、充满血的肿胀身体,血从口流出。狂犬病是人畜共患病,处于狂躁期会引发不可预测的暴力与攻击行为。某些边缘系统疾病也影响性行为。埋在冷湿处的尸体往往延迟分解,使皮下组织呈蜡样。关于"液态血",某些疾病延长了液体期;分解最终发生、内脏溶解时产生的气体使各部分(如生殖器、面部)膨胀、舌头伸出、血沫从口流出。

人类狂犬病症状:多数为"狂躁型"而非麻痹型——失眠、不自主激动、恐水(hydrophobia 旧名)、肌肉痉挛、怕照镜子及其他恐怖怪异表现。面部肌肉痉挛使嘴唇缩回露齿、发出不可辨识之声。有狂躁病人(如同狂犬)冲向人攻击、咬人。间歇期患者口滴血。也有极度亢进——多日持续勃起、频繁性交与猛烈强奸企图(参 Warrell, 1977)。狂犬病可经多种方式人传人:动物/人咬伤、生殖器黏膜等(Warrell, 1977)。众多理论试图解释这一传说,从单纯迷信到精神分裂。疫情期间尸体有时埋得浅被狗/狼挖出,使"吸血鬼爬出坟墓"的想法合理化。19 世纪前的流浪狂犬患者其攻击与亢进可支撑吸血鬼传说。Gómez-Alonso (1998) 文献详实、令人信服地论述了狂犬病产生吸血鬼传说的可能性。狂犬病史与人们对它的态度既迷人、引人入胜、又常恐怖、有时滑稽。

18–19 世纪英国狂犬病:Ritvo (1987) 关于维多利亚时代英国动物–人社会学的史学著作给出一些惊人事实与看法。19 世纪英国疫情(数字其实很小)造成浩劫并催生滑稽的法律与观点。杰出兽医 George Fleming 在其论文中写道:"可能有根据假设强烈的性兴奋会产生狂犬病……"。性与罪、狂犬病之间的关联可被含蓄暗示——维多利亚时代的人喜闻乐见。低下阶层的"curs and lurchers"(偷猎犬)被认为特别易得狂犬病。19 世纪狂犬病地理范围与发生频率扩大,但死亡率从未很高:1877 年 79 人死亡,1879 年 35 人,1875 年 47 人(约 200 万人口对应 2 例)。19 世纪后期的英国公民被谋杀的几率是死于狂犬病的 10 倍以上。Ritvo (1987) 引用当时《Kennel Review》的有趣评论:定义 hydrophobia(狂犬病)为"一种抓住人、驱使人毁灭狗的奇特疯狂"。Surgeon General Charles Alexander Gordon 在上议院狂犬病委员会作证时表达了类似观点。

这些 20 世纪下半叶的遗迹当然仍存在。狂犬病现状:仍是严重疾病,存在于除英国、爱尔兰、瑞典、澳大利亚、新西兰等少数国家外的几乎所有国家。WHO (1998) 是数据来源,1994 年全球估计:人病例约 35,000–45,000(欧洲 10–20、北美 4–8、拉美 200–400、非洲 500–5000、亚洲 35,000–45,000、印度 30,000–40,000)。1994 年孟加拉国 3000 例,美国 6 例(1995、96、97 各 4 例),法国 1994 年 1 例、1995 与 1996 各 3 例。动物病例远高:1994 年美国 8224 例确诊,孟加拉国 960 例实验室确诊加 3500 例非实验室确诊。

法国自 1986 年起连续两年在春夏两季用疫苗诱饵(口服)开展疫苗接种,Aubert (1997) 详述其效果。1989–1996 年间治疗地区动物狂犬病几乎根除。法国建立了从英吉利海峡到瑞士的"免疫屏障"阻止疫情南扩。过去 20 年的所有犬狂犬病均见于进口动物,1995 年最后一例,可能本来可以靠更严的边境管理防止。Pastoret (1998) 讨论比利时疫情,根除后 1994 年再现。Barrat 与 Aubert (1993) 评论法国 1989 年峰值后的下降,部分归因于模型显示的振荡。在全球人流物流背景下,无病国家难以避免再度引入。英国因英吉利海峡隧道与蝙蝠可能带病而"妄想"未减。比利时部分地区发现受感染蝙蝠。Teulières 与 Saliou (1995) 指出 1970–1993 年法国无本土病例,但 14 人在流行区被传染后回国死于狂犬病。人用疫苗方案自 1988 年基于肌注:第 0 日两剂两部位、第 7 与 21 日加强;无失败报告。美国 CDC 推荐暴露者第 0、7、28、365 日接种;保护期 3 年。

吸血蝙蝠在墨西哥、拉美是牛狂犬病爆发的源头。亚洲、拉美、非洲主要为地方性犬狂犬病。多数人通过狂犬动物咬伤或抓伤传染,但与受感染蝙蝠共处洞穴中也可经气溶胶传播。美国狂犬病罕见但几乎都是蝙蝠咬伤引起。1980–1999 年间 25 例中除 3 例外均源自蝙蝠。睡时被咬因针状牙难以察觉,可能是多数案例的成因。CDC 建议醒后若在房间发现蝙蝠须接种(现为 4 周 5 针)。有些传播途径离奇悲惨:14 岁女孩因被患狂犬病的狗舔外阴而传染;人传人最恐怖的是一名妇女通过角膜移植(来自一名感染男子)而传染(Houff et al. 1979),双方死于麻痹型后才发现其眼中含狂犬病毒。角膜移植也涉及克雅氏病(Creutzfeldt–Jakob)人传人案例(Duffy et al. 1974),后者被视为 BSE(疯牛病)的人类形式,由食用受感染牛的牛肉感染。

13.4 狐狸狂犬病空间传播 I:背景与简单模型(The Spatial Spread of Rabies Among Foxes I: Background and Simple Model)

过去几百年中欧洲反复经历狂犬病疫情;当前疫情开始前 50 年左右狂犬病在欧洲消失的真正原因未知;本节及后续分析所给出的模型提供一种可能解释。当前欧洲兽疫约 1939 年始于波兰,稳步西移 30–60 公里/年;河流、高山、高速公路仅能暂时减缓。赤狐是当前欧洲疫情的主要携带者与受害者。狂犬病扩散似行波(图 13.3)。

狂犬病是中枢神经系统的病毒感染,经直接接触传播;狗是向人传播的主要媒介。在欧美人感染罕见,每年仅个别死亡,但欠发达国家多得多。狂犬病对家畜与野生动物影响严重:1980 年法国一年报告 314 例家畜、1280 例野生动物。狂犬病足以引起关切并需深入研究控制策略(13.6 节)。图 13.3 给出狂犬病在法国传播的进程(数据来自 Centre National d'Études sur la Rage,每两年一取 1969–1977 年法国东北部)。Macdonald (1980) 详述当时法国状况,描述了疫苗控制的效果及停止后的反弹。之后多国通过疫苗诱饵成功控制。

美国东海岸也有狂犬病快速传播:主要携带者是浣熊。弗吉尼亚狩猎俱乐部从佐治亚与佛罗里达进口受感染浣熊,使该疫情大幅加速。若再回到图 13.1,恰如卷 I 第 10 章讨论的均匀流行病系统,疫情过后一部分易感者存活。能在空间语境下解析地估计这个存活比例是值得做的。作者用一个非常简单的狂犬病空间传播模型做这件事。

西欧红狐占记录病例约 70%。英国自 1900 年起实际无狂犬病,但可能因非法宠物进口或大陆受感染蝙蝠而短期内再现。英国形势特别严峻,城乡狐、猫、狗密度高。布里斯托尔狐密度约 12 狐/公里²(农村 2–4/公里²)。Macdonald (1980) 的狐与狂犬病专著提供英国大量事实与数据。欧洲狂犬病数据可查法国 Centre National d'Études sur la Rage。Kaplan (1977) 与 Bacon (1985) 编辑的书专门讨论狂犬病种群动力学,给出生物学与生态学背景与一些数据。

理解狂犬病兽疫波前如何向未感染地区推进、控制方法如何阻断、参数如何影响,是重要问题。本章余下各节专门讨论这些空间问题。材料主要来自 Murray et al. (1986) 模型与 Källén et al. (1985) 更简单但不太逼真的模型。13.6 与 13.9 节的控制策略具体针对欧洲当前的狐兽疫,模型类型适用于许多其他空间传播疫情。

空间传播常是非常复杂过程,狂犬病也不例外。建模可尽量包含所有事实(必含很多参数,估计难);或者以尽可能简单的模型抓住关键要素,使得较少参数可以被确定。本章采取后者策略——这些模型虽简单却能提出高度实际的问题并给出重要特征的有用估计。本节描述与分析一个特别简单的模型,但它能给出有用的解析结果。

虽然涉及多种动物,基本且合理的假设是狐(主要媒介)的生态决定疫情扩散动态。进一步假设兽疫的空间传播主要由狂犬病狐的随机无规律游走造成;未感染狐不远游(Macdonald 1980)。把狐种群分两类:易感 S 与染病 I(狂犬病 + 潜伏期)。此模型虽捕捉了兽疫波前的某些方面,但漏掉了狂犬病的基本特征——长达 12–150 天的潜伏期(从被咬到出现临床症状)。13.5 节的更逼真模型会纳入此因素。

为评估控制策略必须先理解狂犬病如何传播。本节首先研究一个特别简化的流行病模型版本 (13.1) 的变种,以捕捉狐群中狂犬病传播的关键要素。然后用其估计兽疫波的关键事实。

考虑狐分两类:感染者 I 与易感者 S;感染者包括狂犬病狐与潜伏期狐。主要假设:(i) 病毒经唾液从感染狐传给易感狐,每头感染狐对易感狐的平均感染率 rI,r 是传播系数,衡量两类间接触率;(ii) 狂犬病必然致死,狐按人均率 a 死亡,即感染狐寿命 1/a;(iii) 狐具领地性,把乡野分成不重叠的领地;(iv) 病毒入中枢神经系统改变行为,入脊髓则瘫痪、入边缘系统则短期攻击性行为,失去领地感随机游走——假设是感染者以扩散系数 D 公里²/年扩散。模型即 (13.1) 但易感狐不扩散(也不考虑幼狐迁移)。一维模型为 (13.12):∂S/∂t = −rIS、∂I/∂t = rIS − aI + D ∂²I/∂x²。预期系统有行波解,其速度依赖参数值。

用 (13.2) 的无量纲化把 (13.12) 化为 (13.13):∂S/∂t = −IS、∂I/∂t = I(S − λ) + ∂²I/∂x²。λ = a/(rS₀) 度量相对死亡率与接触率。预期阈值仍为 λ = 1(参见练习 2)。

行波前沿解取 (13.14) 形式 S(x, t) = S(z)、I(x, t) = I(z)、z = x − ct,c 为波速。边界条件 (13.15):S(∞) = 1、S′(−∞) = 0、I(∞) = I(−∞) = 0。S′(−∞) = 0(而非 S(−∞) = 1)反映了波过之后剩余的未确定数目的存活易感狐。代入 (13.13) 得 (13.16):cS′ = IS、I″ + cI′ + I(S − λ) = 0。在 I = 0、S = 1 附近线性化,并要求 I 恒非负,得到 c ≥ 2√(1 − λ)、λ < 1((13.17))。

本模型能更进一步,求出疫情过后存活的易感者比例。由 (13.16) 第一式 I = cS′/S,代入第二式得 I″ + cI′ + cS′(S − λ)/S = 0。积分得 I′ + cI + cS − cλ ln S = constant。用 z → ∞ 时 S = 1、I = 0、I′ = 0 定常数为 c。令 z → −∞ 用 I = I′ = 0 得超越方程 (13.18):σ − λ ln σ = 1、λ < 1、σ = S(−∞),即存活易感比例 σ 与 c 无关。改写为 (σ − 1)/ln σ = λ < 1 ⇒ 0 < σ < λ < 1。λ 度量疫情严重程度;λ 越小则存活越少(疫情越重)。例如 λ = 0.4 ⇒ σ = 0.1,λ = 0.7 ⇒ σ = 0.5。图 13.4 给出 σ(λ)。

临界分支值 λ = 1(维度 a/(rS₀) = 1)。若 λ > 1,无兽疫波可传播——这在物理上对应 a > rS₀(死亡率高于新感染者的招募率)。给定 r 与 a 存在临界最小狐密度 Sc = a/r,低于此则狂犬病不能在狐群中持续。

λ < 1 时,计算出的兽疫波速为最小可能速度 c = 2√(1 − λ),维度形式 (13.20):c = 2√(D(rS₀ − a))。图 13.5 给出 λ = 0.5 时 S 与 I 的行波前沿数值解;由图 13.4 在 λ = 0.5 时 σ ≈ 0.2。

将图 13.5 中易感狐的定性形式与欧洲数据图 13.6 比较:两图在波后行为有明显差异。(13.13) 只覆盖疫情波前的通过;显然波前通过后易感狐因发现新承载容量更大的环境而开始增加。即 (13.13) 的时间尺度远短于图 13.6 振荡的时间尺度。建模应在模型中加入狐繁殖项。把易感者的第一个方程加 logistic 项得 (13.21):∂S/∂t = −rIS + BS(1 − S/S₀),B 为线性增长率。用 (13.2) 的无量纲化得 (13.22):∂S/∂t = −IS + bS(1 − S)、∂I/∂t = I(S − λ) + ∂²I/∂x²,b = B/(rS₀)。图 13.7 给出由 (13.22) 数值解得到的易感者与感染者的疫情波:初始波后接有反复但较小的疫情;振荡衰减,最终 S → λ、I → b(1 − λ)——(13.22) 的稳态解。这与图 13.6 数据符合良好。

可以用 (13.22) 的维度形式 (13.23) 给出有用的解析结果。引入无量纲量 (13.24) 把 (13.23) 化为 (13.25):U_t = U(1 − U − V)、V_t = αV(U − λ) + V_xx。此即第 1 章 (1.3) 方程(无被捕食者扩散,α、λ 替换 a、b)。稳态为 (0, 0)、(1, 0)、(λ, 1 − λ),最后一个只在 λ < 1 时存在正象限。相空间 (U, V, W),W = V′ 的分析(1.2 节)表明 λ < 1 时存在 (1, 0) 与 (λ, 1 − λ) 之间的行波解,并存在阈值 α:α > α 趋向 (λ, 1 − λ) 的方式为振荡,α < α 为单调(图 1.3)。图 13.7 即 α > α 的例子。回到 13.2 节关于黑死病后续爆发的观察:把 (13.1) 的易感者方程修改为纳入疫情后人口恢复,可以得到主疫情后接续的周期性小疫情,与图 13.6、13.7 类似。这解释稍嫌轻巧,因黑死病涉及人、跳蚤、鼠等多群体互动。尽管简单,结果能定性反映主要现象。

13.5 狐狸狂犬病空间传播 II:三物种 (SIR) 模型(The Spatial Spread of Rabies Among Foxes II: Three-Species (SIR) Model)

为实用地发展控制策略以遏制疫情空间传播,需考虑更逼真且更复杂的模型,以便与已知数据定量比较并更可信地做实际预测。13.4 节模型过于初步不能定量。关键遗漏是 12–150 天的长潜伏期(到发病前),本节考虑更逼真的模型:用它可对流行病学与公共卫生重要的时间与距离给予定量估计。

三物种模型仍以狂犬病狐为空间传播主因。野生狂犬病狐移动数据虽少但非零,部分将用于估计狂犬病狐关键扩散系数。模型虽仍较简单,但部分参数难从现有数据估计。Murray et al. (1986) 是该工作基础,扩展了 Anderson et al. (1981)(仅时间相关情形)的工作,加入空间效应,特别是狂犬病狐的关键空间扩散,13.9 节再加入所有狐的扩散。

三物种 SIR 模型把狐群分为:易感 S、感染但非感染性(潜伏期)I、感染性(狂犬病)R。需要至少三物种的关键是狂犬病病毒在感染动物体内有 12–150 天(有时更长)的长潜伏期,期间行为正常不传播疾病;之后是 1–10 天的临床期。

基本模型假设同 13.4 节(记号略异)但在此重述便于参考:

(i) 无狂犬病时狐种群动态近似为简单 logistic:dS/dt = (a − b)S(1 − S/K),a 为线性出生率、b 为内禀死亡率、K 为环境承载容量。a、b、K 可随栖息地变化,此处取常数;13.8 节对英格兰"实验"令 K 变化。

(ii) 狂犬病由狂犬病狐经直接接触(通常咬伤)传给易感狐;易感狐按平均人均率 βR 感染(β 为接触率常数)。

(iii) 感染狐按人均率 σ 转入感染(狂犬病)期,1/σ 为平均潜伏期。

(iv) 狂犬病必然致死,狂犬病狐按人均率 α 死亡,1/α 为临床病平均时长。

(v) 狂犬病与感染狐继续消耗环境,因非狂犬病原因死亡但有可忽略的健康后代(为完整性而加入)。

空间假设:(vi) 狐具领地性,把乡野分成不重叠范围。(vii) 约一半感染狐得"狂躁型"狂犬病——表现出狂犬病的凶猛症状;其余病毒侵脊髓导致瘫痪。狂躁型狐可能攻击性强、失去方向与领地行为、随机游走——它们是空间扩散主因。

模型 (13.26):

∂S/∂T = aS − bS − (a − b)NS/K − βRS

∂I/∂T = −bI − (a − b)NI/K + βRS − σI

∂R/∂T = −bR − (a − b)NR/K + σI − αR + D ∂²R/∂X²

总种群 N = S + I + R((13.27))。易感狐出生是唯一源项,所有狐都自然死亡(期望寿命 1/b 年)。(a − b)N/K 项代表所有狐对食物供应的消耗。从易感到感染靠 βRS,从感染到感染期靠 σI。狂犬病狐还死于狂犬病以 αR 表示,期望寿命 1/α,并按 D 扩散。除 D 外的典型参数值由 Table 13.1 给出。空间均匀时把 (13.26) 三式相加得 (13.28):dN/dt = aS − bN − (a − b)N²/K − αR,即总种群的等效 logistic 形式。模型写成一维形式但 13.8 节用于英格兰(从假设疫情开始)的二维实际情形。

此模型忽略年轻游荡狐的空间扩散——它们可能在寻领地时被咬、将狂犬病带到新地。本节给出的依据是狂犬病在幼狐中比成狐少见得多(Artois 与 Aubert 1982, Macdonald 1980)。

Table 13.1(参数值,Anderson et al. 1981):平均出生率 a = 1 狐·年⁻¹;内禀死亡率 b = 0.5 狐·年⁻¹;临床病平均时长 1/α = 5 天;平均潜伏期 1/σ = 28 天;临界承载容量 K_T = 1 狐·公里⁻²;传播系数 β = 80 公里²·年⁻¹;承载容量 K = 0.25–4.0 狐·公里⁻²。

空间均匀稳态 (13.29) 经代数运算得到。S₀、I₀、R₀ 非负要求 (13.30):K_T = (σ + a)(α + a)/(σβ),即 K > K_T 时存在非零稳态。

空间均匀 (D = 0) 时若把狂犬病引入稳定健康狐群有三种可能行为,由 K 相对于 K_T 的大小决定。若 K < K_T,疫情最终消失(R → 0、I → 0),种群回到 S = K。K > K_T 时种群绕稳态振荡。对 (S₀, I₀, R₀) 做线性稳定性分析可证:若 K 不太大则稳态稳定、扰动以振荡方式衰减;若 K 足够大则出现极限环。故 K 有两个分支值:K_T 与分隔极限环与稳定稳态的临界 K。

从流行病学证据看,狂犬病在承载容量 0.2–1.0 狐/公里² 之间某处会消失(WHO Report 1973, Macdonald 1980, Steck 与 Wandeler 1980, Anderson et al. 1981, Boegel et al. 1981)。β 不能直接估计——接触难以观察。Anderson et al. (1981) 用 (13.30) 间接估计 β,因为 K_T 与除 β 外的参数已知。Murray et al. (1986) 详细讨论它们对狂犬病空间传播的影响;模型在许多参数估计值附近对变动是稳健的。Bentil 与 Murray (1991) 给出在观察受限情况下获取参数信息的另一种方法。

K > K_T、Table 13.1 的参数选择下,振荡周期 3–5 年、狂犬病平衡持续率 p = (R₀ + I₀)/(S₀ + I₀ + R₀) 为 0–4%((13.31))。与现有流行病学证据吻合(Toma 与 Andral 1977, Macdonald 1980, Steck 与 Wandeler 1980, Jackson 与 Schneider 1984)。

行波兽疫波前与传播速度:引入无量纲量 (13.32),模型化为 (13.33):

∂s/∂t = ε(1 − n)s − rs

∂q/∂t = rs − (μ + δ + εn)q

∂r/∂t = μq − (d + εn)r + ∂²r/∂x²

n = s + q + r

正均匀稳态 (s₀, q₀, r₀) 由 (13.29)/K 得到。疫情发生的条件 K > K_T 等价于 (13.34):0 < d < (1 + δ + ε/μ)⁻¹ − ε。系统依赖 4 个无量纲参数 ε、δ、μ、d,原维度 7 个。取代表承载容量 K = 2 狐/公里²得 ε = δ = 0.003、μ = 0.08、d = 0.46。ε、δ 远小于 1、μ、δ、1 − d 这一点可用于简化 (13.33) 的分析并获得有用的解析结果。

无量纲化再次让参数组合显示对实际场参数变化的等效效果——例如 ε、δ 小意味着疫情期间传染率远大于非狂犬病引起的出生与死亡率。

现在求兽疫波以常速 v 传向未受扰动、无狂犬病区域的解 (13.35):

vs′ = ε(1 − n)s − rs

vq′ = rs − (μ + δ + εn)q

vr′ = μq − (d + εn)r + r″

n = s + q + r

s → 1、q → 0、r → 0 当 ξ → −∞(波前远处)。以下利用 ε ≪ 1、δ ≪ 1。

(13.35) 在正象限有三个稳态:(1, 0, 0)、(0, 0, 0) 与 (s₀, q₀, r₀)。ε、δ 小时到一阶展开为 (13.36):s₀ = d + (ε + εd + δ/μ)d,q₀ = εd(1 − d)/μ,r₀ = ε(1 − d)。s₀、q₀、r₀ 非负要求 (13.34)。

(13.35) 的行波解是 4 维相空间 (s, q, r, r′) 中从 s = 1、q = r = 0 出发的轨迹。线性化 (13.35) 在临界点 (1, 0, 0, 0) 附近给出 4 维一阶系统,解为线性组合 xᵢ exp(λᵢξ)。特征值 (13.37) 之一为 λ = −ε/v < 0,另三个是 f(λ) = λ³ + (μ + δ + ε/v − v)λ² − (d + μ + δ + 2ε)λ + μ(1 − d − ε) − (δ + ε)(d + ε)/v = 0 的根。

f(λ) → ∞ 当 λ → ∞、f(λ) → −∞ 当 λ → −∞。(13.36) 成立时 f(0) > 0 且在 λ = 0 斜率为负。随 v 变化 f(λ) 依次呈现图 13.8(a)、(b)、(c) 三种形状。阈值条件 (13.34) 成立时 f 有一个负实根,随 v 值另有两个正实根或两个复根。复根对应振荡解即负值种群,物理上不现实。临界速度 v_c 是 f = 0 与 df/dλ = 0 同时成立时消去 λ 所得的双重根速度(图 13.8(c))。对 v_c 的大量代数得 (13.38):g(z) = [4μ + (d − μ)²]z³ + 2[3μ(1 − d)(3d + μ) + (d + μ)²(2d + μ)]z² + μ²[(d + μ)² − 6(1 − d)(3d + μ) − 27(1 − d)²]z − 4μ⁴(1 − d)。阈值条件 (13.34) 成立时 g(z) 在 z = 0 负且 d²g/dz² 正;g(z) 略图显示它有唯一正根对应兽疫波的最小可能速度。

接下来证不可能有从 s = 1、q = r = 0 到原点的轨迹。线性化 (13.35) 在原点附近的特征解 (13.39) 中分别取指数项;轨迹可表示为这些特征解的线性组合。充分接近原点时接近原点的轨迹是负指数特征解的线性组合,因此在 s = 0 平面内接近原点。把 (13.35) 中的 ξ 替换为 τ = −ξ(即可逆),令 s = 0 初始则 s = 0 对所有正 τ 成立——不论 r、q 初值如何。这意味着任何 ξ 处 s = 0 的轨迹其过去与未来所有 ξ 上 s = 0;故从 s = 1 进入 s = 0 平面再接近原点的轨迹不存在。

因此行波仅当存在从 s = 1 到 (s₀, q₀, r₀) 的轨迹时可能——这要求 (13.34) 成立。在 (s₀, q₀, r₀) 附近线性化 (13.35)(大量代数后)得一阶展开的特征值 (13.40):

λ₁, λ₂ = (1/2)[(v − μ)/v ± √{(v − μ)/v}² + 4(μ + d)]

二阶展开为 (13.41):

λ₃, λ₄ = ±(i/v)√{εμd(1 − d)/(μ + d)} − (εd/[2v(μ + d)²])(μ(1 − d)(μ/v² − 1) + (μ + d)²)

λ₁ 正;任意接近 (s₀, q₀, r₀) 的解为对应 λ₂、λ₃、λ₄ 的特征解线性组合。由于 |λ₂| ≫ |Re(λ₃, λ₄)|,其特征解的振幅比复特征值对应的衰减快得多。故波尾(ξ 大)解由复特征值主导:

s − s₀ ∼ [A cos(ωξ/v) + B sin(ωξ/v)] exp(−λξ/v)

q − q₀ ∼ (ω/μ)[A sin(ωξ/v) − B cos(ωξ/v)] exp(−λξ/v)

r − r₀ ∼ (ω/d)[A sin(ωξ/v) − B cos(ωξ/v)] exp(−λξ/v)

ω 是波长(复特征值虚部除以 v),λ 是衰减率(实部除以 v)。A、B 由轨迹接近方式决定,不可由线性分析得到。

利用 ε、δ 小作渐近近似((13.36)):r₀ = μq₀/d。(13.41) 同样显示波尾 r − r₀ ∼ μ(q − q₀)/d——感染与狂犬狐密度在全波上比例相似,仅尺度不同。图 13.9、13.10 的全非线性数值模拟中此相似对整个波都成立,提示在复杂度下三物种模型可以高近似程度简化为二物种模型。即可用二物种(易感 + 狂犬狐)替代三物种 SIR 体系;感染但未发病的狐由狂犬狐简单标度 q(ξ) ∼ dr(ξ)/μ((13.42))。Murray et al. (1986) 给出奇异摄动分析(基于 μ 小于 d 与无量纲波速 v、但大于 ε、δ)。

兽疫行波前沿须数值求解 (13.33)——初始 s = 1(S = K)各处、x = 0 处少量狂犬狐。(13.34) 成立时形成兽疫波,并以近常速自初始集中处向外传播。违反 (13.34) 时狂犬病消失,种群回到环境承载容量。图 13.9 给出按欧洲当前兽疫参数解出的行波前沿——主要由大量狐死亡的疫情前阵与后接的振荡尾组成,后者的每次疫情复发较前次更弱,振荡逐渐衰减到常数非零值(狂犬与感染狐为零)。图 13.10 给出按英国参数解出的行波——易感者密度振荡更剧烈。

Table 13.1 的参数下 (13.34) 满足且 ε、δ ≪ 1、d、μ、1 − d。波后接振荡尾。解析上最小速度 v = z^(1/2) 为 (13.38) 的正根 z。图 13.11 给出 0 ≤ d ≤ 1 范围内该根 v 的等高线图。所有数值模拟解出的波似均以此最小速度传播,由 (13.32) 维度形式 (13.44):V = √(DβK) v。例如 Table 13.1、扩散系数 200 公里²/年、承载容量 2 狐/公里²下,从 (13.32) 算 d 与 μ、从图 13.11 读 v,得维度传播速度 V = 51 公里/年。

在 (s₀, q₀, r₀) 附近的线性分析表明对充分大 t 波趋于 (13.41) 的衰减振荡。原 (x, t) 变量下 (13.45):

s(x, t) = s₀ + A cos[ω(t + x/v) + ψ] exp[−λ(t + x/v)]

q(x, t) = q₀ + (1/μ)(s − s₀)′

r(x, t) = r₀ + (μ/d)(q − q₀)

一阶 ε、δ 中,无量纲波数 (13.46):ω = √(εμd(1 − d)/(μ + d)) + O(ε^(3/2)),衰减率 (13.47):λ = (εd/[2(μ + d)²])(μ(μ/v² − 1)(1 − d) + (μ + d)²)。A、ψ 为常数。易感者振荡与感染及狂犬者均相位差 90°(对称性在下阶 ε、δ 中打破)。r − q 比例关系 (13.42) 在合理参数下数值模拟中普遍成立。

Murray et al. (1986) 的奇异摄动分析给出几个有用的渐近估计。例如首次疫情中感染与狂犬狐最大密度 (13.48):

r_max ≈ μ(ln d + 1 − d)/d

q_max ≈ d(ln d + 1 − d)/d

维度形式 (13.49):

R_max ≈ (σK_T/α)(ln(K_T/K) + K/K_T − 1)

Q_max ≈ K_T(ln(K_T/K) + K/K_T − 1)

K ≤ K_T 时不发生疫情,R_max = Q_max = 0。两者随 K > K_T 增大。

由 (13.38) 得无量纲波速 v = z^(1/2),再由 (13.47) 得衰减率 λ。λ 总为正。这意味着无扩散版本 (13.26)(即 D = 0)对足够大 K > K_T 能表现的极限环行为在加入扩散时消失:振荡总衰减到常数态 (s₀, q₀, r₀)。维度衰减率为 βKλ。维度周期 (13.50):

T = 2π/√{((α + σ + b)[(a − b)(α + b)σ]^(1/2)(1 − (α + b)/(βK))}

T 随 K 增大而减小——一般地疫情前狐密度越大,疫情后波后周期性疫情间隔越长,与某些观察一致(Macdonald 1980)。但数值发现前阵附近非线性重要处,疫情间隔可能随 K 增大而增加(图 13.9、13.10)。维度波长 L = VT。

扩散系数 D 的估计与波速、疫情波长对 D 变化的灵敏度:要计算兽疫维度速度 V,从而算出周期与波长,须估计 D——它是狂犬狐覆盖地面速率的度量。野生狂犬狐行为所知甚少,使 D 难估。Andral et al. (1982) 在野外跟踪三只成年狂犬狐——给捕获狐接种狂犬病毒、配发信号项圈、原捕获地释放。他们在潜伏期跟踪狐以确定其家园范围与正常行为,在狂犬期跟踪以观察疾病引起的行为变化。狐变狂犬后每日活动模式改变。画图显示三只狐的潜伏期范围与狂犬期主要位移——所有三只都在狂犬期某一刻离开家园但都未走得太远。Murray et al. (1986) 用 Andral et al. (1982) 结果以"原始"方式估扩散系数:

D ≈ (1/N) Σⱼ (起点到死亡点的直线距离)² / (4 × 起始时间)

用狂犬期起点到死亡点距离与狂犬期近似时长,得 D ≈ 50 公里²/年。因三只狐中有两只死亡时距起点比其平均距离更近,这很可能是 D 的下界。粗略上界可从任一只狐离起点最远距离估——一只狐在狂犬期过半时离起点最远达 2.7 公里,给出 D 的上界 330 公里²/年。

扩散系数还有别的估计方法。Källén et al. (1985) 估二物种模型时假设感染狐在一个月潜伏期末才离家(即在它们被假定变狂犬时),取平均领地 5 公里²,得 D = 60 公里²/年。给三物种模型确定 D 需要估计临床病发病后狐离开领地的平均率。若 N 只感染狐中第 j 只在变狂犬后 tⱼ 时间离开领地,可由

(1/N) Σⱼ 1/tⱼ

估 k。由于约一半感染狐发展为麻痹型狂犬病可假定永不离开家园,对约 N/2 狐 tⱼ 无穷。对狂躁狐,假设其中也有一半永不离开,其余在其后 6 天(病程)内均匀离开,则

k ≈ (1/N) Σⱼ^(N/4) 1/tⱼ = (1/24) Σⱼ^6 1/j 天 ≈ 40 年⁻¹

取平均领地 5 公里²(Toma 与 Andral 1977;Macdonald 1980),得 D = 190 公里²/年。另一方法:估计狂犬狐的平均自由程与速率。Andral et al. (1982) 观察狐在狂犬期平均每日共走 9 公里。假设这具有代表性,例如狂犬狐连续走 100 米后被干扰转向,则 D = 速率 × 步长给出 330 公里²/年——与前面估计的上界相同。所有这些方法在有足够观察时原则上一致;目前野生狐行为观察数据不足,无法获得更佳估计。由于波速正比于 D^(1/2),把 D 从 50 改到 330 公里²/年使 V 增为 2.6 倍。Table 13.2 给出给定 D = 200 公里²/年时波速与波长对承载容量的灵敏度。

另一难估参数是传播系数 β。如前所述,可由 (13.30) 反向估计。狐种群密度的绝对值实际上难获——通常通过报告死亡/射杀/毒杀的狐数估计并对占总数百分比作假设,或通过与已知狐密度地区地形的比较。K_T 估计特别难,文献中估为 0.2–1.2 狐/公里²(WHO Report 1973, Steck 与 Wandeler 1980, Macdonald et al. 1981, Gurtler 与 Zimen 1982)。由于 K/K_T 只涉及种群大小的比较,这个比值可能比 K 与 K_T 各自更易获得。一个相关问题:定量结果对参数不确定性的灵敏度如何。Murray et al. (1986) 详细讨论。

13.6 基于波传播入无疫情区的控制策略:狂犬病屏障宽度估计(Control Strategy Based on Wave Propagation into a Nonepidemic Region: Estimate of Width of a Rabies Barrier)

Murray et al. (1986) 提出一种可能的控制策略——在行进的波前前方的易感狐种群中建立"保护屏障",把易感狐密度降到临界密度 K_T 以下。例如在丹麦的日德兰半岛成功实施;意大利、瑞士部分地区也曾以各种强度执行,结果各异(Macdonald 1980, Westergaard 1982)。屏障可由扑杀或疫苗接种产生。扑杀释放了领地,可能让年轻狐快速占领反而促进疾病扩散。疫苗接种对生态破坏少、更有效也更经济。

要使"阻断"有效需对其宽度与可允许的易感狐密度有合理估计。本节解析推导保护阻断区所需宽度,给出 (13.33) 全方程系统数值模拟的部分结果。以下"感染狐"指所有带狂犬病的狐,无论是否感染性。

观察狂犬病兽疫波在固定地点的通过:每次疫情后接长静止期(图 13.9、13.10)。时空尺度上二次疫情波距主波足够远——当二次疫情到达时第一次波要么已过阻断、要么已基本消失。每一次后续疫情较前次更弱。因此合理假设能消除第一次疫情的减员方案对后续疫情也有效——故只需考虑使第一次疫情停止所需的阻断宽度。阻断宽度依赖于阻断内易感狐种群密度。

由于空间扩散由确定性扩散机制建模,从严格数学角度感染狐密度不可能在任何处为零——把狐密度视为空间时间连续且用经典扩散模拟狂犬狐扩散。因此不能让兽疫波进入有限宽度的阻断区后确定另一侧感染狐密度是否仍为零——它总为正,只是指数小。故无论多宽的阻断最终也有足够感染狐泄漏使另一边出现疫情。因此必须考虑的是"已足够小以使感染狐到达阻断对侧的概率可接受"的时间。

由于控制方案旨在保持狐密度小,把阻断区视为承载容量低于 K_T 的区,并假设疫情前阻断内狐密度已降到此值。为估计阻断宽度,把易感狐密度降低区设在 x = 0 到 x → ∞ 处。先给出全系统 (13.33) 的数值模拟结果,再于本节后半得到近似解析结果。

图 13.12 与 13.13 显示当兽疫波自左进入阻断区时发生什么。兽疫波在承载容量低于临界值 K_T 时不能传播。感染狐最大密度在 x = 0。感染波进入 x > 0 区后扩散、振幅衰减、感染狐总数减少。最终剩余感染狐数少于 p 狐/公里²,p 为某小数。tc(p) 为发生此情况的时间。选 p 足够小使狂犬狐在临界时间后遇到易感狐的概率可忽略。波在阻断区不能传播故只是衰减——所有时间感染狐密度都在阻断边缘最大并随 x 指数衰减(实际为 x² 衰减)。选阻断宽度为感染狐密度为原点值的给定(小)分数 m 的点 xc:

I(x_c, t_c) + R(x_c, t_c) = m[I(0, t_c) + R(0, t_c)] (13.51)

现有证据表明不可能把狐从某区完全消除——70% 减员是能达到的最好(Macdonald 1980)。图 13.14 给出对不同临床病平均时长 1/α 情况下,阻断宽度关于减员百分比的关系。

图 13.14 的数值模拟中:阻断外 βK 保持 160 年⁻¹;临界时间感染狐数取 p = 0.5 狐/公里²;m 取 10⁻⁴;除 α 外的参数取自 Table 13.1。由此 d = (α + 0.5 年⁻¹)/(160 年⁻¹)、(13.30) 给出阻断外 K = 149/(α + 0.5 年⁻¹) 狐·公里⁻²·年⁻¹。例如狂犬期平均 3.8 天时 d = 0.6、阻断外 K = 1.5 狐/公里²。减员方案能令疫情到达前阻断内承载容量降到 0.4 狐/公里²则 s_b = 0.26,图 13.14 给 x_b = 15。取 D = 200 公里²/年,(13.32) 给出预测的阻断宽度 17 公里。p 与 m 的选择取决于多保守——Murray et al. (1986) 讨论模型对它们的灵敏度。tc 时所有计算中 I + R 最大值小于 0.15 狐/公里²;m = 10⁻² 时阻断保护侧感染狐少于 0.0015 狐/公里²。

13.7 狂犬病控制阻断宽度的解析近似(Analytic Approximation for the Width of the Rabies Control Break)

可以解析地确定阻断宽度对参数的近似函数依赖。兽疫波到达后阻断区内各种狐种群密度的行为应与下述理想化情形相似:在 t = 0 时刻 x = 0 处引入集中于一点的感染与狂犬狐局部密度(总 I、R 与兽疫波相同),并在承载容量到处等于阻断初始狐密度的域内。令 t = 0 时 r = r₀δ(x)、q = q₀δ(x) 作为 Dirac delta 函数初值——所有 r₀ 狂犬狐初始集中于 x = 0(即假设承载容量为零、s = 0)。

先假设对 x ≥ 0 全部易感狐被消除(例如通过免疫或扑杀)。以下进一步近似认为感染与狂犬狐方程中的非线性项可忽略。由于 ε、δ 是小参数,这是合理近似。数值计算阻断宽度时也发现忽略这些项宽度不变——这是又一佐证。假设下 (13.33) 化为线性 (13.52):

∂q/∂t = −μq

∂r/∂t = μq − dr + ∂²r/∂x²

由对称性,可把 x = 0、t = 0 处 δ 源向 x ≥ 0 区移动的问题替换为初值 (13.53):q(x, 0) = 2q₀δ(x)、r(x, 0) = 2r₀δ(x)(δ 函数幅值加倍以补偿只考虑半轴),考虑全实轴 −∞ < x < ∞。感染狐向阻断的传播由 (13.52) + (13.53) 描述。所关心量:阻断种群衰减到给定水平 p 的时间 t_c,由

√(KD/β) ∫₀^∞ [q(x, t_c) + r(x, t_c)] dx = p (13.54)

隐式定义;阻断宽度 x_c 由 (13.55):q(x_c, t_c) + r(x_c, t_c) = m[q(0, t_c) + r(0, t_c)] 隐式定义(m = 10⁻⁴ 为任意选定的相对衰减阈值)。

先估 t_c。对 (13.52) 沿 x 从 0 到 ∞ 积分得 ODE (13.56):

dQ/dt = −μQ

dF/dt = −dF + dQ*

其中 Q(t) = ∫₀^∞ q(x, t) dx(潜伏期狐总数)、F(t) = ∫₀^∞ [q(x, t) + r(x, t)] dx(总感染狐数)。初值 F(0) = q₀ + r₀、Q(0) = q₀。第一个方程解出 Q(t) = q₀ exp(−μt) 代入第二个得 (13.57):dF/dt = −dF* + dq₀ exp(−μt)。给定初值解为 (13.58):

F*(t) = [q₀ + r₀ − dq₀/(d − μ)] exp(−dt) + [dq₀/(d − μ)] exp(−μt)

临界时间 t_c 由 (13.54) 解 F*(t_c) = p√(β/(KD)) 给出。(13.58) 右端两项都含指数。合理场参数下 d > μ 且 d − μ = o(1/t_c),故第一项相对第二项可忽略(t_c 充分大时)。验证后忽略第一项,代数方程解出 (13.59):

t_c ≈ (1/μ) ln[d√(KD/β) q₀ / (p(d − μ))]

典型值 d = 0.46、μ = 0.08。√(KD/β) q₀ 可由图 13.13 与 q ≈ dr/μ 估——总感染狐数为

∫_{-∞}^{∞} (I + R) dX = √(KD/β) (1 + μ/d) q₀

图 13.13 中 ∫ ≈ 6.9 狐/公里,给出 √(KD/β) q₀ ≈ 5.9 狐/公里。p = 0.5 狐/公里时 (13.59) 估 t_c ≈ 33。exp(−dt_c) 与 exp(−μt_c) 之比约 3 × 10⁻⁶,证明在 (13.58) 中可忽略较小指数。

接下来估阻断宽度 x_c。求解 (13.52) + (13.53)。(13.52) 第一个给出 (13.60):q(x, t) = 2q₀δ(x) exp(−μt)。代入第二个得 (13.61):

∂r/∂t = −dr + ∂²r/∂x² + 2q₀μδ(x) exp(−μt)

带初值 (13.53) 解为 r(x, t) = 2r₀(1/√(πt)) exp(−x²/(4t) − dt) + exp(−μt) r(x, t),其中 r(x, t) 是 (13.62):∂r/∂t = (μ − d)r + ∂²r/∂x² + 2q₀μδ(x) 带齐次初值的解。解此用 Laplace 变换。记 r 的 Laplace 变换为 ρ(x, s) = ∫₀^∞ r(x, t) exp(−st) dt,Re s > 0。则 ρ 满足 (13.63):d²ρ/dx² + (μ − d − s)ρ = −2q₀μδ(x)/s,−∞ < x < ∞,Re s > 0。只对 x > 0 解为 (13.64):ρ(x, s) = μq₀ exp[−(s + d − μ)^(1/2) x] / [s(s + d − μ)^(1/2)]。反演得 r(x, t) = μq₀ (1/2πi) ∫_C exp[−(s + d − μ)^(1/2) x] exp(st) / [s(s + d − μ)^(1/2)] ds,C 为 Bromwich 围道。被积函数奇点为 s = 0 极点与 s = −(d − μ) 分支点。分支割可沿负实轴左到分支点左;故围道可变形到负实轴上下。t = t_c 时只需求 r*(x, t),可设 t ≫ 1 在 (13.64) 积分中。用最陡下降法(参 Murray 1984 第 6 章)知积分主贡献来自 s = 0 留数;分支割贡献指数小比较——若

(x/(2t))² ≪ d − μ (13.65)

则可如此近似。r(x, t) 的渐近解为 (13.66):

r(x, t) ∼ r₀(1/√(πt)) exp(−x²/(4t) − dt) + (μq₀/√(d − μ)) exp(−μt − (d − μ)^(1/2) x)

估阻断宽度:(13.55) 不能直接用——(13.60) 给出 q(x, t) 总是含 δ 函数。改为 (13.67):r(x_c, t_c) = m r(0, t_c)。用 (13.65) 与 t ≫ 1 可再忽略 (13.66) 第一项相对第二项(前者因 x²/(4t) ≪ 1 即小 x 区域也指数小)。则 (13.67) + (13.66) 估阻断宽度 (13.68):

x_c ∼ √(d − μ) ln(1/m)

取 m = 10⁻⁴ 与之前估 t_c 的参数:(13.65) 在 t = t_c、x = x_c 时易验证——(x_c/(2t_c))² ≈ 0.05、d − μ ≈ 0.38,满足 ≪ 关系。x_c 表达式到主阶与临界时间 t_c 无关——t_c 的计算只是为了验证 t 大假设。

维度形式 (13.69):

X_c ∼ (1/√(βK)) √(D/(α + b − σ)) (−ln m)

典型参数取 Table 13.1。(13.68) 中 x_c 对 d 与 m 的依赖大致符合图 13.14。x_c 应不很敏感于 p——Murray et al. (1986) 表明当阻断内承载容量不太接近临界值时确实如此。

13.8 二维兽疫前沿与可变狐密度的效应:对英格兰疫情的定量预测(Two-Dimensional Epizootic Fronts and Effects of Variable Fox Densities: Quantitative Predictions for a Rabies Outbreak in England)

一般狐种群不均匀而随局部环境的宜人度与承载容量变化。英格兰就极典型——有趣的是其中部分最高密度(高出 2–3 倍)在城市如布里斯托尔。

先给出一维模型在二维情形下的结果:把 (13.26) 方程中 R 的扩散项改为 D∇²R。承载容量 K 与初始易感狐密度在方形区域上为均匀值,仅中央有一小块值不同。把狂犬狐均匀分布在方形一边上,使一维疫情前阵贯穿方形,数值求解。图 13.15 给出中央小块具较高初始易感狐密度的情形。

由图 13.15(b) 可见前阵在高承载容量区移动更快。首疫情过后剩余狐群在 K 较高的"口袋"内比其周围略低。对低密度口袋情形则反。低密度口袋对紧邻区有某种保护——其外圈不会有那么多狂犬病案例且最终易感狐密度较高。13.6 节的阻断区也呈此特征,源于低密度区不向邻区扩散那么多狂犬狐——存在偏好的扩散方向。高密度口袋则相反。疫情可从高密度口袋前移出前阵(图 13.15(c) 中心图)。此聚焦效应可解释某些疫情在前阵之前出现的情况。这些效应也可能是图 13.3 中兽疫前沿曲折形状的成因。

英国自 1900 年起实际无狂犬病(一战后曾有小型疫情),主要靠严格检疫法与公众意识。海峡对面法国北部狂犬病近、私人船运增多使疾病短期内在英国出现几乎不可避免。英国若出现狂犬病将特别严重——英国城乡狐密度都高。另一令人担忧点是城市狐与猫的相容性(Macdonald 1980)。若无控制措施——实际必不会如此——疫情会在英国快速移动。可用模型估计狐群若被引入时疫情前阵的位置。

Macdonald (1980) 给出英国狐密度估计图(不含城市高密度口袋)。Murray et al. (1986) 把英格兰下半部覆盖网格,按 Macdonald 图给各方格赋密度。归一化到 [0, 1] 的等高线(密度 1 对应春季 2.4 成年狐/公里²)作图 13.16。模型基于全年度平均狐密度。引入疫情前种群在繁殖后达年度高点再回落到春季成年数;平均约为繁殖前后均值。公母比 1.2:1、母狐平均年产 3.7–4.2 仔(Lloyd et al. 1976)。故平均种群约为春季成年数 1.9 倍,归一化 1 对应承载容量 4.6 狐/公里²(最深阴影)。

用图 13.16 的承载容量与初始狐密度,假设疫情始发于南安普顿附近,对二维形式 (13.33) 数值求解。Table 13.1 的参数,扩散系数 200 公里²/年。在 Los Alamos National Laboratory 的 CRAY XMP-48 上数值模拟约 120 分钟。结果见图 13.17、13.18。图 13.17 给出每 120 天前阵位置。如此高狐密度下疫情迅速覆盖研究大部。4 年内前阵已实际到达曼彻斯特。图 13.18 序列显示与均匀情形一样,狂犬病例主要集中在前阵窄带内;易感狐群被疫情扫除并部分恢复后才开始下一波。图 13.18 显示约 7 年后南安普顿开始第二次疫情。

这些定量预测当然只是粗略估计。Macdonald (1980) 强调其狐密度图只是基于狐生态学知识的"有根据的猜测"。如前述,对狂犬狐行为所知不足使扩散系数无法精估;前阵速度可能为计算结果的一半到三分之四。模型忽略河流等地理因素——河流会成为疫情通道,沿岸顺行加速、垂直方向暂阻。但此相对简单的 SIR 模型为英格兰不受控疫情的进展提供了看似合理的定量首估,也为估计现实阻断宽度提供了方法,至少能严重阻碍疾病传播。

模型纳入了疾病与狐生态的许多显著特征。简单到能为除扩散系数外的所有参数给出相当可靠的估计(扩散系数给出一系列可能值)。模型对不同环境下疫情波行为的预测提供了对疫情空间传播与传播机制的某些定量洞见。例如:疫情空间传播的最初始原因到底是感染狐的扩散还是带病健康年轻狐的迁移,或是两者同等重要,这些未知。把某机制孤立出来可确定若该机制为主时疫情波的行为,再与欧洲大陆观察比较以判断是否为主导因素。结果显示感染狐的混乱游走足以解释当前疫情的大多数行为。研究以年轻狐迁移为主要原因的模型将是有趣的。已知一部分狐对狂犬病免疫。13.9 节把这种效应纳入模型框架。

模型与现有流行病学证据吻合较好,尽管扩散系数大小有不确定性。初始狐密度 2 狐/公里²下(与多数欧洲大陆报告密度相似),任何合理 D 选值下模型给出疫情前阵速度 25–65 公里/年,覆盖观察的 30–60 公里/年范围。前阵速度随狐密度增大,到临界密度时降为零。模型还预测第一次疫情后狂犬病约 5 年内基本消失然后再现,第二次较第一次弱——与欧洲多地现象吻合。模型另一有趣特征是高密度区疫情前移加强——可能解释"在主前阵之前"出现的疫情。

降低易感狐密度的条带可以阻挡疫情前阵、保护阻断前方未感染区。要高效应用此控制方法须有有效阻断区宽度的指示。图 13.14 给出无量纲估计。2 狐/公里²初始、减员 80% 有效时图 13.14 给阻断宽度 10–25 公里(依 D 而定)——与丹麦、瑞士部分有效阻断同量级。丹麦以 20 公里宽的密集控制条带加相邻 20 公里较弱条带。

应当用什么方法遏制疫情值得讨论。模型显示疫苗接种比毒气或毒杀更有效——前者限制感染狐扩散,后者促进扩散。Ontario 用浸有疫苗的鸡头证明相当有效——依赖于狐的清道夫行为。城市狐不必然如此(Stephen Harris 1988 私人通讯)。疫苗接种一般问题是某物种的疫苗接种水平可能诱发另一物种的疫情,红狐与灰狐间似乎如此。

狂犬病到达英国及其他无病区的概率不小。疾病传播与传播方式研究在它到达前就显得很重要。英国许多地区狐密度远高于大陆,疫情在英国可能进程不同。图 13.17、13.18 总结了对特定 D 与英格兰南部狐种群估计的模型预测。最令人不安的是疫情会以约 100 公里/年的速度迅速扫过中部。同样令人不安的是疫情前阵过去数年后再次出现——相对无疫情的时期会催生自满。

本章模型类型有更广适用性:害虫、killer bees(南美传播数据见 Taylor 1977)、动物、植物等的空间传播。

13.9 狐免疫对狂犬病空间传播的影响(Effect of Fox Immunity on the Spatial Spread of Rabies)

已知一定比例的狐天然免疫于狂犬病。Murray 与 Seward (1992) 量化其对空间传播的影响——本节简述他们对前述模型的修改(见原文献详情与更全的结果)。结果用现实免疫规模估计:免疫对初始疫情波传播速度影响小,但影响振荡尾的周期性疫情行为。模型也研究阻断宽度对疫情的影响并加入易感与感染狐的空间扩散。最后研究"阻断宽度是否依赖于减员或疫苗"这一假设。结果:除非空间扩散率大、免疫类显著增加,否则阻断宽度变化不显著。本节讨论其模型,因方法不限于狂犬病空间传播。也包括所有物种的扩散。

前述 SIR 模型假设所有狂犬狐死亡。事实是部分狐确实康复,康复者中一定比例形成免疫。Steck 与 Wandeler (1980) 多个实验研究表明约 2% 感染红狐形成免疫。野外狐免疫状态难评。Steck 与 Wandeler (1980) 数据提示疫情前阵过后存活狐中免疫的不到 8%。但美国灰狐(Urocyon cinereoargenteus)与红狐同时传播时,估计免疫狐比例可高达 20%——混合物种下模型不直接适用,可通过参数选取扩展。Wandeler (1987) 指出,临床病存活的可靠记录罕见与野外动物血清中常见狂犬病中和抗体形成对比——尚不清楚后者是否证明临床病存活。整体看免疫发展可能但对疾病传播影响小。研究把免疫类引入 13.5 节模型是值得的——预期小免疫比例下修正模型结果与原模型差异不大。发展此模型不仅教学上有意义,也确认上述信念、定量空间传播效应(随免疫类增加)以及免疫向仔狐传递的效应。

现把狐群分四类:易感 S、潜伏期 I、狂犬 R、免疫 Z。仍基于狂犬病毒 12–135 天的长潜伏期(临床表现正常不传播),临床期 1–10 天。

关键假设含 13.5 节 (i)–(vi) 并加:(i) 狂犬病并非必然致死。狂犬狐按人均率 α 死亡、按人均率 γ 康复形成免疫。γ 值由狂犬狐康复并形成免疫的百分比 p = γ/(α + γ) 给出。不形成免疫的康复者视为易感类,不单独处理。(ii) 狂犬与感染狐继续消耗环境、因非狂犬病原因死亡但健康后代可忽略。(iii) 免疫狐可有易感或免疫后代——研究两种极端情形:全为易感或全为免疫。

修正模型 (13.70):

∂S/∂T = (a − b)(1 − N/K) + a*Z − βRS

∂I/∂T = βRS − σI − [b + (a − b)N/K] I

∂R/∂T = σI − αR − γR − [b + (a − b)N/K] R + D_R ∂²R/∂X²

∂Z/∂T = γR + (a − a*)Z − [b + (a − b)N/K] Z

总种群 N = S + I + R + Z。免疫狐若全为易感后代则 a = a;若全为免疫后代则 a = 0。仅研究一维问题,重点考察加免疫类对 (i) 兽疫波速度、(ii) 主疫情后周期性疫情行为、(iii) 控制措施的影响。

兽疫传播速度:用同 (13.29) 的无量纲变量但加一组度量免疫的 (13.71),得到 (13.72):

∂s/∂t = ε(1 − n)s + (ε + δ)* z − rs

∂q/∂t = rs − (μ + δ + εn)q

∂r/∂t = μq − (v + d + εn)r + ∂²r/∂x²

∂z/∂t = vr + [(ε + δ) − (ε + δ)*] z − (δ + εn) z

n = s + q + r + z,(ε + δ)* = (ε + δ) 若免疫类全为易感后代,或 = 0 若全为免疫后代。p = γ/(α + γ)。其他维度参数值取 Table 13.1 与 D_R = 200 公里²·年⁻¹。

先看免疫狐有易感后代的情形。Murray 与 Seward (1992) 对 p = 2%、5%、10%、15%、20% 数值求解四类模型。三类与四类模型的初始波与周期性疫情形状变化很小(似图 13.9、13.10)。引入免疫群的效果为:(i) 初始波速度减小;(ii) 初始疫情中感染与狂犬群水平较低;(iii) 疫情发生时易感群减少不严重;(iv) 周期性疫情间隔缩短。前三点是直觉的;第四点由前三推出。这些效应随免疫百分比增加而更显著。13.5 节三类的渐近与数值结果:K = 2 狐/公里² 时初始波速约 51 公里/年、K = 4.6 狐/公里² 时约 103 公里/年。四类模型波速见表 13.3。

加免疫群(带易感后代)的主要效应在初始波的尾。在三类模型中,K = 2 狐/公里² 时易感狐需约 5 年才能恢复至足以发生二次疫情;K = 4.6 狐/公里² 时 11 年。四类模型中首次周期性疫情出现显著更早(表 13.4)。Murray 与 Seward (1992) 发现四类模型中初始疫情对易感群的减少不如三类模型严重,种群水平回升更快——这解释了周期性疫情间隔的缩短。

四类模型另一显著变化是波尾振荡的衰减增加。空间均匀情形((13.70) D_R = 0)有非平凡稳态由 (13.70) a* = a 定义。三类模型 ((13.70) Z = 0、γ = 0) 中 (13.30) 给出临界承载容量 K_T。Murray 与 Seward (1992) 数值估计四类模型中稳态与 K_T。免疫百分比增加时,数值解趋稳态值的时间缩短。

免疫狐全为免疫后代情形:取相同 K、p 数值求解。免疫后代对初始疫情传播无影响——波速同表 13.2(数值验证)。原因:免疫群到初始疫情过后才存在。主要效应仍在初始波的尾,与"易感后代"四类与原三类模型差异显著。初始疫情过后,免疫狐占种群相当大比例,对后续疫情有显著衰减效应。K = 2 狐/公里² 时第二次疫情出现时间大于三类模型且随免疫百分比增加(表 13.5)。K = 4.6 狐/公里²、2% 免疫时二次疫情在 18 年后;5% 免疫时 21 年后。K = 4.6 狐/公里² 时 p > 5% 的所有图都很类似。

这些差异可由 (13.70) a* = 0 时的稳态解释。此时无物理上合理的正 I、R 稳态;模型需疫情消失后才有稳态,此时 (13.70) 化为剩余总狐 S + Z 的 logistic 增长律。稳态 S₀ + Z₀ = K,S₀、Z₀ 相对值由初值决定。系统接近稳态时间依赖 K 与免疫百分比。K 或 p 增大时更快接近稳态。

无论是易感或免疫后代,四类模型对初始疫情传播影响小。波速只在假设高免疫百分比时显著变化。四类模型主要效应在初始波的尾;效应在"易感后代"与"免疫后代"间差异很大。

只设全为免疫或全为易感后代是简化。实际中免疫与易感狐杂交,后代免疫比例取决于免疫与易感群的相对大小。模型可加入遗传方程但比三、四类复杂得多。模型预测低自然免疫率(2%–5%)对波速影响小,这与"缺乏免疫是狂犬病在狐间传播的助因之一"的观察一致(Blancou 1988)。K = 4.6 狐/公里² 时低免疫率最显著效应是减小二次疫情前时间。

免疫对与狂犬病"阻断"相关的控制措施的影响:传染病模型的主要用途之一是评估各种控制策略遏制疾病。如 13.6 节所述,一可能方法是在初始波前引入"阻断"——易感狐减员至承载容量以下的区,使兽疫波不能传播。可用模型再估所需阻断宽度。实践中阻断可由减员产生(丹麦用此法——生育季加强狩猎与狐穴毒气 Wandeler 1987)或疫苗接种(瑞士成功 Wandeler et al. 1987)。可用模型比较这两种方法。如前所述减员的潜在困难是降低的种群密度可能鼓励狐扩散入该区从而降低阻断效果。可在三类模型前两方程加入扩散项 D_S ∂²S/∂X² 与 D_I ∂²I/∂X² 来模拟此扩散。扩散项平均阻断内外种群。Yachi et al. (1989) 研究三类模型加易感与感染狐扩散(同一扩散率),发现这可使兽疫波速度显著增大。Murray 与 Seward (1992) 取易感与感染狐的扩散率小于狂犬狐的。疫苗接种可在四类模型中以初始免疫群 Z 非零表示。仍可通过数值求解(阻断设 x = 0 至 x → ∞、承载容量降、初值减员)估所需阻断宽度;感染狐总数 < F 狐·公里⁻¹(F 远小于临界承载容量)的时间 t_c(F);阻断宽度 x_c 为感染狐密度为 x = 0 处小分数 m 的点(同 (13.51))。

降低狐种群密度有两种方式:减 K 或减初值 S。减 K 对应持续控制方案,13.6 节采用。阻断宽度关于阻断内 K、参数 d、F、m 的依赖见图 13.14。仅减初值 S 是"一次性"建立阻断(而非持续控制)。直观上后者不如前者有效,但计算结果示两种情形下阻断宽度相近。Murray 与 Seward (1992) 先考虑用三类模型(含所有狐的扩散、不含免疫)效果。(13.70) Z = γ = 0 加 D_S、D_I 项得 (13.73)。D 估法同 13.5 节:领地大小 5 公里² × 出生率 1 狐/年 = D_S = D_I = 5 公里²/年。最可能低估——Garnerin et al. (1986) 基于离散模型估幼狐扩散距离限 8 公里。数值结果对 5、20、50、200 公里²/年计算。

Murray 与 Seward (1992) 用 13.6 节方法(F = 0.5 狐·公里⁻¹、m = 10⁻⁴)重算无扩散的阻断宽度作为参照,再设计两倍于算出的阻断宽度作数据区,置于两个正常密度区之间,观察兽疫波流入阻断。扩散下狂犬病最终仍会穿越阻断但比无阻断时慢得多。两种控制情景:K = 2 狐·公里⁻²,一次性减阻断内易感群;持续减阻断内承载容量。两种情形下阻断宽度相似;D_S = D_I = 200 公里²/年时减初值易感群似更有效。两种策略同样有效——兽疫波穿越阻断的时间大致相同。积分时间约一年,阻断在疫情前约三个月建立。一次性策略若在疫情到达更早前建立效果较弱。一次性策略下当阻断内易感群为 1.2 狐·公里⁻²时不能把总感染狐减至 F = 0.5 狐·公里⁻¹。Murray 与 Seward (1992) 给出各种 D 与 K 的多种阻断场景图。

疫苗接种效果以四类模型(含免疫狐 Z)模拟——Z 代表疫苗接种狐。仍用 13.6 节方法估阻断宽度。解模型方程在 S(X, 0) 减的初值下传播,初始疫苗狐 Z(X, 0) = K − S(X, 0)。若疫苗狐有易感后代,阻断实质由"一次性"方式建立——疫苗接种仅一次性用于初始群。疫苗狐有"免疫"后代时可粗略模拟持续疫苗方案。图 13.19 比较三类(减员建阻断)与四类(疫苗接种建阻断)模型 K = 2 狐·公里⁻² 下的基本阻断宽度。四类模型阻断宽度普遍较大。疫苗狐有易感后代情形与有易感狐扩散入阻断类似。如图 13.20 所示,三类模型 D 在 20–50 公里²/年时阻断宽度类似四类模型。持续疫苗方案效果为何不如减员不太清楚。如上所述,四类模型初始疫情中感染与狂犬群较低;但这些群在该情形下衰减比三类模型慢。结果是总感染群降到 F 以下所需时间更长,疫情在减员易感群中稍多扩散,算出的阻断宽度更大。难判断这种表观扩散多少来自微分方程性质、多少来自免疫狐对环境的实际效果。持续疫苗方案有两优势:阻断宽度稍小;可在更高狐密度下建阻断。

K = 4.6 狐·公里⁻² 时建阻断更难。疫苗狐有易感后代时不能在易感狐密度高于 0.4 狐·公里⁻² 建阻断。设"免疫"后代、减员密度 1.38 狐·公里⁻²、F = 1.5 狐·公里⁻¹ 时阻断宽度 36 公里;F = 0.5 公里⁻¹ 时此密度下建不了阻断。

四类模型中加天然免疫狐,阻断宽度减小但 2–5% 比例下变化小。20% 天然免疫时(无论易感或"免疫"后代)阻断窄 5–10 公里。图 13.20 给出"免疫"后代情形(持续控制方案)。

考虑这些数学模型仅给出阻断宽度的粗略估计,作者给出以下观察:建阻断最有效策略是执行减员方案。只有当减员导致相邻狐显著扩散入阻断区时疫苗接种才更有效。在高承载容量区建阻断困难。持续控制方案一般给出稍小的阻断宽度且可在比"一次性"建阻断更高的狐密度下建阻断。

本章主要关注疫情的空间传播。疫情已在某社区时存在的重要且有趣的控制策略问题。一有趣且实用的模型是 Frerichs 与 Prawda (1975) 提出的——处理哥伦比亚某城区的城市狂犬病。卷 I 第 10 章讨论了他们的方法在牛结核病与牛、獾交互作用(提供疾病储存宿主)的修正应用。

本章模型类型有更广适用性——害虫、killer bees(南美传播数据见 Taylor 1977)、动物、植物等的空间传播。

若干注意事项:13.6 节关于"阻断区种群灭绝意味着什么"的讨论在本节建模中同样适用。连续模型有助于理解疾病传播与空间动态并给出有用的定性常是定量预测,但仍可被合理批评。例如 Mollison (1991) 对确定性连续与随机模型提出相关观点。即使这里讨论的相对简单模型若加入随机性将复杂数量级。多数情况下两种方法的区别是"能否对调查情况做点什么";研究随机性效应无疑会很有启发。离散模型可更逼真(虽然参数估计更困难)。狐以家族群生活于独立领地,感染更可能在整个家族内传播。繁殖离散非连续等。连续空间离散时间模型可纳入这些效应。

建模总有"简洁与可估参数"与"含更多方面但参数难估难解难释"之间的权衡。一个明显缺陷是疫情后波后存在疾病"储存宿主"而计算种群基本在感染群灭绝水平——离散模型也有此问题。相关问题是假设平均潜伏期 1/σ 而实际为分布的潜伏期(变化大)。英格兰犬只检疫中有 6 个月后才发病的极罕见案例。最近 Fowler (2000) 重新检视 13.6 节模型并研究这两个方面(潜伏期分布与灭绝),证明潜伏期分布的引入可解释灭绝为何不发生。他进一步得出感染狐最低密度的渐近估计。故即使连续模型也能纳入灭绝。

本章个人批注

本章是卷 I 第 10 章(非空间流行病模型)向真实空间流行病学的"应用层"扩展。但 Murray 的写法非常克制——13.1 节先给一个最简单的双物种扩散 SIS 模型,把 13.1 那个 λ < 1 的阈值与卷 I 10.2 节"无空间时"的阈值在形式上对照、复述 Fisher–Kolmogoroff 的最小波速结论,这种"先简单后复杂"的处理是 Murray 一贯的思路。13.4 节是这一思路的"中段"——把扩散仅给 I、S 不动,得到一个"易感者存活比例 σ 由 σ − λ ln σ = 1 隐式决定"的结果,是少见的从行波解中解出"过波后剩余"这一可观测量。13.5 节引入长潜伏期形成的三物种 SIR 模型后,行波前沿的解析重点变成"特征值 λ₃、λ₄ 的虚部对应振荡尾波长"——这是把分析特征值的实部→波速、虚部→波长的标准技巧做了一次。读到这里我反复感受到的是:作者把 13.5 节给出的"三物种可约化为二物种 + 标度 q(ξ) ∼ dr(ξ)/μ"作为一个非平凡的观察而不仅仅是技术细节——它意味着在合理参数范围内"感染但未发病"这一类的具体动力学对宏观行波形态并不重要。13.6–13.7 节把模型用于实际控制("狂犬病屏障"),其解析核心是把"阻断区内感染狐随时间衰减到 F 狐/公里"的时间与"在阻断内按扩散衰减的指数宽度"分开估计,然后给出 x_c ∼ √(d − μ) ln(1/m) 的简洁形式((13.68)),其中 m 是阻断对侧相对源侧的衰减比——参数 d、μ 都是无量纲的,d − μ 出现是因为感染狐既按 σ 流出(转狂犬)也按 α 死亡(致死)。13.8 节"英格兰实验"是本书很特别的"定量预测"段落——Murray 直接给出一张南英格兰 0–1 归一化的狐密度图并把方程在 CRAY XMP-48 上跑了 120 分钟得到 4 年内前阵到达曼彻斯特的图;他警告这不是预测而是"若无控制"情形下的估计,但他同时强调英国狐密度(K = 4.6 狐/公里²)远高于欧洲(K = 2 狐/公里²),导致周期性疫情间隔从 5 年降到 11 年——这一对比是 13.5 节表 13.2 中"周期随 K 增大而减小"的具体兑现。13.9 节把 Murray 与 Seward (1992) 的四类模型(加免疫 Z)作为对前三类模型"在合理小免疫比例下结论不变"的检验——这是少见的"先做模型再检验自己的预测"的做法。其中 20% 免疫的极端情形作为上限是为了说明模型在哪些参数下结论失效。结尾的"若干注意事项"是 Murray 风格的诚实标注:他承认 Fowler (2000) 把分布潜伏期纳入后可以解释"灭绝"这一长期被忽略的难题——这是少见的作者在自己书的章节中给"自己的工作遗留问题"做正面的处理。整章读完后最深的印象是:作者把"控制策略"作为建模的目标之一,但从不混淆"模型的预测"与"应做什么"——他反复指出"模型预测不能替代防控决策"。

与上下章的衔接(一段话)

上一章(第 12 章)把"激活–抑制"作为神经图样形成的主线,从卷 I 第 2 章已建立的色散关系技术出发,在眼优势条纹、视幻觉、贝壳图样三个看似毫不相关的生理场景中复用同一数学骨架。本章(第 13 章)则把"反应扩散 + 阈值"作为流行病空间传播的主线,其数学基础(即 13.1 节的双物种扩散 SIS 模型与最小波速)正是第 1 章 Fisher–Kolmogoroff 方程与卷 I 第 10 章非空间 SIS 模型的直接合成;与第 12 章不同的是,本章把模型从"产生空间图样"的视角切换到"产生空间行波"——同一套线性化→特征值→实部→波速、虚部→波长的技巧在这里用得最为直接。13.4 节用 I 扩散、S 不动的简化 SIS 模型得到行波解中"过波后存活比例 σ 由 σ − λ ln σ = 1 隐式决定",是第 12 章所讨论的激活–抑制卷积方程(其解的形态由核结构决定)所不具有的"显式可解"特征。13.5 节进一步把模型升级到三物种 SIR,得到的"振荡尾"是第 1 章 activator–inhibitor 振荡与卷 I 第 10 章 SIR 振荡的合并:基本传染率 R₀ = 1/λ 与卷 I 10.2 节、卷 II 13.1 节完全同源;K_T 表达式 (13.30) 与第 1 章 (1.3) 的 steady-state 分析形式上对应。13.8 节的"英格兰实验"是本书很独特的"二维异质介质中的行波模拟"段落——把 Murray et al. (1986) 在 CRAY XMP-48 上跑出的图 13.17、13.18 直接嵌入本章。这一节事实上把本章的"二维问题"明确提出,但作者并未系统讨论二维行波的选择性稳定性(卷 I 第 2 章 2.5 节讨论过平面波的横向稳定性),只是把它当作"工具"使用。13.9 节是本书中少有的"加一类、把模型变复杂、再与原模型比较"的研究范式——把免疫 Z 加入并检验结论在合理小免疫比例下不变,这与第 12 章反复强调的"机制不唯一"态度相呼应。下一章(第 14 章)则把注意力从"动物/人之间的疾病"转向"动物之间的社会行为"——狼的领地性作为空间模式形成的另一种机制,把第 1 章的捕食者–猎物模型与本章的 SIS 模型都置入到"群体动力学 + 空间"的更大框架中。