我近半年至少被三种传染病仿真跑出来像直线的情况折磨过模型参数设置不合理、潜伏期逻辑写错、还有一种是压根没想清楚SEIR和SIR到底差在哪。NetLogo这类基于个体的建模工具其实特别适合把SEIR从一组公式变成一群小人的动态行为但很多人一上手就被海龟变量和感染半径搞晕了。这篇内容就是我基于NetLogo实现SEIR传染病仿真程序的全过程总结主线包括为什么SEIR这类仓室模型适合用NetLogo搭、四个仓室状态怎么用代码驱动流转、如何设计对照组跑参数扫描、以及仿真曲线到底该怎么解读。适合正在做人群行为模拟、流行病学入门建模、或者准备用NetLogo交仿真作业/做课程设计的读者。如果你只想看某个部分的代码逻辑或参数表可以直接跳到对应章节。1. 为什么偏偏用NetLogo搭建SEIR从仓室模型到个体交互的跨越很多人在接触传染病模型时第一反应就是用Python或者Matlab写一个SEIR微分方程组四行代码解完画出四条曲线。这当然是最快的路径但它把传染病在人群中的传播这个复杂过程压缩成了纯粹的数学产物没有空间概念没有个体差异没有接触网络更没有某个人为什么被感染而另一个人没有的这种微观机制。NetLogo的定位恰恰相反。它天生就是为基于个体建模设计的每一个turtle海龟/个体都可以被赋予状态、位置、移动规则和交互规则。你可以让易感者在空间里自由游荡然后设定感染者路过时按一定概率传染它们。这种建模方式从底层逻辑上就更接近真实的传播过程。用SEIR模型本身作为例子来谈S、E、I、R四个仓室在微分方程里只是四个变量但到了NetLogo里它们是每个个体自身携带的状态标签。同样是从E变成I这一件事微分方程里发生的是一瞬间的速率变化而NetLogo里则是一个个体在潜伏期计时到达终点后的状态转变。这个区别非常关键——它决定了你能在模型里问出什么样的问题。我实际整体比较下来NetLogo做SEIR的优势有三个空间接触关系明确每个个体在二维空间里移动被感染只能发生在有效接触半径内天然带有空间聚集效应。参数可以个体化不同个体的移动速度、潜伏周期、恢复时间都可以设置成随机分布而不是单一常数。实验做批处理方便BehaviorSpace可以在不修改代码的前提下批量跑几十组参数组合输出对比数据。当然NetLogo也有明显的短板它在数值精度上远不如专业数值计算工具大规模模拟比如几十万个体时运行速度会明显变慢而且随机性较强单次运行结果不具备代表性。所以我的建议是你如果真的只是要算一个R0精确值用Python求解方程更合适但如果你想观察空间隔离、人员流动模式对疫情曲线的影响NetLogo就是很好的选择——这也正是它进入仿真课程体系的原因。2. SEIR的四个仓室到底在流动什么数学基础与代码概念对齐在动笔写NetLogo代码之前我强烈建议先花半小时把SEIR的数学结构搞明白否则后面写代码时你会反复遇到这段逻辑该放在哪个状态里的困惑。2.1 四个仓室的含义与流转方向SEIR模型把人群划分为四个状态状态英文全称含义在本仿真中的表现SSusceptible易感者尚未接触病原体绿色个体可被感染EExposed潜伏者已接触病原体但尚不具备传染性黄色个体带潜伏计时IInfectious感染者具备传染能力红色个体可传染其他人RRecovered恢复者康复且获得免疫力蓝色个体不再参与传播状态流转方向只有一个S → E → I → R。模型的核心假设是感染者在恢复后就获得抗体不会再次感染个体在潜伏期内不传染他人。这两个假设让模型变得简洁清晰但与现实中的很多传染病并不完全吻合这点我在后面第6章再展开谈。从微分方程的角度看SEIR可以用下面四个方程描述dS/dt -β·S·I/NdE/dt β·S·I/N - σ·EdI/dt σ·E - γ·IdR/dt γ·I其中β是有效传播率σ是潜伏者转化为感染者的速率σ 1/潜伏期时长γ是感染者的恢复速率γ 1/传染期时长N是总人口数。2.2 从方程到NetLogo的概念映射很多人在这一步就卡住了方程是连续的、基于总人口比例的而NetLogo是个体离散的、基于交互的。这两者怎么统一我的做法是做一个概念映射表微分方程元素在NetLogo中的对应实现β有效传播率由 传播概率 × 有效接触次数 复合而成σ潜伏转感染速率每个E状态个体的潜伏期倒计时γ感染恢复速率每个I状态个体的感染期倒计时人口总量N模型中创建的海龟总数S/I/R的比例用全局变量实时统计四类状态的数量也就是说在NetLogo中我们不再显式声明β和γ这两个速率常数而是把它们拆解为更底层的个体行为一只有传染性的海龟在移动中进入另一个易感海龟的感染半径的概率有多大、一个潜伏者经过多少天转为感染者、一个感染者经过多少天康复。这是ABM和ODE最本质的差别——宏观参数被微观行为替代了。我最开始写代码的时候犯过的一个错误是把潜伏期和恢复期直接写成了全局变量所有个体步调一致地转换状态。跑出来的结果是曲线抖动严重因为真实世界里每个个体的潜伏期和恢复期应该存在个体差异。后来改成从正态分布里为每只海龟抽取一个具体时长曲线立刻平滑了很多。别小看这个改动它体现的是ABM建模的精髓差异是世界的常态。2.3 基本再生数R0在NetLogo场景下的理解在传统SEIR模型中R0 ≈ β/γ含义是一个感染者在完全易感人群里平均能传染多少人。但在NetLogo这种基于空间接触的模型里β不是一个直接设定的数值而是一个涌现结果——它由传播概率、感染半径、个体移动速度、空间密度共同决定。我建议在NetLogo版本中用这样一个估算公式R0 ≈ 传播概率 × 平均有效接触人数 × 感染持续天数其中平均有效接触人数可以粗略地用感染半径内的平均海龟数量来近似。为什么说是近似因为NetLogo里的接触是动态的、非均匀的海龟在空间中的分布不可能是完全均匀的这导致真实R0会比估算值略高或略低这是个正常现象也是基于个体建模的固有特性。3. 从空白界面到能跑的模型核心代码逻辑与构建顺序NetLogo建模最忌讳一上来就写go函数。我自己的经验是先把setup的部分打磨好因为setup决定了初始状态是否正确而几乎所有仿真跑出奇怪结果的问题都出在初始状态上。3.1 全局参数与个体变量的声明这一步相当于给模型搭好基础设施。我在NetLogo的Code标签页里写的第一段代码就是变量的声明globals [ S-count E-count I-count R-count total-infected ] turtles-own [ state ; 取值 S E I R latent-period ; 该个体的潜伏期时长 infectious-period ; 该个体的感染期时长 latent-timer infectious-timer ]这段代码的含义非常直接四个全局变量用来统计各状态个体数量海龟自身的四个属性用来跟踪个体状态和个体差异。latent-period和infectious-period在设置时用正态分布抽样生成后续不会再改变经纬两个timer则是随着时间步进递增的计数器。3.2 初始化世界海龟初始状态设置to setup clear-all set-default-shape turtles person create-turtles initial-population [ setxy random-xcor random-ycor set state S set color green set latent-period 0 set infectious-period 0 set latent-timer 0 set infectious-timer 0 ] ask n-of initial-infected-number turtles [ set state I set color red set infectious-period random-normal infectious-mean infectious-sd ] reset-ticks update-counts end这里有一个值得注意的细节我在很早期建模时踩过坑**如果直接把所有初始感染者设置为潜伏很多天之后爆发你会发现最开始几天一个红色个体都没有曲线像是死了一样。**更合理的初始条件是先放几个已经在传染状态的个体这样才能让传播过程立刻启动。如果你确实想模拟零号病人刚被感染的场景那么初始感染者应该设置为E状态而不是I状态。这是两个不同的实验场景对应不同的初始条件不要混在一起。在初始化中还有一种常被忽视的细节——海龟之间的空间重叠。NetLogo默认随机放海龟完全可能出现两只海龟同时占用同一个坐标点。对一般空间模型来说这无所谓但在接触传播模型里两个个体重叠意味着它们必然接触传播概率会被无形放大。我通常会在setup后加一行ask turtles [ fd 1 ]之类的轻微位移或者直接用setxy random-xcor random-ycor加while any? other turtles-here [ setxy random-xcor random-ycor ]来消除重叠。这个细节对结果的影响显而易见。3.3 主循环移动、接触传播、状态转换go过程是整个模型的心脏它每一tick要做三件事让个体移动、让感染者传染接触者、让所有个体推进自己的状态计时器。to go move-turtles transmit advance-state update-counts tick if (I-count E-count) 0 [ stop ] end to move-turtles ask turtles [ rt random 90 - 45 fd move-speed ; 四周环绕世界不产生边界聚集 if xcor max-pxcor [ set xcor min-pxcor ] if xcor min-pxcor [ set xcor max-pxcor ] if ycor max-pycor [ set ycor min-pycor ] if ycor min-pycor [ set ycor max-pycor ] ] end后续逻辑to transmit ask turtles with [ state I ] [ let neighbors-in-radius turtles in-radius infection-radius with [ state S ] ask neighbors-in-radius [ if random-float 100 transmission-probability [ set state E set color yellow set latent-period random-normal latent-mean latent-sd set latent-timer 0 ] ] ] end to advance-state ask turtles [ if state E [ set latent-timer latent-timer 1 if latent-timer latent-period [ set state I set color red set infectious-period random-normal infectious-mean infectious-sd set infectious-timer 0 ] ] if state I [ set infectious-timer infectious-timer 1 if infectious-timer infectious-period [ set state R set color blue ] ] ] end这段代码有几个关键逻辑需要注意移动是随机的每只海龟每tick随机转向-45到45度然后前进一个固定步长。这种随机游走模型既简单又满足大部分人群移动场景如果要模拟居家隔离或交通管制只需要把move-speed在特定条件下设为0即可。传播发生在感染者的角度代码里是用ask turtles with [ state I ]做主体然后查询它们半径内的易感者这样做效率高而且逻辑清晰。不要把逻辑反过来说成易感者找感染者虽然结果一样但代码可读性差很多。同一只海龟单次tick只能被感染一次如果一个感染者周围有两只感染者它会被第一步感染然后在第二个感染者的询问中它已经变成了E状态不再符合with [ state S ]条件所以不会二次感染。这符合SEIR假设。用random-float 100 transmission-probability表示概率NetLogo中random-float 100返回0到100之间的随机数把它和传播概率比如5到20之间比较就相当于按百分比概率做判断。当然你也可以用random-float 1 传播概率/100但前者更直观一些。3.4 统计与绘图to update-counts set S-count count turtles with [ state S ] set E-count count turtles with [ state E ] set I-count count turtles with [ state I ] set R-count count turtles with [ state R ] set total-infected R-count I-count E-count end这一段没啥技术含量但必须写在go循环里否则绘图和输出的数据不会实时更新。BehaviorSpace的输出变量也是直接引用这几个全局变量。4. 用BehaviorSpace跑参数扫描一次完整的实验设计NetLogo有一个被很多人低估的功能叫BehaviorSpace它允许你在不改代码的前提下批量执行不同参数组合。对SEIR模型来说这几乎是必备技能——因为单次仿真结果受随机性影响太大不看多组参数的平均趋势完全无法结论。4.1 一个可执行的实验方案假设我们要研究潜伏期长短对疫情峰值的影响可以在BehaviorSpace中这样设置参数扫描值latent-mean3, 5, 7infectious-mean7transmission-probability10infection-radius1.2initial-population800initial-infected-number5move-speed0.4每个参数组合重复10次然后取平均值。这样一个实验跑下来就是3×1030次仿真看起来不多但如果每次跑到曲线稳定需要3000多个tick总计算量也很可观。所以我在跑BehaviorSpace之前通常会先单次运行一遍确定需要的tick数上限然后在BehaviorSpace中把终止条件设为I-count E-count 0这样仿真会在疫情结束后自动停止省掉了大量浪费的计算时间。4.2 BehaviorSpace的结果导出与预处理BehaviorSpace运行完成后结果会以表格形式呈现每一行是一次仿真每一列是一个输出变量。默认输出包含了[run number]、[step]和你选定的参数及统计量。我一般会勾选上这些指标S-count、E-count、I-count、R-count在每个时间步的快照total-infected的最终值用于对比整体爆发规模最高I-count发生的时间步用于衡量疫情峰值的到达时刻对于导出后的CSV文件我的处理习惯是在Python里用pandas按run number分组提取每个仿真的峰值I-count、峰值时间、最终累计感染数等聚合指标然后再用箱线图展示不同潜伏期参数下的差异分布。这一步才是真正产出结论的关键——没有统计上处理的多组对比几组原始曲线放在一起只是一堆乱线。4.3 参数敏感性怎么看扫描结果做敏感性分析时我总结出了三个最常用的观察维度峰值的出现时间随参数变早还是变晚这是最直观的效应。潜伏期越短E→I转化越快I峰值出现得就越早。在仿真里通常可以看到latent-mean3的那组曲线在感染数上冲到200人时latent-mean7的组才刚刚进入指数增长段。累积感染总数随参数变化是否明显如果潜伏期参数从3天改到7天最终感染总数几乎没有变化那说明系统对潜伏期不敏感瓶颈在传播阶段。如果变化明显说明提前隔离潜伏者的干预策略对累计感染数有实质性影响。曲线上锯齿的幅度锯齿多、抖动大通常意味着感染半径或个体移动速度设置过大接触机会不均匀锯齿少说明模型更接近混合均匀假设宏观曲线更光滑。5. 仿真曲线能告诉我们什么SEIR输出解读与常见误区很多人在仿真跑完之后只会把四条曲线截图发出去完全不解释这些线的形状意味着什么这其实浪费了模型一半的价值。SEIR曲线里面蕴含着大量的信息关键在于你怎么读它。5.1 四条典型曲线段的含义如果你完整跑一次SEIR通常会得到下面这样的宏观特征潜伏期平坦段前几个tickE曲线上升但I曲线几乎没有变化因为所有感染者都处于潜伏阶段还没有传染能力。这个阶段的长度大致就对应初始感染者的平均潜伏期。如果你看到I曲线在开始一段时间内完全是一个平台不要以为程序坏了这是SEIR的典型特征叫潜伏延滞。指数增长段一旦第一批潜伏者转为感染者传播进入正反馈过程I曲线呈近似指数增长。这段曲线越陡说明传播概率、感染半径或者人口密度越高。如果你在日志里记录每个tick的I-count用指数拟合这段曲线算出的增长率可以用来估算实际的有效传播率。这个拟合操作只需要几行Python代码就能做强烈推荐一试。峰值与下降段当易感者数量下降到一定程度新增感染数逐渐少于恢复数I曲线开始下降。注意I曲线的峰值比新增感染数的峰值晚出现几个tick因为感染人数是存量而不是流量——这个时间差恰好等于平均传染期的长度。长期尾段曲线末端会拖出一个长尾这是因为部分个体潜伏期较长或感染期较长它们会在别人都已恢复后再完成自己的病程。长尾越明显说明个体差异标准差设置得越大。5.2 一个很多人犯的解读错误把潜伏期0和潜伏期0的两条I曲线放在一起看的时候大部分人觉得区别只是峰值晚来几天。但实际上潜伏期的影响远不止是时间平移更长的潜伏期意味着感染者进入I状态更晚、每个I状态个体面对的人口密度可能更小因为空间中的易感者被之前几代感染者消耗了一部分最终累积感染总数会有可测量的差异。这个效应在传染病防控中的实际含义就是压低潜伏期的干预措施比如对密接者提前隔离不仅推迟峰值还会压低峰值高度和最终感染规模。另外一个常见误区是直接比较两条看似相同参数的曲线。NetLogo默认的随机数种子不同每次运行结果都会有波动如果你不通过BehaviorSpace多次重复取平均或者至少固定随机种子再对比看到的两条曲线差异可能完全来自随机噪声而不是参数不同导致的。所以我在每次对比实验前都先跑5次同一参数看波动范围再决定是不是真的需要保存好几轮BehaviorSpace的数据。5.3 从仿真中识别合理与不合理的曲线我刚调试模型的时候屡次遇到这类问题某个参数组合下I曲线在短时间内冲到一个极高值又迅速归零曲线形状像一根针。这种尖峰型曲线通常说明模型处于过度接触状态——感染半径太大或移动速度太快导致每个感染者几乎每tick都在感染大半个地图上的个体传染源很快就耗尽了。对应的解决办法是缩小感染半径、调低传播概率或者增大地图面积。反过来还有一种情形I曲线以极低的速率爬升跑了几千tick还没有形成明显的峰值E曲线却一直有累积。这通常是因为初始感染者太少加上传播概率太低传播网络基本处于准消亡状态。这种模型设置如果是为了研究疫情在低接触条件下能否自愈是有意义的但如果你的目的是观察典型SEIR爆发过程那就需要适当增加初始感染人数或传播概率。6. 模型的短板在哪里三个必须知道的边界条件NetLogo版SEIR模型在课堂演示和概念验证层面非常好用但如果把它当作严谨的流行病学预测工具那就超出了它的能力范围。我在这里坦诚分享一下我在实际使用中观察到的三个最明显的短板。6.1 无差别随机接触忽略了现实的社交网络NetLogo的基本模型中所有个体在整个空间内随机游走任何两个个体只要距离够近就会发生传播。这显然与现实不符在实际人群中每个人的接触对象有很强的异质性——家庭成员、同事、同学、陌生人接触频率完全不同。我在NetLogo中尝试过改进接触结构比如用网络扩展库NW为个体创建社交关系边让传播更容易发生在熟人之间把陌生人接触作为低概率事件。加入这个设定之后仿真结果明显更合理指数增长段更平滑峰值更低疫情存续时间更长。如果你要做进一步研究建议关注带网络结构的ABM传染病模型这是比纯随机接触更贴近实际的方向。6.2 空间移动模型过于简单无法模拟物理隔离等干预代码里的移动规则就是随机游走世界环绕这其实隐含了一个假设人群完全自由流动没有任何物理边界或干预策略。如果你想模拟封闭管理、学校停课、居家办公这些调节措施就得显式修改move-speed或者引入特定区域禁止进入的空间限制。这个改动本身不复杂但我建议你想清楚一件事**你是在研究干预策略本身的机制还是只想验证这个策略有效**后者用ODE模型加几行参数修改就能完成没必要动用空间仿真。只有当你关心干预措施在不同居住密度、不同空间分布下的效果差异时ABM的空间能力才真正派上用场。6.3 免疫假设过于绝对R状态并非永久免疫SEIR模型里R终态恢复者不再参与传播。但在现实中许多病原体的免疫保护会随时间衰减或者因为变异株产生免疫逃逸导致恢复者重新变回易感者。NetLogo实现这一点其实很简单给R状态的个体加一个免疫持续时间到时自动重置为S状态并恢复绿色或者通过一个全局变量waning-immunity-duration来控制。加上这个机制后模型会变成SEIRS呈现周期性的多波流行特征——如果你研究的是流感这类容易反复感染的疾病这一步几乎是必需的而原始SEIR没法处理这个问题。最后再说一个我总结出来但经常被忽略的小提醒模型验证永远要比模型参数本身多花时间。你在NetLogo中随意改一个数字也许就会跑出完全不同的疫情所以归根结底与其问模型跑得漂不漂亮不如先问我设置的参数在真实世界中有没有对应的测量依据。参数来源可靠模型结论才站得住脚。这是我在做了一整轮SEIR仿真之后最想压在文末的一句实在话。
