英国新冠与变种病毒的传播及传染力融入数学建模的探讨（网页版）

本文已参与「新人创作礼」活动，一起开启掘金创作之路。

前言

首先我想先感谢我的老师以及我的队友，没有他们的辛勤付出，这次数学建模比赛论文估计也写不出来😂🙏另外因为刚刚学习两个月，再加上比赛从开始到结束只有三天时间，所以时间实在是非常紧，所以论文写的也很仓促，即使比赛后又修改了一些，肯定还是不可避免的有些错误，还请各位大佬指正🙏因为是网页版，为方便网友测试程序也会直接放在正文里，可能不是特别严谨，切勿过于较真。

问题

新冠病毒（SARS-CoV-2）被检测出来后，在传播的过程中，发生了包括Alpha、Delta和Omicron在内的多个变种，它们陆续成为全球传播的主要变种病毒。以英国疫情为例，英国依次经历了分别以原始病毒（SARS-CoV-2）、Alpha、Delta和Omicron变种病毒为主的疫情阶段。为了更好地预判出原始病毒和变种病毒的传播情况，根据其在不同的气候、温度、地理位置等不同的条件下，结合实际，建立数学模型研究下列问题：

根据英国日新增的病例数据，讨论划分出原始病毒和变种病毒的传播阶段，结合感染病毒人数，建立出英国的原始型、Alpha和Delta三种病毒各自的传播速率及比较不同病毒的传染力和持续时长的数学模型。
根据天气数据，比较多种天气因素，分析对英国379个地区的原始病毒、Alpha和Delta三种病毒传播的影响。

假设与符号说明

假设

假设感染人数数据连续
假设病毒突变时间不重合
假设两毒株间有“竞争关系”
英国人口数量较为稳定，且病毒感染人数会有一个最大值
环境条件对种群的增长有阻滞作用，并随种群密度增加而按比例增加
人类免疫系统对不同毒株间的特异性很敏感
假设后期alpha型病毒感染人数为上限
假设只有阳性、阴性、死亡以及痊愈四种情况，并且阳性均住院
假设原始型的死亡人数为0，全部确诊病例均去医院且疫苗接种率为0

符号说明

问题一

解决思路

我们认为解决本问题的核心是要依据短期数据预测中期数据，并且要合理推测出不同毒株主要的传播时间段。

首先我们将病毒的传播看作为一个种群的繁衍。通过这样一个抽象化的过程，以及考虑到英国的总人口是有限的，所以我们想到了“S”型模型，所谓“S”型模型也就是logistics模型。

然而logistics模型通常是用来进行中长期数据的预测，而目前整理数据后发现只有delta型毒株有较为完整的感染人数数据。而alpha型毒株和原始型毒株情况则分别为只有后期数据以及无具体数据。所以根据假设3，我们应当先尽量完美的拟合，找到数学特征，再与已知中期数据的毒株感染人数进行对比，从而确认时间节点。整体流程如下流程图：

流程图.jpg 通过查阅和整理英国新冠病毒每周基因测序的抽测结果^[1]分别得到两种毒株不同时期的感染人数数据（见附件）其中，alpha毒株以及delta毒株感染人数散点图分别如下：

alpha&delta.jpg 绘制代码如下：

import pandas as pd
import matplotlib.pyplot as plt

# 设置中文字体，不要删，不设置就都是框框
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']
# 读取csv文件
df = pd.read_csv('alpha_delta.csv')
# 先画alpha的
plt.subplot(2, 1, 1)
# 下面color参数设置成这样是为了画出来好看点 lol
plt.scatter(df.x, df.alpha, color = plt.cm.get_cmap('Set3')(0))
plt.title('alpha感染人数散点图')
plt.xlabel('时间')
plt.ylabel('人数')
# 再画delta的
plt.subplot(2, 1, 2)
plt.scatter(df.x, df.delta, color = plt.cm.get_cmap('Set3')(2))
plt.title('delta感染人数散点图')
plt.xlabel('时间')
plt.ylabel('人数')
# 这个不要删，要不然字会重叠上
plt.tight_layout()
# 设置分辨率，防止画质糊化
plt.rcParams['savefig.dpi'] = 600
# 保存图片
plt.savefig('alpha&delta.jpg')

由散点图可观察到两数据有一定的线性关系。并且alpha感染人数属于短期数据，而delta感染人数数据属于中期数据。通过观察alpha毒株数据形如反比例函数且又呈线性关系，因此建立数学模型：

将（1-1）反比例函数式带入一次函数中：

通过SPSS进行回归，得到结果如下：

效果图如下：

alpha逆函数.png 在回归结果中，调整R²=0.861，较为接近1，说明拟合优度良好。F检验中，F并不大于临界值，对应概率P值小于0.001，回归方程整体并不显著，但是通过拟合图仍能看出一定趋势。系数β₀、β₁对应的t检验概率P值均小于0.001，说明该系数比较显著。

根据delta感染人数散点图可以看出，前期感染人数为0，而后期数据则在>1.5×10⁶处感染人数增长缓慢或停止，说明我们得到的是中期数据，因此我们依据假设4和假设5建立logistics模型。首先建立微分方程：

通过SPSS进行回归，得到结果如下：

效果图如下：

在回归结果中，调整R²=0.931，接近1，说明拟合优度很好。F检验中 F值较大，对应概率P值小于0.001，回归方程整体较为显著。 N对应的t检验概率P值均小于0.001，说明该系数比较显著。虽然常量对应的t检验概率P值为0.006，但仍小于0.05，所以我们仍认为常量仍具有较好的统计性差异。

由alpha感染人数图可以看出在x=10以前感染人数出现骤减，而结合delta感染人数图可以看到在x=10时正是delta感染人数骤增的时候，所以我们认为delta的感染人数增加对alpha感染人数增加有抑制作用，验证了假设3的合理。因此我们推测delta的开始传播也在x=10之前。而结合新冠病毒有一定的潜伏期无症状这一点来看，我们推测很有可能是在x=3，也就是2021年6月3日开始传播。我们有通过查询新闻获得了发现delta病例的第一天，为2021年6月15日，这之间的时间间隔与新冠的潜伏期十分接近。

alpha型的中期数据预测

根据上文得到的经验，我们根据假设7利用logistics模型以及6月之前的总感染人数以及alpha发现时间进行预测。

通过SPSS进行回归，得到结果如下：

在回归结果中，调整R²=0.848，说明拟合优度较好。F检验中 F值较大，对应概率P值小于0.001，回归方程整体较为显著。 N的系数以及常量对应的t检验概率P值均小于0.001，说明均比较显著。

时间段的划分以及原始型的中期数据预测

接着我们尝试对不同毒株开始传播时间进行预测。首先我们先根据实际alpha第一例病例发现时间划分出了原始型以及alpha型时期实际数据，如下图。

时间点.jpg 绘制代码如下：

import pandas as pd
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']
df = pd.read_csv('new_cases.csv')
raw = df[df.x < 32]
alpha = df[(df.x >= 32) & (df.x < 513)]
delta = df[(df.x >= 513) & (df.x < 716)]
other = df[df.x >= 716]

plt.subplot(2, 2, 1)
plt.plot(raw.x, raw.cumCasesBySpecimenDate, color = plt.cm.get_cmap('Set3')(0), label = '原始型')
plt.plot(alpha.x, alpha.cumCasesBySpecimenDate, color = plt.cm.get_cmap('Set3')(2), label = 'alpha')
plt.plot(delta.x, delta.cumCasesBySpecimenDate, color = plt.cm.get_cmap('Set3')(3), label = 'delta')
plt.plot(other.x, other.cumCasesBySpecimenDate, color = plt.cm.get_cmap('Set3')(4), label = '其他')
plt.fill_between(raw.x, raw.cumCasesBySpecimenDate, facecolor = plt.cm.get_cmap('Set3')(0), alpha = 0.25)
plt.fill_between(alpha.x, alpha.cumCasesBySpecimenDate, facecolor = plt.cm.get_cmap('Set3')(2), alpha = 0.25)
plt.fill_between(delta.x, delta.cumCasesBySpecimenDate, facecolor = plt.cm.get_cmap('Set3')(3), alpha = 0.25)
plt.fill_between(other.x, other.cumCasesBySpecimenDate, facecolor = plt.cm.get_cmap('Set3')(4), alpha = 0.25)
plt.vlines(32, 0, 2e+7, linestyles = 'dashed', colors = 'red', label = '发现时间')
plt.vlines(15, 0, 2e+7, linestyles = 'dashed', colors = 'blue')
plt.vlines(513, 0, 2e+7, linestyles = 'dashed', colors = 'blue', label = '推测时间')
plt.vlines(716, 0, 2e+7, linestyles = 'dashed', colors = 'blue')
plt.vlines(690, 0, 2e+7, linestyles = 'dashed', colors = 'black', label = '混合时期')
plt.legend(loc = 2, bbox_to_anchor = (1.05,1.0), borderaxespad = 0.)
plt.title('不同毒株传播时间点')
plt.subplot(2, 2, 2)
plt.plot(raw.x, raw.cumCasesBySpecimenDate, color = plt.cm.get_cmap('Set3')(0))
plt.vlines(15, 0, 60, linestyles = 'dashed', colors = 'blue')
plt.fill_between(raw.x, raw.cumCasesBySpecimenDate, facecolor = plt.cm.get_cmap('Set3')(0), alpha = 0.25)
plt.title('原始型时期')
plt.subplot(2, 2, 3)
plt.plot(alpha.x, alpha.cumCasesBySpecimenDate, color = plt.cm.get_cmap('Set3')(2), label = 'alpha')
plt.title('alpha时期')
plt.fill_between(alpha.x, alpha.cumCasesBySpecimenDate, facecolor = plt.cm.get_cmap('Set3')(2), alpha = 0.25)
plt.subplot(2, 2, 4)
plt.plot(delta.x, delta.cumCasesBySpecimenDate, color = plt.cm.get_cmap('Set3')(3), label = 'delta')
plt.vlines(690, 0, 2e+7, linestyles = 'dashed', colors = 'black', label = '混合时期')
plt.fill_between(delta.x, delta.cumCasesBySpecimenDate, facecolor = plt.cm.get_cmap('Set3')(3), alpha = 0.25)
plt.title('delta时期')
plt.tight_layout()
plt.rcParams['savefig.dpi'] = 600
plt.savefig('时间点.jpg')

我们看到原始型时期在x=20～30之间有一段处于上升状态，而又考虑到病毒有一定的潜伏期且期间无症状、发现需要时间以及当时英国政府可能因为当时确诊病例较少而放松警惕，所以我们决定将推测时间向前15天。

接下来我们又对原始型数据进行了同样的预测，得到结果如下：

在回归结果中，调整R²=0.859，说明拟合优度较好。F检验中 F并不大于，可能受样本量限制，但影响并不大，对应概率P值小于0.001，回归方程整体较为显著。N的系数对应的t检验概率P值均小于0.001，较为显著，而常量检验P值为0.055，不太显著。

拟合效果如下图：

原始型logistics.png

传播速率

通过阅读文献，我们发现logistics模型的导数具有对称性，也就是瞬时增长量会在一点达到最大值，而这一点则正好是k/2。^[2]我们已知logistics模型的上限，可以直接求出。 matlab求解代码：

a = readmatrix('/Users/yanghaolin/Downloads/2022北京高校数学建模校际联赛赛题/B题/k.csv');

display(a(:, end) / [2; 2; 2]);

传染力

通过了解，我们将用流行病学上的基本传染数代表传染力。通过查询文献^[3]，我们获得以下公式：

为了简化模型方便我们计算，我们将时间依赖性传播率代替为最大瞬时增长量与当时的人口总量之比，并根据假设8得到下式：

在此模型的基础上，我们有考虑了易感情况。最终考虑因素思维导图如下：

传染力考虑因素.png 根据基本生物学知识、生活经验以及简单的生物信息学知识，易感程度越高应当使病毒的传染力增加，所以需要将原本R₀乘易感程度，得到下式：

通过检索遗传学文献^[4]，我们利用python计算欧洲人与世界平均多基因风险评分之比，结果为0.98667，小于1，说明易感风险略降低。根据查到的数据^[5]可知易感年龄占比为37% 接下来经过计算即可算出感染力。

PRS python计算代码如下：

import pandas as pd
import numpy as np

df = pd.read_csv('genetics.csv')
freq_arr = np.array(df.iloc[:, -2:])
risk_arr = np.array([df.OR]).T
prs = np.prod(np.power(risk_arr, 2 * freq_arr), axis = 0)
print(prs[0]/prs[-1])

问题二

解决思路

根据题目所给的附件我们有了英国各地区有关温度、露点、湿度等气候数据。我们的目标是利用所得到的数据，去比较分析温度、露点、湿度等气候因素对三种不同病毒传播的影响。对于上述资料与题目。我们通过使用典型相关性分析，将具有显著统计意义的因素留下，剩下因素忽略不计。

以牛津的数据为例

我们将6种天气因素（气温、露点、湿度、风速、气压、降水）以及牛津每日新增阳性人数利用SPSS进行典型相关性分析，并通过对比其皮尔逊系数以及其P值来确定两者之间是否相关以及相关程度。

天气与每日新增阳性人数相关

截屏2022-05-03 下午4.21.59.png 整体天气因素与每日新增阳性人数相关性为0.382，接近于0.4，说明接近于中等相关的水平，而其P值趋近于0，说明统计意义显著。最终得出结论天气与每日新增阳性人数相关

温度、湿度对每日新增阳性人数影响最大

截屏2022-05-03 下午4.22.11.png 我们通过python进行数据可视化，方便我们分析。

python代码如下：

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from seaborn import *

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']
df = pd.read_csv('pearson.csv', index_col = 0)
plt.subplot(1, 2, 1)
heatmap(abs(df))
plt.title('各因素间皮尔逊系数')
plt.subplot(1, 2, 2)
stem = plt.stem(df.iloc[:, 0], basefmt = '--')
plt.xticks(np.arange(len(df.index)), np.array(df.index), rotation = 45)
plt.hlines([0.2, -0.2], [0, 0], [6, 6], color = 'pink', linestyles = 'dashed', label = '弱相关')
plt.legend()
plt.setp(obj = stem, color = plt.cm.get_cmap('Set3')(0))
plt.title('新冠新增人数与天气因素间皮尔逊系数')
plt.tight_layout()
plt.rcParams['savefig.dpi'] = 600
plt.savefig('heatmap_pearson.jpg')

每日新增阳性人数我们可以初步看出气温、露点、湿度三个温度对每日新增阳性人数影响最大且三因素P值均小于0.01，说明三因素具有统计显著性。然而当我们观察到露点与气温这两个因素的相关性时发现皮尔森系数高达0.862且P值小于0.01，说明露点与气温呈强相关。但因与每日新增阳性人数相关的因素已经有了气温这一因素，且其相关性高于露点，所以我们决定舍弃露点这一因素，直接得出结论：温度、湿度对每日新增阳性人数影响最大。

参考文献

[1]王玉龙.“J”型与“S”型曲线的数学推导及相关概念论证[J].生物学教学,2021,46(06):76-79.

[2]胡晓冬,董辰辉.MATLAB从入门到精通[M].北京:人民邮电出版社,2018

[3]Tianxiang Yue, Bin Fan, Yapeng Zhao, John P. Wilson, Zhengping Du, Qing Wang, Xiaozhe Yin, Xiaonan Duan, Na Zhao, Zemeng Fan, Hui Lin, Chenghu Zhou. Dynamics of the COVID-19 basic reproduction numbers in different countries[D]. Science Bulletin, 2021, 66(03):229-232.

[4]Roberts, G. H. L. et al. AncestryDNA COVID-19 host genetic study identifies three novel loci. Preprint at doi.org/10.1101/202… (2020).

[5]华经产业研究院.2010-2020年英国人口数量及人口性别、年龄、城乡结构分析[R].华经情报网,2021