基于Holling II型功能反应的疟原虫感染动力学建模与稳定性分析

Dynamic Modeling and Stability Analysis of Plasmodium Infection Based on Holling Type II Functional Response

摘要: 疟疾是由按蚊传播的典型媒介传染病，一直以来是全球公共卫生领域的重大挑战。本文基于生物数学建模理论，构建了带有Holling第II型功能反应函数模型，描述红细胞–染虫红细胞–免疫效应(NETs)三者的动力学行为。我们定性分析了当$R_0\le 1$ 时，无病平衡点$E_0$ 全局渐近稳定；当$\text{1}<R_0<1+\frac{\eta }{\sigma }\left(\alpha +\frac{\beta }{d}\right)$ 时，边界平衡点$E_1$ 局部渐近稳定；当$R_{\text{0}}>1+\frac{\eta }{\sigma }\left(\alpha +\frac{\beta }{d}\right)$ 时，正平衡点$E^{\ast }$ 局部渐近稳定。本文研究在一定程度上为疟疾防控提供理论支持和实践指导。

Abstract: Malaria, a vector-borne disease transmitted by Anopheles mosquitoes, remains a major global public health concern. Based on the theory of biomathematical modeling, in this study, we develop a dynamic model incorporating Holling type II functional response to characterize interactions among erythrocytes, malaria-infected erythrocytes, and neutrophil extracellular traps (NETs) as immune effectors. We qualitatively analyzed that when $R_0\le 1$ , the disease-free equilibrium point $E_0$ is globally asymptotically stable; when $\text{1}<R_0<1+\frac{\eta }{\sigma }\left(\alpha +\frac{\beta }{d}\right)$ , the boundary equilibrium point $E_1$ is locally asymptotically stable; when $R_{\text{0}}>1+\frac{\eta }{\sigma }\left(\alpha +\frac{\beta }{d}\right)$ , the positive equilibrium point $E^{\ast }$ is locally asymptotically stable. This study aims to provide theoretical support and practical guidance for malaria control.

1. 引言

疟疾是经过雌性按蚊叮咬或输血传播的一类严重危害人类健康的血液寄生虫疾病，它由疟原虫感染引起的，世界卫生组织(WHO)将它列为世界三大传染病之一，在全球范围内造成了巨大的健康负担，至今仍然造成每年近2.5亿感染和60万死亡病例，且输入性疟疾在中国部分地区仍呈上升趋势。长期以来一直是全球公共卫生领域面临的一大挑战。虽然前人在生物数学领域里做了许多研究工作和建立了诸多模型，但随着研究的深入，以及对疟原虫的生物学特性和致病机制有了更深的认识，发现仍有许多关键问题亟待解决。疟原虫在宿主内的生存和繁殖过程十分复杂，涉及多种代谢途径和与宿主细胞的相互作用，想要全面揭示疟原虫致病的本质，还需要进一步探索疟原虫与宿主相互作用的动力学模式，即构建合适的疟原虫致病模型并研究其动力学行为。

在本研究中，我们聚焦的是疟原虫侵染人体内红细胞的行为来建立相应的疟原虫致病数学模型，并深入研究模型的平衡点存在性和稳定性等动力学行为。通过对模型的分析，我们期望能够揭示疟原虫与宿主相互作用过程中关键因子和关键环节的动力学特征，显示TatD-like DNase是疟疾寄生虫生存于人体内的一种手段，为新型抗疟药物的产生以及疟疾疫苗的研制提供新的视角和策略。

2. 模型的提出

人类感染疟疾后，疟原虫会为了防卫释放一种DNA酶(TatD-like DNase)降解中性粒细胞胞外的NETs去抵抗免疫清除。2016年，文献[1]研究的是恶性疟原虫抵抗宿主天然免疫机理方面的重要工作。因此，定义$x\left(t\right)$ 表示正常红细胞的数量，$y\left(t\right)$ 表示感染红细胞的数量，$z\left(t\right)$ 表示免疫效应NETs的数量，我们考虑宿主体内三者之间的相互关系，建立含功能性反应的动力学模型[2]：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=\lambda -dx-\frac{\beta xy}{1+\alpha y},\hfill \\ \frac{\text{d}y}{\text{d}t}=\frac{\beta xy}{1+\alpha y}-ry-\delta yz,\hfill \\ \frac{\text{d}z}{\text{d}t}=\sigma yz-\eta z.\hfill \end{array}$ (1)

其中，$\lambda$ 表示正常红细胞的生成速率，$d$ 表示正常红细胞的死亡率，$\beta$ 表示疟原虫对红细胞的侵染率，$r$ 表示感染红细胞的死亡率，$\delta$ 表示NETs捕杀感染红细胞的能力，$\sigma$ 表示NETs的生成率，$\eta$ 表示Nets

的分解率，新加入的参数$\alpha$ ，$\frac{\text{1}}{\alpha }$ 表示的是一个饱和常数，模拟感染红细胞的生长速度是其最大速率的一半。

3. 动力学分析

3.1. 解的非负性与有界性

引理1 (非负性)：在初始条件$x\left(0\right)>0,y\left(0\right)\ge 0$ 和$z\left(0\right)\ge 0$ ，系统(1)的解是非负的。

证明：反证法，假设存在最小时间$t_0>0$ ，有$x\left(t_0\right)=0$ 。把$x\left(t_0\right)$ 代入(1)中的一式有$\frac{\text{d}x}{\text{d}t}=\lambda \ge 0$ ，则

对$\forall \epsilon >0$ ，当$t\in \left(t_0-\epsilon ,t_0\right)$ 时，有$x\left(t\right)<0$ ，由于$t_0$ 是$x\left(t_0\right)=0$ 的最小时间，且$x\left(0\right)>0$ 。由解对初值的连续依赖性，则必有$x\left(t\right)>0$ 。这与之前的假设矛盾，故对$\forall t>0$ ，都有$x\left(t\right)>0$ 。同理可证，$y\left(t\right)\ge 0$ ，$z\left(t\right)\ge 0$ 。

引理2 (有界性)：系统(1)的解是有界的。

证明：由(1)有

$\frac{\text{d}\left(x+y+\frac{\delta }{\sigma }z\right)}{\text{d}t}=\frac{\text{d}x}{\text{d}t}+\frac{\text{d}y}{\text{d}t}+\frac{\delta }{\sigma }\frac{\text{d}z}{\text{d}t}=\lambda -dx-ry-\frac{\delta \eta }{\sigma }z\le \lambda -m\left(x+y+\frac{\delta }{\sigma }z\right),$

其中$m=\max \left\{d,r,\eta \right\}$ 。令$v=x+y+\frac{\delta }{\sigma }z$ ，则$\frac{\text{d}v}{\text{d}t}\le \lambda -mv$ ，从而有

$\underset{t\to \infty }{\lim }\sup v\left(t\right)=\underset{t\to \infty }{\lim }\sup \left(x\left(t\right)+y\left(t\right)+\frac{\delta }{\sigma }z\left(t\right)\right)\le \frac{\lambda }{m},$

故(1)的解有界。

由引理1，2知(1)的正不变集为$\Gamma =\left\{\left(x,y,z\right)|x+y+\frac{\delta }{\sigma }z\le \frac{\lambda }{m},x>0,y\ge 0,z\ge 0\right\}$ 。

3.2. 平衡点的稳定性与基本再生数

令系统(1)中等号左边式子等于0，有

$\{\begin{array}{l}\lambda -dx-\frac{\beta xy}{1+\alpha y}=0\\ \frac{\beta xy}{1+\alpha y}-ry-\delta yz=0\\ \sigma yz-\eta z=0\end{array}$ (2)

我们容易求出有无病平衡点$E_{\text{0}}=\left(\frac{\lambda }{d},0,0\right)$ 。对于(1)，基本再生数是描述了一个感染细胞在其生命周期内产生的新感染细胞的平均数量，定义正常红细胞的基本再生数为

$R_{\text{0}}=\frac{\lambda \beta }{dr}.$

引理3 以下是成立的。

1) 当$R_{\text{0}}\le \text{1}$ 时，系统(1)存在有无病平衡点$E_{\text{0}}$ ；

2) 当$\text{1}<R_0<1+\frac{\eta }{\sigma }\left(\alpha +\frac{\beta }{d}\right)$ 时，系统(1)存在有无病平衡点$E_{\text{0}}$ 和边界平衡点：

$E_1=\left(x_1,y_1,z_1\right)=\left(\frac{r\left(1+\alpha y_1\right)}{\beta },\frac{d\left(R_0-1\right)}{\beta +d\alpha },0\right);$

3) 当$R_{\text{0}}>1+\frac{\eta }{\sigma }\left(\alpha +\frac{\beta }{d}\right)$ 时，系统(1)存在无病平衡点$E_{\text{0}}$ ，边界平衡点$E_1$ 和正平衡点：

$E^{\ast }=\left(x^{\ast },y^{\ast },z^{\ast }\right)=\left(\frac{\lambda }{d+\frac{\beta y^{\ast }}{1+\alpha y^{\ast }}},\frac{\eta }{\sigma },\frac{r}{\delta }\left(\frac{R_0}{1+\alpha \frac{\eta }{\sigma }+\frac{\beta }{d}\cdot \frac{\eta }{\sigma }}-1\right)\right).$

证明：我们容易证明1，2成立。

现证3，由(2)中的三式，有$y^{\ast }=\frac{\eta }{\sigma }$ ，把(2)整理有

$\{\begin{array}{l}\lambda -d{x}^{\ast }-\frac{\beta x^{\ast }{y}^{\ast }}{1+\alpha y^{\ast }}=0,\\ \left(\frac{\beta x^{\ast }}{1+\alpha y^{\ast }}-r-\delta z^{\ast }\right)y^{\ast }=0,\end{array}$

从而有

$x^{\ast }=\frac{\lambda }{d+\frac{\beta y^{\ast }}{1+\alpha y^{\ast }}},$

$z^{\ast }=\frac{1}{\delta }\left(\frac{\beta x^{\ast }}{1+\alpha y^{\ast }}-r\right)=\frac{1}{\delta }\left(\frac{\frac{\lambda \beta }{d+\frac{\beta y^{\ast }}{1+\alpha y^{\ast }}}}{1+\alpha y^{\ast }}-r\right)=\frac{1}{\delta }\left(\frac{\lambda \beta }{d\left(1+\alpha y^{\ast }\right)+\beta y^{\ast }}-r\right)=\frac{r}{\delta }\left(\frac{R_0}{1+\left(\alpha +\frac{\beta }{d}\right)y^{\ast }}-1\right).$

定理4 当$R_0\le 1$ 时，系统(1)无病平衡点$E_0$ 是全局渐近稳定的。

证明：定义Lyapunov函数

$V_0=x-x_0-x_0\ln \frac{x}{x_0}+y+\frac{\delta }{\sigma }z$

对Lyapunov函数沿系统(1)的解求全导数，有

$\frac{\text{d}{V}_0}{\text{d}t}=\left(1-\frac{x_0}{x}\right)\frac{\text{d}x}{\text{d}t}+\frac{\text{d}y}{\text{d}t}+\frac{\delta }{\sigma }\frac{\text{d}z}{\text{d}t}$

$\begin{array}{c}=\left(1-\frac{x_0}{x}\right)\left(\lambda -dx-\frac{\beta xy}{1+\alpha y}\right)+\left(\frac{\beta xy}{1+\alpha y}-ry-\delta yz\right)+\frac{\delta }{\sigma }\left(\sigma yz-\eta z\right)\\ =\left(1-\frac{x_0}{x}\right)\left(d{x}_0-dx\right)-\frac{\beta xy}{1+\alpha y}+\frac{\beta x_{0}y}{1+\alpha y}+\frac{\beta xy}{1+\alpha y}-ry-\frac{\delta \eta }{\sigma }z\\ =-\frac{d}{x}{\left(x-x_0\right)}^2+y\left(\frac{\beta x_0}{1+\alpha y}-r\right)-\frac{\delta \eta }{\sigma }z\\ \le -\frac{d}{x}{\left(x-x_0\right)}^2+y\left(\beta x_0-r\right)-\frac{\delta \eta }{\sigma }z\\ =-\frac{d}{x}{\left(x-x_0\right)}^2+ry\left(R_0-1\right)-\frac{\delta \eta }{\sigma }z.\end{array}$

由于$R_0\le 1$ ，有$\frac{\text{d}{V}_0}{\text{d}t}\le 0$ ，且$\frac{\text{d}{V}_0}{\text{d}t}=0$ 当且仅当$x=\frac{\lambda }{d},y=0,z=0$ ，故$\left\{\left(x,y,z\right)|\frac{\text{d}{V}_0}{\text{d}t}=0\right\}$ 的最大不变集只含无病平衡点$E_0$ 。由LaSalle不变集原理，可知，当$R_0\le 1$ 时，$E_0$ 全局稳定。

下面，我们考虑边界平衡点$E_{\text{1}}$ 与正平衡点$E^{\ast }$ 的稳定性，使系统(1)在点$E=\left(x,y,z\right)$ 处线性化，得雅可比矩阵，记作$J\left(E\right)$ ，有

$J\left(E\right)=\left(\begin{array}{ccc}-\alpha -\frac{\beta y}{1+\alpha y}& -\frac{\beta x}{{\left(1+\alpha y\right)}^2}& 0\\ \frac{\beta y}{1+\alpha y}& \frac{\beta x}{{\left(1+\alpha y\right)}^2}-r-\delta z& \delta y\\ 0& \sigma z& \sigma y-\eta \end{array}\right).$ (3)

定理5 当$\text{1}<R_0<1+\frac{\eta }{\sigma }\left(\alpha +\frac{\beta }{d}\right)$ 时，系统(1)边界平衡点$E_1$ 是局部渐近稳定的。

证明：由(3)可知

$J\left(E_{\text{1}}\right)=\left(\begin{array}{ccc}-d-\frac{\beta y_1}{1+\alpha y_1}& -\frac{\beta x_1}{{\left(1+\alpha y_1\right)}^2}& 0\\ \frac{\beta y_1}{1+\alpha y_1}& \frac{\beta x_1}{{\left(1+\alpha y_1\right)}^2}-r& \delta y_1\\ 0& 0& \sigma y_1-\eta \end{array}\right),$

则$J\left(E_{\text{1}}\right)$ 的特征方程为

$\left(\wedge -\sigma y_1+\eta \right)\left[\left(\wedge +d+\frac{\beta y_1}{1+\alpha y_1}\right)\left(\wedge -\frac{\beta x_1}{{\left(1+\alpha y_1\right)}^2}+r\right)+\frac{{\beta }^2{x}_1{y}_1}{{\left(1+\alpha y_1\right)}^3}\right]=0.$

显然$\wedge -\sigma y_1+\eta =\text{0}$ ，有

$\begin{array}{c}{\wedge }_{\text{1}}=\sigma y_1-\eta \\ =\frac{\sigma d}{\beta +d\alpha }\left[R_0-1-\frac{\eta }{\sigma \alpha }\left(\beta +d\alpha \right)\right]\\ =\frac{\sigma d}{\beta +d\alpha }\left[R_0-\left(1+\frac{\left(\frac{\beta }{d}+\alpha \right)\eta }{\sigma }\right)\right]<0,\end{array}$

故${\wedge }_{\text{1}}<0$ ，考虑

$\left(\wedge +d+\frac{\beta y_1}{1+\alpha y_1}\right)\left(\wedge -\frac{\beta x_1}{{\left(1+\alpha y_1\right)}^2}+r\right)+\frac{{\beta }^2{x}_1{y}_1}{{\left(1+\alpha y_1\right)}^3}=0$ ，

化简，有

${\wedge }^{\text{2}}+\left(-\frac{\beta x_1}{1+\alpha y_1}+r+d+\frac{\beta y_1}{1+\alpha y_1}\right)\wedge -\frac{d\beta x_1}{{\left(1+\alpha y_1\right)}^2}+r\left(d+\frac{\beta y_1}{1+\alpha y_1}\right)=0.$

令

$A=-\frac{\beta x_1}{1+\alpha y_1}+r+d+\frac{\beta y_1}{1+\alpha y_1}，B=-\frac{d\beta x_1}{{\left(1+\alpha y_1\right)}^2}+r\left(d+\frac{\beta y_1}{1+\alpha y_1}\right),$

由于A，B > 0，所以特征根具有负实部，从而$E_1$ 是局部渐近稳定的。

由(3)知

$J\left(E^{\ast }\right)=\left(\begin{array}{ccc}-d-\frac{\beta y^{\ast }}{1+\alpha y^{\ast }}& -\frac{\beta x^{\ast }}{{\left(1+\alpha y^{\ast }\right)}^2}& 0\\ \frac{\beta y^{\ast }}{1+\alpha y^{\ast }}& \frac{\beta x^{\ast }}{{\left(1+\alpha y^{\ast }\right)}^2}-r-\delta z^{\ast }& \delta y^{\ast }\\ 0& \sigma z^{\ast }& \sigma y^{\ast }-\eta \end{array}\right),$

则$J\left(E^{\ast }\right)$ 的特征方程是

${\wedge }^{\text{3}}+b_1{\wedge }^2+b_2\wedge +b_3=0,$ (4)

其中

$\begin{array}{l}{b}_1=d+\frac{\beta y^{\ast }}{1+\alpha y^{\ast }},\\ b_2=dr+\frac{\beta }{1+\alpha y^{\ast }}\left(r{y}^{\ast }-d{x}^{\ast }\right)+\left(\frac{\beta x^{\ast }}{{\left(1+\alpha y^{\ast }\right)}^2}-r\right)\left(\frac{\beta y^{\ast }}{1+\alpha y^{\ast }}+d+\eta \right),\\ b_3=\eta \left(d+\frac{\beta y^{\ast }}{1+\alpha y^{\ast }}\right)\left(\frac{\beta x^{\ast }}{{\left(1+\alpha y^{\ast }\right)}^2}-r\right).\end{array}$

由Routh-Hurwitz判据，(4)所有根具有负实部当且仅当

$b_1>0,b_1{b}_2-b_3>0,b_3\left(b_1{b}_2-b_3\right)>0.$ (5)

因此，我们有

定理6 当$R_{\text{0}}>1+\frac{\eta }{\sigma }\left(\alpha +\frac{\beta }{d}\right)$ 时，若系统(1)满足(5)，则正平衡点$E^{\ast }$ 是局部渐近稳定的。

4. 数值模拟

通过数学软件Matlab从一个量化的角度验证模型的理论结果，直观地展示模型的动力学行为。系统中各变量和参数的初值是参考文献[1] [3]-[6]，令$\lambda =4.15\times {10}^4$ ，$d=8.3\times {10}^{-3}$ ，$\beta =2\times {10}^{-9}$ ，$\alpha =5\times {10}^{-4}$ ，$r=1$ ，$\delta =1\times {10}^{-8}$ ，$\sigma ={2.5}^{-5}$ ，$\eta =0.05$ 。取定初始值$x_0=5\times {10}^6$ ，$y_0=1\times {10}^4$ ，$z_0=1\times {10}^{-4}$ 。

Figure 1. Dynamics of complete Plasmodium clearance mediated by the host immune system

图1. 宿主免疫系统完全清除疟原虫的动态过程

Figure 2. Dynamics of the critical steady state

图2. 临界稳态的动态过程

图1中，$R_0=0.01<1$ ，数值模拟结果看出无病平衡点$E_0=\left(5\times {10}^6,0,0\right)$ 是全局稳定的，表明当$R_0\le 1$ 时宿主感染疟疾后能恢复到健康状态，即消除了疟原虫。

图2中取$\beta =2.5\times {10}^{-7}$ ，$\sigma =5\times {10}^{-5}$ ，其他参数固定不变，算出$R_0=1.25$ ，$\frac{\eta }{\sigma }\left(\alpha +\frac{\beta }{d}\right)\approx 0.5301$ ，$\text{1}<R_0<1+\frac{\eta }{\sigma }\left(\alpha +\frac{\beta }{d}\right)$ ，边界平衡点$E_1\approx \left(4943181.81,471.59,0\right)$ 是稳定的。从图观察到NETs的数量显著

减少并最终走向0，即体内生成的免疫因子随即被疟原虫分泌的一种DNA酶(TatD-like DNase)降解，使得疟原虫成功在宿主体内建立感染。

图3中取$\beta =2.5\times {10}^{-7}$ ，$\sigma =5\times {10}^{-4}$ ，其他参数为固定值，有$R_0=1.25$ ，计算得到$\frac{\eta }{\sigma }\left(\alpha +\frac{\beta }{d}\right)\approx 0.053$ ，从而$R_0>1+\frac{\eta }{\sigma }\left(\alpha +\frac{\beta }{d}\right)$ ，如图3所示的数值模拟结果，结果表明存在稳定的正平衡点

$E^{\ast }\approx \left(4985697.94,100,18707093.82\right)$ 。该动力学特征揭示宿主免疫效应因子NETs未能实现疟原虫的完全清除，致使疟原虫种群与NETs浓度最终收敛于动态平衡状态。此现象表明疟原虫成功突破宿主免疫防御机制，在宿主体内形成感染。

Figure 3. Numerical simulation of Plasmodium overcoming host immune defenses and establishing infection

图3. 疟原虫突破宿主免疫防御并建立感染的数值模拟结果

5. 结论

疟疾作为一种蚊媒传染病，其病原体与宿主的相互作用机制对公共卫生安全具有重要影响。本研究构建了基于Holling第II类功能反应的ODE模型，系统刻画红内期感染过程中正常红细胞、感染红细胞与中性粒细胞胞外诱捕网(NETs)的动态互作关系。

通过对模型动力学行为的分析，我们揭示了正常红细胞的基本再生数$R_{\text{0}}=\frac{\lambda \beta }{dr}$ 对系统动力学特征的决定性作用：当$R_0\le 1$ 时，无病平衡点$E_0$ 是全局渐近稳定的；当$\text{1}<R_0<1+\frac{\eta }{\sigma }\left(\alpha +\frac{\beta }{d}\right)$ 时，边界平衡点$E_1$ 是局部渐近稳定的；当$R_{\text{0}}>1+\frac{\eta }{\sigma }\left(\alpha +\frac{\beta }{d}\right)$ 时，正平衡点$E^{\ast }$ 是局部渐近稳定的，这个模型在正平衡点

处并没有出现Hopf分支现象。理论结果一定程度上有助于揭示疟疾的发病机制，为开发新的治疗策略提供一定的理论依据。

结合流行病学数据和文献[7]指出，我国恶性疟疾分布已随防控措施完善显著缩小至桂、琼、滇局部地区，但奎宁类药物广泛使用引发的疟原虫抗药性正成为重要防控障碍。本模型通过研究病原体–宿主互作机制，为揭示疟疾致病机理提供了理论框架，但需指出当前研究侧重某方面的微观动力学过程，对蚊媒种群动态、气候变化等宏观环境因子及宿主个体差异以及其他的免疫机制如T细胞、B细胞介导的适应性免疫和细胞因子等的影响尚未充分纳入模型体系。后续研究将在多尺度建模与临床验证方面持续深入，以提升模型对复杂传播环境的解释力。

基金项目

湘南学院科学研究项目青年项目“蚊媒传染病防控的数学建模和动力学分析”(校发【2022】6号附件)。

NOTES

^*通讯作者。