引子

对传染病传播规律的研究有着十分悠久的历史，一般认为最早该学科最早始创于1760年数学家丹尼尔·伯努利(Daniel Bernoulli)的一项对接种预防天花的研究，他最早基于概率论提出了第一个具有微分方程数学形式的数学模型，用来定量分析天花接种后的预期效果，论证了全民接种将显著提高人均寿命的观点，奠定了传染病动力学的学科基础。

第二个在该领域作出标志性贡献的人是英国医生罗纳德·罗斯(Ronald Ross)，在他的一项关于疟疾传播规律的研究中，提出了人-蚊双宿主模型，这是传染病动力学中的第一个媒介传播模型，他本人也因为揭露了疟疾传播的核心规律而斩获了1902年的诺贝尔医学奖。

随后，两位英国数学家威廉(William Ogilvy Kermack)和安德森(Anderson Gray McKendrick)在他们合署的一篇论文A Contribution to the Mathematical Theory of Epidemics中，开创性地提出了SIR模型(Susceptible-Infected-Recovered Model)，这是传染病动力学中一个十分经典的数学模型，迄今为止仍被广泛使用。

SIR模型

在现今的传染病动力学中，目前主要沿用的方法仍是SIR模型，该模型将总人群(N)分为易感者(S, Suspectible)、感染者(I, Infected)和康复者(R, Recovered)三类，研究传染病传播的数学规律。该模型的基本假设是规模不变的封闭群体(Closed Population)，即：

不考虑出生和死亡；
不考虑感染致死；
不考虑迁入和迁出。

即：

$N = S + I + R$ $\frac{dS}{dt} + \frac{dI}{dt} + \frac{dR}{dt} = 0$

并在此基础上假设存在一个恒定的传染速率$\beta$和康复速率$\gamma$，建立了如下微分方程组：

$\left\{\begin{align*} &\frac{dS}{dt} = -\beta SI \\ &\frac{dI}{dt} = \beta SI - \gamma I \\ &\frac{dR}{dt} = \gamma I \end{align*}\right.$

在上述假设下，他们最终证明：

存在一个可预测的流行阈值(即最大感染人数和发生的时间)；
并非所有人都会感染；
流行最终会结束。

如今传染病动力学的大多数模型都是它的改良和拓展。

例如在不考虑感染后获得免疫的假设下，模型简化为SIS模型，如普通感冒和淋病就属于这一类传染病；在免疫期有限的假设下，则发展为SIRS模型，新型冠状病毒(COVID-19)就是属于这一类传染病；加入感染后潜伏期则发展为SEIR模型，适用于像麻疹和埃博拉等具有这一特征的传染病；考虑感染者死亡则发展为SIRD模型等等。

此外，还可以将人群按照年龄结构、空间分布等特征作进一步细分，提出更复杂的年龄结构模型和空间迁移模型，适用于更复杂的传染病动力学研究。

仿真

下文仅以经典SIR模型和SIRS模型为例，提供一个简单可运行的Matlab仿真代码：

% SIR_simulation.m
clear; close all; clc;

%% 参数
N  = 1e4;
I0 = 10;
R0_init = 0;
S0 = N - I0 - R0_init;

b0    = 0.30;
gamma = 0.10;
tspan = [0 200];

R0_basic = b0 / gamma;
fprintf('基本传染数 R0 = %.2f\n', R0_basic);

y0 = [S0; I0; R0_init];

%% RHS
rhs = @(t, y) [
    - b0 * y(1) * y(2) / N ;                  % dS/dt
      b0 * y(1) * y(2) / N - gamma * y(2);    % dI/dt
      gamma * y(2)                            % dR/dt
];

opts = odeset('RelTol',1e-8,'AbsTol',1e-8);
[t, y] = ode45(rhs, tspan, y0, opts);

S = y(:,1); I = y(:,2); R = y(:,3);

%% 求解
[t,y] = ode45(rhs, tspan, y0);
S = y(:,1); I = y(:,2); R = y(:,3);

%% 峰值
[maxI, idxMax] = max(I);
tPeak = t(idxMax);
fprintf('峰值感染人数 = %.0f，发生在 t = %.2f 天（占 %.2f%%）。\n', maxI, tPeak, 100*maxI/N);

%% 绘图
figure('Units','normalized','Position',[0.1 0.1 0.7 0.45]);
plot(t, S, 'b-', 'LineWidth', 1.6); hold on;
plot(t, I, 'r-', 'LineWidth', 1.6);
plot(t, R, 'g-', 'LineWidth', 1.6);
xlabel('时间（天）'); ylabel('人数');
title(sprintf('SIRS 模拟（N=%d, \\beta=%.2f, \\gamma=%.2f）', N, b0, gamma));
legend('S','I','R','Location','Best');
grid on;
% 标注峰值
plot(tPeak, maxI, 'ko', 'MarkerFaceColor','y');
text(tPeak, maxI, sprintf('  感染人数峰值 (t=%.1f, I=%.0f)', tPeak, maxI));

%% S - I相位图
figure('Units','normalized','Position',[0.15 0.15 0.35 0.35]);
plot(S, I, 'LineWidth', 1.4);
xlabel('S'); ylabel('I');
title('相图');
grid on;

输出结果如下图所示：

如果引入一个免疫消退率$\omega$，表示康复者将在一定时间后重新成为易感者，则SIR模型改良为SIRS模型，其微分方程形式如下：

$\left\{\begin{align*} &\frac{dS}{dt} = -\beta SI + \omega R \\ &\frac{dI}{dt} = \beta SI - \gamma I \\ &\frac{dR}{dt} = \gamma I - \omega R \end{align*}\right.$

对应的Matlab仿真代码如下：

% SIRS_simulation.m
clear; close all; clc;

%% 参数
N       = 1e4;     % 总人口
I0      = 10;      % 初始感染者
R0_init = 0;   % 初始康复者
S0      = N - I0 - R0_init;

b0    = 0.30;    % 传播率
gamma = 0.10;    % 康复率
omega = 1 / 30;    % 免疫消退率，分母为免疫天数

tspan = [0 200]; % 模拟时间

R0_basic = b0 / gamma;
fprintf('参数： beta=%.3f, gamma=%.3f, omega=1/%.0f=%.5f\n', b0, gamma, 1 / omega, omega);
fprintf('基本传染数 R0 = %.2f\n', b0/gamma);

y0 = [S0; I0; R0_init];

%% RHS
rhs = @(t, y) [
    - b0 * y(1) * y(2) / N + omega * y(3);    % dS/dt
      b0 * y(1) * y(2) / N - gamma * y(2);    % dI/dt
      gamma * y(2) - omega * y(3)             % dR/dt
];

opts = odeset('RelTol',1e-8,'AbsTol',1e-8);
[t, y] = ode45(rhs, tspan, y0, opts);

S = y(:,1); I = y(:,2); R = y(:,3);

%% 峰值与长期行为
[maxI, idxMax] = max(I);
tPeak = t(idxMax);
fprintf('峰值感染人数 = %.0f，发生在 t = %.2f 天（占 %.2f%%）。\n', maxI, tPeak, 100*maxI/N);

finalI = I(end);
finalS = S(end);
finalR = R(end);
fprintf('模拟结束时（t=%.0f 天）： I=%.2f, S=%.2f, R=%.2f\n', t(end), finalI, finalS, finalR);

%% 绘图
figure('Units','normalized','Position',[0.1 0.1 0.7 0.45]);
plot(t, S, 'b-', 'LineWidth', 1.6); hold on;
plot(t, I, 'r-', 'LineWidth', 1.6);
plot(t, R, 'g-', 'LineWidth', 1.6);
xlabel('时间（天）'); ylabel('人数');
title(sprintf('SIRS 模拟（N=%d, \\beta=%.2f, \\gamma=%.2f, \\omega=1/%.0f）', N, b0, gamma, 1/omega));
legend('S','I','R','Location','Best');
grid on;
% 标注峰值
plot(tPeak, maxI, 'ko', 'MarkerFaceColor','y');
text(tPeak, maxI, sprintf('  感染人数峰值 (t=%.1f, I=%.0f)', tPeak, maxI));

%% S - I相图
figure('Units','normalized','Position',[0.15 0.15 0.35 0.35]);
plot(S, I, 'LineWidth', 1.4);
xlabel('S'); ylabel('I');
title('相图');
grid on;