考虑家、工作场所、学校等地点的拓扑结构的个体随机模型

使用matlab实现考虑家、工作场所、学校等地点的拓扑结构的个体随机模型，模拟流行病传播过程

汤悦

日期

考虑家、工作场所、学校等地点的拓扑结构的个体随机模型

下载讲解视频

1 模型假设与符号定义

1.1 模型假设

1.2 符号定义

1.3 数据储存

  1. 个体状态（stage）
    数据类型是大小为总人口数 ( n ) 的向量，用于存储每个个体的健康状态。其中，0表示易感者，1表示感染者，2表示康复者。通过这个状态变量，模型能够跟踪每个个体的当前健康状态，并据此更新其后续的传播和康复情况。

  2. 个体感染持续时间（durationSinceEnteredTheStage）
    数据类型是大小为总人口数 ( n ) 的向量，用于记录每个感染者从进入感染状态（stage = 1）以来已经经历的时间，单位为小时。从而计算感染个体是否满足康复条件（即感染持续时间是否超过康复期 ( invGamma )。

  3. 易感者、感染者和康复者的数量（S_record、I_record、R_record）
    数据类型是大小为 ( m * 24 ) 的向量，即每小时记录一次，( m ) 是模拟的天数。其中，S_record(i)、I_record(i)、R_record(i) 分别记录第 ( i ) 个小时时的易感者、感染者、康复者数量。这些变量用于记录在模拟过程中每个小时内的人群状态变化，可以绘制出易感者、感染者和康复者随时间变化的曲线。

  4. 新感染和新康复数（new_infected_record、new_recovered_record）
    数据类型是大小为 ( m * 24 ) 的向量，即每小时记录一次，( m ) 是模拟的天数。new_infected_record、new_recovered_record 分别记录第 ( i ) 个小时中的新感染者、新康复者数量。这些变量帮助跟踪传染病的传播速率和康复速率，能够反映出每个时间步的动态变化。

  5. 活动场所（location）
    数据类型是大小为总人口数 ( n ) 的向量，用于存储每个个体在某一时刻所处的活动场所。其中，0表示在家庭中，1表示在办公室中，2表示在学校中，-1表示外出（不参与传播）。该变量与个体的健康状态 stage 结合使用，来计算传播概率，决定易感者是否可能被感染。

  6. 个体是否为学生（is_student）
    数据类型是大小为总人口数 ( n ) 的向量，用于标记每个个体是否为学生，值为1表示是学生，0表示不是学生。它用于决定每个个体在工作日内的活动场所，学生在学校，非学生在办公室。

  7. 家庭和办公室编号（family_ids、office_ids）
    数据类型是大小为总人口数 ( n ) 的向量，分别存储每个个体所属的家庭编号和办公室编号。确保每个个体有家庭，且确保每个非学生在办公室工作。

  8. 标记新感染和新康复者（newly_infected、newly_recovered）
    数据类型是大小为总人口数 ( n ) 的逻辑数组。newly_infected 用于标记每个时间步内哪些易感者（stage == 0）变成了感染者（stage == 1）；newly_recovered 用于标记每个时间步内哪些感染者（stage == 1）经过了康复期后转为康复者（stage == 2）。

  9. 家庭和办公室的成员数（family_sizes、office_sizes）
    数据类型为向量，大小分别为 n_homes 和 n_offices。family_sizes、office_sizes 分别记录每个家庭的成员数量。这些变量是根据 home_PMF 和 office_PMF 概率分布随机生成的，用于模拟家庭和办公室中个体的分布。

  10. 场所的传播概率（prob_home、prob_office、prob_school）
    数据类型为标量，分别用于计算家庭、办公室、学校班级中的传播概率。这些值根据场所的接触率和每个时间步中的感染者数量进行更新，传递给 transmission_prob 变量，后者会用来决定个体是否会被感染。

2 个体属性

2.1 属性定义

每个个体在模型中具有以下属性：

  1. 家庭分配
    每个个体在模型初始化时随机分配到5000个家庭，确保每个家庭至少有一个个体。家庭中的个体数遵循经验分布，分布规则（home_PMF）由用户指定。

  2. 学生与非学生
    个体随机分配为学生（30%）或非学生（70%）。

    • 学生：工作日的上午和下午在学校，早晚时间段在家。
    • 非学生：工作日的上午和下午在办公室，早晚时间段在家。

  3. 办公室分配
    非学生随机分配到1300个办公室，办公室中的个体数遵循经验分布（office_PMF），由用户指定。

  4. 班级分配
    学生随机分配到80个班级，每个班级大约50个学生。

2.2 人群状态初始化

在模型初始化时，个体被随机分配到易感者（S）、感染者（I）和康复者（R）三个状态：

3 个体行为

3.1 各场所之间的运动

个体根据时间段和场所规定进行流动：

1. 工作日：

2. 周末（周六和周日）：

所有个体在家中活动，但有一定概率选择外出（例如去商场、餐厅、公园等）。外出时，不涉及传播。

3.2 各场所内的接触传播规则

在每个场所内，接触发生时有一定概率（ q ）发生传播感染：

个体在每个场所内均匀分布，且接触行为随机，固定数量的接触机会在每个时间段内发生。

3.3 患病状态之间的转化

4 伪代码与流程图

4.1 伪代码

算法1：考虑家、工作场所、学校等地点拓扑结构的个体随机模型模拟

输入：个体数N，每个个体的自然史参数康复期invGamma，个体中的学生比例p_student，家庭、办公室、学校班级数量n_homes、n_offices、n_classes，家庭、办公室成员人数分布概率home_PMF 、office_PMF，家庭、办公室、学校班级中个体每小时的接触频率c_home、c_office、c_school，个体在各特殊时段外出的概率p_out、p_lunchout、p_weekendout

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

初始化第0天的所有个体的状态(是否是学生 x.isStudent, 自然史状态 x.stage, 自发展为感染者状态后至今的时长 x.durationSinceEnteredTheStage);
根据 R_0 解出单次接触感染概率 q;
保存当前处于各自然史状态的人数，保存当日新增的病例数;

while 迭代步数未达到最大模拟天数 M do
    确定当前时间段是工作日还是周末，并按时间段分配个体活动场所;
    
    更新所有个体的新感染状态转变(易感到感染);
    更新所有个体的新康复状态转变(感染到康复);
    
    保存当前处于各自然史状态的人数，保存当日新增的病例数;
    
    if 当日无新增病例 then
        停止迭代并退出循环;
    end
end

绘制模拟结果：各自然史状态的人数随模拟天数的变化曲线，新感染者和新康复者的数量随时间变化曲线

算法1-1：确定当前时间段属性，并按时间段分配个体活动场所

输入：当前迭代步的时间点

输出：当前时段每个个体的活动场所

for 当前时间段 do
    确定当前时间段为工作日还是周末;
    
    if 是周末 (is_weekend) then
        更新 x.location = home;
        对每个个体生成随机值小于周末外出概率，则标记为外出 x.location = out;
    else 是工作日, 按不同时间段分配 then
        if 是早晨 (0:00 - 8:00) then
            更新所有个体 x.location = home;
        end
        
        if 是上午 (8:00 - 12:00) then
            更新学生 x.location = school;
            更新非学生 x.location = office;
        end
        
        if 是中午 (12:00 - 14:00) then
            更新所有个体 x.location = home;
            一定概率外出 x.location = out;
        end
        
        if 是下午 (14:00 - 18:00) then
            更新学生 x.location = school;
            更新非学生 x.location = office;
        end
        
        if 是晚上 (18:00 - 24:00) then
            更新所有个体 x.location = home;
            一定概率外出 x.location = out;
        end
    end
end

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

输入：上一时刻所有个体的状态

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

for each 个体 x do
    随机生成个体 x 在 stepSize 这段时间内的密切接触者;
    
    根据不同场所的接触率计算家庭、办公室和学校的传播概率 p = 1 - (1 - q)^(接触人数 * 时间单位);
    
    对每个个体 x，计算该个体所在场所的传播概率;
    
    if Uniform(0, 1) < x.p AND x.stage = 易感者 (S) then
        更新 x.stage = I; x.durationSinceEnteredTheStage = 0;
    end
    
    记录新感染者数量，记录新感染更新后的所有个体的状态;
end

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

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

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

for each 个体 x do
    if x.stage == 1 且 x.durationSinceEnteredTheStage > invGamma then
        更新 x.stage = R;
    end
    
    对所有感染者，更新其感染时间 x.durationSinceEnteredTheStage = x.durationSinceEnteredTheStage + stepSize;
end

4.2 流程图

流程图

5 模拟设置

5.1 场景一

5.2 场景二

5.3 场景三

6 模拟结果与解释

6.1 场景一

场景一SIR 该图展示了模拟过程中，易感者（S）、感染者（I）、康复者（R）的数量随时间变化的动态趋势。随着时间推移，易感者的数量逐渐减少，一些易感者通过接触感染者而变为感染者；感染者的数量随时间波动，在大约第34天时达到流行高峰；随着感染者逐渐恢复，康复者数量也逐步增加。最后，易感者、感染者和康复者的数量都趋于稳定，模拟达到平衡状态。 场景一Rate 该图展示了每天（以小时为单位累计）新增感染者和新增康复者的数量。感染速率先上升后下降，在第27天时达到最大29人每小时；康复速率先上升后下降，在35天左右达到最大30人每小时。最终感染速率和康复速率都趋于0，模拟达到平衡状态。

6.2 场景二

场景二SIR 场景二中易感者（S）、感染者（I）、康复者（R）的数量随时间变化的动态趋势如图。随着时间推移，易感者的数量逐渐减少，一些易感者通过接触感染者而变为感染者；感染者的数量随时间波动，在大约第6天时达到流行高峰6650人；随着感染者逐渐恢复，康复者数量也逐步增加。最后，易感者、感染者和康复者的数量都趋于稳定，模拟达到平衡状态。 场景二Rate 场景二每天（以小时为单位累计）新增感染者和新增康复者的数量如图。感染速率先上升后下降，在第7天时达到最大131人每小时；康复速率先上升后下降，在9天左右达到最大131人每小时。最终感染速率和康复速率都趋于0，模拟达到平衡状态。

6.3 场景三

场景三SIR 场景三中易感者（S）、感染者（I）、康复者（R）的数量随时间变化的动态趋势如图。随着时间推移，易感者的数量逐渐减少，一些易感者通过接触感染者而变为感染者；感染者的数量随时间波动，在大约第57天时达到流行高峰52300人；随着感染者逐渐恢复，康复者数量也逐步增加。最后，易感者、感染者和康复者的数量都趋于稳定，模拟达到平衡状态。 场景三Rate 场景三每天（以小时为单位累计）新增感染者和新增康复者的数量如图。感染速率先上升后下降，在第45天时达到最大245人每小时；康复速率先上升后下降，在59天左右达到最大245人每小时。最终感染速率和康复速率都趋于0，模拟达到平衡状态。

6.4 参数分析与讨论

参数的设置反映了流行病传播模型中多个因素的复杂性。传染期决定了个体从感染到康复的时间，而时间步长决定了模拟的精细度。总人口数和学生比例则影响着疫情传播的规模和传播速度。接触频率则是决定传播速度的关键因素之一，尤其是在不同场所的群体活动密度较高时，传播风险显著增加。

通过调节这些参数，模型能够在不同的社会场景下，模拟出适合的疫情传播趋势，帮助研究人员和决策者制定有效的公共卫生策略。在实际的传染病防控中，了解这些参数如何相互作用，是制定精准防控策略的基础。

1.传染期（invGamma）

2.时间步长（Ts）

3.总人口数（n）

4.学生比例（p_student）

5.外出概率（p_out）

6.各地点接触频率（c_home、c_office、c_class）

7 源代码

模拟用时1.663秒

% 参数设置
m = 100; % 模拟天数
n = 1e4; % 总人口数
Ts = 1; % 每步时间单位：小时
c_home = 1; % 家庭接触率 (每小时接触单位)
c_office = 2; % 办公室接触率 (每小时接触单位)
c_school = 5; % 学校接触率 (每小时接触单位)
q = 0.01; % 传播率 (每次接触传播概率)
invGamma = 7 * 24; % 康复期 (小时)，7天转为小时

% 场所分配
p_student = 0.3; % 学生比例30%
p_non_student = 1 - p_student; % 非学生比例70%
n_students = round(n * p_student); % 学生人数
n_non_students = n - n_students; % 非学生人数
n_homes = 5000; % 家庭数量
n_offices = 1300; % 办公室数量
n_classes = 80; % 班级数量

home_PMF = [0.2, 0.3, 0.5]; % 1, 2, 3人的概率
office_PMF = [0.05, 0.05, 0.1, 0.15, 0.2, 0.2, 0.05, 0.1, 0.05, 0.05];  % 1, 2, 3, 4, 5, 6, 7, 8, 9, 10人的概率

% 1. 向量化生成家庭成员数
family_sizes = randsample(1:3, n_homes, true, home_PMF); % 为3000个家庭分配人数，依据home_PMF
family_members = repelem(1:n_homes, family_sizes); % 每个家庭成员编号

% 2. 向量化生成办公室成员数
office_sizes = randsample(1:10, n_offices, true, office_PMF); % 为500个办公室分配人数，依据office_PMF
office_members = repelem(1:n_offices, office_sizes); % 每个办公室成员编号

% 3. 确保总人口数匹配
if length(family_members) < n
    family_members = repmat(family_members, ceil(n / length(family_members)), 1); % 重复家庭成员编号，直到其长度大于等于n
end
family_members = family_members(1:n); % 截取至总人数

if length(office_members) < n
    office_members = repmat(office_members, ceil(n / length(office_members)), 1); % 重复办公室成员编号，直到其长度大于等于n
end
office_members = office_members(1:n); % 截取至总人数

% % 4. 为每个个体分配家庭编号和办公室编号
% individual_family_id = family_members; % 每个个体的家庭编号
% individual_office_id = office_members; % 每个个体的办公室编号

% % 显示分配结果
% disp('家庭分配样本：');
% disp(family_sizes(1:100)); % 显示前100个家庭的家庭人数
% disp('办公室分配样本：');
% disp(office_sizes(1:100)); % 显示前100个办公室的人数

% 初始化状态
stage = zeros(n, 1); % 0 - 易感者; 1 - 感染者; 2 - 康复者
durationSinceEnteredTheStage = zeros(n, 1); % 自进入感染阶段以来的时间

% 随机挑选10个初始感染者
initial_infected = randperm(n, 10); % 随机选择10个个体作为初始感染者
stage(initial_infected) = 1; % 将这10个个体的状态设为感染者

% 每个时间步的S, I, R数量记录
S_record = zeros(m * 24, 1); % 每小时记录
I_record = zeros(m * 24, 1);
R_record = zeros(m * 24, 1);
new_infected_record = zeros(m * 24, 1); % 新感染人数
new_recovered_record = zeros(m * 24, 1); % 新康复人数

% 每个个体的场所分配
location = zeros(n, 1); % 0 - 家庭; 1 - 办公室; 2 - 学校

% 随机为每个个体分配学生/非学生标志
is_student = rand(n, 1) < p_student; % 每个个体是否为学生（逻辑数组）

% 非学生分配到办公室，学生分配到学校
location(is_student) = 2; % 学生在学校
location(~is_student) = 1; % 非学生在办公室

% 模拟过程
for i = 1:m * 24  % 总小时数
    % 确定当前时间段
    current_day = mod(i - 1, 24 * 7); % 当前是第几天
    is_weekend = current_day == 0 || current_day == 120; % 周六或周日为周末

    % 每个时间段的个体行为
    hour_of_day = mod(i - 1, 24); % 当前是第几个小时（0 到 23）
    
    % 使用逻辑索引来分配活动场所
    location(:) = 0;  % 默认全体个体在家中

    if is_weekend
        % 周末：所有个体在家中活动，有50%概率选择外出
        p_weekendout = 0.5; % 外出概率
        weekendoutside_choice = rand(n, 1) < p_weekendout;
        location(weekendoutside_choice) = -1; % 外出的个体标记为 -1 (不涉及传播)
    else
        % 工作日：按时间段分配个体到各场所
        % 0:00 AM - 8:00 AM (早晨)
        is_morning = hour_of_day < 8; % 早晨的判断，返回逻辑数组，表示每个个体是否在早晨时段
        location(is_morning) = 0; % 所有个体在家中

        % 8:00 AM - 12:00 PM (上午)
        is_am = hour_of_day >= 8 & hour_of_day < 12;
        location(is_am & is_student) = 2; % 学生在学校
        location(is_am & ~is_student) = 1; % 非学生在办公室

        % 12:00 PM - 2:00 PM (中午)
        is_lunch = hour_of_day >= 12 & hour_of_day < 14;
        location(is_lunch) = 0; % 所有个体在家中
        p_lunchout = 0.05; % 外出概率
        lunch_choice = rand(n, 1) < p_lunchout;
        location(is_lunch & lunch_choice) = -1; % 外出的个体标记为 -1

        % 2:00 PM - 6:00 PM (下午)
        is_afternoon = hour_of_day >= 14 & hour_of_day < 18;
        location(is_afternoon & is_student) = 2; % 学生在学校
        location(is_afternoon & ~is_student) = 1; % 非学生在办公室

        % 6:00 PM - 12:00 AM (晚上)
        is_evening = hour_of_day >= 18 & hour_of_day < 24;
        location(is_evening) = 0; % 所有个体在家中
        p_out = 0.1; % 外出概率
        outside_choice = rand(n, 1) < p_out;
        location(is_evening & outside_choice) = -1; % 外出的个体标记为 -1
    end
    
    % 计算传播概率，根据每个时段的不同时间长度进行调整
    prob_home = 1 - (1 - q).^(c_home * Ts * sum(stage == 1 & location == 0) / n); % 家庭传播概率
    prob_office = 1 - (1 - q).^(c_office * Ts * sum(stage == 1 & location == 1) / n); % 办公室传播概率
    prob_school = 1 - (1 - q).^(c_school * Ts * sum(stage == 1 & location == 2) / n); % 学校传播概率

    % 各场所传播更新 (向量化)
    transmission_prob = zeros(n, 1);
    transmission_prob(location == 0) = prob_home; % 家庭传播概率
    transmission_prob(location == 1) = prob_office; % 办公室传播概率
    transmission_prob(location == 2) = prob_school; % 学校传播概率
    
    % 随机传播：易感者接触感染者的传播
    newly_infected = rand(n, 1) < transmission_prob & stage == 0; % 易感者转为感染者
    stage(newly_infected) = 1; % 更新为感染者
    durationSinceEnteredTheStage(newly_infected) = 0; % 对新感染者重置感染时长为0
    new_infected_record(i) = sum(newly_infected); % 记录新感染者数量
    
    % 状态转变更新
    % 更新康复状态
    newly_recovered = stage == 1 & durationSinceEnteredTheStage > invGamma; % 计算新康复者
    stage(newly_recovered) = 2; % 康复者状态
    new_recovered_record(i) = sum(newly_recovered); % 记录新康复者数量
    
    % 更新感染者的感染时长
    durationSinceEnteredTheStage(stage == 1) = durationSinceEnteredTheStage(stage == 1) + Ts; % 更新感染者的感染时长
    
    % 记录S, I, R的数量
    S_record(i) = sum(stage == 0); % 易感者
    I_record(i) = sum(stage == 1); % 感染者
    R_record(i) = sum(stage == 2); % 康复者
end

% 绘制图表
figure;
hold on;
% 横坐标改为以天为单位
time_in_days = (1:m*24) / 24; 

plot(time_in_days, S_record, 'g', 'LineWidth', 2); % 绘制易感者曲线，绿色
plot(time_in_days, I_record, 'r', 'LineWidth', 2); % 绘制感染者曲线，红色
plot(time_in_days, R_record, 'b', 'LineWidth', 2); % 绘制康复者曲线，蓝色
title('SIR模型 (易感者、感染者、康复者)');
xlabel('时间 (天)');
ylabel('人数');
legend('易感者 (S)', '感染者 (I)', '康复者 (R)');
grid on;

% 绘制新感染和新康复人数
figure;
hold on;
plot(time_in_days, new_infected_record, 'm', 'LineWidth', 1); % 新感染者，品红色
plot(time_in_days, new_recovered_record, 'c', 'LineWidth', 1); % 新康复者，青色
title('每日新感染和新康复人数');
xlabel('时间 (天)');
ylabel('人数');
legend('新感染者', '新康复者');
grid on;
汤悦
男神GG Bond全球后援会亚太地区分会长、四川邓紫棋