最新技術為你播報
轉自大數據文摘
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)
相關報道：
更多案例源碼私信小編01領取完整項目源代碼~
閱讀更多 地表嘴強程序員 的文章
關鍵字: 想象偉大的一平方公里 模擬 推測
