导言
前段时间英国采取的“群体免疫”政策,似乎引起了很多关注,当然同时也引起了很多争议。我们希望从应用数学的角度出发,带大家看一下数学上是如何对传染病进行建模的,他背后的原理是什么。
数学难们别被应用数学吓到我们会从最易懂的角度去解析这样一个问题。
首先,一个与“群体免疫”息息相关的概念,即媒体不断提到的R0值。虽然看上去这两个概念并无明显的联系,但我们会在下文为大家解释R0的含义以及它是如何与群体免疫这一概念联系起来的。
注:图片来自网络
01、关于R0与群体免疫
R0,即基本传染数,是量化一种流行病的标准,通俗来说表示了一个得了某种流行病的人平均能传播给几个人。简单地说,R0小于1,意味着得病的人会越来越少,流行病会逐渐消失;R0大于1,流行病则会扩大传播;R0等于1或者在1左右,总有人治愈的同时就有人患病,流行病会成为一种地区性传染病一直存在。一般通过控制流行病传播使R0值降低到1以下就意味着控制疫情的成功。
我们再来看群体免疫的概念。其实群体免疫的概念也很易理解。整个人口群体中应有一部分人免疫病毒(无论是由于感染后痊愈而携带抗体的还是天生免疫),也有一部分人不免疫。那么对于R0值来说,这个传染数应该只会影响不免疫的人群。我们假设免疫人群占总体人群的比例为a,则不免疫的人群比例就是1-a,所以实际可能被感染的人群是R0*(1-a),他其实就代表了一个新的R0,考虑免疫人群下的R0。按照上面的推理,实际上我们把这个数字控制在1以下就达到了胜利,病毒会越传越少。即我们可以得到以下不等式:R0*(1-a)<1。
大家会发现这个不等式很好解,结果是a>1-1/R0. 这就是我们能达成群体免疫所需人群比例的下限。由于现在关于新冠病毒的数据还不完整,我们先假设假设这个R0的数值大概在3-4之间,所以a的数值应该在66.7%到75%之间。因此我们可以得出这样的结论,如果英国人群当中有66.7%到75%之间的人免疫新冠,则R0*(1-a)能够被控制到1以下,也就是成功了。当然实际情况要复杂得多,比如人口的出生率死亡率等等,这里只是最简单的模型预测。
下面的图像也方便大家更直观的理解群体免疫的作用。
注:图片来自网络
以上就是关于R0以及群体免疫最直观的解释了,但实际上其背后的理论层次是要深得多的。为了研究疫情的发展,我们通常会使用一种研究传染病大致走向的动力学模型----SEIR。
02、SEIR模型你需要知道什么?
“SEIR”其实是一组英文单词的缩写,表示的是易感人群(susceptible),处于潜伏期的人群(exposed),感染者(infectious)以及治愈者(recovered)。不难发现,这已正是一个个体从没有生病到痊愈的全部过程。当然这是一个相对基本的数学模型,但是现实生活往往是更加复杂,有更多变量的。
在数学模型中,我们设定某个群体的总人数为N(不考虑出生率和死亡率,因此N是一个常数)。而S,E,I,R,也不难发现,是四个关于时间的函数。那么我们可以列出这样的一个等式,即在任何时刻,都有N=S(t)+E(t)+I(t)+R(t)。同时我们假定所有的S人群的变化都会转化成E人群,而E人群的会全部转化为I人群,以此类推。
首先来看S的函数。
易感人群每一天相较于前一天的改变人数是多少呢?按照我们上面所说,S减少的人,即是E增多的人。为了计算这些被感染的人的数量,首先,需要定义一个系数β,即每人每天接触的人数*被传染的概率(这里的概率指的就是一个未感染的人接触传染者后被传染的概率),比如平均每天一个人会接触到十个人,而每个被感染者接触过的人都有10%的概率被感染,则这个数值就是1。
由于β定义中的概率是针对未被感染的人群,所以我们这里的β*I就表示的是假设全体人群都是易感染群的情况下所有感染者一天能传染多少人,因此需要乘上S/N(由于所有人不全是易感人群也有E和R人群,所以我们乘的是S人群在人群总体中占的比例)以及时间长度t,这便是新一天的变化量了。当然再加上原有的易感人群,我们就能得出关于易感人群的方程(注意易感人群在疾病传染过程中一定是减少的,所以有负号),即:
而为了方便研究每一天的变化,我们把S0挪到左边去,就得到了每天的变化。
我们希望的是研究S关于t的函数,然而如果t以天计算我们得到的就只是个散点图,我们必须要拟合出一个函数来才行。所以我们把时间从天缩小到一个极其微小的单位来让t在每一个点都能有一个对应的函数值。我们可以用一个微小量dt来表示时间非常少,而同样的因为时间很短S的变化也是很小的,所以用dS。那么再把dt除到等式左边来,我们就得到了这样的方程,即:
初中物理都教过我们,路程除以时间是速率,S关于t的函数也是一样的。上面的方程其实就表示了S减少的速率。
解决了S(t)函数的问题,我们再来看其他几个函数,由于已经解释了这一套推导的过程,接下来会根据上面的解释进行分析。
我们已经知道了易感人群关于时间的变化速率的函数,下面我们来看潜伏期人群(exposed)的。我们通过模型可以发现,所有的S的减少人数最后都变成了E的新增人数,所以可以先把上面的式子抄下来,当然不能加负号了。同时我们还要减去E人群减少的人数,因为E人群同时在转化为感染者。(注意我们这里假设潜伏期的人没有传染的能力)
下面我们来处理E转化为I的部分。我们知道潜伏期的人在潜伏期过后会变成感染者,因此我们需要新定义另一个系数α,它代表了单位时间内从潜伏期过渡到感染者的概率,而1/α则表示了潜伏期的长度。举例来说,如果潜伏期为14天,即1/α=14天,那么α就等于1/14,也正好是每一天一个人从潜伏期到被确诊的概率。
那么再使用这个系数α乘上现有E人群的总体,就应该是E人群转化为I人群的速率。所以整体方程为:
我们以同样的思路想一下下面的函数
为了观察感染者是如何变化的,首先我们确定I关于时间t的函数由三部分组成,I0是最开始的感染者,即为一个常数,加上一个潜伏者转成感染者的人(我们提到过这是一个动力学模型,所以上一个的减少量即是这一个的增加量),还要减去治愈的人,这里的γ就是被治愈成功的概率。
为了得到I的改变速率,我们将I0放到等式的左边得到了δI即I的该变量,再把t除过去,这个时候为了函数的完整性,我们将t设为微小量,I的改变量相应的也变得无限小,就可以得到以下式子。
R也就很简单了,同理:
(这里的R0不是我们上面所提到的基本传染数,是recover的初始值,当然在我们从疫情最开始考虑的时候,他默认是0)
R的改变速度自然就是上面I转移过来的速度,即是治愈的概率。
在我们完成了对几个基本方程的列举之后,我们也就可以回答R0是怎么构成的这个问题。
大家如果还记得我们之前定义的几个系数的话,这个就很好解释了。R0其实就等于前面的β/γ(那两个系数)。这也从数学上给我们上了重要的一课:由于β在分子上,所以当每个人每天接触的人较少时,β就会更小,相应地R0也会减小。这体现了隔离的重要性,因为隔离能显著地减少每个人每天接触的人数。
03、如何建模
了解了这些函数之后可能会有一个问题,我们如何去建模呢?
我们会用到一些建模软件完成对这个方程的求解(如matlab,geogebra)
首先我们为了建模要完成数据的收集。方程中很多的常数都需要我们收集生活中实际的数据进行求算,如β,α等。数据来源自选即可,对于武汉的疫情数据可以在国家卫建委的网站上收集。对于代码编写这个就会涉及到一些matlab等软件中的算法,我们不会详细展开,但是我们可以找到网上的一些开源的代码,他们都是建构好的模板,如果不需要添加更复杂的结构可以就按模板来。然后在相应的地方填入数据。
这里提供一个在github上扒的开源代码,有兴趣的可以进去看看。
网页链接:https://github.com/ECheynet/SEIR
04、帝国理工预测
下面我们来看一下帝国理工Covid-19团队的建模预测:
相信大家在听到英国群体免疫提出之后,也听说了帝国理工发布了一个模型预测实行群体免疫那种政策将会有数十万人死亡。那我们带大家一起看一下这样的一个报告。
(以下都是官方公开的文件)
原报告标题为:以非药物干预方法减少covid-19死亡率和医疗需求的影响
他们在官方发布的报告中声明的是使用了transmission model,并没有声明具体使用了哪样一种模型。但是他对各项系数的假设以及模型的准确性是很有参考性的。
他将潜伏期设为了5.1天,也就是说如果我们按他的数据来的话我们式中的α约等于0.2;他将R0假设为2.4。模型上不同的是,他考虑的模型中,假设所有已有症状的人比潜伏期的人有50%的概率更有可能感染别人,而我们的模型中E是百分之百不传染人的。
接下来,他又通过中国的数据预估了死亡率,如下图:
下面这两条曲线展示了如果 do nothing需要的重症监护病床位是极大的。
这是他们根据数据预测出来的死亡曲线:
结语
SEIR模型看似只是用建模软件解个几微分方程式就可以让我们预测到未来的传染病走势,但是实际上数学预测模型中还有很多精细却又复杂难以确认的东西需要注意,对于实际生活中的疫情,疫情的较精确模型远远要比SEIR复杂,一些国家需要考虑到隔离执行力问题,如武汉的隔离执行力就很强,中国还有大量的自我隔离者,即会在易感人群的数量上与其他国家产生很大的偏差,这种情况有些学者就会用到SIQR的模型,以Q来计算到隔离人数;这个模型是没有计算出生与死亡率的,但实际上对于一些严重的病毒,死亡率肯定要计入其中;以及RNA的改变特性,可能会导致一些recover的人重新回到S变成一个循环,像这些情况会用到一些SEIR的变式,一些更复杂的传染病模型,真实的建模当中永远是灵活使用,应用中没有死板的公式。
关于群体免疫,显然它的可行性是有待考证的,仍然需要大量数据的支撑,我们在这里讨论的只是相关概念和理论。
声明:文章收集于网络,版权归原作者所有,为传播信息而发,如有侵权,请联系小编删除,谢谢!
欢迎加入本站公开兴趣群
高性能计算群
兴趣范围包括:并行计算,GPU计算,CUDA,MPI,OpenMP等各种流行计算框架,超级计算机,超级计算在气象,军事,航空,汽车设计,科学探索,生物,医药等各个领域里的应用
QQ群:326600878