原文出處 Cheng-Chih Richard Hsu (所有圖文版權皆屬原著作權人所有)
「老師,請問武漢肺炎疫情什麼時候才會結束?」
這是昨天我學生問我的一句話。
我當下回應是,我不是公衛和流行病學專業,所以我無法回答。不過學生問的這個問題,其實從化學課本中教過的知識,或許就可以嘗試來回答。
時間先拉回到小年夜武漢封城那天,當天官方公布的武漢肺炎 (2019-nCoV) 單日確診人數還只有兩百多人(包括中國國家健康衛生委員會和 WHO 的報告),轉眼到初三在一天內的單日確診數已經達到一千七百多人,而在中國境內的總確診病例也在前天 (1/28) 突破 2003 年 SARS 的數字(註一)。
因為確診數字跳動非常快,而且不斷創新高,相信所有人都想知道:
「到底有多少人被感染,以及疫情何時會結束?」
對於還沒發生的事情,除了上帝本人以外沒有人知道正確答案。但是針對這問題,流行病學家參考過往的經驗,早已提出各種數學模型,來去描述流行病傳播時,隨著時間變化,被感染的人數可能是多少。以 SARS 為例,2003 年那年就有至少十篇的統計模型被提出來,隔年還有人寫了文獻回顧整理(註二)。這些模型令人眼花撩亂,非相關專業背景的人很難知其究竟。
剛好我學生也隨口問了這個問題,因此在這邊,我從化學反應的碰撞理論作為起點,提出一個有高中化學程度就能大致理解的動力學模型,以較直觀的方式來去描述傳染病的傳播方式,以及被感染人數隨時間的變化。接著再嘗試透過過去的經驗和其他專家的意見,來去預測這次武漢肺炎可能的走向(為了避免我被疾管局找去喝咖啡,請先詳閱文章結尾的聲明)。
*以下容我先上一點反應動力學的課,沒興趣的人可以直接跳到後面疫情預測的部分。
【動力學模型】
首先我們要先想像,在一個空間裡面,有兩種人,一種是健康的人(用A表示),另一種是不幸被感染的人(用B表示),A和B在這個空間中會隨機亂走,因此有一定的機率會彼此相碰。當A碰A或B碰B時,並不會有反應,但是當A碰到B的時候,則有一定的機率,被感染的B會將健康的A也變成B(將病原傳給A)。
我們可以把上述的情況用一條化學反應式來代表(反應不可逆):
A(一位健康者)+B(一位被感染者)→ 2 B(兩位被感染者)
如此,這個反應的速率,會跟空間中A和B的數量都分別呈正相關,因此當一群A裡面忽然有了B,隨著時間演進,整個群體會逐漸都變成B。而被感染者B的數量增加最快的時候,並不是最剛開始出現B時,而是A和B差不多是一樣多的時候,因為A跟B相碰的頻率會最高。因此你可以想像,起初被感染的B人數不多時,人數增加得並不快,接著逐漸到邁向高峰,而當跨過最高峰後因為可以被感染的人變少了,被感染人數B所增加的速率也就會迅速降低。
理論上,在不考慮健康者A自然死亡的情況下,A+B的數量會等於區域內的總人口數。但是若在疫情剛開始的時候進行隔離防疫,相當於限制B可以觸及A的機會和範圍,因此A+B的總數量會隨著防疫的落實而大幅減少,形成一個初始的邊界條件,決定了最後的總感染人數。
這樣的一個極為簡化的模型,它的每日「新增」的被感染者(delta [B])就是化學上的二級反應式,寫成數學式子就是:
delta [B] = 「當日患者數」 x 「當日健康者數」 x 常數 ([A][B] Const)
用這樣的式子就可以去計算,隨著疫情演進,在不同時間點被感染者的增加速率。
但是有幾點需要特別注意。第一是這個模型假設一旦被感染就具有傳染力,即使是在潛伏期也是。第二則是因為從感染後進入潛伏期到出症狀,以及檢體送樣後到確診,都需要一點時間,因此確診日期與實際被感染日中間會有誤差,以這次 2019-nCov 武漢肺炎來說,大概抓個 10-15 天的延遲。第三,為了方便起見,我們假定每個人的延遲時間相同,因此直接以官方統計的「確診日期」來套用到模型上面,用來估算每日新增「確診數」。第四,先前的 SARS 或 MERS 都證實存在有少數的超級傳染者,對疫情傳播有關鍵作用,但是在這個簡易的模型中我們假設每個人的傳播能力是一樣的。
理論建立好了,再來就是要去驗證。
【過去案例探討】
本來這個假期給自己訂下的工作是要趁著難得的長假,把手邊幾篇學生已經寫好的初稿改一改,但是從小年夜武漢封城開始,思緒就被疫情牽動,索性回去把自己以前的一些手稿和官方公布資料調出來讀。找了幾個大家比較知道的,像是 2003 年的 SARS、2014-2016 的台灣登革熱疫情、前幾年的伊波拉,以及台灣近期的幾次流感大爆發,大致上都能 fit 進去統計數字中(左圖紅色的校正曲線就是用模型預估的理論值,藍色點為實際值)。
左圖選了幾個代表性的疫情貼給大家看。數據來源包括 WHO 網站和台灣傳染病統計資料查詢系統(註三);這邊要特別稱讚台灣傳染病統計資料查詢系統,數據非常詳盡,而且資料一抓就有,比較可惜的是它的資料多是以周作統計,不像 WHO 很多疫情可以做到逐日更新。
從圖上來看,不管是登革熱或是 SARS,每周/每日的新增人數都是兩側低,中間高,然後右翼(疫情減緩時)數字下降的會比較陡峭。SARS 疫情因為是逐日登載,所以上下跳動的幅度比較大;而登革熱的數據因為是每周紀錄,相當於做了七點平均,因此平滑許多。在登革熱的數據 R^2 值甚至都可以超過 0.95,可以很大程度解釋傳播機制(關於登革熱疫情我在 2015 年就曾討論過,當時估算時沒考慮不同區疫情的長短以及防疫工作的影響,因此那次估出來的數字與實際的數值差了好幾倍,這次的預測會一併考慮此效應)。
為了方便大家理解和比較不同傳染病間疫情的差別,我定義校正曲線上半高寬所橫跨的期間,為疫情的「高峰期」。在數據點半高寬碰觸的兩邊,恰好是確診數增加速率為快的時候(切線斜率最大),而同時也可作為一個比較疫期長短的基準。
譬如 2015 年台南登革熱高峰期(約六周)就遠比 2014 年的高屏區登革熱疫情(約十一周)要短;而同樣是 SARS,在香港疫情高峰期約 40天(3/18~4/27),但在台灣就大概只有 18 天(從圖上也可以清楚看出來,儘管當年 4/24 和平醫院就封院了,但是從事後數據看來,真正的高峰要到 5/10 才來,而此時香港的疫情已經結束的差不多了)。這個高峰期我們在後面預測本次武漢肺炎未來走勢時也相當重要,請大家先牢記。
至於總確診數,就是整個區線底下的線下面積,因此曲線涵蓋的範圍愈大,感染人數愈多。特別注意的是,圖上的縱軸代表的是單位時間(天/周)內「新增」的個案數,因此這裡說的疫情結束是指感染人數不再新增,而不是感染人數歸零。即使感染人數不再新增,可能還是存在大量感染人數,和新增死亡人數。
【武漢肺炎疫情走勢預判】
回頭再來看看此次的 2019-nCoV 武漢肺炎。
在用這個模型作數據推估前,我們先看一下其他流行病學專家怎麼說。
各家統計的 2019-nCoV 感染人數,從幾萬人到上億都有人估計,估最高的是 Eric Toner 在去年十月就預測的新型冠狀病毒會大流行並且殺死六千五百萬人(註四),但是這個數據我一直找不到原始論文,加上這個數字看起來是在毫無防疫措施的前題下推估的,因此先暫且擱一邊。
另一個比較可信的,是 Lancaster 大學的 Jonathan Read 等人在 1/24 放上medRxiv 的一篇論文,在它的第一版(註五)中提到,在沒有有效的防疫措施下,到 2/04 那天之前光是武漢地區就會有約 20 萬人被感染,並且隨著春運,疫情會在中國各大城市遍地開花。不過在他們 1/27 上傳的第二版論文中(註六),已經刪去了該段 20 萬人染病的預測,並在修正了參數之後,更新了其中一項數字,推估 1/22 之前在武漢的總感染約是兩萬人。這篇因為是極少數的預測中有公開他們的原始論文內容的,因此後面的估計會有很大比重是參考這篇。不過須注意的是這篇並沒有討論何時會是疫情高峰。
所有的預測模型都有一個特性,就是數據點愈多,模型的預測力會愈強,因此愈是後期的預測,儘管數學工具不同,但是預測的數字會愈接近真實狀況。這兩天有兩位醫界專家不約而同分別作出預測,第一位是當年的抗 SARS 大將鍾南山,他認為疫情大概在 2/09 達頂峰,但是沒有預估最終的確診人數(註七)。第二位是香港大學醫學院院長梁卓偉,他推估在 1/27 日武漢已經有四萬多人感染,而且因為中國其他大城市也會出現大規模傳染,因此整體疫情要到四月底才封頂(註八)。
綜合上述預測內容,搭配手邊現有的確診數據,我嘗試以這個簡單的動力學模型來推估可能的確診數。我這邊給了三個劇本,分別以感染人數從多到少、每日確診人數高峰從遠到近來排序。(因為篇幅的關係,模型參數的設定、優化等細節,就不在此贅述)
第一個狀況,比較接近 J. Read 在第一版的論文中提到的,防疫效果低,感染人數在 2/04 破 20 萬的情況(注意可被觀察到的確診數要推遲到 2/18 以後),另外參考當年 SARS 在台灣高峰期大約 18 天,可以得到右圖的橘色線。在這個劇本中,確診數高峰要到 2/09 左右之後才開始,並且持續約三周,最高峰時每天新增的確診病例會達到一萬人,最後總人數約 24 萬上下(請自行抓正負 50%)。而這個劇本也會很類似 2003 年台灣的 SARS,當年和平醫院封院,距離高峰期還有兩個禮拜;而這次武漢封城 (1/23) 距離 2/09 也大約是兩周多一點的時間。不過這個狀況的前題是前面封城的防疫效果不佳,到二月初才停止人傳人的情況,而 J. Read 那篇第二版自己也拿掉了這個預測,因此是否走這個劇本,我個人比較抱持懷疑的態度。
狀況二(綠線),就是武漢等城市封城隔離後收到不錯的防疫效果,因此在 1/23 之後人傳人的情況大幅收斂,往後推遲兩周被確診,高峰會落在二月上旬,疫情高峰期縮短為兩週,並且從 1/31(也就是明天)開始進入高峰期。在高峰期每天會有 2000 ~ 4000 個新增確診病例,並持續到 2/13。
這個劇本很接近鍾南山的預測,但在套用模型後則多了一個確診人數總數,這邊推估約六萬人左右。因此狀況一與狀況二要區分,就是觀察從最高峰的 2/09 之後是否新增數有下降趨勢,如果下降幅度明顯,到了 2/13 以後低於每日兩千人,那可能就代表高峰期已經結束,大概可以判斷這波疫情有機會在二月下旬走完。現在單從人數要判斷是狀況一還二,尚言之過早,但是以這次武漢封到街上跟《我是傳奇》電影場景差不多的情況,我比較傾向是未來會走這個劇本。
狀況三(紅線),最佳劇本,現在就是高點!但這個劇本連中國官方都沒人這樣樂觀,絕大部分的官員都是說防疫措施在元宵節 (2/08) 後見效,因此機會應該不高。但不排除一個特殊的情況,就是若這次 2019-nCoV 是以輕症甚至無症為主,而大家基於「不想去醫院作檢測以降低感染風險」的心態,反而造成很多病例無法被官方確認。若真的這樣走,這個特殊情況是我目前能想得到的可能性。畢竟這次武漢封的這樣徹底,可以說是人類史上最大規模的一個社會學實驗,這樣大規模的隔離如何影響人類的行為,也只能事後諸葛了。另外特別提一下,這個劇本的總確診人數是兩萬人,差不多就是 J. Read 推估 1/22 前的感染數,至於結果如何,再過一周見分曉。(2/02 更新:這個版本從 1/30 ~ 2/02 的新增確診人數來看已經破功,這個劇本可以丟掉了!)
在這邊特別強調一下,圖上的預測都是在武漢肺炎疫情沒有在其他中國境內城市大規模群聚感染的情況下才成立,因為若從其他城市開始新的散播,會爆發出新的獨立疫情,並且和舊的曲線交疊,會讓數據判讀相對複雜。若要考慮這樣的情況,那就會有一個圖上沒畫出來的「狀況四」,也就是港大醫梁卓偉講的,其他大城市遍地開花,疫情就真的要落到三四月才能善了(註九)。(2/02 更新:國際醫學權威期刊《The Lancet》在 2/01 發表了一篇論文,根據他們的流行病學模型推估,疫情要到四月之後才會達到高峰)。
【結論與聲明】
在這篇文章中我試著以簡化後的反應動力學模型來描述傳染病的傳播模式,試著讓具有高中以上化學基礎的人以直觀的方式理解傳染性疾病的爆發與疫情特徵。此外,若是在武漢以外的城市沒有大規模發生疫情傳播的情況下,此次中國的 2019-nCoV 武漢肺炎疫情有機會在二月底就能收尾,視防疫成效高低,最終的確診人數可能會落在六萬到二十四萬人之間。
另外本篇主要是作為教育用途,嘗試以理科教科書中教過的觀點來幫助大家理解疫情的可能進展。本篇提出的模型非常簡陋,有非常多的假設,與實際狀況差異甚遠,不能當作真正的流行病學預測,也未經同儕審查 (peer review),未經同儕審查,未經同儕審查(很重要所以要說三次),亦不能作為正式的學術報告。
再次強調,本篇的目的是要讓有高中化學背景的人,以比較直觀的方式去理解疾病傳染是怎麼回事。實際上,流行病學中經典的 SIR model,就是在A和B之後,多了一個R(代表康復者)。若對 SIR model 有了解的人,會發現本篇的模型只有 S 和 I(對應到A和B),所以沒有辦法算 reproductive number (R0) 這項(相當於被感染者一直具有傳染力,這點很不符合真實世界的情況)。
https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SIR_model
作者:台大化學系助理教授 徐丞志
【數據更新 2/01 】
本文發表於 1/30,圖上呈現的數據只統計到 1/29。此後的中國境內每日新增確診人數補充如下:
1/30 新增 1982 人(劇本二估 1755 人)
1/31 新增 2102 人(劇本二估 2116 人)
2/01 新增 2590 人(劇本二估 2512 人)
2/02 新增 2829 人(劇本二估 2929 人)
2/03 新增 3235 人(劇本二估 3340 人)
【參考文獻】
[註一]武漢肺炎確診人數數據來源
WHO – Novel Coronavirus (2019-nCoV) situation reports
https://www.who.int/emergencies/diseases/novel-coronavirus-2019/situation-reports
中華人民共和國國家衛生健康委員會
http://www.nhc.gov.cn/xcs/yqfkdt/gzbd_index.shtml
鳳凰新聞數據即時更新
https://news.ifeng.com/c/special/7tPlDSzDgVk
[註二]
Chris T. Bauch, et al. (2005) “Dynamically Modeling SARS and Other Newly Emerging Respiratory Illnesses: Past, Present, and Future.” Epidemiology, 16, pp. 791-801.
[註三]
WHO – Cumulative Number of Reported Probable Cases of SARS
https://www.who.int/csr/sars/country/en/
台灣傳染病統計資料查詢系統
https://nidss.cdc.gov.tw/ch/Default.aspx
[註五]
Novel coronavirus 2019-nCoV: early estimation of epidemiological parameters and epidemic predictions (1/24 ver. 1.0)
https://www.medrxiv.org/content/10.1101/2020.01.23.20018549v1.full.pdf
专家:武汉病毒感染人数14天内将突破25万 (1/26)
https://www.voachinese.com/a/experts-predict-quick-spread-of-coronavirus-01252020/5260321.html
[註六]
Novel coronavirus 2019-nCoV: early estimation of epidemiological parameters and epidemic predictions (1/27 updated. ver. 2.0)
https://www.medrxiv.org/content/10.1101/2020.01.23.20018549v2
[註七]
大陸專家鍾南山預估 武漢肺炎約再十天達到高峰 (1/29)
https://udn.com/news/story/7333/4309518
[註八]
武漢肺炎:香港專家估計逾4萬人感染 5月「見頂」 (1/27)
https://www.bbc.com/zhongwen/trad/chinese-news-51267519
[註九]
何清漣專欄:武漢肺炎引發的全球恐慌才剛開始 (1/28)
https://www.upmedia.mg/news_info.php?SerialNo=80236
Joseph T. Wu, et al. (2020) “Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study.” The Lancet, in press. https://doi.org/10.1016/S0140-6736(20)30260-9