如果新冠病毒是在亞美尼亞爆發會發生什麼？程式設計師用Python進行模擬

• 2020 年 2 月 24 日
• 筆記

大數據文摘出品

來源：towardsdatascience

編譯：大萌、趙吉克、武帥、錢天培

2019年末，在中國武漢爆發的冠狀病毒疫情衝擊了整個金融市場和實體經濟。這座總人口超過千萬，春運期間流動人口超過500萬的巨型城市的災難在世界各地引發了一連串蝴蝶效應，也在全球普通民眾中引發恐慌。

2020年1月30日，2019-nCoV已被世界衛生組織列為國際關注的突發公共衛生事件（PHEIC）。在撰寫本文時，尚未發現經過醫學研究驗證的具體治療方法。

此外，一些關鍵的流行病學指標，如基本再生數(一個感染者傳染的平均人數，即R0值)仍然未知。當今時代，全球互聯，這類流行病更是借勢成為全球範圍內的重大威脅之一。

可以推測，如果2020年全球發生災難性事件(大致定義為導致不少於1億人傷亡的事件)，最有可能的原因恰恰是某種流行病——而不是核災難，也不是氣候災難，等等。全球範圍內的快速城市化進一步加劇了這一問題，人口密集且頻繁流動的城市儼然變成了疾病擴散網路的傳播節點，使得防疫系統變得極為脆弱。

武漢的這場災難也引發了全球對於城市規劃和防疫政策的思考。如果疾病不是在武漢，而是在另一座人口規模更小、人口流動也更弱的城市，它的傳染性和感染人數，又會是怎樣的故事？

在這篇文章中，我們將討論當流行病襲擊一個城市時會發生什麼，應該立即採取什麼措施，以及這對城市規劃、政策制定和管理有什麼影響。

我們將以亞美尼亞首都，一個人口剛過百萬、以鋼鐵和葡萄酒著名的城市埃里溫市為例進行研究，建立數學模型並模擬冠狀病毒在該市的傳播，研究城市流動模式如何影響疾病的傳播。

城市流動性

有效、高效和可持續的城市流動性對現代城市的運作至關重要。它已經被證明會直接影響城市的宜居性和經濟產出(GDP)。然而，一旦發生疫情，它就會火上澆油，擴大疾病的傳播。

那麼，讓我們先來看看埃里溫市在一個平面坐標繫上的聚合OD流動網路(Origin-Destination)，以了解城市流動模式的空間結構：

接著，如果我們觀察網格的總流入量，我們會看到或多或少的單中心空間組織，其中一些網格的日流入量較高但位於中心之外:

現在，假設一種流行病在城市的任意地點爆發。它將如何傳播？我們能做些什麼來控制它？

流行病學建模

為了回答這些問題，我們將建立一個簡單的房室模型來模擬傳染病在城市中的傳播。

當一種流行病爆發時，其傳播動力會有顯著變化，這取決於最初感染的地理位置及其與城市其他地區的連接性。這是最近關於城市人口流行病的數據驅動研究得出的最重要見解之一。然而，正如我們將在下文進一步看到的，各種結果都要求採取類似的措施來控制疫情，並在規劃和管理城市時考慮到這種可能性。

註：compartmental model，房室模型，也稱SIR模型，是一種簡化的傳染病數學模型。

由於我們的目標是展示城市流行病傳播的一般原理，而不是建立一個實時的精確模型，我將參考《自然》雜誌的一篇文章，通過簡單地修改經典SIR模型即可滿足我們的需求。

相關鏈接：

https://www.nature.com/articles/s41467-017-02064-4

這個模型把人口分成三部分。在時間t的每個位置i，其三個房室如下：

• Si,t：尚未感染或易感的人數；
• Ii,t：感染疾病並有能力將疾病傳播給易感人群的人數；
• Ri,t: 由於康復或者死亡，在被感染後從受感染組中移除的人數。這個群體中的個體沒有能力再次感染該疾病或將感染病傳染給他人。

在我們的模擬中，時間將是一個離散變數，因為系統的狀態是以日為單位進行建模的。在t時刻j點的完全易感人群中，感染病出現的概率為：

βt是t時刻的傳輸率; mj,k是從k地到j地的流動性, xk,t 和yk,t 代表t時刻k地和j地的感染和易感人群數，其中xk,t = Ik,t / Nk，yj,t = Sj,t / Nj，Nk和Nj代表k地和j地的人口總數。然後，我們繼續模擬一個隨機過程，將這種疾病引入完全易感人群所在的地區,其中Ij,t+1 是概率為h(t,j)的伯努利隨機變數。

一旦感染在隨機地點出現，疾病不僅會在該地點傳播，還會通過攜帶者傳播到他處。這就是以OD流矩陣為特徵的城市流動模式發揮關鍵作用的地方。

此外，為了確定疾病是如何通過感染者傳播的，我們需要考慮其R0值。此處，其中y表示的是治癒率，也可以認為是二次感染率。在撰寫本文時，新型冠狀病毒的基本再生數估計值在1.4到4之間。凡事做最壞的準備，因此我們假設R0值為4。需要注意的是，R0值是一個有期望值的隨機變數。為了讓事情更有趣一點，我們在每個地區採用不同的R0值進行模擬，其中R0值服從均值為4的伽瑪分布:

現在我們來討論所建立的模型：

βk,t是t時刻k地的轉移率，α是刻畫出行方式傾向的參數。

上述的模型十分簡潔：為了求得t + 1時刻的j地的尚未感染或易感的人數，我們需要從t時刻的j地的尚未感染或易感的人數中減去j地本地感染的人數（第一個方程的第二部分），還要減去從其他地方來到j地的感染者數，這些外來的感染者通過其傳輸率進行加權計算（第一個方程的第三部分）。由於總人口數Nj = Sj + Ij + Rj，我們需要將減去的部分移至感染組，同時將治癒的部分移至Rj,t+1(第二個和第三個方程)。

模擬建立

在此分析中，我們將使用由當地共乘公司gg提供的GPS數據獲得的一個典型日的總OD流量矩陣作為埃里溫市交通模式的代表。

接下來，我們需要每個250×250m網格單元的人口計數，通過按比例縮放提取的流量計數來近似計算，從而使不同位置的總流入量之和接近埃里溫市110萬人口的一半。這是一個大膽的假設，但對結果影響不大。

減少公共交通？

第一次模擬，我們將模擬背景設定為一個高度依賴公共交通的未來城市，設定流動率α=0.9：

可以看到，經過大約8-10天左右的時間感染人數比例迅速增加至70%，達到峰值，但此時僅有小部分（約10%）的人康復。至100天時，疫情逐漸緩解，康復人數比例達到了驚人的90%！現在，我們再來看一下如果將公共交通強度α降低至0.2時，是否有利於緩解傳染病的傳播。這可以解釋為採取嚴厲措施來降低城市流動性（例如實施宵禁），或者增加私家車出行比例，以減少人們出行期間感染的機率。

可以發現在這種假設下，疫情在16至20天左右到達頂峰，峰值感染人數比例明顯降低（約45%），並且此時康復人數為之前的兩倍（約20%）。在疫情結行將結束時，易感染人群比例也是之前的兩倍（約24% vs. 約12%），這意味著更多的人躲過了這場疫情。正如人們所期望的，通過實施嚴厲的管控措施來臨時降低城市的流動性對於減少傳染病傳播有明顯作用。

隔離熱門區域？

接著，再來看另一個直觀想法——隔離一些關鍵區域能否得到預期的效果。為了測試這一想法，先挑選人流量位於前1%的區域：

接著完全限制這些區域的進出，建立有效的隔離制度。從這張圖我們可以看出，在埃里溫市這些位置主要位於市中心，另兩個位置是兩家最大的購物商場。將α取中間值，即0.5，我們得到如下結果：

感染人數比例的峰值更小了（約35%），並且更重要的是，在疫情行將結束時，大約一半的人未被感染，說明該種方法能夠幫助人們有效的降低感染風險！

如下動圖顯示了高度依賴公共交通場景下的結果：

結論

該實驗絕不是說我們已經構建了準確的傳染病模型（甚至模型中不涉及任何傳染病學的基礎知識），我們的目標是在傳染病爆發時能夠即時了解城市交通網路對傳染病傳播的影響。

隨著人口密度、流動性和互動性的增強，我們的城市更容易發生「黑天鵝」事件，並且變得更加脆弱。例如，我們從這個模型可以發現在關鍵地區實施隔離制度或者採取嚴苛的措施來控制人員流動能夠在疫情期間發揮巨大作用。

但還有一個十分重要的問題就是，如何在執行這些措施期間，使得城市功能和經濟的損壞最小化？

此外，傳染病傳播機制也是一個活躍的研究領域，該領域的成果必須要滲透並整合到城市規劃、政策制定和城市管理當中，以使我們的城市更安全更抗打擊。

上述模型程式碼如下：

import numpy as np    # initialize the population vector from the origin-destination flow matrix    N_k = np.abs(np.diagonal(OD) + OD.sum(axis=0) - OD.sum(axis=1))    locs_len = len(N_k)                 # number of locations    SIR = np.zeros(shape=(locs_len, 3)) # make a numpy array with 3 columns for keeping track of the S, I, R groups    SIR[:,0] = N_k                      # initialize the S group with the respective populations      first_infections = np.where(SIR[:, 0]<=thresh, SIR[:, 0]//20, 0)   # for demo purposes, randomly introduce infections    SIR[:, 0] = SIR[:, 0] - first_infections    SIR[:, 1] = SIR[:, 1] + first_infections                           # move infections to the I group      # row normalize the SIR matrix for keeping track of group proportions    row_sums = SIR.sum(axis=1)    SIR_n = SIR / row_sums[:, np.newaxis]      # initialize parameters    beta = 1.6    gamma = 0.04    public_trans = 0.5                                 # alpha    R0 = beta/gamma    beta_vec = np.random.gamma(1.6, 2, locs_len)    gamma_vec = np.full(locs_len, gamma)    public_trans_vec = np.full(locs_len, public_trans)      # make copy of the SIR matrices    SIR_sim = SIR.copy()    SIR_nsim = SIR_n.copy()      # run model    print(SIR_sim.sum(axis=0).sum() == N_k.sum())    from tqdm import tqdm_notebook    infected_pop_norm = []    susceptible_pop_norm = []    recovered_pop_norm = []    for time_step in tqdm_notebook(range(100)):        infected_mat = np.array([SIR_nsim[:,1],]*locs_len).transpose()        OD_infected = np.round(OD*infected_mat)        inflow_infected = OD_infected.sum(axis=0)        inflow_infected = np.round(inflow_infected*public_trans_vec)        print('total infected inflow: ', inflow_infected.sum())        new_infect = beta_vec*SIR_sim[:, 0]*inflow_infected/(N_k + OD.sum(axis=0))        new_recovered = gamma_vec*SIR_sim[:, 1]        new_infect = np.where(new_infect>SIR_sim[:, 0], SIR_sim[:, 0], new_infect)        SIR_sim[:, 0] = SIR_sim[:, 0] - new_infect        SIR_sim[:, 1] = SIR_sim[:, 1] + new_infect - new_recovered        SIR_sim[:, 2] = SIR_sim[:, 2] + new_recovered        SIR_sim = np.where(SIR_sim<0,0,SIR_sim)        # recompute the normalized SIR matrix        row_sums = SIR_sim.sum(axis=1)        SIR_nsim = SIR_sim / row_sums[:, np.newaxis]        S = SIR_sim[:,0].sum()/N_k.sum()        I = SIR_sim[:,1].sum()/N_k.sum()        R = SIR_sim[:,2].sum()/N_k.sum()        print(S, I, R, (S+I+R)*N_k.sum(), N_k.sum())        print('n')        infected_pop_norm.append(I)        susceptible_pop_norm.append(S)        recovered_pop_norm.append(R)

相關報道：

https://towardsdatascience.com/modelling-the-coronavirus-epidemic-spreading-in-a-city-with-python-babd14d82fa2