基于数据驱动的一类具有信息效应的传染病模型研究

Research on a Class of Infectious Disease Model with Information Effects Based on Data Driven

摘要: 本文我们提出了一类具有信息效应的SS_aEAIRD传染病动力学模型来揭示非药物干预措施对传染病传播的影响。首先，给出了模型的基本再生数，研究了无病平衡点的局部渐近稳定性；其次，基于我国COVID-19疫情传播的病例数据，利用bootstrap抽样、极大似然估计以及非线性最小二乘拟合对该模型的参数进行了估计；通过95%置信区间可以看出模型拟合的可靠性；最后，为了确定关键参数对疫情发展的影响，还进行了不确定性和敏感性分析。研究结果表明，加大信息宣传效应会使个人的行为发生改变，有意识的人群增多，对疾病的防护意识增强；还会降低感染者对无意识易感者的有效传染率，减少感染的人群，从而更好地控制疾病的传播。

Abstract: In this paper, we propose a class of SSaEAIRD infectious disease dynamics models with information effects to reveal the impact of non pharmacological interventions on the spread of infectious diseases. Firstly, the basic reproduction number of the model is given, and the local asymptotic stability of the disease-free equilibrium point is studied; secondly, based on the case data of COVID-19 transmission in China, the parameters of the model are estimated by using bootstrap sampling, maximum likelihood estimation, and nonlinear least squares fitting. The accuracy of model fitting can be seen through the 95% confidence interval; finally, in order to determine the influence of key parameters on the development of the disease, uncertainty and sensitivity analyses are also conducted. The research results indicate that increasing the effect of information dissemination can lead to changes in human behavior, an increase in conscious population, and an improvement in disease prevention. It will also reduce the effective transmission rate of infected individuals to unconscious susceptible individuals, reduce the number of infected individuals, and better control the spread of the disease.

1. 引言

自人类发展以来，传染病一直是危害人类身心健康和社会发展的一个重要因素，全世界每年有四分之一的死亡是由传染病造成的，其中绝大多数发生在低收入和中等收入国家 [1] 。有一些疾病，如SARS，流感以及COVID-19，它们都有着相似的特征，并且会表现出一些明显的特征，如空间传播迅速和症状明显 [2] 。这些与全球化趋势日益加剧和信息技术发展有关的特点，预计也会出现在新发的传染病中。因此，重要的是要完善经典的数学模型，通过增加大规模新闻报道的维度来反映这些特征，这些维度不仅对个人行为有很大的影响，而且对公共干预和控制政策的形成和实施也有很大的影响 [2] 。

在新冠肺炎疫情抗击过程中，经济和信息技术的迅猛发展使得媒体能够迅速传播有关传染病的信息。这种信息传播在潜移默化中影响公众的心态，引起了一定的恐慌情绪，促使人们改变日常行为，比如减少社交接触，以降低感染风险。这一现象表明，在传染病的爆发和蔓延过程中信息传递的变得越来越重要，对疫情的流行趋势有着显著地影响。据此，许多学者结合新冠病毒的传播特性，对具有信息影响的传染病模型进行了研究，并对疫情的预测、防控策略进行了深入的探讨。Cui等 [3] 提出了一类人群具有Logistic自然增长和媒体报道影响的SEI传染病模型建立了一个传染率函数描述媒体报道对疾病传播的影响，研究表明m值较小时，模型可能会出现一个Hopf分支。Xiao等 [4] 考虑了受到媒体影响，传染病动力学模型成为一个切换系统，其中媒体影响可能多次打开和关闭，研究了2003年SARS疫情和2009年流感大流行的多波流行病。Gui等 [5] 提出了一个改进的SEIRS模型，考虑了具体信息导致大流行发展的时间延迟以及对时间免疫反应的复苏性。Leonardo等 [6] 在SEIR模型的基础上，引入了不同比例的遏制效应，考虑潜伏期感染的传播，以评估流行病爆发期间发病率和死亡率的比率。王琪等 [7] 研究了在信息作用的影响下心理效应的媒介传染病SEIRS模型，并对改变自身行为，染病者的数量进行评估。史学伟等 [8] 研究了一类具有信息变量，饱和发生率和等级治愈率的SIR传染病模型，并探讨了无病平衡点和地方病平衡点的存在性和稳定性。张杰豪等 [9] 研究了一类疾病发生率与潜伏期人群移出率受信息因素影响的SEIR传染病模型，通过数学分析得到了模型平衡点的存在性。Manuel等 [10] 探讨了引入时变接触率降低社区传播所需的行为改变的效果，以及疾病在非必要活动暂停的人群中传播的后果。魏永越等 [11] 考虑到新冠肺炎的感染人群类型以及隔离措施，建立了一类SEIR^+CAQ传播动力学模型，并对全国和湖北省的新冠肺炎疫情的发展趋势进行了预测。Bavel等 [12] 探讨了在社会信息的影响下，行为变化对传染病的多个阶段的影响，并为政策制定者提出了多种见解。Weitz等 [13] 对美国的新冠肺炎疫情进行分析，建立了一个死亡驱动意识的模型。探讨了在信息宣传的影响下，个人保护措施随着死亡率的增加而增加。这些研究工作对我们认识COVID-19的传播规律，进而为预防和控制疫情提供了重要的参考。

然而，反映信息效应在突发传染病的疫情影响中，考虑有意识和无意识易感者之间可以相互转化的动力学模型尚不多见，本文通过对新冠肺炎传播规律和特点的分析，并考虑到COVID-19的潜伏期和无症状均具有传染性和信息影响的大小，建立了一类SS_aEAIRD传染病传播模型来研究信息效应对疫情传播的影响。通过观察受信息效应影响后有效传播率和因病死亡率，分析调节信息效应对易感者的转化率预测累计确诊病例数的变化情况。

2. 模型建立

2.1. 数据来源

本文使用的COVID-19数据来源于中华人民共和国国家卫生健康委员会，数据包括2022年2月20日至2022年5月1日中国累计确诊病例数、累计治愈病例数和现有确诊病例数 [14] ，如图1所示。

Figure 1. The number of cumulative confirmed cases, cured cases, and existing confirmed cases in China from February 20 to May 1, 2022

图1. 2022年2月20日至5月1日中国累计确诊、累计治愈和现有确诊病例

2.2. 模型建立

根据COVID-19传染病的传播特点，随着互联网技术的不断进步，人们获取消息的途径也越来越丰富且快速，在传染病扩散期间，相关信息宣传会对COVID-19的传播产生影响，因此本文主要研究了信息效应对新型冠状病毒传播的影响。在SEIR型(易感–潜伏–感染–恢复)流行病学模型的基础上考虑加入无症状感染者和受信息影响后有防护意识的易感者。现对模型作出如下假设：

1) 不考虑人口的出生率和死亡率；

2) 由于受信息影响后有意识的易感者防护意识较强，被感染的可能性很低，因此不考虑有意识的易感者转化为潜伏者；

3) 潜伏者和无症状感染者均具有传染性。

根据上述的假设，现将总人口(N)分为七个群体：无意识的易感者(S)、有意识的易感者(S_a)、潜伏者(E)、无症状感染者(A)、有症状感染者(I)、恢复者(R)和死亡者(D)。本文考虑受信息影响后有意识的易感者和无意识易感者之间的转变，其中，m是受信息宣传影响后无意识的易感者变为有意识的易感者的转化率，ω是失去防护意识后有意识的易感者变为无意识的易感者的转化率。下面给出各仓室之间的传播仓室图，如图2所示，相应的传染病模型见(1.1)，同时模型(1.1)中的参数和变量的含义见表1。

Figure 2. Flow chart of COVID-19 transmission in various compartments

图2. COVID-19在各仓室的传播仓室图

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=-\frac{\beta \left(I+\theta A+vE\right)}{N}S-mS+\omega S_a\\ \frac{\text{d}{S}_a}{\text{d}t}=mS-\omega S_a\\ \frac{\text{d}E}{\text{d}t}=\frac{\beta \left(I+\theta A+vE\right)}{N}S-\sigma E\\ \frac{\text{d}A}{\text{d}t}=\sigma \left(1-\rho \right)E-\alpha A-{\gamma }_{1}A\\ \frac{\text{d}I}{\text{d}t}=\sigma \rho E+\alpha A-{\gamma }_{2}I-\delta I\\ \frac{\text{d}R}{\text{d}t}={\gamma }_{1}A+{\gamma }_{2}I\\ \frac{\text{d}D}{\text{d}t}=\delta I\end{array}$ (1.1)

3. 无病平衡点的稳定性

研究模型(1.1)的动力学性质首先要确定平衡点，故本节将讨论无病平衡点的稳定性。在初始时刻，人群中全是易感者，故模型(1.1)对应的疾病消亡的无病平衡点 $P^0\left(S^0,S_a^0,0,0,0\right)$ ，根据文献 [15] 中再生矩阵的方法，得到模型(1.1)的基本再生数为

${\mathcal{R}}_{\text{0}}=\frac{\beta \nu }{\sigma N}{S}^0+\frac{\beta \theta \left(1-\rho \right)}{M_{1}N}{S}^0+\frac{\beta \left(\alpha +{\gamma }_1\rho \right)}{M_1{M}_{2}N}{S}^0=\frac{\beta (\nu M_1{M}_2+\theta \left(1-\rho \right)\sigma M_2+\left(\alpha +{\gamma }_1\rho \right)\sigma )}{\sigma M_1{M}_{2}N}{S}^0$

其中， $M_1=\alpha +{\gamma }_1$ ， $M_2={\gamma }_2+\delta$

通过模型(1.1)，我们发现 $E,A,I$ 三个方程与其他方程无关，所以我们可以在 $\left(E,A,I\right)$ 组成的子系统上讨论疾病的灭绝性。

定理1. 当 ${\mathcal{R}}_0<1$ 时，模型(1.1)的 $\left(E,A,I\right)$ 组成的子系统的零解局部渐近稳定；当 ${\mathcal{R}}_0>1$ 时，该零解不稳定。

证明： $\left(E,A,I\right)$ 组成的子系统的系数矩阵为

$J=\left[\begin{array}{l}\lambda -\beta \nu \frac{S^0}{N}\text{+}\sigma -\beta \theta \frac{S^0}{N}-\beta \frac{S^0}{N}\\ \sigma \left(1-\rho \right)-\left(\alpha +{\gamma }_1\right)\text{0}\\ \sigma \rho \alpha -\left({\gamma }_2+\delta \right)\end{array}\right]$

显然， $J$ 有三个特征根 ${\lambda }_i\left(i=1,2,3\right)$ ，且 ${\lambda }_i\left(i=1,2,3\right)$ 满足特征方程 ${\lambda }^3+a_1{\lambda }^2+a_2\lambda +a_3=0$ ，其中

$a_1=M_1+M_2-\beta \nu \frac{S^0}{N}+\sigma$ $a_2=M_1{M}_2-\beta \nu \frac{S^0}{N}\left(M_1+M_2\right)+\sigma \left(M_1+M_2\right)-\theta \beta \sigma \left(1-\rho \right)\frac{S^0}{N}-\beta \sigma \rho \frac{S^0}{N}$

$a_3=-\beta \nu M_1{M}_2\frac{S^0}{N}+\sigma M_1{M}_2-\theta \beta \sigma \left(1-\rho \right)M_2\frac{S^0}{N}-\beta \sigma \alpha \left(1-\rho \right)\frac{S^0}{N}-\beta \sigma \rho M_1\frac{S^0}{N}$

当 ${\mathcal{R}}_0<1$ 时，

$\frac{\beta \nu }{\sigma N}{S}^0<1\iff \frac{\beta \nu }{N}{S}^0<\sigma \Rightarrow a_1>0$ ,

$a_3=\sigma M_1{M}_2\left(1-{\mathcal{R}}_0\right)>0$ ,

$\beta \nu \frac{S^0}{N}{M}_1{M}_2+\sigma \theta \beta \left(1-\rho \right)\frac{S^0}{N}{M}_2+\sigma \beta \frac{S^0}{N}\left(\alpha +\rho {\gamma }_1\right)<\sigma M_1{M}_2\iff$

$\beta \nu \frac{S^0}{N}{M}_1{M}_2+\sigma \theta \beta \left(1-\rho \right)\frac{S^0}{N}{M}_2+\beta \sigma \frac{S^0}{N}\left[\rho M_1+\alpha \left(1-\rho \right)\right]<\sigma M_1{M}_2$

有

$\beta \nu \frac{S^0}{N}{M}_1{M}_2+\sigma \theta \beta \left(1-\rho \right)\frac{S^0}{N}{M}_2<\sigma M_1{M}_2\Rightarrow \beta \nu \frac{S^0}{N}{M}_1+\sigma \theta \beta \left(1-\rho \right)\frac{S^0}{N}<\sigma M_1$

$\beta \nu \frac{S^0}{N}{M}_1{M}_2+\beta \sigma \rho \frac{S^0}{N}{M}_1<\sigma M_1{M}_2\Rightarrow \beta \nu \frac{S^0}{N}{M}_2+\beta \sigma \rho \frac{S^0}{N}<\sigma M_2$

$a_2=M_1{M}_2-\beta \nu \frac{S^0}{N}\left(M_1+M_2\right)+\sigma \left(M_1+M_2\right)-\theta \beta \sigma \left(1-\rho \right)\frac{S^0}{N}-\beta \sigma \rho \frac{S^0}{N}>0$

$\begin{array}{l}{a}_1{a}_2-a_3=\left[\left(M_1+M_2\right)-\beta \nu \frac{S^0}{N}+\sigma \right]\left[M_1{M}_2-\beta \nu \frac{S^0}{N}\left(M_1+M_2\right)+\sigma \left(M_1+M_2\right)-\theta \beta \sigma \left(1-\rho \right)\frac{S^0}{N}-\beta \sigma \rho \frac{S^0}{N}\right]\\ +\beta \nu \frac{S^0}{N}{M}_1{M}_2-\sigma M_1{M}_2+\theta \beta \sigma \left(1-\rho \right)\frac{S^0}{N}{M}_2+\beta \sigma \alpha \left(1-\rho \right)\frac{S^0}{N}+\beta \sigma \rho \frac{S^0}{N}{M}_1\\ =a_1\left(a_2-M_1{M}_2\right)+a_1{M}_1{M}_2+\beta \nu \frac{S^0}{N}{M}_1{M}_2-\sigma M_1{M}_2+\theta \beta \sigma \left(1-\rho \right)\frac{S^0}{N}{M}_2+\beta \sigma \alpha \left(1-\rho \right)\frac{S^0}{N}+\beta \sigma \rho \frac{S^0}{N}{M}_1\end{array}$

$\begin{array}{l}=a_1\left(a_2-M_1{M}_2\right)+\left(M_1+M_2\right)M_1{M}_2+\sigma M_1{M}_2-\beta \nu \frac{S^0}{N}{M}_1{M}_2-\sigma M_1{M}_2+\beta \nu \frac{S^0}{N}{M}_1{M}_2+\beta \sigma \alpha \left(1-\rho \right)\frac{S^0}{N}\\ +\beta \sigma \rho \frac{S^0}{N}{M}_1+\theta \beta \sigma \left(1-\rho \right)\frac{S^0}{N}{M}_2\\ =a_1\left(a_2-M_1{M}_2\right)+\left(M_1+M_2\right)M_1{M}_2+\beta \sigma \alpha \left(1-\rho \right)\frac{S^0}{N}+\beta \sigma \rho \frac{S^0}{N}{M}_1+\theta \beta \sigma \left(1-\rho \right)\frac{S^0}{N}{M}_2\\ \text{>0}.\end{array}$

根据Routh-Hurwitz判据，当 ${\mathcal{R}}_0<1$ 时，矩阵 $J$ 的特征方程的所有特征值均具有负实部，因此零解局部渐近稳定。当 ${\mathcal{R}}_0>1$ 时， $a_3<0$ ，该特征方程有一个正根，因此零解不稳定。

4. 参数估计与数据拟合

根据世界卫生组织(World Health Organization (WHO), 2020)的相关报告，新冠肺炎疫情的潜伏期大约是5.2天，因此发病率 $\sigma =1/5.2$ 。

我们根据泊松分布生成500个每日报告病例数和累计确诊病例数的时间序列数据的自举样本。也就是说，对于每个数据点，我们假设每日报告病例数遵循泊松分布，并假设平均值为观察到的计数。基于这一假设，我们在预设的分布中抽取500个的日报告病例数，从而得到500个时间序列数据样本。首先使用最小二乘法对500个日报告病例数和累计确诊病例的bootstrap样本的每个时间序列进行拟合，我们使用MATLAB软件对部分参数采用先验分布的最小二乘法对数据进行模型拟合，其中ODE系统由“ODE45”函数求解，使用“fmincon”函数搜索目标函数的最优解。在进行拟合过程中，对于未知变量的每个初始值提供了先验分布。对于 $\omega$ 、 $m$ 、 ${\gamma }_1$ 、 ${\gamma }_2$ 和 $\delta$ 通过建立似然函数，使用极大似然估计得到最优的参数值。最佳的拟合结果如图3所示，具体参数估计值见表1。根据拟合出的值预测了2022年4月21到5月1日的累计确诊病例和实际数据进行了对比验证，如图4所示，我们可以发现，该拟合值可以更好地描述疾病的传播趋势。

然后，我们将模型(1.1)拟合到每日报告病例的500个bootstrap样本中的每个样本。并且，500次拟合可以产生500个参数集，因此产生500个累计报告病例数的解。然后，我们计算500个解的2.5%和97.5%

Table 1. The meaning and values of parameter in the model (1.1)

表1. 模型(1.1)中参数的含义及估计值

Figure 3. Fitting chart of the cumulative number of confirmed cases

图3. 累计确诊病例数拟合图

Figure 4. Comparison chart of fitting validation of the cumulative number of confirmed cases

图4. 累计确诊病例数的拟合验证对比图

分位数，以生成拟合结果分位数的95%置信区间。结果如图5所示，我们可以看出拟合值是落在置信区间内的。

Figure 5. 95% confidence interval for the cumulative number of confirmed cases

图5. 累计确诊病例数的95%置信区间

5. 误差分析

均方根误差(RMSE)是衡量预测值与实际观测值之间差异的常用指标。它是预测误差的平方和的平均值的平方根。在统计学中，它用于量化预测或估计的误差大小。RMSE越小，意味着模型的预测越接近实际值。其公式为

$\text{RMSE}=\sqrt{\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(y_i-{\hat{y}}_i\right)}^2}}$

其中，n是观测值的数量， $y_i$ 是第i个观测值， ${\hat{y}}_i$ 是第i个预测值。

在文中所建立的传染病动力学模型中，根据公式计算得实际累计确诊病例数和所拟合的累计确诊病例数之间的均方根误差为598.396。由于数据最大值达到了十几万以上，所以用该指标得到的误差值是合理的，并且相对于数据较大的模型拟合和预测有一定的准确性。

除此之外，我们还通过决定系数(R²)来评估模型拟合的好坏。R²接近于1，表明模型对数据的拟合度很高，说明模型对数据的拟合越好，能够很好地解释数据的变异性。其定义为：

$R^2=1-\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(y_i-{\hat{y}}_i\right)}^2}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(y_i-\bar{y}\right)}^2}}$

其中， $y_i$ 是第i个观测值， ${\hat{y}}_i$ 是第i个预测值， $\bar{y}$ 是实际观测值的平均值。

根据收集到的累计病例数，再结合拟合得到的数据，计算得到R² = 0.9068，意味着模型在预测和解释目标数据是有效的，该模型对于之后新发的传染病有一定的参考价值。

6. 敏感性分析

理论分析表明，基本再生数 ${\mathcal{R}}_0$ 是决定疾病是否流行的重要参数。为了更好地分析各参数与 ${\mathcal{R}}_0$ 的关系，我们先假设相关的参数服从正态分布，然后利用拉丁超立方抽样对参数空间进行抽样，最后计算PRCCs (偏秩相关系数)值来对模型(1)的 ${\mathcal{R}}_0$ 进行不确定性和敏感性分析并确定各参数对阈值的影响大小。若PRCCs的值为正，则表示该参数的增加会导致 ${\mathcal{R}}_0$ 的增加；若PRCCs的值为负，则表示该参数的增加会导致 ${\mathcal{R}}_0$ 的减小。对于模型(1.1)本文选取8个参数 $\beta$ 、 $\nu$ 、 $\theta$ 、 ${\gamma }_1$ 、 ${\gamma }_2$ 、 $\alpha$ 、 $\delta$ 、 $\rho$ 进行敏感性分析，分析结果如图6所示。从图6中可以看出，在这些参数中， $\beta$ 、 $\nu$ 、 $\theta$ 、 $\alpha$ 、 $\delta$ 、 $\rho$ 对 ${\mathcal{R}}_0$ 有正影响， ${\gamma }_1$ 、 ${\gamma }_2$ 对 ${\mathcal{R}}_0$ 有负影响。其中，与 ${\mathcal{R}}_0$ 有强相关性的参数为：无意识的易感者接触有症状感染者的有效传播率 $\beta$ 和无意识的易感者接触无症状感染者的调节因子 $\theta$ 。因此，在初始阶段，不受任何因素影响的情况下，对于减少和控制疾病传播：可以通过控制有症状感染者感染一个易感者的传播率 $\beta$ 使其降低，控制无症状感染者的接触率调节因子 $\theta$ 使其降低，即国家采取措施，执行严格的居家隔离政策，使无意识的易感人群尽量避免接触染病者。随着疾病的蔓延和扩散，媒体宣传、人类行为及社会反映等信息效应会对疾病的传播产生影响，人们会有一定的防护意识，因此，在建立及模型的基础上，改变和信息效应相关的参数 $m$ 和 $\omega$ ，观察累计确诊病例数的变化情况，具体的变化如图7和图8所示。我们发现，失去防护意识后有意识的易感者到无意识的易感者的转化率越大，意味着无意识易感者的人数增多，这些人群对于接触染病者没有一定的防护意识，那么被感染的概率就会增大，累计感染的人数就会大幅度上升；此外，对于受信息效应影响后无意识易感者变为有意识易感者的比率变大，进一步说明人们会采取各种防护措施来降低自己会被传染的可能性，从图8中，也呈现出了当 $m$ 变大时，累计确诊的病例数会减少。

Figure 6. Results of Partial Rank Correlation Coefficients (PRCCs)

图6. 部分相关系数(PRCCs)的结果

7. 结论

本文使用模型分析了信息效应对COVID-19传播的影响，我们发现信息效应在控制新冠肺炎疫情蔓延中发挥着重要的作用，可以让公众意识到对于疾病的控制应该做到预防为主，治疗为辅。而且还得到了如下结论：

1) 根据人群接触信息从而加强自身防范意识这一特点，有较强防护意识的易感人群更容易对传染病进行规避，会采取一定的措施，比如：减少外出降低与感染者的接触率、外出回去之后及时消毒、及时接种疫苗等。

Figure 7. The effect of $\omega$ on the number of cumulative confirmed cases

图7. $\omega$ 对累计确诊病例数的影响

Figure 8. The effect of $m$ on the number of cumulative confirmed cases

图8. $m$ 对累计确诊病例数的影响

2) 根据95%置信区间的图，发现所拟合出的累计确诊病例数不仅是在合理的范围而且误差是很小的，说明该模型能较好地反映该疾病的趋势。

3) 为了让有防护意识的人群逐渐增多，要通过网络、报纸、电视、广播等媒体平台对人们进行传染病危害的宣传教育，使他们对疾病做出积极有效的预防，这样做会降低感染者的比率。并且通过m的值对累计病例数的结果发现，如果增大受信息影响后无意识易感者变为有意识的比率，那么累计感染者的数量会有所下降。

4) 观察估计出来的参数，我们也可以发现，在受到信息效应影响后，有症状感染者的比例和人们的因病死亡率都很低，说明人们的自我保护意识很强，即使无意识的易感人群被感染之后，病情也不是很严重，通过一般的治疗都是可以恢复的。

最后，为了让新型冠状病毒得到及时的控制，阻止蔓延。政府要加大信息宣传的力度，让人们及时获得关于疾病的发展动态，进而从个人行为上发生一定的变化，这样做既能减少感染者的数量，而且还会降低死亡的风险。通过图4的拟合和验证图，发现该模型拟合预测的数据和实际数据是比较吻合的，说明在未来的新发传染病的控制中，该模型可以为其提供一定的理论参考。

该研究还有很多改进之处：第一，由于信息获取的时间可能滞后，可以考虑加入时滞后进行分析；第二，对于新型冠状病毒的病发和个人体质和免疫力有关，因此可以考虑加入年龄结构分析。