keywords: SIR model, Random walk, Matlab simulation
(Article was published on author’s personal Facebook in February 3rd)
前幾天看到一位台大化學系教授發表了一篇關於武漢肺炎的傳染病模型(註1)，並利用分子動力學來解釋感染患者經由接觸後擴散疾病，引起許多有興趣的老師們討論以及大眾們的好奇心。這種生活化的現象，透過科學的建模及模擬分析也讓我覺得有趣。在除夕前(武漢肺炎爆發之際)，我就開始關注每日的感染人數，對於感染人數的增長除了擔心以外更好奇疫情隨者時間演化的情形。總之，我稍微閱讀了徐教授的原文，並參考目前廣泛使用的流行傳染病模型，最後我也提出幾個我的看法及修正的模型。最後也寫了個程式，在所有人類random walk的作用下，觀察疫情在空間和數量對時間的變化。(疾管署也做過類似研究 註2)
本人沒有專業醫學和公共衛生知識，因此所提出的模型非常粗淺，只能在某種程度去逼近真實現象，歡迎更多醫學相關的朋友、或是物理化學統計相關的朋友參與討論。
(以下分析並無參考研究結果或是相關文獻，為本人興趣分析，if there are any problem feel free to let me know)
Model 1
首先，先簡單回顧一下台大徐教授所提的”利用化學動力學來解釋傳染病模型”。所描述的是假設在人群當中，有分為受感染的傳染者A以及尚未受感染仍然健康的未感染者B，傳染的模式是經由A和B發生接觸(或者理解成化學的發生分子碰撞)，並以一定的機率(此機率在化學當中就是平衡常數，受到Arehnius eqn中溫度影響)將疾病傳染給B。
簡單來說就可以理解成 A+B->2A，其化學速率定律式即為
傳染的機率正比於傳染者和未感染者接觸的次數([A][B])，就是常見的二級反應，而這個類似的微分方程式在掠食-獵物模型當中也曾見過(Lotka–Volterra equations/predator–prey equations)，描述的掠食者和獵物相遇的機率。
假設現在我們有感染者對時間變化的方程式後，理論上我們就可以計算出感染者隨時間演化的情形。但在分析前，我們在仔細思考並強調一下這個模型的假設：
- 感染的人是均勻分布在健康人群中
- 所有感染的病患不會被隔離亦不會接受到治療，而是持續地在傳染疾病
這個微分方程式應該是原文(註1)所使用的方程式(我不太確定，仍須確認)，但我認為在實際上模擬可能有問題。(內文有提到當感染人數~未感染人數，則感染速率最快)
因為從微分方程式當中(eq1)，可以發現這個疫情的擴散最後會收斂在[A]=1，也就是說代表所有人類最終都會感染到疾病，正常來說這應該不是最終的結果，應該只有部分的人類(3~10%?? or less)感染。
可以觀察到感染人數最終收斂在[A]=1，且每日感染人數呈現中型曲線，在中間時刻感染速率最快
但是實際上感染人數佔總體比例低(~5%)，所以我們先修正一下原本的Model 1
Model 2
由於感染人數[A]極少，而健康人數[B]=(1-[A])~1，因此 [B]→1， 因此速率定律式可以近似成
則
從解可以知道感染人數是隨著時間呈現指數性的成長，但沒有收斂的情形。因此此解僅能解釋爆發期間短期的改感染人數變化，無法解釋整個階段長期的變化。
因此從第一個建模的結果，和對現實生活的比較可以發現這第一個模型有許多問題
- 感染的人並不是均勻的分佈在健康人群中，感染源通常是從第一個感染者開始傳給附近的人，在逐漸的輻射擴散散布出去，進而感染到更多人(換言之，就是說感染者並不是像化學溶液有攪拌石，每個時刻都是均勻分布的，而是漸漸地擴散出去)也就是說第一個化學動力學的模型，在空間上就和現實有點落差。
- 感染的病患事實上會被篩檢出來，而且被驗出會馬上遭受隔離，而失去持續傳染的能力。但是在第一個動力學模型，我們卻假設所有的病患都持續在人群當中傳染疾病(除非在衛生署毫無作為的情形!?)與實際現況有明顯落差。因此事實上，每個時刻擁有傳染能力的只有在潛伏期當中(已經感染，但尚未被檢出隔離者)的患者才會在人群當中傳染
- 從疾病爆發到開始蔓延，在不同的疫情階段感染率會有改變，因為透過政策的宣導(提醒要戴口罩、洗手等以及疫苗的開發)，因此在潛伏期有傳染能力的病患在初期十分容易感染健康的民眾，但隨者疫情來到末期，即便感染者和未感染者接觸其被傳染的機率也會逐漸降低。但在第一個模型當中，我們也沒有考慮到這方面的問題，使得疫情爆發到蔓延皆是以反應速率k(固定機率)傳染。
首先我們先來改善第三個問題
Model 3
假設自從傳染病被診斷出來後，因為衛生署和政府機關的宣導之下，民眾提高防疫意識，使得感染率會下降，我們可以假設傳染機率從原本的r(t)=r_0變成隨時間遞減的r(t)=r_0-beta*t，隨者時間增加被傳染率遞減，因此我們可以修正我們的分子動力學傳染模型為
同樣的利用[B]~1條件，且僅考慮r_0-beta*>0的狀況
整理以後可以發現[A(t)]感染累積人數，其實就是以 t=r_0/β，半高寬(標準差)為1/√β的高斯曲線呢！而我們僅考慮r_0-beta*t>0，因此感染累積人數就是半邊的高斯曲線，如圖(二)中紫色曲線
若是把感染累積人數對時間微分，即可以得到每日新增感染人數r_A
可以發先他是有點歪掉的高斯曲線(多了r_0-beta*t 這一項)，可以發現疫情的擴散一開始平緩，接者會逐漸陡峭，然後上升逐漸遲緩並漸漸來到高峰後，每日新增感染人數便會快速的降低，即疫情結束。
小結：可以發現多了一點的修正，我們的模型更接近真實現象，也更加符合常理了
接者我們來改善第一個問題
Model 4
在Model 1當中，化學分子動力論是假設所有傳染的人群在空間上均勻分布，但是實際上疫情是從第一個傳染者開始往外擴散。因此我們這邊必須修正在每個瞬間傳染者所接觸到的健康人群是個時間的函數，亦即在疫情剛爆發時，感染者比較集中，且會隨著時間慢慢擴散，所接觸的人群會逐漸增加。從擴散方程式(diffusion equation)，或是隨機漫步(random walk)，我們可以知道疫情的擴散半徑是正比於時間的根號(Infectious Radius ∝ t^1/2)0。因此在二維平面當中，擴散半徑所涵蓋的面積即為 A=πR²∝ t。所以我們感染者所接觸到健康的人可以修正成
因此我們可以進一步的再修正我們的分子動力學傳染模型為
而這個微分方程式可以經由簡單的變數分離法積分得到
最後我們修正第二個問題(潛伏期)
Model 5
在原本的分子動力學模型當中，我們是考慮所有感染的病患都在持續在人群當中持續傳染而未受到隔離。但在現實生活中，只有在潛伏期當中的人才不會遭受隔離持續感染
因此我們要修正我們的方程式，並假設潛伏期為 d。
由於這個方程式我不知道怎麼解出解析解(若是有人知道辦法麻煩指教一下!!)，所以我只能用matlab去得到數值解(用歐拉法跑回圈)
只有在(t, t-d)之內遭受感染的人才會持續在人群中感染，其餘在t-d之前感染的會遭到隔離並接受治療(痊癒或死亡?)
從圖(三)左邊可以發現每日新增感染數會先有一個平緩的潛伏，然後進入爆發期直到疫情最嚴重的巔峰期，然後就快速降低受到控制。而右圖就是感染累積人數隨者時間的變化。
以上是透過微分方程式並做了一些簡單的修正來描述疫情擴散的情形。
Epidemic expand simulation
最後我也利用程式去模擬整個疫情擴散的過程，並用動畫去呈現，就可以看到疫情在空間和數量對時間的變化情形。
假設在一個固定的空間之中(無限邊境)，所有人都是處於隨機走動(Random walk)的情形之下，令第一個感染者在原點遭受感染，並開始傳染給附近的人(Infectious radius)在感染前到潛伏期(Incubation time)結束。若是健康的人距離感染者在傳染半徑之內，健康的人會有一定的機率(Transmission Rate) 遭到感染。而防疫部門會開發疫苗和投藥，使得未感染者受感染率的機率會隨著時間線性下降(Prevention rate)。
最後模擬的結果在n=20000人在350x350單位的空間當z中，傳染半徑11單位，隨機走動速率3單位，傳染率0.3%，之下經過2000單位時間。模擬結果可以發現第一個感染者(x,y)=(0,0)出現以後，疫情先進入一段潛伏期，接者開始快速擴張，直到病情最高峰後，進入控制期後就快速消滅了。經過多次模擬統計結果給出，在這樣的結果下約有10%得人口遭受感染，且感染區域呈現輻射的樣子。此外累積人數和每日新增感染人數的圖形也和前面利用微分方程式所算出的結果相近，所以可以確認前面所建立的微分方程式模型在某程度上是可信的。
最後，再次強調此模擬僅考慮隨機漫步 (Random walk)，因此考慮人群無使用交通工具的情形(更好的模擬可以使用Lévy flight)。此外所使用的傳染半徑、傳染率、空間大小、潛伏期皆沒有參考任何醫學數據(也找不太到)，總之就是一個沒有參考任何數據所建立的定性定量模型（因此在量的部分無法估計），僅希望能對疫情擴張的現象做出一些科學描述。
願疫情早日平復，世界重回健康和平。
後記
P.s 本文無公布對未來的任何預測
第一、過多的預測只是導致恐慌
第二、隨意預估疫情走向會被疾管局開罰300萬
第三、本文只做作為學術、教學、交流討論用途
如有任何興趣或是問題，歡迎討論
P.P.s 本文於 2/3 發表於個人 Facebook， 內容包含理論推導以及動畫模擬，內容並無過多更新
待疫情結束後，希望我還會回來回顧這篇文章。
Reference
1. 武漢疫情預測 Cheng-Chih Richard Hsu ,
https://www.facebook.com/photo.php?fbid=10217029362669127
一篇被廣為盛傳的文章，國內外媒體和記者們爭相報導2.利用整合式傳染病監測軟體以協助流行性感冒的監測,行政院衛生署疾病管制局九十五年度科技研究發展計畫
https://www.cdc.gov.tw/uploads/files/24309cc6-96a4-42b7-84ab-9f10794081f7.pdf?fbclid=IwAR35tDNacEkiyl1zuT2WTZBeyeoY2MlrRmw_e8KIzaV58-h8lB4Th1xiaC0
一篇疾管署和交大資工合作的研究計畫3.Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study
https://www.thelancet.com/action/showPdf?pii=S0140-6736%2820%2930260-9&fbclid=IwAR38f0MK9ORdIeeKCz04Up4I4jkMHsZulfu6uUCpzmA1uE8Ta7afrXdAxa8
一篇由香港科大醫學系探討在交通下疫情的擴散情形研究，發表於 Lancet