SIR模型

SIR 模型把一个封闭人群分为易感者、感染者和移出者，用组间转移建立微分方程。它研究一个传播过程怎样因接触而增长，又怎样因易感人群减少而转弱。经典模型的作用是解释这些机制之间的关系；具体地区、具体疾病的判断还需要独立的数据、适用性检验和更多机制。

一天内有哪些人在改变分组

设一个虚构人群共有 $N=1000$ 人。初始有 990 名易感者、10 名具有传染性的感染者，移出组人数为 0。全篇参数与图都是教学假设，没有使用现实疫情数据。

记 $S(t)$ 为时刻 $t$ 的易感人数，$I(t)$ 为当时具有传染性的感染人数，$R(t)$ 为已经移出传播过程的人数。这里移出可以代表康复后在观察期内不再易感，或其他使人不再参与传播的状态；模型不进一步区分这些途径。时间单位取日，观察期间没有出生、外来输入和组外迁移，并假设移出者不会重新回到易感组。

再作两个机制假设。首先，人群均匀混合：一个感染者遇到易感者的比例为 $S/N$。把有效接触频率记为 $\beta$，则一个感染者每单位时间产生新感染的速率为 $\beta S/N$，全部感染者对应的转移率为 $\beta SI/N$。其次，每个感染者以恒定速率 $\gamma$ 移出，所以总移出率为 $\gamma I$。这两项总转移率的单位均为人/日；参数 β、γ 的单位均为每日。

箭头上的量是每单位时间的转移人数，不是已经转移的人数。

每组的变化等于流入减流出，于是 $$S^{\prime }=-\beta \frac{SI}{N},I^{\prime }=\beta \frac{SI}{N}-\gamma I,R^{\prime }=\gamma I.$$ 例如取 $\beta =0.3{\mathrm{d}}^{-1}$、$\gamma =0.1{\mathrm{d}}^{-1}$。在初始时刻，新感染速率是 $0.3\times 990\times 10/1000=2.97$ 人/日，移出速率是 $0.1\times 10=1$ 人/日，所以 $I^{\prime }(0)=1.97$ 人/日。这个小数是连续模型的瞬时变化率，不表示现实中存在“0.97 个人”。

SIR 的分组思路可追溯到 Kermack 与 McKendrick 1927 年的流行病数学研究。此处使用的是常系数、无潜伏期的基本形式。Jeffrey R. Chasnov《Mathematical Biology》§4.3，pp.51–53

守恒与非负性先保证模型有意义

把三条方程相加，两次转移都在一组中减掉、另一组中加回，得到 $$(S+I+R{)}^{\prime }=0.$$ 因此 $S+I+R=N$ 始终成立。

从 $S^{\prime }=-(\beta I/N)S$ 出发，将它看成关于 $S$ 的一阶线性方程，可以写成 $$S(t)=S(0)\exp (-\frac{\beta }{N}{\int }_0^{t}I(q)dq).$$ 同理 $$I(t)=I(0)\exp \left({\int }_0^t(\frac{\beta S(q)}{N}-\gamma )dq\right).$$ 所以非负初值不会产生负人数；正的 $S(0),I(0)$ 在每个有限时刻仍为正。又因为 $R^{\prime }=\gamma I\ge 0$，移出人数单调增加。三者非负且总和固定，也就都不会超过 $N$。右侧是光滑函数，状态始终留在这个有界区域，因此解可以持续向前延拓。

为了减少数字，可改用比例 $s=S/N,i=I/N,r=R/N$。除以 $N$ 后方程为 $$s^{\prime }=-\beta si,i^{\prime }=\beta si-\gamma i,r^{\prime }=\gamma i,s+i+r=1.$$ 这套写法中 $\beta ,\gamma$ 仍以每日为单位。有的教材在人数字式里写 $S^{\prime }=-\tilde{\beta }SI$，它的参数是 $\tilde{\beta }=\beta /N$，单位为每人每日。两种约定均可使用，但不能把它们的参数数值直接交换。

为什么感染者先增长，后来又减少

把感染方程提取 $I$： $$I^{\prime }=\gamma I(\frac{\beta }{\gamma }\frac{S}{N}-1).$$ 定义基本再生数 ${ℛ}_0=\beta /\gamma$。它与人数函数 $R(t)$ 是不同对象。在恒定移出风险的解释下，一个感染者保持感染状态的概率按 $e^{-\gamma t}$ 衰减，平均感染持续时间为 $${\int }_0^{\mathrm{\infty }}{e}^{-\gamma t}dt=\frac{1}{\gamma }.$$ 当周围几乎全是易感者时，有效传播速率约为 $\beta$，相乘得到 $\beta /\gamma$。随着易感比例下降，有效再生数变为 ${ℛ}_{0}s(t)$。

本例 ${ℛ}_0=3$，初始有效再生数为 $3\times 0.99=2.97>1$，感染人数起初增长。当感染人数为正时，增长或下降只取决于 $s$ 是否大于 $1/3$。如果开始时就有 ${ℛ}_{0}s(0)\le 1$，由于 $s$ 只会下降，感染人数不会在之后自行转为增长。

初期若易感人数尚未明显改变，可暂用 $S\approx 990$，得到近似式 $I^{\prime }\approx 0.197I$，即 $I\approx 10e^{0.197t}$。这只描述易感耗竭很小的阶段。继续用同一个指数预测很长时间，会忽略方程中正在下降的 $S$。

参数 β=0.3/日、γ=0.1/日。虚线标出感染人数峰值时刻，约为第 26.63 日。

图上 $I$ 的峰值不是每天新感染数的峰值。新感染速率是 $J=\beta SI/N$；当 $S,I>0$ 时， $$\frac{J^{\prime }}{J}=\frac{S^{\prime }}{S}+\frac{I^{\prime }}{I}=\beta (s-i)-\gamma .$$ 而 $I^{\prime }=0$ 的条件是 $\beta s=\gamma$。此刻 $J^{\prime }/J=-\beta i<0$，新感染速率已经在下降。明确纵轴是哪一个量，才能正确谈论“峰值”。

不先求时间函数，也能算出峰值高度

将 $i^{\prime }$ 除以 $s^{\prime }$，在 $s,i>0$ 处消去时间： $$\frac{di}{ds}=\frac{\beta si-\gamma i}{-\beta si}=-1+\frac{1}{{ℛ}_{0}s}.$$ 从初值 $(s_0,i_0)$ 积分到 $(s,i)$，得到 $$i-i_0=-(s-s_0)+\frac{1}{{ℛ}_0}\ln \frac{s}{s_0}.$$ 也就是说，每条轨迹上 $i+s-\ln s/{ℛ}_0$ 为常数。这个关系能画出 $(s,i)$ 相轨迹，并从状态间的关系求峰值。

本例峰值处 $s=1/3$，所以 $$i_{max}=0.01+0.99-\frac{1}{3}+\frac{1}{3}\ln \frac{1/3}{0.99}=\frac{2}{3}-\frac{1}{3}\ln 2.97\approx \mathrm{0.3038126824.}$$ 乘以 1000，感染人数峰值约为 303.81 人；当时移出比例为 $1-1/3-i_{max}\approx 0.3628539843$。峰值时间不能仅从这条相轨迹关系直接读取；对时间方程作数值积分，并定位 $s(t)=1/3$ 的穿越事件，才得到约 26.63 日。

箭头沿时间增加方向，因此 s 从右往左减少。相轨迹给出峰值高度；时间曲线给出到达峰值的时间。

传播结束后，为什么仍有易感者

由于 $r^{\prime }=\gamma i$ 且 $r\le 1$，总积分 ${\int }_0^{\mathrm{\infty }}i(t)dt$ 有限。又因系统右侧在有界状态区域内有界，$i$ 不会形成越来越窄却保持固定高度的尖峰；有限积分结合这种平滑性推出 $i(t)\to 0$。这里的结束是趋近意义，并非感染人数在某个有限时刻突然精确变成零。

另一方面，利用 $s^{\prime }=-\beta si$ 和 $r^{\prime }=\gamma i$，可得 $$\ln \frac{s(t)}{s_0}=-{ℛ}_0(r(t)-r_0).$$ 因为 $r(t)-r_0\le 1$，初始易感比例为正时有 $s(t)\ge s_0{e}^{-{ℛ}_0}>0$。令时间趋于无穷，并用 $r_{\mathrm{\infty }}=1-s_{\mathrm{\infty }}$，得到终值方程 $$\ln \frac{s_{\mathrm{\infty }}}{s_0}=-{ℛ}_0(1-s_{\mathrm{\infty }}-r_0).$$ 本例就是 $s_{\mathrm{\infty }}=0.99e^{-3(1-s_{\mathrm{\infty }})}$。在 $0<s<1/3$ 内，函数 $\ln (s/0.99)+3(1-s)$ 的导数 $1/s-3$ 为正；两端符号不同，所以该区间恰有一个根。数值求根得 $$s_{\mathrm{\infty }}\approx 0.0587973648,r_{\mathrm{\infty }}\approx \mathrm{0.9412026352.}$$ 即最终约剩 58.80 名易感者，约 941.20 人曾进入感染并随后移出。这 941.20 包含初始的 10 名感染者；初始时刻之后的新感染累计数是 $990-58.80\approx 931.20$。

剩余易感者没有继续全部感染，是因为感染源在减少。当感染者变少时，易感者面对的接触风险随之下降。守恒不等于“所有人最终都会进入移出组”。如果一开始 $I(0)=0$，且没有外来感染输入，则方程始终给出 $I=0$，不会自行启动传播。

改一个机制参数，曲线怎样变

保持初值和 $\gamma =0.1$ 不变，把 $\beta$ 依次设为 0.08、0.15、0.3。它们的初始有效再生数分别为 0.792、1.485、2.97。第一个情形从一开始就下降；后两个情形先上升，但峰值高度与发生时间不同。

这些曲线只比较同一模型中的假设参数，β 的改变在现实中由什么因素造成，需要另作研究。

模型的无病状态是一整条平衡状态集合：$I=0,S+R=N$。沿感染方向的小扰动满足近似增长率 $\beta S/N-\gamma$；沿这条平衡线移动，却仍是另一个平衡点。因此不能把整个阈值讨论简写为“某一个无病点吸引所有邻近状态”。一般稳定性的语言见平衡点与稳定性。

如果传播存在明显的潜伏期，应增加尚未具有传染性的分组；如果接触网络差异很大，均匀混合假设需要调整；如果有再次易感、人口更新或持续输入，$R$ 单调增加与上述终值关系也可能改变。有限人数下，真实个体转移还存在随机性，特别是感染者很少时；这时可进一步用随机过程描述。数值方法可以把给定方程求得更精确，却不能替代对这些机制的判断。

参考资料

• Jeffrey R. Chasnov，Mathematical Biology，§4.2–4.3，pp.50–53，恒定移出率、SIR 方程与阈值；其人数字式使用未除以 N 的传播系数，须按正文换算。
• W. O. Kermack、A. G. McKendrick，A Contribution to the Mathematical Theory of Epidemics，Proceedings of the Royal Society A，115（1927），700–721，早期一般传播模型研究。