引言
歷史上的大瘟疫不僅是人類(lèi)文明的重大考驗(yàn)，更深刻影響了社會(huì)結(jié)構(gòu)、經(jīng)濟(jì)發(fā)展和政治格局的演變。公元前 430 年雅典瘟疫的暴發(fā)削弱了雅典的政治和軍事力量，直接影響了伯羅奔尼撒戰(zhàn)爭(zhēng)的走向。中世紀(jì)的黑死病導(dǎo)致歐洲人口銳減，封建制度逐漸崩潰，加速了資本主義的興起。
每一次瘟疫的暴發(fā)都推動(dòng)了疾病管理、社會(huì)組織和科學(xué)技術(shù)的進(jìn)步，促進(jìn)了公共衛(wèi)生體系和醫(yī)療保健的發(fā)展。研究歷史瘟疫的傳播模式，不僅有助于我們更好地理解疾病本身的傳播規(guī)律，還能為未來(lái)潛在疫情的防控提供寶貴的經(jīng)驗(yàn)和科學(xué)依據(jù)。本文將通過(guò)經(jīng)典的 SIR 模型分析傳染病傳播的基本機(jī)制，給出流行病學(xué)關(guān)鍵參數(shù)的計(jì)算方法，并以倫敦大瘟疫期間亞姆村的鼠疫為例，模擬鼠疫的傳播過(guò)程，探討其傳播規(guī)律及防控策略。
數(shù)據(jù)
亞姆村（Eyam），位于英格蘭德比郡，因 1665-1666 年鼠疫期間的自我隔離而被稱(chēng)為“鼠疫之村”。為防止瘟疫擴(kuò)散至鄰近地區(qū)，村民自愿隔離，付出巨大代價(jià)：全村 350 人中僅 83 人幸存（圖 2）。
亞姆村的疫情源于 1664-1666 年倫敦的大瘟疫。1665 年 9 月初，村里的裁縫收到了一包從倫敦寄來(lái)的含鼠疫跳蚤的布料，不久便去世，此后疫情迅速在村內(nèi)蔓延。疫情在 1666 年 5 月似乎有所減緩，但到了 6 月再次暴發(fā)。以 1666 年 6 月 18 日為時(shí)間起點(diǎn)，亞姆村健康和感染人數(shù)見(jiàn)表 1。
模型
本文要介紹的 SIR 模型[2]是流行病學(xué)中最基礎(chǔ)的傳染病數(shù)學(xué)模型之一，用于預(yù)測(cè)傳染病在人群中的傳播動(dòng)態(tài)。該模型將人群劃分為三類(lèi)：易感者（Susceptible）、感染者（Infectious）和移除者（Removed），三類(lèi)人群的定義分別如下：
易感者（S）：指未得病者，但缺乏免疫能力，與感染者接觸后容易受到感染。
感染者（I）：當(dāng)前正攜帶并能夠傳播病菌或病毒的人群，他們能將疾病傳給易感者。
移除者（R）：因康復(fù)獲得免疫力或因病死亡而從感染鏈中被移除的人群。這部分人不再參與感染和被感染過(guò)程。
如圖 3 所示，用 、 和 分別表示三類(lèi)人群的數(shù)量，SIR 模型涉及兩個(gè)主要過(guò)程：傳播和移除。
圖 3: SIR 傳染病模型示意圖
傳播過(guò)程會(huì)使感染者數(shù)量 增加，易感者數(shù)量 減少，而移除過(guò)程則會(huì)使移除者數(shù)量 增加，感染者數(shù)量 減少。具體過(guò)程如下：
傳播：假設(shè)在所有個(gè)體均為易感的人群中，一個(gè)感染者平均每天造成 個(gè)新病例（傳染率）。那么 個(gè)感染者每天會(huì)造成 個(gè)新病例。如果人群并非完全易感，則預(yù)期的新病例數(shù)應(yīng)為 ，這里的 代表接觸到的易感者的比例。新增感染者數(shù) = 新減少易感者數(shù) = 傳染率（ ） 感染者數(shù)（ ） 接觸者中易感者的比例（ ）。
移除：感染者經(jīng)過(guò)一段時(shí)間后會(huì)因康復(fù)或死亡而變?yōu)橐瞥?，這個(gè)轉(zhuǎn)變的速率依賴于感染者的恢復(fù)率或死亡率。假設(shè) 是感染者處于感染狀態(tài)的平均天數(shù)，則每天有 = 1/ 的感染者轉(zhuǎn)變?yōu)橐瞥摺Ｐ略鲆瞥邤?shù) 新減少感染者數(shù) 感染者移除率（ ） 感染者數(shù)（ ）。
如果用 表示時(shí)間（單位：天），則上述兩個(gè)過(guò)程造成系統(tǒng)中三種人群數(shù)量隨時(shí)間的演化可用以下微分方程組表示： 其中， 和 分別表示傳染率和感染者移除率。將以上三式相加可知，在 SIR 系統(tǒng)中，總?cè)丝跀?shù)是一個(gè)不隨時(shí)間變化的常數(shù)：為了求解這個(gè)微分方程組，需要給定初始條件。通常假設(shè)在初始時(shí)刻 = 0，有一部分易感者 和感染者 ，而移除者 為零。因此，初始條件可以表示為： 如果已知 、 ，以及傳染率 和感染者移除率 ，就可以通過(guò)對(duì)以上微分方程組進(jìn)行數(shù)值積分來(lái)得到三種人群數(shù)量隨時(shí)間的變化。
圖 4: 基本傳染數(shù)為 2 的示意圖
基本傳染數(shù)（又稱(chēng)基本再生數(shù)） 是流行病學(xué)中重要的參數(shù)，它表示在全易感人群中（沒(méi)有任何預(yù)防手段介入并且所有人對(duì)此病原體沒(méi)有免疫力的情況下），一個(gè)感染者在其感染周期 內(nèi)平均能傳染的人數(shù)（圖 4），即 注意，為避免與第 0 天的移除者數(shù)量 混淆，本文采用花體 表示基本傳染數(shù)。通常 的值越大，表明疾病的傳播潛力越強(qiáng)，控制疾病越困難。當(dāng) < 1 時(shí)，每個(gè)感染者傳染給不到一個(gè)人，傳染病將逐漸消退。如果 > 1，傳染病將以指數(shù)級(jí)速度蔓延，可能發(fā)展為大規(guī)模流行病。例如，COVID-19 奧密克戎變異株的 約為 7[3]，遠(yuǎn)大于 1，顯示其具有高度傳染性。然而，這種情況通常不會(huì)無(wú)限持續(xù)，因?yàn)橐赘腥丝趯⒁蚋腥竞笏劳龌颢@得免疫而逐漸減少。若 = 1，傳染病將在人群中持續(xù)存在，形成地方性流行。
在 SIR 模型中，判斷一種傳染病能否持續(xù)流行不僅取決于基本傳染數(shù) ，還與當(dāng)前易感者人數(shù)有直接關(guān)系。可以將 SIR 模型中描述感染者變化的方程重寫(xiě)為如下形式： 其中， 為臨界易感者數(shù)。上式表明，當(dāng)易感者人數(shù)大于臨界易感者數(shù)（ > ）時(shí)， > 0，感染人數(shù)將繼續(xù)增加；相反，當(dāng) < 時(shí)， < 0，新增的感染者數(shù)小于移除的感染者數(shù)，感染人數(shù)將不斷減少。感染者的最大數(shù)量 出現(xiàn)在 = 處，此時(shí) = 0。
還可以通過(guò)分析相軌跡來(lái)研究 SIR 模型的特征。由于模型中的三個(gè)微分方程是解耦的，只需要考慮微分方程組中的前兩式。將兩式相除可以消去 ，得到： 對(duì)上述方程積分，可以得到 - 相平面中的軌跡： 常 數(shù) 相平面中的軌跡如圖 5 所示，相軌跡描述了易感者 和感染者 之間的關(guān)系。初始時(shí)刻， 所有初始值 和 滿足 + = 。對(duì)于 > 0，則有 + < 。隨著時(shí)間推移， 逐漸減少， 先增加后減少，最終趨于零。相軌跡的形狀取決于初始條件 和 。當(dāng) > 時(shí)，感染人數(shù) 會(huì)先增加到一個(gè)峰值，然后逐漸減少，表明發(fā)生了疫情；而當(dāng) < 時(shí)，感染人數(shù) 單調(diào)減少，表明疫情不會(huì)發(fā)生。
圖 5: 在 相平面中的相軌跡
以上相軌跡公式對(duì)任意時(shí)間 都成立，因此還可以用來(lái)估計(jì)基本傳染數(shù) 或傳染率 的數(shù)值： 特別地，當(dāng) （疫情結(jié)束）時(shí)，顯然有 = 0，令 = ，則上式變?yōu)閇4]：以上兩式都可以估計(jì)基本傳染數(shù) 或傳染率 ，不過(guò)后一式對(duì)于正在發(fā)生或剛剛開(kāi)始的疫情并不適用，因?yàn)樗枰?的值，而前一式則沒(méi)有這個(gè)限制。
如果發(fā)生了疫情，我們可能想知道疫情的嚴(yán)重程度。在前文的分析中，已經(jīng)知道感染者的最大數(shù)量 出現(xiàn)在 = 處。由相軌跡公式可知，當(dāng) = 時(shí)， 對(duì)于任何初始值 和 > ， 相軌跡從 > 開(kāi)始，從圖 5 中可以看到 從 開(kāi)始增加，直到達(dá)到 ，這意味著會(huì)出現(xiàn)疫情。相反，如果 < ，則 從 開(kāi)始減少，不會(huì)暴發(fā)疫情。
結(jié)果
在當(dāng)前亞姆村的鼠疫問(wèn)題中，初始三種人群的數(shù)量分別為 = 254、 = 7 和 = 0（見(jiàn)表 1），人群總數(shù) = 261。研究表明，人類(lèi)鼠疫的潛伏期最長(zhǎng)為 6 天，病程約為 5.5 天。假設(shè)整個(gè)感染期（從感染到死亡）平均為 11 天[5]，因此感染者移除率 = = 1/11。為了求解 SIR 模型的微分方程組，還需要估算出傳染率 的值，圖 6 是應(yīng)用表 1 中不同時(shí)間點(diǎn)的數(shù)據(jù)估算的結(jié)果。從圖中不難看出，從第 31 天的數(shù)據(jù)開(kāi)始，估算的 值都穩(wěn)定在 0.145 附近。接下來(lái)，本文將分別從疫情正在發(fā)生和已經(jīng)結(jié)束兩種情景出發(fā)，模擬三種人群數(shù)量隨時(shí)間的變化。
圖 6: 不同時(shí)間點(diǎn)數(shù)據(jù)估算的 值疫情正在發(fā)生的情況
假設(shè)疫情正在發(fā)生，且僅有早期的數(shù)據(jù)。例如，應(yīng)用第 46 天的數(shù)據(jù)，可以估算出亞姆村鼠疫的傳染率、基本傳染數(shù)和臨界易感者數(shù)分別為： 由于 > 1 且 = 254 > ，這意味著疫情將繼續(xù)擴(kuò)大。在獲得微分方程模型所需的參數(shù)后，可以對(duì)模型的微分方程組進(jìn)行數(shù)值積分，得到三種人群數(shù)量隨時(shí)間的變化，結(jié)果如圖 7 所示。
圖 7: 三種人群數(shù)量隨時(shí)間的變化（ ）
圖中實(shí)心點(diǎn)表示實(shí)際數(shù)據(jù)點(diǎn)，曲線表示模型預(yù)測(cè)，空心點(diǎn)表示模型預(yù)測(cè)的第 108 天三種人群的數(shù)量。模型結(jié)果顯示，當(dāng)易感者數(shù) 下降到約 155.9 以下時(shí)，感染者數(shù) 才開(kāi)始降低。從圖中可以看出，模型預(yù)測(cè)與實(shí)際數(shù)據(jù)非常吻合：在疫情早期，易感者人數(shù)呈現(xiàn)先緩慢減少，然后加速下降，最終逐漸趨于穩(wěn)定；感染者數(shù)則是先逐漸增加，達(dá)到高峰后開(kāi)始減少，并最終下降為 0；移除者數(shù)則從 0 開(kāi)始，隨著疫情的發(fā)展不斷增加，增速先加快后減緩，最終也趨于穩(wěn)定。
疫情已經(jīng)結(jié)束的情況
對(duì)于已經(jīng)結(jié)束的疫情，亞姆村最終的易感者人數(shù)剩余 = 83。據(jù)此可估算出傳染率、基本傳染數(shù)和臨界易感者數(shù)分別為： 以上參數(shù)估計(jì)值與基于第 46 天數(shù)據(jù)給出的結(jié)果差異較小。由此求解得到三種人群數(shù)量隨時(shí)間的變化如圖 8 所示。
圖 8: 三種人群數(shù)量隨時(shí)間的變化（ ）
與應(yīng)用疫情早期數(shù)據(jù)估計(jì)參數(shù)得到的結(jié)果（圖 7）對(duì)比，可以發(fā)現(xiàn)，應(yīng)用疫情結(jié)束時(shí)的數(shù)據(jù)給出的結(jié)果與實(shí)際情況更加吻合，但兩者差異并不顯著。這表明，即便是在疫情早期，也能通過(guò)數(shù)據(jù)分析得出較為準(zhǔn)確的參數(shù)估計(jì)。
結(jié)論
通過(guò)對(duì)亞姆村鼠疫的歷史數(shù)據(jù)進(jìn)行分析，本文應(yīng)用 SIR 模型成功模擬了瘟疫的傳播過(guò)程。通過(guò)微分方程模型有效地估算了亞姆村鼠疫的基本傳染數(shù) 和臨界易感者數(shù) ，并模擬了易感者、感染者和移除者隨時(shí)間的動(dòng)態(tài)變化。模型結(jié)果與歷史數(shù)據(jù)具有較高的一致性，揭示了傳染病在相對(duì)封閉環(huán)境中的傳播規(guī)律。
本文的研究不僅驗(yàn)證了 SIR 模型在歷史疫情分析中的有效性，還為未來(lái)類(lèi)似疫情的防控提供了理論依據(jù)。通過(guò)量化傳染病的傳播參數(shù)，可以更好地理解疾病的傳播機(jī)制，為公共衛(wèi)生政策的制定提供科學(xué)支持。
參考文獻(xiàn)
[1]D Sulsky. Using real data in an sir model. The University of New Mexico, Albuquerque, USA, 2012.
[2]Wikipedia contributors. Compartmental models in epidemiology, 2020.: https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology
[3]Wikipedia contributors. Basic reproduction number, 2024.: https://en.wikipedia.org/wiki/Basic_reproduction_number
[4]Fred Brauer. Compartmental models for epidemics. Department of Mathematics, University of British Columbia Vancouver, BC V6T 1Z2 Canada, 2008.
[5]Graham F Raggett. Modelling the eyam plague. Bull. Inst. Math. and its Applic, 18(221-226):530, 1982.
來(lái)源：數(shù)學(xué)模型
編輯：未
轉(zhuǎn)載內(nèi)容僅代表作者觀點(diǎn)
不代表中科院物理所立場(chǎng)
如需轉(zhuǎn)載請(qǐng)聯(lián)系原公眾號(hào)
熱門(mén)跟貼