自去年12月以来，新冠病毒引发的疫情对城市活动造成了很大的影响。如何准确了解病毒的传播方式，从而帮助城市更好地提出措施？使用建模方法也可以发挥作用。本文是一篇Python教程，教你如何模拟流行病在家中的传播。本文以亚美尼亚共和国首都埃里温市为例，对冠状病毒在城市中的传播进行数学建模和模拟，观察城市流动模式对病毒传播的影响。读者也可以根据文末示例代码自行使用。目前尚未发现经医学研究标准确认的特效药。此外，还有很多重要的流行病学指标仍然未知，比如基本感染数R0（在流行病学中，基本感染数是指在没有外界干预、没有免疫力的情况下感染某种疾病的人数。其他人的平均人数患有传染病的人感染了该疾病）。新型冠状病毒还有很多不确定性。当今的全球社会比以往任何时候都更加互联和流动。由于小世界效应，此类流行病构成了重大的全球威胁。有人推测，如果2020年发生全球性灾难事件（大致定义为超过1亿人伤亡），那么最有可能的原因是某种流行病，而不是核灾难、气候灾难等。城市化，因为人口稠密的城市成为疾病传播网络中的传播节点，因此变得极其脆弱。本文探讨了当一个城市遭受大流行病袭击时会发生什么，应立即采取哪些措施，以及对城市规划、政策制定和管理的影响。本文以亚美尼亚共和国首都埃里温市为例，对冠状病毒在城市中的传播进行数学建模和模拟，观察城市流动模式对病毒传播的影响。城市交通高效和可持续的城市交通对于现代城市的运作至关重要，并直接影响城市的宜居性和经济产出(GDP)。但当流行病爆发时，城市流动性却火上浇油，助长了疾病的传播。让我们首先看一下埃里温的起点-终点(OD)交通流在统一笛卡尔网格上以聚合方式形成的网络，以了解城市流模式的空间结构：仔细观察整体网格单元对于流入，您会看到一个单中心的空间组织，中心的单元格每天流入量很高：现在，假设流行病在城市的随机位置爆发。它将如何传播？可以做些什么来控制它？流行病建模为了回答上述问题，我们首先构建一个简单的分区模型来模拟传染病在城市中的传播。在传染病爆发期间，其传播动态存在显着差异，具体取决于首次感染发生的地点及其与城市其他地区的连通性。这是近期城市人口流行病学数据驱动研究的最重要结论之一。但是，尽管病毒传播情况不同，但城市需要采取类似措施来控制传染病，并加强城市规划和管理。个人很难运行一个流行病模型，所以本文旨在展示流行病在城市传播的一般规律，而不是建立一个精确准确的流行病模型。本文将沿用Nature文章《Multinational patterns of seasonal asymmetry in human movement influence infectious disease dynamics》中介绍的方法，根据具体需求对经典的SIR模型进行修改。该模型将城市人口分为三个部分。对于时间点t的每个位置i，三个隔间是：S_i,t：未感染或易感人群（易感人群）；I_i,t：感染该疾病的传染人群（infectedgroup）；R_i,t：感染该病后，因痊愈或死亡而脱离感染组的种群。这群人不再感染病毒，也没有传染给其他人。在此模拟中，时间是一个离散变量，因为系统状态是按日建模的。假设第j个地点在第t天没有人感染该病（易感人群），则此处出现病例的概率为：其中β_t为第t天的感染率，m_j,k表示第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为例。然而，我们应该注意，这实际上是一个随机变量，报告的数字只是估计值。本文尝试在每个位置用从均值为4的候选伽马分布采样的不同R0值进行模拟：现在我们可以得到模型动力学：其中β_k,t表示（随机）感染率，α系数表示城市公共交通工具利用率。上述等式所描述的模型的动力学非常简单：以t+1日地点j为例，我们需要从易感人群中去除地点j的感染者比例（第一个公式的第二项）人口S_j,t和来自城市其他地方的感染者比例由感染率β_k,t（第一个公式的第三项）加权。由于总人口N_j=S_j+I_j+R_j，我们需要将刚刚移除的部分移至感染组（第二个公式），将治愈的部分移至R_j,t+1（第三个公式）。模拟设置为了分析方便，本文以聚合方式使用来自当地拼车公司ggGPS数据的每日OD交通流矩阵，并用它来表示埃里温市的交通模式。接下来，我们需要每个250×250m网格单元中的人数。我们通过缩放提取的OD交通流量数字来近似这个数字，以便不同位置的总流入量接近埃里温市总人口的一半（110万）。这实际上是一个大胆的假设，但由于改变这部分仍然给出非常相似的结果，本文坚持下去。减少公共交通出行？在第一个模拟中，假设可持续的公共交通主导城市未来的流动性，即α=0.9：从上图我们可以看到感染者的比例立即攀升，大约在8-10天达到峰值，近70%的人口被感染，只有少数人（约10%）康复。到病毒消退的第100天，康复率达到了惊人的90%！我们再来看看，将公共交通密度降低到α=0.2，对减缓疫情传播有没有效果。这也可以理解为：采取强有力的措施减少城市交通出行（如戒严），或者增加私家车的比例，以减少人口流动过程中的感染机会。我们在第16-20天看到一个高峰，感染者的比例要小得多（约45%），而治愈者的比例是后者的两倍（约20%）。到疫情结束时，易感人群增加了一倍（约24%对12%），这意味着更多人逃脱了疫情。正如预期的那样，暂时减少城市交通出行的强有力措施对疾病的传播产生了巨大影响。隔离热门场所？现在，我们再来看另一个直观的想法：完全隔离几个重点热点，就能达到疾病控制的预期。为此，本文选择了城市交通流量前1%的位置，并完全阻断这些位置的人员流入和流出，有效地创建隔离区。从上图可以看出，埃里温的这些地点大部分都位于市中心，其中有两个是最大的购物中心。选择α=0.5，我们得到：我们可以看到高峰时期的感染者比例较低（约35%），最重要的是，在疫情结束时，大约有一半人口未被感染！下图显示了公共交通比例高时病毒传播的状态：结语本文并没有进行精确的流行病建模（即使对流行病学有基本的了解），而是粗略地了解小世界网络效应的影响在城市。随着人口密度、流动性和活力的增长，城市更容易发生“黑天鹅事件”，也更加脆弱。没有高效的危机处理能力和机制，再聪明、再可持续的城市也毫无意义。例如，我们看到，关键地区的隔离或控制人口流动的强有力措施可以在健康危机中发挥重要作用。当然，传染病传播的确切机制仍然是一个活跃的研究领域，未来我们可以看到更多使用数据进行模拟的方法。这项研究的结果可以与城市规划、政策制定和管理相结合，使城市更安全、更有活力。上面描述的模拟的代码如下所示：importnumpyasnp#initializethepopulationvectorfromtheorigin-destinationflowmatrixN_k=np.abs(np.diagonal(OD)+OD.sum(axis=0)-OD.sum(axis=1))locs_len=len(N_k)#numberoflocationsSIR=np.zeros(shape=(locs_len,3))#makeanumpyarraywith3columnsforkeepingtrackoftheS,I,RgroupsSIR[:,0]=N_k#initializetheSgroupwiththerespectivepopulationsfirst_infections=np.where(SIR[:,0]<=thresh,SIR[:,0]//20,0)#fordemopurposes,randomlyintroduceinfectionsSIR[:,0]=SIR[:,0]-first_infectionsSIR[:,1]=SIR[:,1]+first_infections#moveinfectionstotheIgroup#rownormalizetheSIRmatrixforkeepingtrackofgroupproportionsrow_sums=SIR.sum(axis=1)SIR_n=SIR/row_sums[:,np.newaxis]#initializeparametersbeta=1.6gamma=0.04public_trans=0.5#alphaR0=beta/gammabeta_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)#makecopyoftheSIRmatricesSIR_sim=SIR.copy()SIR_nsim=SIR_n.copy()#runmodelprint(SIR_sim.sum(axis=0).sum()==N_k.sum())fromtqdmimporttqdm_notebookinfected_pop_norm=[]susceptible_pop_norm=[]recovered_pop_norm=[]fortime_stepintqdm_notebook(范围(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('totalinfectedinflow:',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_infectSIR_sim[:,1]=SIR_sim[:,1]+new_infect-new_recoveredSIR_sim[:,2]=SIR_sim[:,2]+new_recoveredSIR_sim=np.where(SIR_sim<0,0,SIR_sim)#recomputethenormalizedSIRmatrixrow_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')感染ed_pop_norm.append(I)susceptible_pop_norm.append(S)recovered_pop_norm.append(R)
