[[315177]]

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

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

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

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

城市流動

高效持續(xù)的城市流動對現(xiàn)代城市的運轉(zhuǎn)舉足輕重，并直接影響城市的宜居性和和經(jīng)濟輸出（GDP）。但是，當流行病爆發(fā)時，城市流動卻火上澆油，推動了疾病的傳播。

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

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

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

建模流行病

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

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

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

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

在此次模擬中，時間是離散變量，因為系統(tǒng)狀態(tài)是在每天的基礎(chǔ)上建模的。假設(shè)第 t 天地點 j 處的所有人均未感染該疾?。ㄒ赘腥巳海?，則此處出現(xiàn)病例的概率為：

其中 β_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 處的人口規(guī)模。

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

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

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

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

現(xiàn)在我們可以得到模型動態(tài)：

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

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

模擬設(shè)置

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

減少公共交通出行？

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

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

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

我們看到峰值在第 16-20 天到來，感染者比例要小得多（約 45%），治愈者比例是之前的 2 倍（約 20%）。到流行病結(jié)束時，易感人群比例是之前的 2 倍（大約 24% vs. 12%），這意味著更多人從該流行病的魔爪下逃脫。和預計的一樣，采取有力措施暫時減少城市交通出行對疾病擴散情況有巨大影響。

隔離熱門場所？

現(xiàn)在，我們來看另一個直觀思路：完全隔離少數(shù)關(guān)鍵熱門地區(qū)可以達到疾病控制預期。為此，本文選擇了交通流量處于城市前 1% 的地點，并且完全阻止這些地點的人口流入和流出，有效地建立起隔離區(qū)。

從上圖中可以看到，埃里溫市的這些地點大部分位于市中心，還有兩個是最大的購物中心。選擇 α = 0.5，得到：

我們可以看到，峰值處的感染者比例更低（約 35%），最重要的是，在疫情結(jié)束時，大約一半的人口沒有被感染！

下圖展示了公共交通比例較高時的病毒傳播態(tài)：

結(jié)論

本文并未執(zhí)行精確的流行病建模（甚至對流行病學也僅具備基礎(chǔ)的了解），而是為了粗略了解在傳染病爆發(fā)時小世界網(wǎng)絡(luò)效應對城市的影響。隨著人口密度、流動性和活力的增長，城市更容易遭遇「黑天鵝事件」，也更加脆弱。

如果沒有高效的危機處理能力和機制，城市再智能再可持續(xù)也都沒有了意義。例如，我們看到在重點地區(qū)設(shè)置隔離區(qū)，或者采取強力措施控制人口流動，在衛(wèi)生危機中能夠起到重要作用。

當然，傳染病的確切擴散機制仍是活躍的研究領(lǐng)域，未來我們可以看到更多利用數(shù)據(jù)進行模擬的方法。這種研究成果可以和城市規(guī)劃、政策制定和管理整合，從而使城市更安全、更穩(wěn)健。

上述模擬的代碼如下所示：

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)