數學模型論新冠清零及 「與疫共存」
之兩難困境
陳博彥
先前 《就新冠疫情數學建模談起》 一文 [1] 提出台灣在嚴格邊境管制下, 所產生之零確診 (zero confirmed) 結果進行初步相關評估分析。 新近由於 COVID-19 具高傳染力之 Omicron 等病毒株高速傳播擴散, 時中部長更言及 「社區傳染」 或 「與疫共存」 可能已是在所難免。 因此 學生令人訝異地問筆者: 「當零確診已無可能時, 零確診的模型是否不再適用, 而須適時修正?」
再者, 誠如季刊編者對此文之註解 [2], 指出所有疫情速率常數並非定值, 其動態時間趨勢無法 知曉。 此點確實是不爭之事實, 但是最終平衡穩態行為之預測分析, 可做為疫情智庫決策及全 民防疫作為判斷上極重要之關鍵參考, 此分析提出各種應對防疫策略下, 判斷在何種引導原則 下疫情才可走向安全無虞之疫情控制條件之下。 又者, 「零確診」 若已然不再是可能選項, 又要 如何 「與病毒共存」? 以下就現今所見之疫情演進來 「滾動式修正」 補敘前文未竟之相關內容 [1], 並且澄清前述幾個問題以了解目前究竟應採 「堅壁清野」(清零) 或是該要面對 「與病毒共 存」(不再清零) 來看待此類傳播力強化, 但可能重症或致死機會較低的病毒變異株, 以進行深入 評估。 首先針對新冠疫情提出之流行病學 (Epidemiology) [3], 假設若未知感染源之破口隱藏 存在條件(零確診不再成立: Y (t) > 0) 下, 來修正 SIR 型流行病學模型 (SIR-type epidemic model) 如下。
(A) 反應機制 (I)—針對一般疫情傳播狀況而言:
X−→βY Y −→ν′ Z, Y −→Ø X,
其中 ν = ν′+ Ø > Ø, 事實上, 目前台灣因 COVID-19 有效控制疫情, 死亡人數相對於總人 口數甚微情況下, 可設定系統內總人口數為定值: X + Y + Z = N0, 如此三聯立正微分方程 可由此關係式降階為 2 × 2 之正微分方程組, 因此可建立之疫情流行系統方程式如下:
dX
dt = − βXY + ØY + m = f1(X, Y ), X(t = 0) = Xi, dY
dt = βXY − νY + m′ = g1(X, Y ), Y(t = 0) = Yi >0,
其中 X, Y, Z 分別代表未受感染者 (Susceptibles)、 已被感染者 (Infectious)、 自免疫系統去 除者 (Recovered) (如終生免疫者跟死亡者)。 而 β, Ø, ν, ν′ 分別代表疾病傳播常數, 疫病復 原回復後可再被感染 (Reinfection) 速率常數, 感染者總去除速率常數, 感染後得以免疫速率 常數。 m, m′ 分別表因防疫破口, 無法有效阻斷, 所引入之未感染者速率及已被感染者速率。 切 記: 所有瘟疫之結局皆是所有人口遠離流行疫情系統, 才能完全解封來過正常生活。 唯一差別 在於疫情控制成功與否, 是 Z 將以不幸亡者身分離世、 抑或是活著恢復健康著離開此流行病學 系統。 換言之, 防疫成功的目標指數是尋求在最少死亡人數且最多免疫者條件下, 脫離此疫情陰 霾。 當 dXdt = dYdt = 0 條件下, 其平衡穩態滿足
dX
dt = − βXY + ØY + m = 0, (1) dY
dt = βXY − νY + m′ = 0, (2) 由 (1) + (2) 可得到 Y0 = m+ m′
ν −Ø , 代入 (2) 可得 X0 = ν
β −m′ ν −Ø
m+ m′ = νm+ Øm′ β(m + m′). 亦即此時平衡穩態 SS (steady state) 是
(X0, Y0) =νm+ Øm′
β(m + m′),m+ m′ ν −Ø
.
依據此系統條件來決定穩態 (SS) 特性及分析其穩定性 (Stability) [4]。 若假設 ∆X = X−X0, 以及 ∆Y = Y −Y0, 則可推得如下線性化 ODEs:
d
dt∆X = − βY0∆X + (Ø − βX0)∆Y, (3) d
dt∆Y = βY0∆X + (βX0−ν)∆Y. (4) 此系統方程組之特徵值(eigenvalues) 可由如下行列式來決定
det(A − λI) =
−βY0−λ Ø − βX0
βY0 βX0−ν − λ
= 0, 亦即
λ2+ (ν + βY0−βX0)λ + βY0(ν − Ø) = 0, (5) 由此可決定出其特徵值 λ1, λ2 必須滿足上式, 並且以下式表之:
λ2−(λ1+ λ2)λ + λ1λ2 = (λ − λ1)(λ − λ2) = 0.
若系統可穩定收斂至穩態 (X0, Y0) 下, 則兩特徵值皆必須滿足具有負實部之先決條件下, 亦即 必須滿足
(i) β(Y0−X0) + ν > 0, (ii) ν − Ø > 0.
事實上, 由於 ν > Ø 恆成立, 因此即使是不再清零的疫情系統中, 若仍要維持於長期疫情控制 收斂穩定之下, 以達到疫情在完全可控制狀態下, 顯然其必要條件是滿足 X0 < Y0+ νβ (此不 等式之確定, 亦即表示病毒溯源傳播疫調有其絕對之必要性。 任何時間未感染人數、 已感染人數 與平衡 X0, Y0 之動態值以及各病毒變異株相關之動力學參數 ν, β 皆須確實全程掌握, 在此定 義: Y0+ νβ 為 「關鍵」 族群數目)。 若無法有效持續追蹤, 清零之可能目標即已無可能。 單由此 點可知, 台灣在一開始即已在吹哨人提醒下, 當機立斷超前部署進行出入境管制戒備 (才有機會 能設定清零 Y0 = 0 為目標), 而能在世界防疫戰中取得可掌可控之領先地位。 反觀世界一些國 家早在疫情爆發初期, 即錯失良機導致國門失守, 甚至因未知傳播鏈過多, 即已然放棄疫調, 從 此就被迫注定走入可能必須 「與病毒共存」 的道路上。 防疫成功與否, 在於是否能控制未被感染 人口數目始終小於 「關鍵」 族群數目: Y0+νβ, 以達到疫情可調可控之狀況下。 由上述訊息更可 指出: 現今若能與有效疫苗覆蓋率合併計算, 即可確實推論出是否疫情系統已在可掌可控的發 展情況下, 漸次步出瘟疫陰霾?
再者, 此系統平衡穩態 (X0, Y0) 先天存在著三個令人感到頗有意義的極限疫情情況, 在此 分別點出來進行分析如下:
(i) 第一個極限穩態: lim
m→0(X0, Y0) = (Øβ,mν′′), 代表即使採取更嚴格之邊境隔離檢疫管制, 若 在人口中存在隱形帶原無症狀感染者存在, 則會發現系統中永遠無法達到零確診之情況, 此亦即是不少國家日日出現上萬人確診的可能原因所在。
(ii) 第二個極限穩態: lim
m′→0(X0, Y0) = (νβ,mν′), 說明在實施完全有效入境邊境篩檢, 來剔除受 到感染者, 卻不完全封城, 仍可自由出入境, 可能亦無以避免在人口群眾中始終有確診者 存在。 台灣目前對此種可能發生之機會仍要格外小心注意。
(iii) 第三個極限穩態: lim
m→0,m′→0(X0, Y0) = (βν,0) 或 (Øβ,0), 代表只有在完全鎖國封城之特 殊條件下, 才能達到所謂不沾鍋的 「完全零確診」 結果。
上述三種極限穩態皆說明若無完全封關、 封城或是鎖國條件下, 「社區感染」 或是 「與病毒 共存」 的發生機率會是始終存在的。 時中部長所言實有所本。 一旦出現有破口未能有效堵絕, 未 知傳播鏈隱形存在下, 如此即造成無法有效隔絕所有可能傳播感染機會, 亦代表疫情走上 「與病 毒共存」 之機會即難以避免。 無可諱言地, 自始至終皆未曾完全封城或封關的台灣疫情而言, 所 言之 「與病毒共存」 的危機事實上是處處存在、 而且亦是在所難免。 但即使如此, 若能時時掌握 疫調, 全民配合防疫生活確實執行, 疫情可掌可控仍是一可達到的目標。 此三種極限穩態更清楚 詮釋前文 [1] 所提現今全球可能極為關注的三個穩態: SS1(0, 0) (SS1 即全民脫離疫情系統,
自由出入往來。 後續反應機制 (II) 將補充說明此穩態之存在性。), SS2(Øβ,0), SS3(νβ,0) 存在。
新近由於各種疫苗施打逐漸覆蓋, 究竟會對流行疫情系統內族群人口產生何種效應? 以下 即是就 「有效疫苗覆蓋」 來探討疫情系統穩態平衡狀況, 更對何時才能 「完全解封」, 「自由出入 國際旅遊」 做出之沙盤推演說明如下:
(B) 反應機制(II)—針對有效疫苗被發展出來進行施打注射覆蓋 (亦即因有效疫苗注射施打而 免疫: X → Z) 之疫情傳播狀況 (施行有效疫苗注射施打, 即可離開此疫情系統) 而言:
X−→βY Y −→ν′ Z, Y −→Ø X, X−→Ψ Z,
其中 ν = ν′+ Ø > Ø, 系統中所有人口數可為 X + Y + Z = N0, 因此在可掌可控之疫情系 統所建立之流行病學方程式下, 當 dXdt = dYdt = 0 之平衡穩態下, 須滿足:
dX
dt = − βXY + ØY − ΨX = 0, X(t = 0) = Xi, dY
dt = βXY − νY = 0, Y(t = 0) = Yi >0, dZ
dt = ν′Y + ΨX.
由上述方程式可解出兩個穩態 (X0, Y0) 分別為 SS1(0, 0), SS2(νβ,(Ø−ν)βΨν )。 但第二個穩態並 不存在 (Y0 <0)。 因此 SS1(0, 0) 為唯一存在之穩態。 代入線性方程組可得
d
dt∆X = (−βY0−Ψ)∆X + (Ø − βX0)∆Y, (3′) d
dt∆Y = βY0∆X + (βX0−ν)∆Y. (4′) 取定唯一存在之穩態 SS1(0, 0) = (X0, Y0) 代入可得
d
dt∆X = − Ψ∆X + Ø∆Y, (3′-1) d
dt∆Y = − ν∆Y. (4′-1) det(A − λI) = 0可得兩個特徵值分別皆為負數的 −Ψ 與 −ν。 解析解為
(∆X, ∆Y ) =
∆Xie−Ψt+ Ø∆Yi
Ψ − ν(e−νt−e−Ψt), ∆Yie−νt , 當達到平衡穩態時 lim
t→∞∆X = lim
t→∞∆Y = 0。 換言之, 在有效疫苗注射覆蓋條件下, 平衡後任 何時間皆維持在遠離此瘟疫的狀況, 至此不再有人存在此疫情系統中 (亦即完全解封正常生活)。
亦即是 X(t) = X0, Y (t) = Y0 = 0。 此結果即是前述的完全脫離瘟疫的平衡穩態 SS1(0, 0)。
然而反觀, 另一完全與此反應機制背道而馳的例子則是在疫情爆發初期, 由於並無有效防疫方
法及策略, 因此除了個人衛教自我防護及加強醫療備援外, 對政府下達之公共衛生之政策又應 當如何因應? 以下即是所討論的狀況分析:
(C) 反應機制 (III)—針對有效疫苗未被發展出來, 因此無法產生有效施打預防疫苗注射覆蓋, 而且即使是受感染後康復, 仍會再次被感染狀況而言:
X−→βY Y −→ν′ Z, Y −→Ø X, Z−→Ψ′ X,
其中 ν = ν′+ Ø > Ø, 系統內總人口數為 X + Y + Z = N0。 由此可推得其疫情系統微分方 程組如下:
dX
dt = − βXY + ØY + Ψ′X, X(t = 0) = Xi, dY
dt = βXY − νY, Y(t = 0) = 0, dZ
dt = ν′Y −Ψ′Z, Z(t = 0) = 0.
將總人口數 N0 代入, 原式可推導為 dX
dt = − βXY + ØY + Ψ′(N0−X − Y) = 0, X(t = 0) = Xi, dY
dt = βXY − νY = 0, Y(t = 0) = 0.
由上述兩個方程式可推導解出兩個穩態 (X0, Y0) 分別為 SS1(N0,0), SS2ν
β,Ψ′(N0− ν
β) Ψ′−Ø + ν
=ν
β,Ψ′(N0− ν
β) Ψ′+ ν′
.
對於 SS1(N0,0), 即是台灣一直希望能將疫情控管在零確診 (清零) 的平衡穩態狀況下。 然而 SS2ν
β,Ψ′(N0− ν
β) Ψ′+ ν′
即是現今世界不少國家所提到之 「與病毒共存」(例如: 英國所提之 「佛系防疫」) 的平衡穩態。
而且在此平衡穩態下, 在整個人口族群中, Y0 : Z0 亦會呈現幾乎固定比例 (亦即 Ψν′′)。 換言之, 當疫情開始爆發, 即已經無法有效掌握 「零號病人」 去向, 無法有效疫調追蹤隱形傳播鏈時, 此 時就算後續疫情控制得宜, 可能 「與病毒共存」 亦是在所難免 (例如: 目前新聞報導提及不少國 家, 日日皆有上萬人確診者, 即是實際發生之案例)。 以下即是證明滿足此兩個平衡穩態
SS1(N0,0), SS2ν
β,Ψ′(N0− ν
β) Ψ′+ ν′
之必要條件之推導過程。 首先再次代入線性方程組可得 d
dt∆X = (−βY0−Ψ′)∆X + (Ø − Ψ′−βX0)∆Y, (3′′) d
dt∆Y = βY0∆X + (βX0−ν)∆Y, (4′′) 取定存在之穩態 SS1(N0,0) = (X0, Y0) 代入可得
d
dt∆X = − Ψ′∆X + (Ø − Ψ′ −βN0)∆Y, (3∗-1) d
dt∆Y =(βN0−ν)∆Y. (4∗-1) 上述兩 ODEs 可解得解析解為
∆X =∆Xie−Ψ′t+ (Ø − Ψ′−βN0)∆Yi
e−(ν−βN0)t−e−Ψ′t ,
∆Y =∆Yie(βN0−ν)t.
若 N0 < βν, 當達到平衡穩態時, 可得 lim
t→∞∆X = lim
t→∞∆Y = 0 疫情可有清零 (零確診) 之 穩定平衡點。 反之, 若 N0 > βν, 則不會在平衡點 SS1(N0,0) 穩定存在 (亦即疫情失控而遠離 SS1)。 針對平衡穩態
SS2ν
β,Ψ′(N0− ν
β) Ψ′+ ν′
取定線性方程組 (3∗), (4∗) 之 det(A − λI) = 0 可得兩個特徵值 λ1, λ2 必須滿足如下方程 式:
λ2+βΨ′(N0− ν
β)
Ψ′+ ν′ + Ψ′
λ+ βΨ′
N0− v β
= 0, 亦即
λ2−(λ1+ λ2)λ + λ1λ2 = (λ − λ1)(λ − λ2) = 0.
若 N0 > νβ 條件下, 可得到兩個具有負實部之特徵值, 因此疫情可達到趨於穩態之平衡點
SS2ν
β,Ψ′(N0− ν
β) Ψ′+ ν′
,
(亦即 「與病毒共存」 之平衡狀態)。 由此可以說明何以在未知瘟疫爆發時, 第一時間就必須要做 到最有效之隔離檢疫, 避免個體間產生不必要之接觸而形成傳播, 才能致使區隔出來之每個小 單位區域中 N0 < νβ, 如此疫情才有逐漸式微, 不易形成流行的機會。 事實上, νβ 於流行病學上 更稱之為 「閥值人口」(threshold population) 數目。 新近新聞指出丹麥成為歐盟首個解除疫 情限制措施國家, 不少國家亦會跟進。 此作法亦即主動採取 「與病毒共存」 之策略。 取定此策略
就代表爾後之現有醫療量能, 已經有足夠把握來應付 Y0 : Z0 幾乎呈現固定比例 Ψ′
ν′ 之確診比
例情況。
(D) 反應機制 (I∗)—針對前文 [1] 所提一般疫情傳播狀況來看:
X−→βY Y −→ν′ Z, Y −→Ø X,
其中 ν = ν′+ Ø > Ø。 若假設總人口數定值: X + Y + Z = N0, 可建立出流行病學疫情系 統方程式如下:
dX
dt = − βXY + ØY = f (X, Y ), X(t = 0) = Xi, dY
dt = βXY − νY = g(X, Y ), Y(t = 0) = Yi >0.
顯然在任何 X0 值條件下, 只要 Y0 = 0 之解, 皆可以滿足上述兩式為平衡穩態之結果。 其線性 化 ODEs 亦是式 (3) 及 (4), 若式 (6) 中設定 Y0 = 0 為解, 則可得
λ2+ (ν − βX0)λ = 0, (6) 其中兩個特徵值分別為 0 與 βX0−ν, 其相對應之特徵矩陣則分別為
"
1 0
# 與
"
Ø − βX0
−βX0 + ν
#
。
若起始條件
"
∆X(t = 0)
∆Y (t = 0)
#
=
"
∆Xi
∆Yi
#
=
"
Xi−X0 Yi−Y0
#
, 其矩陣解則可表為
"
∆X
∆Y
#
= C1
"
1 0
#
+ C2e(βX0−ν)t
"
Ø − βX0
−βX0+ ν
# , 其中 C1 = Xi−X0− Yi−Y0
ν−βX0(Ø − βX0) = C1(Xi, Yi), 而且C2 = ν−βXYi−Y00。 此式子亦指出 在疫情爆發當下, 第一時間若能超前部署, 確實有效執行隔離分區檢疫, 致使起始條件即能滿足 Ni = Xi ≃X0 < νβ, 其中 Ni 及 Xi 分別代表起始系統內所有人口數及未被感染人口數(亦即 起始時間 「零號確診病人」 尚未出現, 並無感染現象發生)。 如此疫情才能逐漸壓制下來, 漸次走 向 「清零」。除此之外, 更可得到關係式如下:
(i) lim
t→∞∆X = Xi−X0− Yi−Y0
ν − βX0(Ø − βX0) 或是
t→∞lim X = X0 + C1(Xi, Yi) = Xi− Yi−Y0
ν − βX0(Ø − βX0), (ii) lim
t→∞∆Y = 0 或是 lim
t→∞Y = Y0.
以上結果至少強調出三個疫情發展重點: (a) 瘟疫發生時, 追蹤感染傳播鏈之疫調工作極為重 要, 因為任何時間之 X(t), Y (t) 以及各變異株傳染相關參數 β, Ø, ν 皆主宰著疫情是否可有
效控制達到可掌可控的情況下。 (b) 最為重要的是, 若在防疫無缺口來掌握疫情條件下, 更要確 認在任何時間, X(t) < νβ 此必要條件務必皆能滿足。 (c) 第一時間起始疫情條件 (Xi, Yi) (Yi
表起始時間之確診人口數目) 主宰著當確診數目歸零時, 仍存留於感染系統中, 可能具有被感染 風險之人口數目 (「吹哨人」 效應) (即 lim
t→∞X = Xi− Yi−Y0
ν−βX0(Ø − βX0))。
致謝
感謝教育部教學實踐研究計畫 (PEE1100793) 及科技部計畫 (MOST 109-2221-E-197- 016-MY3) 之經費補助, 更感謝藉由此計畫執行過程中能得力於國立宜蘭大學教學發展中心的 夥伴們持續協助及鼓勵, 以及校方教師社群計畫支助協助能聯結跨校與人文社會科學專業教授 們 (玄奘大學社會工作學系王敏菱、 中國文化大學心理輔導學系陳柏霖、 玄奘大學藝術設計學 院王振邦、 國立台東大學生命科學系黃祥恩、 國立東華大學教育與潛能開發學系張德勝等位教 授) 有腦力激盪實質交流, 方有如此人文關懷研究的進行。 最後亦致上於心底最誠摯的感謝: 能 在 1991年間於台大公共衛生學院林瑞雄院長手下擔任研究工作, 才能有機會對流行病學有初步 的入門體會。 更影響到後續論文能以創新的病毒疫情感染流行病學模型來理解微生物生產製程 [5], 因為所有科學知識學習皆收斂於一事: 真理。 在此皆一併致上最誠摯的感恩之意。 並將此文 獻給在天上的筆者指導教授台大化工系黃世佑教授以及加州大學爾灣分校的 Professor Henry C. Lim。
參考文獻
1. 陳博彥, 林育秀, 王敏菱。 就新冠疫情數學建模談起。 數學傳播季刊, 45(4), 71-76, 2021。
2. 梁惠禎。 編者的話。 數學傳播季刊, 45(4), 1, 2021。
3. W. O. Kermack and A. G. McKendrick, Contributions to Mathematical Theory of Epidemics-I, Bulletin of Mathematical Biology, 53(1/2), 33-55, 1991.
4. B. Y. Chen, Stability analysis on SARS epidemics, Math. Spectrum: Letter to the
Editor, 40(1), 38-39, 2007/2008.
5. B. Y. Chen, A Study on Temperature Induction of Bacteriophage λ in Escherichia coli.
PhD Dissertation. University of California, Irvine, December 1995.
—本文作者陳博彥任教於國立宜蘭大學化學工程與材料工程學系—
