SEIAQRW模型

使用matlab实现SEIAQRW模型的拟合与干预措施效果评估

王静

日期

SEIAQRW模型

1 概述

1.1 模型假设

1.2 模型描述

传播途径包括人传人及水或食物传人的传染病，若在研究过程中需考虑隔离仓室，可以用SEIAQRW模型描述其传播过程，其中S代表易感个体（Susceptible），E代表潜伏期个体（Exposed），I代表病例/有症状感染者（Symptomatic），A代表无症状感染者（Asymptomatic），Q代表被隔离的病例，即隔离者（Quarantined），R代表恢复个体（Recovered），W代表受污染的水或食物（transmission medium）。SEIAQRW模型的基本原理是，通过一定的传播速率，易感者变为潜伏者，潜伏者按照一定的比例和发展速率分别转化为病例或无症状感染者，病例按照一定的隔离速率和隔离比例变为隔离者，无症状感染者、未被隔离的病例以及隔离者按照一定的恢复速率转变为恢复者，并一定的排泄速率向传播介质排放病毒，最后介质中的病毒含量还需考虑其消亡情况。消毒措施在模型中主要通过降低人传人传染率来影响传染病传播过程，因为病毒类传染病一般极低剂量即可致病，因此认为只有消毒措施达到100%的时候，水或食物传人的传染率才为0。

1.3 模型仓室及参数含义

仓室 含义
易感者(S) 尚未感染但易感的个体
潜伏期者(E) 病原体侵入机体至最早出现临床症状的这段时间内被感染的个体
病例(I) 经过潜伏期后出现症状并能够传播疾病的个体
无症状感染者(A) 经过潜伏期后没有出现症状并能够传播疾病的个体
隔离者(Q) 因为诺如腹泻中，只有病例才会表现出症状被隔离，因此这里指被隔离的病例
恢复者(R) 已经从疾病中恢复并获得免疫力的个体
传播介质(W) 被污染的，能够传播疾病的水或食物
参数 含义
β 人传人传染率，控制感染者向易感者传播病毒的速度
βw 水或食物传人传染率，控制受污染的水或食物向易感者传播病毒的速度
k 无症状感染者的病毒传播能力与病例的比值
p 无症状感染者占所有的感染者的比例
ω 病例的潜伏期
ω' 无症状感染者的潜隐期
q 隔离比例
γ 病例、隔离者的恢复率
γ' 无症状感染者的恢复率
γ'' 病例出现后被送去隔离的速率
μ 病例向外环境中排泄病毒的速率
μ' 无症状感染者向传播介质中排泄病毒的速率
ε 传播介质中病毒的消亡速率

2 模型流程图及方程

![png](SEIAQRW model.png)

SEIAQRW模型的核心方程为： $$ \begin{cases} \frac{dS}{dt} = -\beta S (I + kA) - \beta_w S W \\
\frac{dE}{dt} = \beta S (I + kA) + \beta_w S W - p \omega' E - (1 - p) \omega E \\
\frac{dI}{dt} = (1 - p) \omega E - q \gamma'' I - (1 - q) \gamma I \\
\frac{dA}{dt} = p \omega' E - \gamma' A \\
\frac{dQ}{dt} = q \gamma'' I - \gamma Q \\
\frac{dR}{dt} = (1 - q) \gamma I + \gamma' A + \gamma Q\\
\frac{dW}{dt} = \mu I + \mu' A - \epsilon W \end{cases} $$

3 适用疾病流行条件

1.存在易感人群：即研究人群中需要有足够数量的易感人群。 2.存在有效的传播途径：这个模型涉及的传播途径包括人传人途径和水或食物传播途径。 3.存在潜伏期。 4.存在有效的传染率：传染率一方面决定了易感者（S）被感染者（I）或隐性感染者（A）传染的概率，另一方面也决定了易感者在接触受污染的传播介质后被传染的风险。只有当传染率足够高时，疾病才能在人群中迅速传播。 5.在研究期间，不需要考虑人群免疫，感染者（I）或隐性感染者（A）经过恢复期后会转变为恢复者，不会出现死亡。 6.对于该种传染病的防控措施中，消毒和隔离措施是比较常见的防控手段。

4 代码（matlab）

4.1 条件设定

4.1.1 模拟场景及拟解决问题

某疾控接到报告，某小学多名学生出现呕吐、腹泻症状，疾控工作人员迅速赶往现场进行调查，通过症状检查、询问人员接触史和饮食情况及现场采样检测后，确定这是一起水或食物传人引起的诺如腹泻暴发性疫情，传染源为食堂的一道鱼，工作人员在开展调查的同时，已经与学校领导合作开展防控工作，如要求食堂将与鱼一起采购并放在一起的食材进行有效处置、对病例进行隔离治疗、开展环境消杀以及健康教育。在写总结报告时，某工作人员想对该起疫情的诺如病毒传播特性有更深的了解，并想探究他认为最有效的两种防控措施在疫情传播过程中的作用，因此他构建了该SEIAQRW模型（易感-潜伏期-病例-无症状感染-隔离-恢复-传播介质模型），帮助理解疫情的动态演化过程。首先是对某诺如暴发的流行曲线进行拟合，获取相关参数的取值；然后将消毒措施与消毒+隔离措施合理的纳入到模型中，模拟不同干预强度下的防控效果，并通过图表直观展示模型结果。

4.1.2 参数在模型中的具体设定

​ s: 当前易感者人数 ​ e: 当前潜伏者人数 ​ i: 当前病例人数 ​ a: 当前无症状感染者人数 ​ g：当前隔离者人数 ​ r: 当前恢复者人数 ​ w: 当前传播介质中病毒含量 ​ md.b : 人传人传染率 (感染者每人每日传播的概率) ​ md.bW：传播介质（受污染的水或食物）传人的传染率 ​ md.kappa : 无症状感染者的病毒传播能力与病例的比值 ​ md.p : 无症状感染者占所有的感染者的比例 ​ md.q : 病例被送去隔离的比例 ​ md.omegap : 无症状感染者的潜隐期 ​ md.omega : 病例的潜伏期 ​ md.gamma：病例恢复率 ​ md.gammap : 无症状感染者的恢复率 ​ md.gammapp : 隔离者的被隔离的速率 ​ md.c : 无症状感染者的病毒排泄能力与病例的比值 ​ md.epsilon : 病例向传播介质中排泄病毒速率，另外现有文献资料也将这个参数的取值假设为病毒消亡速率 ​ dt: 时间步长，默认值为1

4.1.3 初始条件设定

    - 几个仓室的比例读取某诺如暴发的流行曲线数据资料来设置
    - `k` 设置为0.05。
    - `p` 设置为0.3。
    - `ω` 设置为1。
    - `ω'` 设置为1。
    - `γ` 设置为1/3。
    - `γ'` 设置为0.03846。
    - `μ` 设置为0.1。
    - `q'` 设置为0.3,0.6,0.9。
    - `γ''` 设置为2，即认为病例会在半天内被送去隔离，1/（1/2）。
    -`ε` 设置为0.1。
    - 其余参数通过模型拟合得到。

4.2 代码

4.2.1 函数及模型定义部分代码

% 定义extractOutbreakData函数
function datai = extractOutbreakData(data, num)

idx1 = find(contains(data.Number, num2str(num)));
idx1 = idx1(1);
idx2 = find(contains(data.Number, num2str(num+1)));
idx2 = idx2(1) - 1;

datai = data(idx1:idx2, 2:3);
datai(isnat(datai.Date), :) = [];
end

% 定义interpolateData函数
function [xRefined, yRefined] = interpolateData(x, y, stepSize, method)
xRefined = min(x) : stepSize : max(x);
yRefined = interp1(x, cumsum(y), xRefined, method);
yRefined = [y(1), diff(yRefined) / stepSize];
end

% 定义SEAIQRW模型的函数
classdef dynamicalModelSEIARQW
  % 定义SEAIQRW模型中会出现的各个参数，并设置为非负的
  properties
        N double {mustBeNonnegative, mustBeFinite}
        %beta double {mustBeNonnegative, mustBeFinite}
        kappa double {mustBeNonnegative, mustBeFinite}
        p double { mustBeNonnegative,mustBeFinite}
        omega double {mustBeNonnegative, mustBeFinite}
        omegap double {mustBeNonnegative, mustBeFinite}
        gamma double {mustBeNonnegative, mustBeFinite}
        gammap double {mustBeNonnegative, mustBeFinite}
        gammapp double {mustBeNonnegative, mustBeFinite}
        q double {mustBeNonnegative, mustBeFinite}
        mu double {mustBeNonnegative, mustBeFinite}
        epsilon double {mustBeNonnegative, mustBeFinite}
        b
        bW
        c
        i0
        e0
        w0
    % 设置用来模拟存储时间和状态变量的参数
        tSpan 
        xInit double {mustBeNonnegative, mustBeFinite}
        xFinal 
        t
        x
    end
    % 定义了一个名为derivatives的函数，它计算模型在时间t和状态变量x下的导数
   methods
        function dxdt = derivatives(md, t, x)
    % 状态向量x中提取各个状态变量的值。
            s = x(1);
            e = x(2);
            i = x(3);
            a = x(4);
            r = x(5);
            g = x(6);
            w = x(7);
     %模型微分方程                                                         
            dsdt = -md.b * s * (i + md.kappa * a) - md.bW * s * w;
            dedt = -dsdt - (md.p * md.omegap + (1-md.p) * md.omega) * e;
            didt = (1-md.p) * md.omega * e - (md.q * md.gammapp + (1-md.q) * md.gamma) * i;
            dadt = md.p * md.omegap * e - md.gammap * a;
            drdt = (1-md.q) * md.gamma * i + md.gamma * g + md.gammap * a;
            dgdt = md.q * md.gammapp * i - md.gamma * g;
            dwdt = (i + a * md.c) * md.epsilon - md.epsilon * w;

    % 将所有状态变量的导数组合成一个向量dxdt
            dxdt = [dsdt; dedt; didt; dadt; drdt; dgdt; dwdt];
        end
     % 这几行定义了一个名为predict的函数，它使用ode23函数来预测模型在给定时间跨度md.tSpan和初始条件md.xInit下的情况。
        function [t, x] = predict(md)
            odefun = @(t,x) derivatives(md, t, x);
            [t, x] = ode23(odefun, md.tSpan, md.xInit(:));
        end
     % 定义了一个名为objectivefunction的函数，模型预测与实际数据之间的误差。
        function loss = objectivefunction(parametersToBeFitted, md, tData, iData)
            md.b = parametersToBeFitted(1); % parameter b
            md.bW = parametersToBeFitted(2); % parameter bW 
            md.c = parametersToBeFitted(3); % parameter c
            md.xInit(2) = parametersToBeFitted(4); % variable e
            md.xInit(3) = parametersToBeFitted(5); % variable i
            md.xInit(7) = parametersToBeFitted(6); % variable w

            [t, x] = predict(md);
            dailyIncidentRate = (1-md.p) .* md.omega .* x(:,2);
            predictedCummulateI = cumtrapz(t, dailyIncidentRate);
            predictedCummulateI = interp1(t, predictedCummulateI, tData);
            observedCummulateI = cumtrapz(tData, iData);% 这几行计算了模型的预测累计感染人数，并使其接近实际累计感染人数。

            loss = norm(predictedCummulateI - observedCummulateI) ^ 2;
        end
       
     % 模型拟合，获取参数。
       function [fittedModel, t, x] = fit(md, tData, iData, parameterNames, initialGuessOfParameters, lowerBoundOfParameters, upperBoundOfParameters)
            % 需要拟合的参数为：'b', 'bW', 'c', 'e0', 'i0', 'w0';
            if ~isempty(md.b),  b = md.b;  else b = 0;  end
            if ~isempty(md.bW),  bW = md.bW;    else bW = 0;    end
            if ~isempty(md.c),  c = md.c;  else c = 0;  end
            if ~isempty(md.e0),  e0 = md.e0;    else e0 = md.xInit(2);  end
            if ~isempty(md.i0),  i0 = md.i0;    else i0 = md.xInit(3);  end
            if ~isempty(md.w0),  w0 = md.w0;    else w0 = 0;  end
            defaultParameterValues = [b, bW, c, e0, i0, w0];
            md.tSpan = [min(tData), max(tData)];
            temp = string({'b', 'bW', 'c', 'e0', 'i0', 'w0'}) == string(parameterNames(:));
            temp = mat2cell(temp, ones(1,numel(parameterNames)), 6);
            idx = cellfun(@(x)find(x), temp);               
            defaultParameterValues(idx) = initialGuessOfParameters;
            S.type = '()';
            S.subs = {idx};
            objfun = @(parametersToBeFitted) objectivefunction(subsasgn(defaultParameterValues, S, parametersToBeFitted), md, tData, iData);
            initialTryOfParameters = defaultParameterValues(idx);
            A = [];
            b = [];
            Aeq = [];
            beq = [];
            lb = lowerBoundOfParameters;
            ub = upperBoundOfParameters;
            fittedParameters = fmincon(objfun, initialTryOfParameters, A, b, Aeq, beq, lb, ub);
            fittedParameters = subsasgn(defaultParameterValues, S, fittedParameters);
            
          %  创建一个新的模型fittedModel，并使用拟合后的参数值更新其属性；
            fittedModel = md;
            fittedModel.b = fittedParameters(1);
            fittedModel.bW = fittedParameters(2);
            fittedModel.c = fittedParameters(3);
            fittedModel.e0 = fittedParameters(4); % fitted initial e(0)
            fittedModel.i0 = fittedParameters(5); % fitted initial i(0)
            fittedModel.w0 = fittedParameters(6); % fitted initial w(0)
            fittedModel.xInit(2) = fittedModel.e0;
            fittedModel.xInit(3) = fittedModel.i0;
            fittedModel.xInit(7) = fittedModel.w0;

            [t, x] = predict(fittedModel);
            fittedModel.t = t;
            fittedModel.x = x;
            fittedModel.xFinal = x(end,:);
        end

        function [modelList, t, x] = piecewiseFit(md, tData, iData, breakPoints, parameterNames, initialGuessOfParameters, lowerBoundOfParameters, upperBoundOfParameters)
            tStart = min(tData);
            tEnd = max(tData);
            breakPoints(breakPoints <= tStart | breakPoints >= tEnd) = [];
            breakPoints = [tStart; breakPoints; tEnd];% 找到实际暴发资料中发病率最高值对应的时间点，即断点。
            segmentCount = numel(breakPoints) - 1;% 计算时间段的数量（断点数量减一）。
            modelList = repmat(md, [segmentCount, 1]);
            observedCummulation = iData(1);
            predictedCummulation = iData(1);
         % 创建一个模型对象的矩阵，用于存储每个时间段的模型，并初始化空的时间向量和状态变量矩阵。
            t = [];
            x = [];
            for i = 1:segmentCount
                % data for this segment
                idx = tData >= breakPoints(i) & tData <= breakPoints(i+1);
                tDatai = tData(idx);
                iDatai = iData(idx);
                % init values for this segment
                model0 = md;
                model0.tSpan = [breakPoints(i), breakPoints(i+1)];
                % 找到断点后，根据其是否存在2个时间段开展模型拟合。
                if i > 1
                    model0.xInit = modelList(i-1).xFinal(1:7);
                    model0.e0 = modelList(i-1).xFinal(2);
                    model0.i0 = modelList(i-1).xFinal(3);
                    model0.w0 = modelList(i-1).xFinal(7);
                end
                if i > 1
                    idx = contains(parameterNames, {'e0', 'i0', 'w0'});
                    parameterNames(idx) = [];
                    initialGuessOfParameters(idx) = [];
                end
                [modeli, ti, xi] = fit(model0, tDatai, iDatai, parameterNames, initialGuessOfParameters, lowerBoundOfParameters, upperBoundOfParameters);
               modelList(i) = modeli;
                if i == 1
                    t = [t; ti];
                    x = [x; xi];
                else
                    t = [t; ti(2:end)];
                    x = [x; xi(2:end, :)];
                end
            end
        end
        
% 计算流行病结束的时间
        function d = last(md)
            % compute the time instance of the end of epidemic (the date that daily new incidence lesser than 1)
            maxDuration = 3e3;
            md.tSpan = [md.tSpan(1), maxDuration];
            [t,x] = predict(md);

            obj = (1-md.p) * md.omega * x(:, 2); % dailyNewIncidence
            d = find(obj > 1 / md.N);

            if numel(d) == 0
                d = inf;
                return;
            end
            d = d(end);
            d = t(d);

            if isempty(d) || d == 0 || d == maxDuration 
                d = inf;
            end
        end  
    end % 
end %

4.2.2 模型模拟及拟合部分代码

close all; clear; clc;
%读取流行曲线数据资料
outbreakNumber = 1;
opts1 = detectImportOptions("allData.xlsx", 'Sheet', '流行曲线');
opts2 = detectImportOptions("allData.xlsx", 'Sheet', '疫情基本信息');
data = readtable('allData.xlsx', opts1, 'Sheet', '流行曲线');
dataInfo = readtable('allData.xlsx', opts2, 'Sheet', '疫情基本信息');

# 初始变量设置
model = dynamicalModelSEIARQW;
model.N = dataInfo.PopulationSize(dataInfo.ID == outbreakNumber);
datai = extractOutbreakData(data, outbreakNumber);
datai(datai.Incidence == 0, :) = [];
iData = datai.Incidence / model.N;
dData = datai.Date;
tData = days(dData - min(dData));
[tRefined, iRefined] = deal(tData, iData);

# 初始参数设置
variableNames = {'s', 'e', 'i', 'a', 'r', 'q', 'w'};
model.kappa = 0.05;
model.omega = 1;
model.omegap = 1;
model.p = 0.3;
model.gamma = 1 / 3;
model.gammap = 0.03846;
model.gammapp = 2;
model.epsilon = 0.1;
model.q = 0;
model.e0 = iData(1) / ((1-model.p) * model.omega);
model.i0 = iData(1);
model.w0 = 0;
model.xInit = [1, model.e0, model.i0, 0, 0, model.w0];

parametersToBeFitted = {'b', 'bW', 'w0', 'c'};
initialGuessOfParameters = [1e-1, 1e-1, 1e-1, 0.3];
lowerBoundOfParameters = [0, 0, 1e-3, 1e-3];
upperBoundOfParameters = [1, 1, 1, 1];  

# 将实际流行曲线中第一次出现发病高峰值的时间设置为断点，将其前面的数据设为用来模型拟合的数据
    breakPoints = tData(iData == max(iData));
    breakPoints = breakPoints(1);
    [modelList, t, x] = piecewiseFit(model, tRefined, iRefined, breakPoints, parametersToBeFitted, initialGuessOfParameters, lowerBoundOfParameters, upperBoundOfParameters);
    dailyNewIncidence = (1-model.p) * model.omega * x(:,2); 
    [model1, t1, x1] = fit(model, tRefined(tRefined <= breakPoints(1)), iRefined(tRefined <= breakPoints(1)), parametersToBeFitted, initialGuessOfParameters,           lowerBoundOfParameters, upperBoundOfParameters);
dailyNewIncidence1 = (1-model.p) *  model1.omega * x1(:,2);

# 设置干预措施进行模型模拟，以及进行可视化
% 设置两种干预措施，消毒、消毒+隔离
interventionNames = {'Sterilization',  'Quarantine+Sterilization'};
for k = 1:numel(interventionNames)
    % 先将实际数据用蓝色圆圈进行绘图，然后是用来拟合的这部分数据用红色实线绘图
    ax = nexttile;
    hold on;
    plot(tData, iData, 'bo');
    plot(t1, dailyNewIncidence1, 'r-');
    % 设置干预措施实施强度，为0.3,0.6,0.9
    intensity = 0:0.3:0.9;
    currentCase = 0; 
     lastRecord = zeros(numel(intensity)+1,1);
  for i = 1:numel(intensity)
        model2 = model1;
        model2.xInit = model1.x(end, 1:7);
        % 因为消毒措施牵扯到通过改变解除后单次接触感染概率来改变人传人传染率，所以我们要根据拟合得到的人传人传染率的大小，判断使用哪个公式计算单次接触感染概率系数，再将消毒措施的强度系数与概率系数结合后计算消毒后的人传人接触感染概率系数。因为只有消毒强度达到100%后，水或食物传人传染率系数才会改变，变为0，但因为我们这里消毒强度未考虑100%，所以这里在模拟过程中不需要改变水或食物传人传染率系数
        if model2.b >= 1
            e = model2.b / 15;
            switch k
                case 1
                    % Sterilization
                    y = 1 - intensity(i);
                    model2.b = y * e * 15;
                     currentCase = 1;
                case 2
                % quarantine+Sterilization
                    model2.q = intensity(i);
                     y = 1 - intensity(i);
                     model2.b = y * e * 15;
                      currentCase = 2;
            end
        else
                 e = 1 - (1-model2.b) ^ (1/15);
            switch k
                case 1
                    % Sterilization
                    y = 1 - intensity(i);
                    model2.b = 1 - (1 - y * e) ^ 15;
                     currentCase = 1;
                case 2
                % quarantine+Sterilization
                    model2.q = intensity(i);
                     y = 1 - intensity(i);
                     model2.b = 1 - (1 - y * e) ^ 15;
                      currentCase = 2;
            end
       end

model2.tSpan = [breakPoints(1), last(model2)];
lastRecord(i)=model2.tSpan(2);
        [t2, x2] = predict(model2);


        dailyNewIncidence2 = (1-model.p) * model2.omega * x2(:,2);
       
%这部分均为可视化
        if i == 1
            plot(t2, dailyNewIncidence2, 'k-.');  % 无干预模拟结果，使用黑色虚线
                elseif i ~= 1
                    if  currentCase  == 1
                        plot(t2, dailyNewIncidence2, 'r-.');  % 消毒措施模拟结果，使用红色虚线
                    else
                        plot(t2, dailyNewIncidence2, 'g-.');  % 消毒+隔离措施模拟结果，使用绿色虚线
                end
        end
  end
% 设置横纵坐标名称、标签、图的题目、X轴和Y轴的最大值最小值
    xlabel(ax, 'Time (in days)');
    ylabel(ax, 'Daily New Incidence');
    legend1 = {'Observations', 'Fitted Incidence', 'Predication Without Intervention'};
    legend2 = "Predication With Intervention Intensity " + intensity(2:end) * 100 + "%";
    legend([legend1, legend2], 'Location', 'northeast');
     title(interventionNames{k});
    ax.XLim(1) = 0;
    ax.XLim(2) = 80;
    ax.YLim(1) = 0;
    ax.YLim(2) = 0.06;

end

4.3 结果

png

结果分析

王静
无效内卷反抗者