武漢肺炎的基礎再生數是甚麼概念？
聲明：我的專長是能源轉型路徑與系統分析，不是傳染病學；如果有公衛專家的意見和本文有衝突，請以他們的說法為準。
最近關於新型冠狀病毒（俗稱武漢肺炎）的新聞裡，有一份研究推估出它的基礎再生數竟然達3.6（後已下修至2.5），並據而推論到2/4光是武漢一地的肺炎感染人數將達到25萬人。為了明白這個數字的意義，我看了說明該份研究方法學的補充文件，原來這是「SEIR傳染病學模型」中一個重要參數。
參考這個網站，SEIR模型將可能會遭受疾病感染的人群分成四種：易感（S）、受曝（E）、感染（I）、恢復（R），各個體以類似催化反應機制的方式，逐步依序從一種族群演進成另外一種，其動力學方程式為如下：
可以發現，如果把四式相加，淨值為零。亦即，簡化版的SEIR模型中，S+E+I+R四種族群的總人數是固定的（有些模型會考慮出生、死亡和施打疫苗等動態因素的影響），式中以N表示。另外從微分方程中也可以看出，族群之間的演變是單向的，亦即個體只能由SEIR這樣的路徑演變，不能倒轉次序 （有些模型會考慮免疫能力喪失的情況，使一部份人口由R族群回復成S族群）。
底下是利用RK-4數值方法跑上述動力系統的結果，各參數數值如下：N=10000人、𝛽=1／小時、𝛼=0.02／小時、𝞬=0.03／小時。其中黑線代表S族群的演變、紅線代表E族群、藍線代表I族群、綠線代表R族群。初始條件為S族群9900人、E族群100人。
有趣的事情發生在第一式除上第四式的時候。此時我們會得到：
第五式告訴我們，一個個體走完SEIR這個過程的同時，會有多少易感族群被曝露；這種邊際替代率的概念亦可簡單想成「平均一個新的被感染個體，會造成多少易感族群受曝」。其中的R_0即為基礎再生數。
（注：現實中一般條件下S遠大於E、I、與R，因此式5可近似成dS/dR=-R_0）
因此，該論文初始推測武漢肺炎的基礎再生數為3.6，即代表每一位新的感染者，平均會造成3.6位易感族群受曝。值得一提的是，一般的流感基礎再生數大約為2，因此這份研究得出的推估確實令人擔憂。
假設用更新過的基礎再生數2.5來看，即代表需要有效隔絕60%的易感族群不受暴露（(2.5–1)/2.5），才能將有效基礎再生數降至1；而只有當有效基礎再生數不超過1，傳染病才得以抑制不再成長（因為平均一個新的被感染個體只會造成一個新的易感個體受曝）。
不過，目前仍在對武漢肺炎研究和防疫的早期，一些和該研究使用的假設不符合的事證，比如該病毒可能一開始就有人傳人跡象、該病毒在潛伏期即有可能有感染力等等，都會讓這份研究的推估結果和現實出現誤差。不管如何，眼下只能盡力做好自己能做的防疫功課（勤洗手、戴口罩、沒事別外出等等），然後祈禱台灣不要發生大規模感染了。
附錄：SEIR模型簡易模擬程式碼（R語言）
####SEIR model
#S_t = -beta * I * S
#E_t = beta * I * S — alpha * E
#I_t = alpha * E — gamma * I
#R_t = gamma * I##Parameters
Population = 10000
beta = 1 / Population
alpha = 0.02
gamma = 0.03
TIME = 500
dt = 1### R-K Solver
## ODE
ODE <- function(r){
dr = c()
dr[1] = -beta * r[3]* r[1]
dr[2] = beta * r[3]* r[1] — alpha * r[2]
dr[3] = alpha * r[2] — gamma * r[3]
dr[4] = gamma * r[3]return(dr)
}## RK algorithm
RK <- function(r, dt){
k_1 = ODE(r)
k_2 = ODE(r + k_1 * dt / 2)
k_3 = ODE(r + k_2 * dt / 2)
k_4 = ODE(r + k_3)r_new = r + (k_1 + 2 * k_2 + 2 * k_3 + k_4) / 6
return(r_new)
}r_0 = c(Population — 100, 100, 0, 0)
r = matrix(ncol = 4, nrow = TIME)r[1, ] = r_0for(t in 2:TIME){
r[t, ] = RK(r[t-1, ], dt)
}plot(c(1:TIME)/24, r[,1],type=”l”,xlab=”Time (Hours)”,ylab=”Population”)
lines(c(1:TIME)/24, r[,2], col=”red”)
lines(c(1:TIME)/24, r[,3], col=”blue”)
lines(c(1:TIME)/24, r[,4], col=”green”)
