高维传染病模型的流形分析

Manifold Analysis of Higher-Dimensional Epidemiological Model

摘要: 本文主要研究高维传染病模型如SEIR模型，其中在某种条件下可获得一维中心流形的存在性结论。

Abstract: This paper focuses on high-dimensional Epidemiological models such as the SEIR model, in which conclusions about the existence of a one-dimensional central manifold can be obtained under certain conditions.

1. 引言

传染病动力学通过构建数学模型，定量解析疾病传播机制与种群动态的相互作用，其核心在于结合社会因素与生物学特性建立动态系统。该方法通过微分方程描述易感者(S)、感染者(I)、康复者(R)等群体的状态转换，实现对疫情发展过程的动态推演。经典SIR模型经扩展后形成包含潜伏期(E)的SEIR模型，并进一步融入空间扩散、社会网络传播等复杂维度[1] [2]。

在经典传染病模型的构建中，通常假定个体传播强度与易感者数量呈线性正比关系。然而，该理论框架在特定边界条件下显现出明显的局限性——当易感者密度突破临界阈值时，模型将推导出单位时间内单个感染者可实现无限次传播的悖论。这与现实传播相矛盾，而非线性传染率$a{I}_1^u{S}_1^v$ 的饱和特性弥补了这一条件，通过量化个体传播能力的边际递减效应来完善模型的生物现实性。

2. 模型分析

具有非线性传染率的高维传染病模型如下：

$\begin{array}{l}{S^{\prime }}_1=-a{I}_1^u{S}_1^v+b-\mu S_1\\ {E^{\prime }}_1=a{I}_1^u{S}_1^v-\left(\beta +\mu \right)E_1\\ {I^{\prime }}_1=\beta E_1-\left(\gamma +\mu \right)I_1\\ {R^{\prime }}_1=\gamma I_1-\mu R_1\end{array}$ (2.1)

其中，参数均为正数，需要进行分析的总人口记为$N_1$ ，可以划分为易感者、潜伏者、染病者和移除者，分别表示为$S_1$ 、$E_1$ 、$I_1$ 、$R_1$ ，即$N_1=S_1+E_1+I_1+R_1$ 。非负常数$\mu$ 表示死亡率，$b$ 为出生率，$\epsilon$ 表示潜伏者到染病者的转化率，$\gamma$ 为康复率。

借助无量纲变换(假设出生率与死亡率相等情况下)

$t=\left(\gamma +\mu \right)t_1$ ，$m=\frac{\mu }{\gamma +\mu }$ ，$n=\frac{\beta }{\gamma +\mu }$ ，$p=\frac{a}{\gamma +\mu }$

注意$n>1$ 时，$\beta >\gamma +\mu$ ；$n<1$ 时，$\beta <\gamma +\mu$ ；并且$p<1$ 。

系统(2.1)可以等价于：

$\begin{array}{l}{S^{\prime }}_2=-p{I}_2^u{S}_2^v+m-m{S}_2\\ {E^{\prime }}_2=p{I}_2^u{S}_2^v-\left(m+n\right)E_2\\ {I^{\prime }}_2=n{E}_2-I_2\end{array}$ (2.2)

取值$v=1$ ，$u=2$ 系统则为：

$\begin{array}{l}{S^{\prime }}_3=-p{I}_3^2{S}_3+m-m{S}_3\\ {E^{\prime }}_3=p{I}_3^2{S}_3-\left(m+n\right)E_3\\ {I^{\prime }}_3=n{E}_3-I_3\end{array}$ (2.3)

由计算可知当${\left(\frac{mn}{m+n}\right)}^2-\frac{4m}{p}=0$ ，系统(2.1)存在一个疾病消除平衡点$P_0\left(1,0,0\right)$ 及一个地方病平衡点$P^{\ast }\left(\frac{1}{2},\frac{m}{2\left(m+n\right)},\frac{mn}{2\left(m+n\right)}\right)$ ，且疾病消除平衡点始终为局部渐近稳定的。

3. 一维中心流形存在性

令$S_3^{\ast }=S_3-\frac{1}{2}$ ，$E_3^{\ast }=E_3-\frac{m}{2\left(m+n\right)}$ ，$I_3^{\ast }=I_3-\frac{mn}{2\left(m+n\right)}$ ，则可得

$\begin{array}{l}{S}_3^{\ast }{}^{\prime }=-p{I}_3^{\ast }{}^2{S}_3^{\ast }-\frac{p}{2}{I}_3^{\ast }{}^2-\frac{pmn}{m+n}I{S}_3^{\ast }-\frac{pmn}{2\left(m+n\right)}{I}_3^{\ast }-2m{S}_3^{\ast }\\ E_3^{\ast }{}^{\prime }=p{I}_3^{\ast }{}^2{S}_3^{\ast }+\frac{p}{2}{I}_3^{\ast }{}^2+\frac{p\alpha \beta }{m+n}{I}_3^{\ast }{S}_3^{\ast }+\frac{p\alpha \beta }{2\left(m+n\right)}{I}_3^{\ast }+m{S}_3^{\ast }-\left(m+n\right)E_3^{\ast }\\ I_3^{\ast }{}^{\prime }=n{E}_3^{\ast }-I_3^{\ast }\end{array}$ (3.1)

这样系统(2.3)的正平衡点$P^{\ast }\left(\frac{1}{2},\frac{m}{2\left(m+n\right)},\frac{mn}{2\left(m+n\right)}\right)$ 就转化为系统(3.1)的平衡点$\left(0,0,0\right)$ ；当参数满足$\frac{2m\left(m+n-1\right)}{m+n}>1$ 时，依据特征值情况以及中心流形定理[3]，可将向量场约化到一维中心流形的向量场上。

情形1 取$m=n=1$ ，由$\Delta ={\left(\frac{mn}{m+n}\right)}^2-\frac{4m}{p}=0$ 可得$p=16$ ，此时系统转化为：

$\begin{array}{l}{S}_3^{\ast }{}^{\prime }=-16I_3^{\ast }{}^2{S}_3^{\ast }-8I_3^{\ast }{}^2-8I_3^{\ast }{S}_3^{\ast }-4I_3^{\ast }-2S_3^{\ast }\\ E_3^{\ast }{}^{\prime }=16I_3^{\ast }{}^2{S}_3^{\ast }+8I_3^{\ast }{}^2+8I_3^{\ast }{S}_3^{\ast }+4I_3^{\ast }+S_3^{\ast }-2E_3^{\ast }\\ I_3^{\ast }{}^{\prime }=E_3^{\ast }-I_3^{\ast }\end{array}$ (3.2)

由特征方程

${\lambda }_{}^{\ast }{}^2+\left(3m+n+1\right){\lambda }_{}^{\ast }+\left[\left(2m-1\right)\left(m+n\right)-2m\right]=0$

得出特征根为：

${\lambda }_1^{\ast }=0,{\lambda }_2^{\ast }=-1<0,{\lambda }_3^{\ast }=-4<0$ ，

利用坐标变换

$\left(\begin{array}{c}{S}_3^{\ast }\\ E_3^{\ast }\\ I_3^{\ast }\end{array}\right)=\left(\begin{array}{ccc}-2& -4& 2\\ 1& 0& -3\\ 1& 1& 1\end{array}\right)\left(\begin{array}{c}{\xi }^{\ast }\\ {\eta }^{\ast }\\ {\varsigma }^{\ast }\end{array}\right)$

或者$\left(\begin{array}{c}{\xi }^{\ast }\\ {\eta }^{\ast }\\ {\varsigma }^{\ast }\end{array}\right)=\left(\begin{array}{ccc}\frac{1}{4}& \frac{1}{2}& 1\\ -\frac{1}{3}& -\frac{1}{3}& -\frac{1}{3}\\ \frac{1}{12}& -\frac{1}{6}& \frac{1}{3}\end{array}\right)\left(\begin{array}{c}{S}_3^{\ast }\\ E_3^{\ast }\\ I_3^{\ast }\end{array}\right)$

系统(3.2)转化为可应用中心流形的形式[4]：

$\left(\begin{array}{c}{\xi }^{\ast }{}^{\prime }\\ {\eta }^{\ast }{}^{\prime }\\ {\varsigma }^{\ast }{}^{\prime }\end{array}\right)=\left(\begin{array}{ccc}0& 0& 0\\ 0& -1& 0\\ 0& 0& -4\end{array}\right)\left(\begin{array}{c}{\xi }^{\ast }\\ {\eta }^{\ast }\\ {\varsigma }^{\ast }\end{array}\right)+\left(\begin{array}{c}g\\ 0\\ -g\end{array}\right)$ ，

其中

$g=4{\left({\xi }^{\ast }+{\eta }^{\ast }+{\varsigma }^{\ast }\right)}^2\left(2{\varsigma }^{\ast }-2{\xi }^{\ast }-4{\eta }^{\ast }\right)-2\left({\xi }^{\ast }+{\eta }^{\ast }+{\varsigma }^{\ast }\right)\left({\xi }^{\ast }+3{\eta }^{\ast }-3{\varsigma }^{\ast }\right)$ ，

对应的中心流形为：${\eta }^{\ast }=g_1\left({\xi }^{\ast }\right),{\varsigma }^{\ast }=g_2\left({\xi }^{\ast }\right)$ ，利用中心流形近似计算方程

$D_{x}h\left(x,\epsilon \right)\left[Ax+f\left(x,h\left(x,\epsilon \right),\epsilon \right)\right]-Bh\left(x,\epsilon \right)-g\left(x,h\left(x,\epsilon \right),\epsilon \right)=0$

其中

$x\equiv S_3^{\ast },y\equiv \left(E_3^{\ast },I_3^{\ast }\right),h=\left(g_1,g_2\right),A=0$

$B=\left(\begin{array}{cc}-1& 0\\ 0& -4\end{array}\right)$ ，

$\begin{array}{l}{g}_1=4{\left({\xi }^{\ast }+{\eta }^{\ast }+{\varsigma }^{\ast }\right)}^2\left(2{\varsigma }^{\ast }-2{\xi }^{\ast }-4{\eta }^{\ast }\right)-2\left({\xi }^{\ast }+{\eta }^{\ast }+{\varsigma }^{\ast }\right)\left({\xi }^{\ast }+3{\eta }^{\ast }-3{\varsigma }^{\ast }\right)\\ g_2=-4{\left({\xi }^{\ast }+{\eta }^{\ast }+{\varsigma }^{\ast }\right)}^2\left(2{\varsigma }^{\ast }-2{\xi }^{\ast }-4{\eta }^{\ast }\right)+2\left({\xi }^{\ast }+{\eta }^{\ast }+{\varsigma }^{\ast }\right)\left({\xi }^{\ast }+3{\eta }^{\ast }-3{\varsigma }^{\ast }\right)\end{array}$

令

$\begin{array}{l}{\eta }^{\ast }=f_1\left({\xi }^{\ast }\right)={\rho }_1{\xi }^{\ast }{}^2+{\rho }_2{\xi }^{\ast }{}^3+{\rho }_3{\xi }^{\ast }{}^4+\cdots \\ {\varsigma }^{\ast }=f_2\left({\xi }^{\ast }\right)={\kappa }_1{\xi }^{\ast }{}^2+{\kappa }_2{\xi }^{\ast }{}^3+{\kappa }_3{\xi }^{\ast }{}^4+\cdots \end{array}$

计算可得：

${\rho }_1=0,{\rho }_2=0,{\kappa }_1=\frac{1}{2},{\kappa }_2=2$ ，

那么

$\begin{array}{l}{\eta }^{\ast }=f_1\left({\xi }^{\ast }\right)=0\\ {\varsigma }^{\ast }=f_2\left({\xi }^{\ast }\right)=\frac{1}{2}{\xi }^{\ast }{}^2+2{\xi }^{\ast }{}^3+\cdots \end{array}$

系统(3.2)在中心流形上的降阶形式可以转化为：

$\begin{array}{c}{\xi }^{\ast }{}^{\prime }=A{\xi }^{\ast }+f\left({\xi }^{\ast },h\left({\xi }^{\ast }\right)\right)\\ =2\left({\xi }^{\ast }+{\varsigma }^{\ast }\right)\left[4\left({\varsigma }^{\ast }{}^2-{\xi }^{\ast }{}^2\right)-{\xi }^{\ast }+3{\varsigma }^{\ast }\right]\\ =2\left({\xi }^{\ast }+\frac{1}{2}{\xi }^{\ast }{}^2+2{\xi }^{\ast }{}^3+\cdots \right)\left[4{\left(2{\xi }^{\ast }{}^2+4{\xi }^{\ast }{}^3+\cdots \right)}^2-4{\xi }^{\ast }{}^2-{\xi }^{\ast }+\left(\frac{3}{2}{\xi }^{\ast }{}^2+6{\xi }^{\ast }{}^3+\cdots \right)\right]\\ =-2{\xi }^{\ast }{}^2-6{\xi }^{\ast }{}^3+\cdots \end{array}$

此时，零解为不稳定的。

情形2 取$m=0$ ，系统(3.1)可转化为：

$\begin{array}{l}{S}_3^{\ast }{}^{\prime }=-p{I}_3^{\ast }{}^2{S}_3^{\ast }-\frac{p}{2}{I}_3^{\ast }{}^2\\ E_3^{\ast }{}^{\prime }=p{I}_3^{\ast }{}^2{S}_3^{\ast }+\frac{p}{2}{I}_3^{\ast }{}^2-n{E}_3^{\ast }\\ I_3^{\ast }{}^{\prime }=n{E}_3^{\ast }-I_3^{\ast }\end{array}$ (3.3)

转化为可应用中心流形的形式：

$\left(\begin{array}{c}{S}_3^{\ast }{}^{\prime }\\ E_3^{\ast }{}^{\prime }\\ I_3^{\ast }{}^{\prime }\end{array}\right)=\left(\begin{array}{ccc}0& 0& 0\\ 0& -n& 0\\ 0& n& -1\end{array}\right)\left(\begin{array}{c}{S}_3^{\ast }\\ E_3^{\ast }\\ I_3^{\ast }\end{array}\right)+p\left(\begin{array}{c}-I_3^{\ast }{}^2{S}_3^{\ast }-\frac{1}{2}{I}_3^{\ast }{}^2\\ I_3^{\ast }{}^2{S}_3^{\ast }+\frac{1}{2}{I}_3^{\ast }{}^2\\ 0\end{array}\right)$ ，$m^{\prime }=0$

中心流形为：$E_3^{\ast }=f_1\left(S_3^{\ast },p\right),I_3^{\ast }=f_2\left(S_3^{\ast },p\right)$ ，利用近似计算方程，

其中

$x\equiv S_3^{\ast },y\equiv \left(E_3^{\ast },I_3^{\ast }\right)\epsilon =p,h=\left(f_1,f_2\right),A=0$ ，

$B=\left(\begin{array}{cc}-n& 0\\ n& -1\end{array}\right)$ ，

$\begin{array}{l}f\left(x,h\left(x,\epsilon \right),\epsilon \right)=-p{I}_3^{\ast }{}^2{S}_3^{\ast }-\frac{p}{2}{I}_3^{\ast }{}^2\\ g\left(x,h\left(x,\epsilon \right),\epsilon \right)=\left(\begin{array}{c}p{I}_3^{\ast }{}^2{S}_3^{\ast }+\frac{p}{2}{I}_3^{\ast }{}^2\\ 0\end{array}\right)\end{array}$

令

$\begin{array}{l}{f}_1\left(S_3^{\ast },p\right)={\rho }_1{S}_3^{\ast }{}^2+{\rho }_2{S}_3^{\ast }p+{\rho }_3{p}^2+\cdots \\ f_2\left(S_3^{\ast },p\right)={\kappa }_1{S}_3^{\ast }{}^2+{\kappa }_2{S}_3^{\ast }p+{\kappa }_3{p}^2+\cdots \end{array}$

计算可得

${\rho }_1=0,{\rho }_2=0,{\kappa }_1=0,{\kappa }_2=0$

可得到一低维控制方程

${\xi }^{\ast }{}^{\prime }=A{\xi }^{\ast }+f\left({\xi }^{\ast },h\left({\xi }^{\ast }\right)\right)$

对系统(3.3)可知：$I_3^{\ast }=0$ 且$E_3^{\ast }=0$ 即$S$ 轴是不变流形，意思是$S$ 轴上每个点都是吸引点。

4. 结论

传染病动力学是一种基于数理建模与系统分析的疾病传播研究方法。该方法通过整合病原体扩散特征、宿主群体动态变化及社会环境变量等多维参数，构建包含微分方程、网络拓扑和随机过程的理论框架，重点解析疫病传播的内在机理与外部驱动因素。

本文主要考虑具有非线性传染率的高维传染病模型，进行定性分析[5]。在参数满足相应条件下，一维中心流形的存在性得以说明，类似的方法也适用于其他传染病模型的全局分析。

当参数在相应范围内变动的时候，系统(2.1)不仅存在边界平衡点还存在一正平衡点，边界平衡点始终为局部渐近稳定的。而当正平衡点满足$\frac{2m\left(m+n-1\right)}{m+n}>1$ 的条件时，系统存在一维中心流形，平衡点不稳定；此外，参数范围内的其他情况，诸如二维中心流形、利用正向不变集分析轨线走向等也将会继续探讨。

这种系统动力学视角不仅可突破横断面数据的局限性，更能通过相空间分析预判流行态势的长期演化轨迹，为应对新发突发传染病提供具有前瞻性的科学依据。对于次数达到四阶及以上的高维系统，由于动力性态非常复杂，具有一定的研究价值；关于高维系统的同宿异宿环问题以及高维系统周期轨道的唯一性问题及分支问题还有待于进一步探讨和研究。

基金项目

本文由2024大学生创新创业训练项目(X202410452210)以及省级青年教师教学示范课程支持。

NOTES

^*通讯作者。