COVID-19分数阶群体时滞模型的稳定性分析

魏庆栋，王 忠，汤家军，陈海燕

（1.火箭军工程大学，西安 710038；2.中国矿业大学（北京），北京 100083）

自新冠疫情防控进入常态化以来，使用合适的数学模型来描述和预测病毒传播情况，评估疫情防控措施效果已经成为学界共识。本文通过研究新冠肺炎传播的动力学特性，引入时滞参数，结合分数阶模型，构建了一种描述新冠肺炎传播状况的分数阶群体时滞模型。

1 模型的建立与分析

1.1 分数阶模型的建立

根据COVID-19的传播特性，在传染病动力学SIR模型中引入潜伏期患者，即已经感染病毒但是尚未表现出症状的患者，也可以称为无症状患者，因此在模型中增加时滞参数，根据现有经验，在经过一段时间的潜伏后，会有一定比例的潜伏患者最终发展成为确诊病例[1]，这部分人群对最终结果的影响离不开潜伏周期这一因素，用参数τ来表示，在完成整数阶模型后，修改得到下列分数阶模型

上述系统当t=0时，会满足下述的初值条件

其中，S（t）为易感人群t时刻的人数，A（t）表示在t时刻无症状感染者的人数，I（t）为在t时刻有症状的感染者的人数，R（t）为恢复健康的人数，λ为正常人群中的自然增长速率下的增长人数，θ为无症状感染者转换为有症状感染者的概率，β1、β2分别为正常人群和无症状感染者与有症状感染者接触后被感染的概率，δA、δI分别为无症状被感染人群和有症状被感染人群的恢复率，d、dA、dI分别为正常人群、无症状感染者和有症状感染者的死亡率。在本系统中引入时滞参数τ来模拟模型的潜伏期参数。

1.2 模型的分析

根据系统1（公式（1）所示，下同），可以得出2个平衡点，一个无病平衡点E0=（x0，0，0，0），其中，一个地方病平衡点E*=（S*，A*，I*，R*），该平衡点满足

在对上述2个平衡点进行分析前，利用下一代生成矩阵法求出该模型的基本再生数为

1.2.1 平衡点的非负性分析

定理1.2.2对于系统1来说，地方病平衡点E*有以下结论，当R0＞1时，地方病平衡点E*=（S*，A*，I*，R*）是非负的。

证明：对地方病平衡点E*=（S*，A*，I*，R*）进行求解，得到

由于参数非负，因此S*≥0。对于A*，I*，R*，根据等式，可以得知主要取决于λ－dS*。当R0≥1时，同时因为

为证明分数阶模型平衡点的稳定性可使用如下引理：

引理1[2]分数阶系统的平衡点满足局部渐进稳定的性质，当且仅当在相应平衡点处的Jacobi矩阵J=的所有特征值λi满足条件。

1.2.2 无病平衡点的局部渐进稳定性

设特征值为sα，则模型1（公式（1）所示，下同）在无病平衡点E0处的特征行列式为

定理1.2.3模型1满足下列结论，当τ≥0时，该模型的无病平衡点E0有如下结论：

如果R0＜1，则E0是局部渐进稳定的。如果R0≥1，则E0是不稳定的。

证明：根据特征矩阵（3），进行推导化简可以得到

根据（sα+d）2，得到特征值因此。

将行列式部分用下列形式进行表示

P1（s）=0存在2个根，一个根为一定有负实部；另一个根，当R0＜1时因此此时它也一定有负实部。

1.2.3 地方病平衡点的局部渐进稳定性

当R0＞1时，模型1存在一个非负的地方病平衡点E*=（S*，A*，I*，R*），设特征值为sα，则地方病平衡点E*处的特征行列式为

先讨论τ=0时，令γ=sα，得到

其中参数d1，d2，d3可以表示为

而m11，m12，m13，m21，m22，m32，m33可以表示为

因此当τ=0时，E*有如下性质：当R0＞1时，如果式（8）满足n=3时的霍尔维茨准则，即式（8）的所有根都具有负实部，则E*是局部渐进稳定的。

当τ≠0时，根据特征行列式（6），可以化解得到如下的等式

将虚部和实部分开得到

其中参数为

对式（10）经过推导化简，令w2α=y，得到了如下的结果

进一步简化为

其中各参数表示意义如下

因此当式（13）满足霍尔维茨准则时，式（13）的根均为负值，则式（10）的解中存在负实部，此时地方病平衡点是局部渐进稳定的，综上给出如下定理。

定理1.2.4对于模型1的地方病平衡点E*，在R0＞1的情况下，若式（8）与式（13）均满足引理1，则可以得出结论，当τ≥0时，E*是局部渐进稳定的。

2 数值模拟与结果

为便于后续计算，首先通过NSFD算法对分数阶模型进行离散化处理[3]，然后进行数据模拟，主要模拟各类人群随时间变化的发展趋势。以此来研究带有时滞的分数阶模型的平衡点的稳定性。分2种情况对模型进行数值模拟，一是对模型（1）在不同的阶数下随时间变化的情况进行模拟，二是对模型在同一分数阶下，取不同时滞时，随时间变化的情况进行模拟。

2.1 对无病平衡点的数值模拟

由于对无病平衡点进行分析时要使得不存在非负的地方病平衡点，因此参数设置为λ=5，d=0.04，为增加结果的可靠性，同时设置了几组不同的阶数进行对比试验，将分数阶阶数分别设置为α=0.7，0.8，0.9，其余参数分别为β1=β2=0.001，d1=0.03，d2=0.1，δA=0.4，δI=0.3，θ=0.6，τ＝10。根据参数设置，此时基本再生数R0≈0.3＜1，且λ-dS*＜0，因此显然不存在非负的地方病平衡点，且E0=（125，0，0，0），根据定理1.2.3，在满足τ≥10的情况下，当R0＜1时，无病平衡点是局部渐进稳定的。利用MATLAB进行数值模拟，结果显示与理论相符合。因此按照这个趋势发展，各部分人群都将趋于稳定，其中只有易感人群还是非零的，其余的都将趋于零。

其余参数不变，将阶数设置为0.9，将时滞分别设定为τ=5，10，30，此时模拟结果各数值都趋于稳定，与定理1.2.3相符合。

2.2 对地方病平衡点的数值模拟

在对地方病平衡点进行模拟时，需要对参数重新进行修正，选择参数λ=10，d=0.02，分数阶阶数分别为α=0.7，0.8，0.9，其余参数保持不变，最终数值模拟的结果如图1所示。

图1 地方病平衡点下模型各变量在取不同阶数时随时间变化的图

此时基本再生数R0=1.16＞1，显然存在非负的地方病平衡点E*=（172，15，23，651）。将假设参数带入式（13）后，根据定理1.2.4，此时地方病平衡点是局部渐进稳定的。显然数值模拟的结果与理论相符合，在该情况下，随着时间的发展，各部分人群都将趋于稳定，但是依然存在患者和无症状的潜伏者，此时疫情防控依然不能松懈。

其余参数不变，将阶数设置为0.9，将时滞分别设定为τ=5，10，30，数值模拟结果如图2所示，此时，τ≥0，地方病平衡点都是局部渐进稳定的。模拟的结果与理论相符合。

图2 地方病平衡点下模型各变量在取不同时滞时随时间变化的图

3 结论

根据构建的分数阶模型，得出了其有2个平衡点，对2个平衡点的局部渐进稳定性分别进行了讨论。通过计算求出了基本再生数R0，证明了在τ≥0的情况下，无病平衡点E0在R0＜1时，是局部渐进稳定的，而当R0≥1时，无病平衡点E0不是局部渐进稳定的；对于地方病平衡点E*来说，当R0＞1时，若式（8）与式（13）中参数均满足引理1，可以得出结论，当τ≥0时，E*是局部渐进稳定的。同时根据2种情况对模型进行了模拟，一种是分数阶阶数不同，另一种是保持分数阶阶数相同但取不同的时滞，通过模拟验证了之前推导的平衡点稳定性的结论是正确的。此外，通过分析模拟，证明了基本再生数R0这一指标的重要性，为了更好地防控疫情，需要采取措施力求将R0保持在1以下。