SIR简单个体随机模型

1 模型假设与符号定义

1.1 模型假设

在SIR简单个体随机模型中，我们有以下假设：

1. 个体行为是随机的：每个个体在不同的健康状态下（易感、感染、康复）会有不同的行为模式，且行为是随机的，受外部环境和接触概率的影响。
2. 状态转化遵循一定规则：个体的状态仅能从易感（S）转为感染（I），从感染（I）转为康复（R），不存在从康复回到易感的情况。转化的速率由传播概率和康复概率决定。
3. 传播速率由接触频率与接触感染概率决定：我们将经典SIR模型中的传播速率β = (S ⋅ I) / N 拆分为接触频率c和单次接触感染概率q，即 β = c ⋅ q。
4. 个体之间的接触传播：个体间的接触传播是随机发生的，且每次接触的感染概率由个体的健康状态和接触频率共同决定。

1.2 符号定义

SIR简单个体随机模型的符号定义如下：

• S：易感者数量（Susceptible）
• I：感染者数量（Infected）
• R：康复者数量（Recovered）
• N：总人口数（N = S + I + R）
• β：传播速率（β = c × q），表示单位时间内每个感染者接触到易感者并成功传播的概率
• c：接触频率，单位时间内个体与他人接触的次数
• q：每次接触发生感染的概率。通过基本再生数 R₀ = (c ⋅ q) / γ 反推每次接触的感染概率q
• γ：康复率，单位时间内从感染状态转为康复状态的比例

1.3 数据储存

1. 个体状态（stage）：每个个体的健康状态，通常用整数值表示。数据类型存储为整数数组。其中：

   • 0 代表易感者
   • 1 代表感染者
   • 2 代表康复者

2. 感染持续时间（durationSinceEnteredTheStage）：表示每个感染者自进入感染状态以来的时间（天数），用于判断是否康复。数据类型存储为数值数组（n个个体，长度为n）。对于感染者来说，随着时间的推移，durationSinceEnteredTheStage 增加，直到达到设定的康复天数（invGamma）为止，个体从感染状态转为康复状态。

3. 汇总数据：用来记录每个时间步的群体整体状态。包括：

   • 易感者数量（S_record）
   • 感染者数量（I_record）
   • 康复者数量（R_record）
   • 实时再生数（Rt_record）

   数据类型均存储为数值数组，长度为模拟的时间步长m。

2 个体属性

2.1 属性定义

1. 自然史状态：每个个体在任意时刻都有一个健康状态，可能是易感、感染或康复。

• 易感者（S）：尚未感染，处于暴露风险中的个体。
• 感染者（I）：正在感染并可传播疾病的个体。
• 康复者（R）：已感染并且康复的个体，假设康复后不再传染。

2. 停留时长：每个个体在不同状态下停留的时间是随机的。我们通常假设状态转化时间服从指数分布，即每个状态的平均停留时间由转化率决定。

2.2 人群状态初始化

• 初始状态：模拟开始时，个体处于易感状态 $\ S_0$ ，感染者 $\ I_0 $ 和康复者 $\ R_0 $ 的数量可以根据初始条件设定。通常情况下，初始感染者较少，且康复者数量为零。
• 时间步长设置：在模型中，时间被离散化为多个小时间步长 $\Delta t $，每个时间步长更新个体的状态。

3 个体行为

3.1 接触传播

接触传播是该模型的核心机制。每个易感者个体与感染者的接触是随机的，而每次接触是否会导致感染取决于接触频率 ( c ) 和单次接触感染概率 ( q )。

• 传播速率：传播速率 $ \beta = c \cdot q $，即接触频率与感染概率的乘积。
• 传播事件：当易感个体 ( S ) 与感染个体 ( I ) 接触时，有概率 $ \beta \cdot \Delta t $ 使得易感个体转为感染者。

3.2 患病状态之间的转化

个体的健康状态随着时间的推移发生转化：

• 易感到感染：这一过程的转化速率通常由感染者数量（ I ）和易感者数量（ S ）共同决定。转化的概率为 $ \lambda = \beta \cdot \frac{I}{N} $，其中 $ \beta $ 是传染率，$ \frac{I}{N} $ 是感染者在总体中的比例。则易感者转为感染者的概率为：

  $$P(S \to I) = \beta \cdot \frac{I}{N} \cdot \Delta t$$

• 感染到康复：感染者转为康复者的速率由康复率 γ 决定。每个时间步长内，感染者有概率 $ \gamma \cdot \Delta t $ 转为康复者（ R ）。这里的 γ 通常被定义为单位时间内的康复率，即感染者的康复概率。则感染者转为康复者的概率为：

  $$P(I \to R) = \gamma \cdot \Delta t$$

4 伪代码与流程图

4.1 伪代码

算法1：均匀混合假设下的简单SIR个体随机模型的模拟

输入：个体数 ( N )，每个个体的自然史参数（康复率 γ），接触频率 ( c )，基本再生数 $\ R_0 $

输出：每日新增感染人数，每日现存的处于各自然史状态的个体数

初始化第0天的所有个体的状态(自然史状态x.stage，每个个体进入感染状态后的持续时间x.durationSinceEnteredTheStage);
根据R_0解出个体单次接触感染概率q;
保存当前处于各自然史状态的人数，保存当日新增的病例数;

while 迭代步数未达到最大模拟天数M do
    更新所有个体的新感染状态转变(易感到感染);
    更新所有个体的新康复状态转变(感染到康复);
    保存当前处于各自然史状态的人数，保存当日新增的病例数;
    
    if 当日无新增病例 then
        停止迭代并退出循环
    end
    
    记录当前时间步的易感者、感染者和康复者数量，计算实时再生数 Rt;
end

绘制易感者、感染者和康复者数量随时间的变化曲线，绘制实时再生数 Rt 变化曲线。

算法1-1：更新所有个体的新感染状态转变（易感到感染）

输入：上一时刻所有个体的状态 输出：新感染更新后的所有个体的状态

For each 个体 x do
    随机生成个体 x 在 stepSize 这段时间内的密切接触者;
    计算个体 x 被感染的概率 Probability = 1 - (1 - q) ^ ((感染者总数 / 总人口数) * (c * Ts));
    
    if 随机数 < Probability 且 stage[i] == 0（易感者） then
        更新 x.stage = I;
        更新所有感染者的感染持续时间 x.durationSinceEnteredTheStage = x.durationSinceEnteredTheStage + stepSize;
    end
end

算法1-2：更新所有个体的新康复状态转变（感染到康复）

输入：新感染更新后所有个体的状态 输出：新康复更新后所有个体的状态

For each 个体 x do
    if x.stage == 1 且 x.durationSinceEnteredTheStage > invGamma then
        x 发展为康复者 R，更新 x.stage = R;
    end
end

4.2 流程图

5 模拟设置

• 总人口数 N：模型中总人口数的大小。假设总人口数 ( N = 10000 )，模拟社区。
• 初始健康状态：( $S_0 = N - I_0 - R_0$ ) 是初始时刻的易感者数量，假设初始感染者较少( $I_0 = 1$ )，初始康复者数量为零( $R_0 = 0$ )。
• 传播速率 β：传播速率由接触频率 ( c ) 和单次接触感染概率 ( q ) 共同决定，设定 ( c = 10 )，( q = 0.05 )。
• 康复率 γ：康复率表示单位时间内从感染状态转为康复状态的概率，通常通过疾病的自然恢复时间来设定，假设 ( $\gamma = \frac{1}{7}$ )，表示平均7天康复一个感染者。
• 时间步长 Δt：时间步长用于离散化模拟过程，通常选择1天或1小时等。这里我们设置为1天。
• 模拟时长 m：模拟的时间总长度，设定为60天( m = 60 )。

6 模拟结果与解释

6.1 模拟结果

模拟结果展示了易感者、感染者和康复者在每个时间步长的数量变化，和有效再生数随时间变化的曲线。 随着感染者( I )数量的增加，易感者( S )数量逐渐减少；感染者( I )数量随着疾病的传播逐渐增加，在第24天达到峰值7000后逐渐下降，最终趋于零；康复者( R )的数量逐渐增加，表明感染者( I )逐渐康复并退出传播链。 模拟结果中，基本再生数 ( $R_0 = 3.5$ )，说明每个感染者在完全易感人群中平均能够将疾病传播给3.5个其他个体。随着感染者的增加，易感者逐渐减少，( $R_t$ )逐渐下降至1以下并趋于0，表明疫情正在渐渐平息。

6.2 参数分析与讨论

• 总人数 ( N )：是指模型中所有个体的总数，包括易感者、感染者和康复者。它通常用于表示目标人群的规模。在大多数流行病模拟中，总人数决定了疾病传播的基础条件，如人群的密度、接触网络等。大规模人群流行病传播更为缓慢，峰值更高，疫情持续时间较长，群体免疫的到达需要更多的感染；小规模人群传播速度更快，峰值较高但衰退也较快。流行病通常较早达到顶峰并迅速结束。
• 传染期（invGamma）：是指个体从感染状态转为康复状态所需的时间，通常以天为单位，传染期的倒数通常称为康复率( $\gamma$ )。如果传染期较长，意味着感染者在传染阶段停留的时间较长，这将增加每个感染者传播病毒的机会，从而可能导致更高的感染人数；较短的传染期则意味着感染者快速康复，减少了传播的机会，可能会导致较低的感染人数，尤其在疾病传播速度较快时，整个流行病会较早达到峰值并迅速消退。
• 接触频率（c）：指每个个体每天与他人接触的平均次数。在此模型中，接触频率决定了易感者和感染者之间发生接触并传播的概率。较高的接触频率（如人群密集的城市环境）将大大增加病毒传播的机会。此时，尽管其他参数保持不变，病毒的传播速度会显著加快，感染者数量的增长也会更为迅速；在接触频率较低的情况下（如社交距离较远或密闭环境中），疾病的传播会变慢，整个流行病的峰值也会降低。 总结：总人数、接触频率、传染期等参数在模拟结果中密切相关。增加总人数 ( N ) 和接触频率 ( c )，会增加流行病传播的机会，而较长的传染期 invGamma 则会延长流行病的传播过程。

7 源代码

模拟用时3.058秒

m = 60; % days for simulation
n = 1e4; % person
Ts = 1; % day 
c = 10; % person per day 
q = 0.05; % R0=(c*q)/Gamma
invGamma = 7; % days

% 初始化状态
stage = zeros(n,1); % 0 - susceptible; 1 - infectious; 2 - removed
durationSinceEnteredTheStage = zeros(n,1); % time duration since entered the stage of infectious
stage(1) = 1; % [0, 1, 2] represent [S, I, R] respectively

% 用于记录每个时间步的S, I, R数量以及Reff
S_record = zeros(m,1); 
I_record = zeros(m,1);
R_record = zeros(m,1);
Rt_record = zeros(m,1); 

% 模拟更新
for i = 1:m
    % 传染更新
    Probability = 1 - (1-q)^((sum(stage) / n) * (c*Ts));
    stage(rand(n,1) < Probability & stage == 0) = 1;

    % 状态转变更新
    durationSinceEnteredTheStage(stage == 1) = durationSinceEnteredTheStage(stage == 1) + Ts;
    stage(stage == 1 & durationSinceEnteredTheStage > invGamma) = 2;

    % 记录S, I, R的数量
    S_record(i) = sum(stage == 0); % 易感者
    I_record(i) = sum(stage == 1); % 感染者
    R_record(i) = sum(stage == 2); % 康复者

    % 计算实时 Rt
    S = S_record(i); % 当前易感者数量
    N = n; % 总人口数量
    R0 = (c*q)/(1/invGamma); % 计算 R0
    Rt_record(i) = R0 * (S / N); % 计算实时 Rt
end

% 绘制S, I, R数量的变化
figure;
hold on;
plot(1:m, S_record, 'g', 'LineWidth', 2); % 绘制易感者曲线，绿色
plot(1:m, I_record, 'r', 'LineWidth', 2); % 绘制感染者曲线，红色
plot(1:m, R_record, 'b', 'LineWidth', 2); % 绘制康复者曲线，蓝色
% 添加标题和标签
title('SIR简单个体随机模型');
xlabel('时间 (天)');
ylabel('人数');
legend('易感者 (S)', '感染者 (I)', '康复者 (R)');
grid on;

% 绘制Rt的变化
figure;
plot(1:m, Rt_record, 'k', 'LineWidth', 2); % 绘制Rt曲线
xlabel('时间 (天)');
ylabel('Rt');
title('实时再生数Rt');
grid on;