[[315177]]

自去年12月以來，新型冠狀病毒所引發的疫情已經給城市活動帶來了很大影響。怎樣確切了解病毒的傳播過程，從而幫助城市更好提出措施？使用建模的方法也能起到一些作用。本文是一篇Python教程，教你在家中也可以建模疫情傳播。本文以亞美尼亞共和國首都埃里溫作為案例，對冠狀病毒在該城市中的蔓延情況進行數學建模和模擬，并觀察城市流動模式對病毒傳播的影響。讀者也可根據文末的示例代碼，自己上手使用。

目前，尚未發現經醫療研究標準確認的特效藥。此外，很多重要的流行病學指標仍屬未知，如基本傳染數 R0（在流行病學上，基本傳染數指在沒有外力介入，同時所有人都沒有免疫力的情況下，一個感染到某種傳染病的人，會把疾病傳染給其他多少個人的平均數）。新型冠狀病毒還存在很多的不確定性。

如今的全球社會，聯通度和流動性空前。由于小世界效應，這類流行病對全球造成重大威脅。有人猜想如果 2020 年發生全球性災難事件（粗略定義為傷亡人數大于 1 億），那么最可能的原因是某種流行病，而不是核災難、氣候災難等。全球城市化的快速發展更加強化了這一點，因為人口密集的城市變成了疾病擴散網絡中的傳播節點，因而變得極度脆弱。

本文將探討當一座城市遭受流行病襲擊時會發生什么，應立即采取哪些措施，以及這些措施對城市規劃、政策制定和管理帶來的影響。本文以亞美尼亞共和國首都埃里溫作為案例，對冠狀病毒在該城市中的蔓延情況進行數學建模和模擬，并觀察城市流動模式對病毒傳播的影響。

城市流動

高效持續的城市流動對現代城市的運轉舉足輕重，并直接影響城市的宜居性和和經濟輸出（GDP）。但是，當流行病爆發時，城市流動卻火上澆油，推動了疾病的傳播。

我們先來看埃里溫市起訖點（origin-destination，OD）交通流以聚集的方式在均勻笛卡爾網格上形成的網絡，以便了解城市流動模式的空間結構：

仔細觀察網格單元的總體流入情況，你會看到一個單中心空間組織，位于中心的單元格具備高日流入量：

現在，假設流行病在這座城市的隨機地點爆發了。它將以怎樣的方式擴散？怎么做才能控制住它？

建模流行病

為了解答以上問題，我們先構建一個簡單的房室模型（compartmental model）來模擬傳染病在城市中的擴散。在傳染病爆發時，根據首次感染發生地點及其與城市其它地區的連通性，其傳播動態存在顯著差異。這是近期對城市人口流行病進行數據驅動研究后得出的最重要結論之一。但是，雖然病毒傳播狀況不同，但城市需要采取類似的措施來控制傳染病，以及執行城市規劃和管理。

個人運行流行病模型難度很大，因此本文旨在展示流行病在城市中的一般傳播原理，而不是構建精細準確的流行病模型。本文將按照 Nature 文章《Multinational patterns of seasonal asymmetry in human movement influence infectious disease dynamics》中介紹的方法進行操作，并根據具體需求修改了經典 SIR 模型。

該模型將城市人口分成三室。對于時間點 t 處的每個地點 i，這三室分別是：

• S_i,t：尚未感染或易感人群（易感組）；
• I_i,t：感染該疾病且具備傳染性的人口（感染組）；
• R_i,t：感染疾病后，由于痊愈或死亡從感染組移出的人口。這組人不再感染病毒，也不會將病毒傳播給他人。

在此次模擬中，時間是離散變量，因為系統狀態是在每天的基礎上建模的。假設第 t 天地點 j 處的所有人均未感染該疾病（易感人群），則此處出現病例的概率為：

其中 β_t 是第 t 天的傳染率，m_j,k 表示從地點 k 到地點 j 的流動情況，x_k,t 和 y_j,t 分別表示第 t 天地點 k 處的感染者比例和地點 j 處的易感人群比例，x_k,t = I_k,t / N_k, y_j,t = S_j,t / N_j，其中 N_k、N_j 分別表示地點 k 和 j 處的人口規模。

接下來我們來模擬將疾病帶入易感人群所在地區的隨機過程，I_j,t+1 是伯努利隨機變量，概率為 h(t,j)。

一旦病毒傳播到隨機地點，則它們不僅在這些地點傳播，還會隨流動人口向其他地點擴散。這就是城市流動模式的重要影響（城市流動模式以 OD 交通流矩陣呈現）。

要想形式化感染者擴散病毒的過程，我們需要基本傳染數 R0。R0 = β_t / γ，其中 γ 表示治愈率。截至本文寫作時，2019-nCoV 的 R0 估計值在 1.4 到 4 之間。本文以最壞的情形為例，即 R0 值為 4。不過我們應該注意，這其實是個隨機變量，報告的數字只是估計值。

本文嘗試在每個地點用不同的 R0 值進行模擬，這些 R0 值采樣自候選伽馬分布，平均值為 4：

現在我們可以得到模型動態：

其中 β_k,t 表示第 t 天地點 k 處的（隨機）傳染率，α 系數表示城市公共交通使用率。

上式描述的模型動態非常簡單：以第 t+1 天的地點 j 為例，我們需要從易感人群 S_j,t 中去除地點 j 處的感染者比例（第一個公式的第二項）以及從這座城市其他地方流入的感染者比例，其權重為傳染率 β_k,t（第一個公式的第三項）。由于總人口 N_j = S_j + I_j + R_j，因此我們需要將剛才去除的部分移到感染組（第二個公式），將治愈者移到 R_j,t+1（第三個公式）。

模擬設置

為方便分析，本文以聚集的方式使用來自當地共乘公司 gg GPS 數據的日常 OD 交通流矩陣，并用它表示埃里溫市的流動模式。接下來，我們需要每個 250×250m 網格單元中的人口數。我們通過按比例縮放提取到的 OD 交通流數量來逼近該數字，使不同地點的總流入量接近埃里溫市總人口（110 萬）的一半。這實際上是一次大膽的假設，但是由于對這部分做出更改仍會得到非常類似的結果，因此本文堅持使用該假設。

減少公共交通出行？

在第一次模擬中，假設可持續的公共交通主導城市未來流動，即 α=0.9：

從上圖中，我們可以看到感染者比例立即攀升，在大約 8–10 天處達到峰值，將近 70% 的人口被感染，只有少量人口（約 10%）痊愈。到第 100 天病毒撤退時，治愈者比例達到令人震驚的 90%！

我們再來看，將公共交通的密度減少到 α = 0.2 對緩解流行病擴散是否有影響。這也可以解釋為：采取有力措施來減少城市交通出行（如實行戒嚴），或增加私家車比例，以降低人口流動過程中的感染幾率。

我們看到峰值在第 16-20 天到來，感染者比例要小得多（約 45%），治愈者比例是之前的 2 倍（約 20%）。到流行病結束時，易感人群比例是之前的 2 倍（大約 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)