包含年龄异质性、疫苗效果的个体随机模型

使用matlab实现包含年龄异质性、疫苗效果的个体随机模型，模拟流行病传播过程

汤悦

日期

包含年龄异质性、疫苗效果的个体随机模型

1 模型假设与符号定义

下载讲解视频

1.1 模型假设

1.2 符号定义

2 个体属性

2.1 属性定义

2.2 人群状态初始化

3 个体行为

3.1 接触传播

个体之间的接触频率由年龄组接触矩阵决定，每个年龄组的个体在每天接触的随机个体数不同。接触矩阵定义了每个年龄组与其他年龄组的接触概率。接触传播的过程中，易感者可能因为接触到感染者而被感染。

3.2 患病状态之间的转化

3.3 接种疫苗

疫苗的效果在防感染和防传播方面有所不同，且会随时间衰退。每个个体的疫苗效果都会在每个时间步长后均匀衰退。

4 伪代码与流程图

4.1 伪代码

算法1：包含年龄异质性、疫苗效果的个体随机模型模拟

输入：个体数 ( N )，每个年龄组个体的自然史参数（康复速率 ( Gamma )），年龄组日均接触频率 ( c )，各年龄组接触偏好矩阵 ( contact_preference_matrix )，疫苗覆盖率 ( vaccine_coverage )，疫苗防感染率 ( vaccine\efficacy\infection )，防传播率 ( vaccine_efficacy_transmission )，衰退率 ( vaccine_decay_rate )
输出：每日新增感染人数、新增康复人数，每日现存的处于各自然史状态的个体数，分年龄组每日现存的处于各自然史状态的个体数，分疫苗组每日现存的处于各自然史状态的个体数

初始化第0天的所有个体的状态(是否接种疫苗 x.isVaccinated，自然史状态 x.stage，年龄组 x.ageGroup，自发展为感染者状态后至今的时长 x.durationSinceEnteredTheStage)；
根据 R_0 解出单次接触感染概率 q；
保存当前处于各自然史状态的人数，保存当日新增的病例数、新增康复人数；
while 迭代步数未达到最大模拟天数 M do
    查找当前易感者 (S) 和感染者 (I) 的索引
    if 没有易感者或感染者 then
        结束模拟；
    end
    基于年龄组接触能力和接触偏好，为易感者创建接触矩阵；
    对每个易感者计算感染概率，更新所有个体的新感染状态转变（易感到感染）；
    更新所有个体的新康复状态转变（感染到康复）；
    保存当前处于各自然史状态的人数，保存当日新增的病例数、新增康复人数；
    更新疫苗效果衰退；
end
绘制模拟结果，各自然史状态的人数随模拟天数的变化曲线，新感染者和新康复者的数量随时间变化曲线

算法1-1：基于年龄组接触能力和接触偏好，为易感者创建接触矩阵

输入：当前迭代步的易感者索引 输出：当前时间每个个体的接触对象

for 当前迭代步的每个易感者 do
    根据易感者的年龄组获取接触能力 c_vector；
    根据接触偏好矩阵生成接触个体的年龄分布 contact_probs；
    for 每个年龄组的易感者 do
        根据偏好矩阵，选择接触的个体的年龄组 contact_age_group；
        随机选择一个年龄组内的个体作为接触者 chosen_contact；
        记录这个易感者接触的每一个个体，生成接触矩阵 contacts_matrix；
    end
end

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

输入：上一时刻所有个体的状态，接触矩阵 contacts_matrix 输出：新感染更新后的所有个体的状态

for each 易感者 do
    读取该易感者实际接触的个体数；
    计算该易感者接触到的感染者数量；
    进一步计算接触到的感染者中接种疫苗的人数 vaccinated_infected_contacts 和未接种疫苗的人数 unvaccinated_infected_contacts；
    根据个体接种疫苗时间和疫苗衰退率，计算个体的剩余疫苗效应；
    根据产生接触的易感者和感染者人数及其各自的剩余疫苗效应，计算传播概率 p；
    if Uniform(0,1) < 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将会发展为康复者，更新 x.stage = R；
        对所有感染者，更新其感染时间 x.durationSinceEnteredTheStage = x.durationSinceEnteredTheStage + stepSize；
    end
end

4.2 流程图

流程图

5 模拟设置

6 模拟结果与解释

6.1 模型结果

全人群SIR 通过模拟结果，绘制出易感者（S）、感染者（I）和康复者（R）数量随时间变化的曲线。随着时间的推移，易感者逐渐转变为感染者，并且最终成为康复者。这反映了一个考虑年龄异质性和疫苗接种的典型SIR模型过程：人群经过一段时间的流行病传播后，易感者数量减少，感染者和康复者数量交替上升。最终，系统达到了一个稳定的状态，其中大部分个体已经恢复或免疫。 新感染新康复 该图显示了在模型中，每一时刻的新增感染者（感染率）和新增康复者（恢复率）数量。最初阶段，感染率大幅上升，因为易感者大量转变为感染者。随着感染者逐渐康复，恢复率逐渐增大。这表明疫情的传播过程和人群的康复过程是同步进行的，在特定时段内，恢复者的数量与感染者数量有较强的相关性。 年龄组SIR 在不同年龄组之间，由于生理和免疫差异，易感者、感染者和康复者的数量呈现出明显不同的趋势。儿童、成人和老年人在接触病毒、感染病毒以及康复的过程中，存在不同的速率和特点。通过比较分析不同年龄组的SIR曲线，可以深入了解年龄对传染病传播过程的影响。 疫苗组SIR 接种疫苗的群体在疫情中的表现明显优于未接种疫苗的群体。接种疫苗后，易感者的数量减少，感染率的峰值也较低，康复率的增长较为平稳。这表明疫苗能够显著减缓病毒的传播并提高个体的免疫力，从而减少感染者和增加康复者的比例。未接种疫苗的群体则因为没有免疫保护，感染者数量较多，康复的速度较慢，最终在较长时间后才逐步恢复。

6.2 参数分析与讨论

参数对模拟结果的影响非常显著。通过调整这些参数的值，可以模拟不同的流行病场景，理解如何通过疫苗接种、年龄组差异、恢复速率等因素来控制疾病的传播。在实际应用中，需要根据不同地区、不同人群的特点来调整相关参数，以获得更精确的模拟结果。

1. 总人口数（n）

当人口规模较小时，个体之间的接触和传播的随机性会更加明显，可能会导致模型的波动较大，尤其在早期阶段（初始感染者较少时）。随着人口规模增大，模型结果趋向稳定。流行病传播的趋势会更加平滑，整体感染曲线会更符合大样本下的真实传播趋势。

人口规模越大，计算量也会增加，特别是在复杂模型或多个变量的情况下，模拟时间会显著增加。

2. 不同年龄组恢复速率（Gamma）

如果恢复速率较高，感染者在较短时间内康复，感染者的数量不会长时间持续，疫情的爆发期较短，传播的时间窗较小。整体感染人数的峰值较低，疫情会较早达到高峰并迅速下降。

恢复慢的情况下，感染者停留在传染状态的时间较长，导致疾病传播速度较快且持久，感染人数的峰值较高，疫情波动较大。

3. 不同年龄组接触能力（c）

如果成人的接触能力较强，感染传播的速率较快，尤其是成人群体之间的接触密切可能导致快速的传播，整体感染人数将增加。

如果群体接触能力较低，传播较为缓慢，尤其是在儿童和老年人群体中，感染者数量的增加速度较慢，疫情的传播更为缓和。

4. 不同年龄组之间的接触偏好概率（contact_preference_matrix）

通常为一个nxn的矩阵，每个元素表示不同年龄组之间接触的概率分布。

如果某些年龄组之间的接触概率较高（如儿童与成人接触概率较高），则这些群体之间的传播会较为迅速，尤其是感染者的传播可能会通过这些频繁接触的群体扩散开来。

如果各年龄组之间的接触概率相对均匀，传播路径较为复杂，可能会导致不同群体的感染率接近，疫情会呈现较为平缓的增长趋势。

5. 疫苗覆盖率（vaccine_coverage）

取值范围为0到1之间，0表示完全没有疫苗接种，1表示所有人都接种疫苗。

如果疫苗接种率低，易感人群较大，疫情爆发的风险较高，感染人数会迅速增加，且难以控制传播。

疫苗覆盖率较高时，大部分群体拥有免疫力，流行病传播会受到抑制，感染人数的峰值显著降低，传播速度也会减缓。

6. 疫苗防感染效果（vaccine_efficacy_infection）

取值范围在0到1之间，0表示疫苗完全无效，1表示疫苗完全有效。

如果疫苗防感染效果较低（如0.2或更低），即使接种了疫苗，个体仍然容易感染。疫情传播速度会受到较小的抑制，感染人数可能仍然较高。

如果疫苗防感染效果较高（如0.8或更高），接种者感染的几率大大降低，疫情传播的总体速度较慢，最终感染人数显著减少。

7. 疫苗防传播效果（vaccine_efficacy_transmission）

取值范围在0到1之间，0表示疫苗完全不能抑制传播，1表示疫苗完全能够阻止传播。

如果疫苗的防传播效果较低，接种者仍然可以传播病毒，尽管感染概率可能较低，但传播路径依然存在，疫情难以控制，整体传播速度较快。

如果防传播效果较高，接种者不仅不易感染，还能够减少感染后传播给他人的概率，整体疫情传播会显著减缓。

8. 疫苗衰退率（vaccine_decay_rate）

取值范围在0到1之间，值越大，表示疫苗效果衰退得越快。

疫苗效果迅速衰退，意味着即使个体在初期受到疫苗保护，随着时间推移，保护效果大幅下降，疫情的爆发可能会延迟，但之后会迅速恶化。

疫苗的效果保持较长时间，能够更长时间

7 源代码

模拟用时75.148秒


m = 80; % 模拟的天数
n = 1e4; % 人口总数
Ts = 1; % 时间步长（天）
q = 0.05; % 单次接触感染概率

% 不同年龄组恢复速率：儿童0.12，成人0.2，老人0.1（人/天）
Gamma = [0.12, 0.2, 0.1];

% 不同年龄组的接触能力（每个年龄组的易感者每天接触的随机个体数）
c = [10, 15, 12]; % 例如，儿童：10，成人：15，老人：12

% 定义接触矩阵（3x3矩阵）,表示不同年龄组个体接触不同年龄组的偏好的概率分布
% 行代表接触者的年龄组，列代表易感者的年龄组（条件概率矩阵）
contact_preference_matrix = [
    0.3, 0.5, 0.2;  % 儿童接触不同年龄组的概率分布
    0.4, 0.4, 0.2;  % 成人接触不同年龄组的概率分布
    0.3, 0.1, 0.6   % 老人接触不同年龄组的概率分布
];

% 疫苗效果参数
vaccine_coverage = 0.3; % 疫苗覆盖率30%
vaccine_efficacy_infection = 0.6; % 防感染率 VEi=60%
vaccine_efficacy_transmission = 0.5; % 防传播率 VEs=50%
vaccine_decay_rate = 0.001; % 疫苗效果每日衰退0.1%

% 初始化个体的状态
stage = zeros(n,1); % 0 - 易感者, 1 - 感染者, 2 - 康复者
age_group = zeros(n,1); % 初始化年龄组：0: 儿童, 1: 成人, 2: 老人
vaccination = rand(n,1) < vaccine_coverage; % 随机分配疫苗接种状态（30%已接种）
vaccination_time = zeros(n, 1);  % 记录每个个体的接种时间
vaccination_time(vaccination == 1) = rand(sum(vaccination == 1), 1) * 300;  % 随机生成模拟开始前疫苗已接种的时间，范围在0至300天之间
remaining_vaccine_efficacy_infection = max(vaccine_efficacy_infection - vaccine_decay_rate * vaccination_time, 0);
remaining_vaccine_efficacy_transmission = max(vaccine_efficacy_transmission - vaccine_decay_rate * vaccination_time, 0);

% 分配年龄组（20%儿童，50%成人，30%老人）
age_group(1:round(0.2*n)) = 0; % 儿童
age_group(round(0.2*n)+1:round(0.7*n)) = 1; % 成人
age_group(round(0.7*n)+1:end) = 2; % 老人

% 记录每个个体自进入感染阶段以来的持续时间
durationSinceEnteredTheStage = zeros(n,1); 

% 初始化感染者
% 随机选择初始感染者，保持人群年龄和疫苗接种分布不变
initial_infected = 5; % 初始感染者数
infected_indices = randperm(n, initial_infected); % 随机选择5个个体作为感染者，不改变全人群的年龄分布和疫苗接种状态
stage(infected_indices) = 1; % 将这5个个体设置为感染者


% 用于记录每个时间步的易感者(S)、感染者(I)、康复者(R)的数量
S_record = zeros(m,1); 
I_record = zeros(m,1);
R_record = zeros(m,1);

% 初始化记录感染率和恢复率的数组
infectionRates = zeros(m, 1);
recoverRates = zeros(m, 1);

% 用于记录每个年龄组的 S、I、R 数量
S_child = zeros(m,1); I_child = zeros(m,1); R_child = zeros(m,1);
S_adult = zeros(m,1); I_adult = zeros(m,1); R_adult = zeros(m,1);
S_elderly = zeros(m,1); I_elderly = zeros(m,1); R_elderly = zeros(m,1);

% 用于记录接种和未接种疫苗的 S、I、R 数量
S_unvaccinated = zeros(m,1); I_unvaccinated = zeros(m,1); R_unvaccinated = zeros(m,1);
S_vaccinated = zeros(m,1); I_vaccinated = zeros(m,1); R_vaccinated = zeros(m,1);

for i = 1:m
    % 找出当前易感者和感染者的索引
    susceptible_indices = find(stage == 0); 
    infectious_indices = find(stage == 1); 
    
    % 如果没有易感者或感染者，终止模拟
    if isempty(susceptible_indices) || isempty(infectious_indices)
        break;
    end
    
   % 创建接触矩阵：每个易感者基于其年龄组的接触能力和偏好
    c_vector = c(age_group(susceptible_indices) + 1);  % 根据易感者的年龄组获取接触能力
    contacts_matrix = zeros(length(susceptible_indices), max(c_vector));  % 初始化接触矩阵
    
    for j = 1:length(susceptible_indices)
        susceptible_person = susceptible_indices(j);
        age = age_group(susceptible_person);
        
        % 根据接触偏好矩阵生成接触个体的年龄分布
        contact_probs = contact_preference_matrix(age + 1, :);  % 获取易感者年龄组的接触概率分布
        
        % 根据概率生成每个易感者的接触个体
        for k = 1:c_vector(j)
            % 根据偏好矩阵，选择接触的个体年龄组
            chosen_age_group = find(rand <= cumsum(contact_probs), 1, 'first') - 1;  % 找到对应的接触年龄组
            contact_age_group = chosen_age_group;  % 记录接触者的年龄组
    
            % 随机选择一个年龄组内的个体作为接触者
            contact_indices = find(age_group == contact_age_group);
            chosen_contact = contact_indices(randi(length(contact_indices)));  % 随机选择接触者
    
            % 记录这个易感者接触的个体
            contacts_matrix(j, k) = chosen_contact;  % 将接触者的索引记录到接触矩阵
        end
    end

    % 记录本次感染和恢复的变化
    new_infected = 0;
    new_recovered = 0;

    % 计算感染概率：考虑接触到的感染者数量、接触率、疫苗效应等
    for j = 1:length(susceptible_indices)
        susceptible_person = susceptible_indices(j);
        age = age_group(susceptible_person);
        
        % 获取接触个体
        contacts = contacts_matrix(j, 1:c_vector(j));  % 只取该易感者实际接触的个体数
        
        % 计算接触到的感染者数量
        infected_contacts = contacts(stage(contacts) == 1);  % 找到所有感染者的接触者
        num_infected_contacts = length(infected_contacts);  % 密切接触者中感染者的数量
        
        % 计算接种疫苗的感染者（01向量，1表示接种，0表示未接种）
        vaccinated_infected_contacts = vaccination(infected_contacts) ==1; 
        % 计算未接种疫苗的感染者（01向量，0表示接种，1表示未接种）
        unvaccinated_infected_contacts = vaccination(infected_contacts) ==0; 

        % 计算个体的剩余疫苗效应
        individual_vaccine_effect_infection = max(vaccine_efficacy_infection - vaccine_decay_rate * vaccination_time(susceptible_person), 0);
        individual_vaccine_effect_transmission = max(vaccine_efficacy_transmission - vaccine_decay_rate * vaccination_time(susceptible_person), 0);
    
        % 计算感染概率：考虑易感者疫苗接种情况（VEi）和感染者疫苗接种情况（VEs）
        infection_probability = 1 - (1 - (1 - individual_vaccine_effect_transmission) * q) .^ ...
            sum(vaccinated_infected_contacts .* (1 - individual_vaccine_effect_infection) + ...
                unvaccinated_infected_contacts .* 1);

        % 如果随机数小于感染概率，则该易感者被感染
        if rand < infection_probability
            stage(susceptible_person) = 1; % 将该易感者感染
            new_infected = new_infected + 1; % 记录感染人数
        end 
    end
    
    % 更新感染者的持续时间
    durationSinceEnteredTheStage(stage == 1) = durationSinceEnteredTheStage(stage == 1) + Ts;
    
    % 计算恢复率：每个感染者是否康复
    for j = 1:n
        if stage(j) == 1 && durationSinceEnteredTheStage(j) > (1 / Gamma(age_group(j) + 1)) % 加1解决索引问题
            stage(j) = 2; % 转为康复者
            new_recovered = new_recovered + 1; % 记录康复人数
        end
    end


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

    % 记录感染率和恢复率
    infectionRates(i) = new_infected;  % 当前迭代中感染人数
    recoverRates(i) = new_recovered;   % 当前迭代中康复人数


    % 分别记录各个年龄组的易感者、感染者、康复者数量
    S_child(i) = sum(stage == 0 & age_group == 0);
    I_child(i) = sum(stage == 1 & age_group == 0);
    R_child(i) = sum(stage == 2 & age_group == 0);

    S_adult(i) = sum(stage == 0 & age_group == 1);
    I_adult(i) = sum(stage == 1 & age_group == 1);
    R_adult(i) = sum(stage == 2 & age_group == 1);

    S_elderly(i) = sum(stage == 0 & age_group == 2);
    I_elderly(i) = sum(stage == 1 & age_group == 2);
    R_elderly(i) = sum(stage == 2 & age_group == 2);

    % 分别记录接种和未接种疫苗的 S、I、R 数量
    S_unvaccinated(i) = sum(stage == 0 & vaccination == 0); 
    I_unvaccinated(i) = sum(stage == 1 & vaccination == 0);
    R_unvaccinated(i) = sum(stage == 2 & vaccination == 0);
   
    S_vaccinated(i) = sum(stage == 0 & vaccination == 1); 
    I_vaccinated(i) = sum(stage == 1 & vaccination == 1);
    R_vaccinated(i) = sum(stage == 2 & vaccination == 1);

    % 每天疫苗效果衰退
    remaining_vaccine_efficacy_infection = max(remaining_vaccine_efficacy_infection - vaccine_decay_rate, 0);
    remaining_vaccine_efficacy_transmission = max(remaining_vaccine_efficacy_transmission - vaccine_decay_rate, 0);
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;

% 绘制感染率和恢复率图
figure;
hold on;
plot(1:m, infectionRates, 'm-', 'LineWidth', 2);
plot(1:m, recoverRates, 'c-', 'LineWidth', 2);
title('模型感染率、恢复率');
xlabel('时间（天）');
ylabel('人数');
legend('新增感染者', '新增康复者');
grid on;

%% 绘制儿童的 SIR
figure;
subplot(1,3,1);
hold on;
plot(1:m, S_child, 'g', 'LineWidth', 2); % 绘制儿童的易感者数量
plot(1:m, I_child, 'r', 'LineWidth', 2); % 绘制儿童的感染者数量
plot(1:m, R_child, 'b', 'LineWidth', 2); % 绘制儿童的康复者数量
title('儿童的 SIR');
xlabel('时间（天）');
ylabel('人数');
legend('易感者(S)', '感染者(I)', '康复者(R)');
grid on;

% 绘制成人的 SIR
subplot(1,3,2);
hold on;
plot(1:m, S_adult, 'g', 'LineWidth', 2); % 绘制成人的易感者数量
plot(1:m, I_adult, 'r', 'LineWidth', 2); % 绘制成人的感染者数量
plot(1:m, R_adult, 'b', 'LineWidth', 2); % 绘制成人的康复者数量
title('成人的 SIR');
xlabel('时间（天）');
ylabel('人数');
legend('易感者(S)', '感染者(I)', '康复者(R)');
grid on;

% 绘制老人的 SIR
subplot(1,3,3);
hold on;
plot(1:m, S_elderly, 'g', 'LineWidth', 2); % 绘制老人的易感者数量
plot(1:m, I_elderly, 'r', 'LineWidth', 2); % 绘制老人的感染者数量
plot(1:m, R_elderly, 'b', 'LineWidth', 2); % 绘制老人的康复者数量
title('老人的 SIR');
xlabel('时间（天）');
ylabel('人数');
legend('易感者(S)', '感染者(I)', '康复者(R)');
grid on;

%% 绘制未接种疫苗的 SIR
figure;
subplot(1,2,1);
hold on;
plot(1:m, S_unvaccinated, 'g', 'LineWidth', 2); % 绘制未接种的易感者数量
plot(1:m, I_unvaccinated, 'r', 'LineWidth', 2); % 绘制未接种的感染者数量
plot(1:m, R_unvaccinated, 'b', 'LineWidth', 2); % 绘制未接种的康复者数量
title('未接种疫苗的 SIR');
xlabel('时间（天）');
ylabel('人数');
legend('易感者(S)', '感染者(I)', '康复者(R)');
grid on;

% 绘制接种疫苗的 SIR
subplot(1,2,2);
hold on;
plot(1:m, S_vaccinated, 'g', 'LineWidth', 2); % 绘制接种的易感者数量
plot(1:m, I_vaccinated, 'r', 'LineWidth', 2); % 绘制接种的感染者数量
plot(1:m, R_vaccinated, 'b', 'LineWidth', 2); % 绘制接种的康复者数量
title('接种疫苗的 SIR');
xlabel('时间（天）');
ylabel('人数');
legend('易感者(S)', '感染者(I)', '康复者(R)');
grid on;
汤悦
男神GG Bond全球后援会亚太地区分会长、四川邓紫棋