菜鸟建模：新冠疫情+疫苗的数学建模(2)
该集对上一集的模型继续完善。上一集https://www.bilibili.com/read/cv9483873
传染病的预防措施分为传染源、传播途径和易感人群三个方面。所以预防传染病的措施有控制传染源、切断传播途径、保护易感人群。即在起点、路上、终点三个方面给它掐断线。
控制传染源的方法主要是隔离、治疗感染者，对接触者进行隔离、筛查、治疗；切断传播途径的方法主要是我们平时遵守的各项防疫要求，比如不得聚众、戴口罩之类的。需要注意的是这并不是“保护易感人群”，——所谓“保护易感人群”是指特异性免疫（比如接种疫苗）和非特异性免疫（比如锻炼身体提高免疫力）。
对于这三大方面的措施，我们对模型进行更细致的设计。
1、对于传染源管理，考虑疫区医疗资源承载量Nc，在未得到外界增援时为人口N的千分之一，如果医疗资源充足，只要确诊就会被立即隔离并接受治疗；若确诊患者溢出承载量，则溢出部分，二者的区别在于治愈率、病亡率。
2、对于传播途径，仍沿用接触率、感染率，不再引入新参数，接触率、感染率的大小，令其与政府防疫力度、对民众的有效宣传教育力度有关。
3、对于保护易感人群，引入新参数疫苗接种率V_rate，V指Vaccination rate.假设疫苗接种后都是有效的，那么疫苗接种率就是人群免疫率，疫苗接种对象为潜伏者或者易感者，只要没发病前接种疫苗我们都认为是有效的（匈奴南下前修好长城就行）。∴(E+S)*V_rate=V,
它与治愈者R一样都是传染病演化的终点。则用如下的状态转移图描述：
对程序进行相应修改：
N=1000000; %N为某一地区人口常量，100万人口的城市
E=0;%E为潜伏者数，初始时潜伏者数为0，E在潜伏期过后可能会变成确诊感染者
I=1;%I为感染者，初始时只出现首例感染者
R=0;%康复者数，假定康复后在考察时间段内具备免疫力
contact_I=20;%为感染者群体每人每日平均接触人数
contact_E=20;%为潜伏者群体每人每日平均接触人数
%每个感染者每日接触到的易感者数为contact_I*易感率=contact_I*S/N
%每个潜伏者每日接触到的易感者数为contact_E*易感率=contact_E*S/N
infect_rate_IS=0.05;%I、S群体间的感染率
infect_rate_ES=0.05;%E、S群体间的感染率
%∴感染者群体每天感染的易感者总数为I*(contact_I*S/N)*infect_rate
%∴潜伏者群体每天感染的易感者总数为E*(contact_E*S/N)*infect_rate
A=1/7;%A为潜伏者转化为感染者的概率(潜伏期的倒数)，潜伏期一般为1~14天，取7天
r1=0.1;%r为每日治愈率，每个患者每日有一定的几率被治好,r1为住院患者治愈率
r2=0.01%r2为非住院患者治愈率，基本指望自愈
d1=0.01;%Death_from_illness_rate为患者每日病亡率，d1为住院患者病亡率
d2=0.1;%d2为非住院患者病亡率
D=0;%D为病亡群体，由感染者群体I转化而来
V=0;%V为接种人群，假设疫苗有效率100%（事实上肯定小于100%）
S=N-E-I-R-D;%S为易感者数，为总人口数-潜伏者数-感染者数-治愈者数-病亡者数
T=1:100;%g考察时间段为T天
NC=N*0.001;%NC为疫区医疗资源承载量为该地区人口的千分之一
for i =1:length(T)-1
if I(i)<=NC
S(i+1)=S(i)- contact_I*infect_rate_IS*S(i)*I(i)/N -contact_E*infect_rate_IS*S(i)*E(i)/N;
E(i+1) = E(i) + contact_I*infect_rate_ES*S(i)*I(i)/N-A*E(i) + contact_E*infect_rate_ES*S(i)*E(i)/N;
I(i+1) = I(i) + A*E(i) - r1*I(i)-d1*I(i);
R(i+1) = R(i) + r1*I(i);
D(i+1)=D(i)+d1*I(i);
end
if(I(i)>NC)
S(i+1)=S(i)- contact_I*infect_rate_IS*S(i)*I(i)/N -contact_E*infect_rate_IS*S(i)*E(i)/N;
E(i+1) = E(i) + contact_I*infect_rate_ES*S(i)*I(i)/N-A*E(i) + contact_E*infect_rate_ES*S(i)*E(i)/N;
I(i+1) = I(i) + A*E(i) - r1*NC-d1*NC-r2*(I(i)-NC)-d2*(I(i)-NC);
R(i+1) = R(i) + r1*NC+r2*NC;
D(i+1)=D(i)+d1*I(i)+d2*(I(i)-NC);
end
if(i>10)
contact_I=1;
contact_E=1;
end
if(i>35&&i<39)%隔离一段时间，民众憋疯了，借机发泄，游行三天
contact_I=2;
contact_E=23;
infect_rate_ES=0.1;
end
if(i>=39)%游行结束，政府加大了戒严力度
contact_I=0.5;
contact_E=0.5;
infect_rate_ES=0.03;
infect_rate_IS=0.015;
d1=0.015;
r1=0.5;
end
end
figure(1)
plot(T,S,T,E,T,I,T,R,T,D);grid on;
xlabel('天'),ylabel('人数'),legend('易感者','潜伏者','传染者','康复者','病亡者'),title('made by 光电面壁人')
figure(2)
plot(T,E,T,I,T,R,T,D);grid on;
xlabel('天'),ylabel('人数'),legend('潜伏者','传染者','康复者','病亡者'),title('made by 光电面壁人')
这个千分之一的承载量是我自己感觉的一个经验数，不同级别的城市其医疗资源承载量有数量级上的差别。100万人的城市在我们国家算是个地级市，没有外界支援的情况下医疗资源全力运转承载千分之一人口数也就是1000人应该差不多。
我家乡15万人的小县城城镇，没有外界支援也就承载个位数的容量，撑死两位数到万分之一的数量级。当然这都是算的城市，对于乡下人口就不好使了。
因为有数量级差别，所以画下边这个图，
100万人的城市，按程序运行后约有1.5%的人病亡。这个结局仍然过于沉重，使结局更美好的措施是外界支援医疗力量、接种疫苗。美国许多城市医疗力量已经崩溃了，援军遥遥无期，能指望的也就是疫苗了。
假设这个游行为多次游行的一个缩影、或者说多次游行的合效果——使潜伏者大增从而摧毁防疫成果。也就好比现在的美国许多城市，挽回成果的唯一举措就是及时接种疫苗，假如说第40天接种疫苗，日接种率是10%S，即10万数量级的疫苗数，接种三天（该城市取得的支援有限），然后在一定时间后，接种人群接受后期疫苗。为简化处理，假设接种人群一次疫苗后就不会感染，相应地修改程序：
N=1000000; %N为某一地区人口常量，100万人口的城市
E=0;%E为潜伏者数，初始时潜伏者数为0，E在潜伏期过后可能会变成确诊感染者
I=1;%I为感染者，初始时只出现首例感染者
R=0;%康复者数，假定康复后在考察时间段内具备免疫力
contact_I=20;%为感染者群体每人每日平均接触人数
contact_E=20;%为潜伏者群体每人每日平均接触人数
%每个感染者每日接触到的易感者数为contact_I*易感率=contact_I*S/N
%每个潜伏者每日接触到的易感者数为contact_E*易感率=contact_E*S/N
infect_rate_IS=0.05;%I、S群体间的感染率
infect_rate_ES=0.05;%E、S群体间的感染率
%∴感染者群体每天感染的易感者总数为I*(contact_I*S/N)*infect_rate
%∴潜伏者群体每天感染的易感者总数为E*(contact_E*S/N)*infect_rate
A=1/7;%A为潜伏者转化为感染者的概率(潜伏期的倒数)，潜伏期一般为1~14天，取7天
r1=0.1;%r为每日治愈率，每个患者每日有一定的几率被治好,r1为住院患者治愈率
r2=0.01;%r2为非住院患者治愈率，基本指望自愈
d1=0.01;%Death_from_illness_rate为患者每日病亡率，d1为住院患者病亡率
d2=0.1;%d2为非住院患者病亡率
D=0;%D为病亡群体，由感染者群体I转化而来
V=0;%V为接种人群，假设疫苗有效率100%（事实上肯定小于100%）
S=N-E-I-R-D-V;%S为易感者数，为总人口数-潜伏者数-感染者数-治愈者数-病亡者数
T=1:100;%g考察时间段为T天
NC=N*0.001;%NC为疫区医疗资源承载量为该地区人口的千分之一
for i =1:length(T)-1
if I(i)<=NC
S(i+1)=S(i)- contact_I*infect_rate_IS*S(i)*I(i)/N -contact_E*infect_rate_IS*S(i)*E(i)/N;
E(i+1) = E(i) + contact_I*infect_rate_ES*S(i)*I(i)/N-A*E(i) + contact_E*infect_rate_ES*S(i)*E(i)/N;
I(i+1) = I(i) + A*E(i) - r1*I(i)-d1*I(i);
R(i+1) = R(i) + r1*I(i);
D(i+1)=D(i)+d1*I(i);
V(i+1)=V(i);
end
if I(i)>NC
S(i+1)=S(i)- contact_I*infect_rate_IS*S(i)*I(i)/N -contact_E*infect_rate_IS*S(i)*E(i)/N;
E(i+1) = E(i) + contact_I*infect_rate_ES*S(i)*I(i)/N-A*E(i) + contact_E*infect_rate_ES*S(i)*E(i)/N;
I(i+1) = I(i) + A*E(i) - r1*NC-d1*NC-r2*(I(i)-NC)-d2*(I(i)-NC);
R(i+1) = R(i) + r1*NC+r2*NC;
D(i+1)=D(i)+d1*I(i)+d2*(I(i)-NC);
V(i+1)=V(i);
end
if i>10
contact_I=1;
contact_E=1;
end
if i>35&&i<39%隔离一段时间，民众憋疯了，借机发泄，游行三天
contact_I=2;
contact_E=23;
infect_rate_ES=0.1;
end
if i>=39%游行结束，政府加大了各方面防疫力度
contact_I=0.5;
contact_E=0.5;
infect_rate_ES=0.03;
infect_rate_IS=0.015;
d1=0.015;
r1=0.5;
if i>=40&&i<=42
S(i+1)=S(i)- contact_I*infect_rate_IS*S(i)*I(i)/N -contact_E*infect_rate_IS*S(i)*E(i)/N-0.1*S(i);
E(i+1) = E(i) + contact_I*infect_rate_ES*S(i)*I(i)/N-A*E(i) + contact_E*infect_rate_ES*S(i)*E(i)/N-0.1*E(i);
I(i+1) = I(i) + A*E(i) - r1*NC-d1*NC-r2*(I(i)-NC)-d2*(I(i)-NC);
R(i+1) = R(i) + r1*NC+r2*NC;
D(i+1)=D(i)+d1*I(i)+d2*(I(i)-NC);
V(i+1)=V(i)+0.1*S(i)+0.1*E(i);
end
end
end
figure(1)
plot(T,S,T,E,T,I,T,R,T,D,T,V);grid on;
xlabel('天'),ylabel('人数'),legend('易感者','潜伏者','传染者','康复者','病亡者','接种者'),title('made by 光电面壁人')
figure(2)
plot(T,E,T,I,T,R,T,D,T,V);grid on;
xlabel('天'),ylabel('人数'),legend('潜伏者','传染者','康复者','病亡者','接种者'),title('made by 光电面壁人')
对比之下，接种疫苗有效地挽回了许多人的生命。
假设第40天时力求全员接种疫苗，疫苗供应充足，日接种率为(0.1+0.002*i)，日接种人数为约为S*(0.1+0.002*i)、E*(0.1+0.002*i)，则
后期疫情平稳了，政府和民众对继续全员疫苗的事不热心了，100万人口的城市在游行后火速集中接种疫苗，然后经政府防控，
约有占疫区0.49%的人病亡，即使游行后火速接种疫苗，虽说多挽回了疫区约0.9%的人员生命，但仍然代价高昂，可见，没事别游行，一游行就会造成不可挽回的巨大损失，疫情期间游行等于集体自残。图上我们可以看到接种疫苗后，疫情很快会被控制住，在接种疫苗过后不久即可以全面恢复经济，当然，要是像欧美某些政客没有疫苗就敢强行恢复经济，无疑是谋财害命。
本模型中所有参数均来自于作者凭生活经验主观臆断，并无严格的事实依据，实际上即使是有专门的统计部门，它们的数据很可能也十分片面，经过各方面汇总后的数据也许会全面一些，但无论如何，新冠疫情的各个参数都是史无前例、起伏不定的，传染病模型细致考虑的话十分复杂，如果你有意见欢迎与我交流。
