考虑社会行为影响的离散传染病模型的建立及分析

Establishment and Analysis of a Discrete Infectious Disease Model Considering the Effects of Social Behavior

摘要

在突发传染病大流行期间， 如何在社会成本损失较小的前提下有效控制疾病传播， 是政府和公共卫生决策者面临的重大挑战。本文通过引入具有信息延迟的效用函数， 刻画了感染的预期成本和风险等社会行为， 建立了接触率受社会行为影响的离散传染病模型并对其进行了动力学性态分析。通过计算得到信息没有延迟时疾病的基本再生数， 给出无病平衡点及正平衡点的存在性条件。通过构造Lyapunov函数的方法证明了无病平衡点是全局渐近稳定的， 并利用Jury判据证明正平衡点是局部渐近稳定的。进一步对在不同时间延迟下的系统的性态进行了数值模拟， 发现时间延迟可能导致系统发生复杂的动力学行为， 包括从单调稳定到振荡稳定以及引发周期性流行病， 甚至在相当长的一段时间内出现多个峰值和激增， 这意味着时间延迟会导致系统性态不确定。

Abstract

During the outbreak of infectious diseases， how to effectively control the disease transmission and small social cost losses， is a major challenge for governments and public health decision-makers. In this paper， we introduced a utility function with information delay， which was described the social behaviors such as the expected cost and risk of infection， established and analyzed a discrete infectious disease model with contact rate that was affected by social behaviors. By calculating， the basic reproductive number of the model was obtained when the information was not delayed， and the existence conditions of the disease-free and positive equilibria were given. Using the method of constructing Lyapunov functions， we proved that the disease-free equilibrium was globally asymptotically stable， and the positive equilibrium was locally asymptotically stable by the Jury criterion. Furthermore， we simulated the behavior of system with the different time delay， and found that the system exhibits complex dynamic behaviors， such as monotonic stability， oscillation stability， periodicity， and even multiple peaks and spikes over a longer period of time with the increase of time delay， which means that time delay may lead to uncertainty of behaviors for the system.

Graphical abstract

关键词

时间延迟 / 离散模型 / 社会行为 / 效用函数 / 稳定性

Key words

time delay / discrete model / social behavior / utility function / stability

0 引　言

在传染病尤其是突发传染病爆发初期， 由于疾病新发且传染性强， 没有有效的药物和治疗方法， 感染人数会急剧增加。随着信息化的迅速发展， 人们通过媒体等方式对疾病有一定了解后， 会选择一些措施去规避风险， 感染人数就会逐渐减少。

传染病的流行会引发人们行为的改变， 行为的改变反过来又会影响传染病的传播， 其相互作用如图 1 所示。行为的改变 （如保持社交距离、 减少出行以及取消聚集性活动等）一方面可以降低疾病的传播， 使许多人免于染病， 并使疾病的基本再生数减小； 另一方面， 如果这些行为改变持续时间长， 虽然传染病得到很好的控制， 但可能会造成经济发展的降速或失业率的增加， 甚至出现社会不稳定因素^［1-2］。因此， 如何在控制疾病传播和社会成本之间找到平衡是值得考虑的问题。

在传染病流行期间， 个体的行为会随着人们感知到的风险而发生变化， 研究者们越来越关注行为变化在传染病预测和控制中所起的作用， 并开展了大量相关研究。Arthur等^［3］和Naidoo等^［4］通过埃博拉、 艾滋病等例子， 强调了将行为变化纳入传染病模型的重要性。Smith等^［5］和Brown等^［6］在一般传染病模型中纳入了保持社会距离等干预措施（实施干预措施会使人群接触率减小）。Fenichel等^［7］建立了一个考虑发生率为c\beta SI/N的传染病模型， 其中c是一个与染病者I有关的接触函数， 结果表明， 公共卫生措施的成本可能超过疾病本身带来的经济负担。此外， Panicker等^［8］不仅讨论了非药物策略（如保持社交距离、 戴口罩）对传染病系统动力学的影响， 还特别关注了时间延迟对实施这些策略的影响， 发现时间延迟可能会加剧传染病的流行， 并使系统出现振荡现象。时间延迟是由于在流行病爆发初期， 人们对传染病信息的获取、 对病毒的识别以及采取行为进行改变是需要时间的， 这会对传染病的传播产生显著影响。因此， 时间延迟也是研究疾病传播需要考虑的一个重要因素。

本文建立了受社会行为影响的离散传染病模型， 通过寻找一个最优接触率c^{\mathrm{*}}来实现在控制疾病传播的同时使成本损失最小。为此， 定义了一个反应接触率对疾病控制和预期成本影响的效用函数U(c)。传染病发生时采取防控措施会导致两方面的影响： 1） 增加政府或医疗机构的成本和个人经济负担， 显然预期成本与接触率c和感染者数量I成正比， 用函数h(c,I)表示； 2） 面对高传染性的突发疾病， 人们会感知风险并主动采取措施， 如减少流动或聚集， 可能导致社会经济发展受阻， 这种经济影响可以通过一个风险厌恶型效用函数g(c)（凹函数）来量化^［9］。结合这两方面的考虑， 可取目标函数为U(c)=h(c,I)+g(c)。另外， 疾病发生时不进行行为改变， 政府和医疗机构也会有一定的医疗成本， 记为{\alpha }_{\mathrm{0}}， 这样， 目标函数U(c)可改写为

U(c)={\alpha }_{\mathrm{0}}+h(c,I)+g(c)。

通过求解\mathrm{d}U/\mathrm{d}c=\mathrm{0}得到最优接触率c^{\mathrm{*}}。因此， 可将易感者S与染病者I接触时感染的概率\beta表示为{\beta }_{\mathrm{0}}（即进行行为改变后， 个体间通过接触被感染的概率）与c^{\mathrm{*}}的组合， 即

\beta ={\beta }_{\mathrm{0}}{c}^{\mathrm{*}}。

2021年， Ronan等^［10］基于经典的离散SIS模型研究了模型

\left\{\begin{array}{l}{S}_{t+\mathrm{1}}=S_t-\beta S_t{I}_t+\gamma I_t，\\ I_{t+\mathrm{1}}=I_t+\beta S_t{I}_t-\gamma I_t，\end{array}\right.（1）

式中： S_t、 I_t分别表示t时刻易感者和染病者的数量； \gamma表示恢复率。为了体现行为的变化， 作者提出具有延迟的目标函数U(c)中g(c)=-{\alpha }_{\mathrm{1}}{(c-\widehat{c})}^{\mathrm{2}}是一个凹函数， 是指偏离最大利益的目标接触率\widehat{c}的惩罚， h(c,I)=-{\alpha }_{\mathrm{2}}\left[\mathrm{1}-{\left(\mathrm{1}-\frac{I_{t-\Delta }}{N}{\beta }_{\mathrm{0}}\right)}^{\mathrm{c}}\right]是感染成本与感染率的乘积， 即感染造成的效用损失， 其中， I_{t-\Delta }是t-\Delta时刻染病者的数量， 时间延迟\Delta表示信息获取和对信息作出反应的延迟。在实际情况中， 效用不会无限制地减少。温利民等^［11］对股票投资策略进行研究时， 在线性函数的基础上结合了指数函数的性质， 给出了凹的混合型指数效用函数g(c)=(\mathrm{1}-{\mathrm{e}}^{-{\alpha }_{\mathrm{1}}c})/{\alpha }_{\mathrm{1}}+c_{\mathrm{0}}。文献［10］与［11］的两个厌恶函数的区别是

\underset{c\to \mathrm{0}}{\mathrm{l}\mathrm{i}\mathrm{m}}g(c)=-{\alpha }_{\mathrm{1}}{(-\widehat{c})}^{\mathrm{2}}，

\underset{c\to \mathrm{0}}{\mathrm{l}\mathrm{i}\mathrm{m}}g(c)=\underset{c\to \mathrm{0}}{\mathrm{l}\mathrm{i}\mathrm{m}}[(\mathrm{1}-{\mathrm{e}}^{-{\alpha }_{\mathrm{1}}c})/{\alpha }_{\mathrm{1}}+c_{\mathrm{0}}]\mathrm{。}

对于突发传染病来说， 人们采取措施（如保持安全距离、 减少出行等）引起行为改变， 个体间的接触趋于0时， 显然会存在效用成本， 因此， 本文中的g(c)将采用文献［11］的形式， h(c,I)将采用文献［10］的形式， 效用函数可表示为

\begin{array}{c}U(c)={\alpha }_{\mathrm{0}}+\frac{\mathrm{1}-{\mathrm{e}}^{-{\alpha }_{\mathrm{1}}c}}{{\alpha }_{\mathrm{1}}}+c_{\mathrm{0}}-\\ {\alpha }_{\mathrm{2}}\left[\mathrm{1}-{\left(\mathrm{1}-\frac{I_{t-\Delta }}{N}{\beta }_{\mathrm{0}})\right)}^c\right]\mathrm{。}\end{array}（2）

为求得最优接触率并使计算简单， 将式（2）在c=\mathrm{0}处泰勒展开到二次幂， 则目标函数近似为

\begin{array}{c}U(c)={\alpha }_{\mathrm{0}}+c_{\mathrm{0}}+\left[\mathrm{1}+{\alpha }_{\mathrm{2}}\mathrm{l}\mathrm{n}\left(\mathrm{1}-\frac{I_{t-\Delta }}{N}{\beta }_{\mathrm{0}}\right)\right]c+\\ \frac{\mathrm{1}}{\mathrm{2}}\left[-{\alpha }_{\mathrm{1}}+{\alpha }_{\mathrm{2}}{\left(\mathrm{l}\mathrm{n}\left(\mathrm{1}-\frac{I_{t-\Delta }}{N}{\beta }_{\mathrm{0}}\right)\right)}^{\mathrm{2}}\right]c^{\mathrm{2}}\mathrm{。}\end{array}（3）

对式（3）两边关于c求导， 得

\begin{array}{l}\frac{\mathrm{d}U}{\mathrm{d}c}=\mathrm{1}+{\alpha }_{\mathrm{2}}\mathrm{l}\mathrm{n}\left(\mathrm{1}-\frac{I_{t-\Delta }}{N}{\beta }_{\mathrm{0}}\right)+\\ \left[-{\alpha }_{\mathrm{1}}+{\alpha }_{\mathrm{2}}{\left(\mathrm{l}\mathrm{n}\left(\mathrm{1}-\frac{I_{t-\Delta }}{N}{\beta }_{\mathrm{0}}\right)\right)}^{\mathrm{2}}\right]c\mathrm{。}\end{array}

求解\mathrm{d}U/\mathrm{d}c=\mathrm{0}可得

c=\frac{\mathrm{1}+{\alpha }_{\mathrm{2}}\mathrm{l}\mathrm{n}\left(\mathrm{1}-\frac{I_{t-\Delta }}{N}{\beta }_{\mathrm{0}}\right)}{{\alpha }_{\mathrm{1}}-{\alpha }_{\mathrm{2}}{\left[\mathrm{l}\mathrm{n}\left(\mathrm{1}-\frac{I_{t-\Delta }}{N}{\beta }_{\mathrm{0}}\right)\right]}^{\mathrm{2}}}\mathrm{。}（4）

假设人们采取戴口罩、 勤洗手及减少出行等措施， 就会将染病者数量控制到很小， 因此， 在分析行为改变后的系统动力学性态时， 为使模型简化， 在I_{t-\Delta }=\mathrm{0}处泰勒展开， 式（4）可写为c=a-b\frac{I_{t-\Delta }}{N}+o(I_{t-\Delta }^{\mathrm{2}})， 其中， a=\mathrm{1}/{\alpha }_{\mathrm{1}}， b={\alpha }_{\mathrm{2}}{\beta }_{\mathrm{0}}/{\alpha }_{\mathrm{1}}。因此取最优接触率为c^{\mathrm{*}}=a-b\frac{I_{t-\Delta }}{N}， 即S与I接触时感染的概率为\beta ={\beta }_{\mathrm{0}}\left(a-b\frac{I_{t-\Delta }}{N}\right)。

在传染病模型中， 饱和发生率能够更好地描述传染病传播中的限制因素， 如保持社交距离、 减少出行等。因此， 本文建立了一个具有时间延迟的， 考虑行为变化影响且为饱和发生率的离散传染病模型， 具体为

\left\{\begin{array}{l}{S}_{t+\mathrm{1}}=S_t-\frac{{\beta }_{\mathrm{0}}\left(a-b\frac{I_{t-\Delta }}{N_t}\right)S_t{I}_t}{\mathrm{1}+\alpha I_t}+\gamma I_t，\\ I_{t+\mathrm{1}}=I_t+\frac{{\beta }_{\mathrm{0}}\left(a-b\frac{I_{t-\Delta }}{N_t}\right)S_t{I}_t}{\mathrm{1}+\alpha I_t}-\gamma I_t，\end{array}\right.（5）

式中： \alpha为饱和率常数； N_t=S_t+I_t是总人口的数量， 假设N_t在封闭种群中是固定的， 且短时间内没有输入输出， 因此， 用N表示恒定的总人口。

1 平衡点的存在性

系统（5）的正向不变可行区域为

D=\{(S_t,I_t)\left|S_t\ge \mathrm{0},\right.I_t\ge \mathrm{0},\left|S_t+\right.I_t\le N\}\mathrm{。}

1.1 无病平衡点

将S_t=N-I_t代入系统（5）的第2个式子， 化简得

\begin{array}{c}{I}_{t+\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{1}+\alpha I_t}\left[I_t^{\mathrm{2}}\left(\alpha -a{\beta }_{\mathrm{0}}+\frac{b{\beta }_{\mathrm{0}}}{N}{I}_{t-\Delta }\right)+\right.\\ \left.\begin{array}{c}\\ \end{array}{I}_t(\mathrm{1}+a{\beta }_{\mathrm{0}}N-b{\beta }_{\mathrm{0}}{I}_{t-\Delta }-\gamma -\alpha \gamma )\right]=f(I_t,I_{t-\Delta })\mathrm{。}\end{array}（6）

当没有时间延迟即\Delta =\mathrm{0}时， f(I_t)是一个I_t的三次多项式，

\begin{array}{c}{I}_{t+\mathrm{1}}=\frac{b{\beta }_{\mathrm{0}}}{N}{I}_t^{\mathrm{3}}-(a{\beta }_{\mathrm{0}}+b{\beta }_{\mathrm{0}}+\alpha \gamma )I_t^{\mathrm{2}}+\\ (\mathrm{1}+a{\beta }_{\mathrm{0}}N-\gamma )I_t =f(I_t)，\end{array}（7）

为分析信息没有延迟时（\Delta =\mathrm{0}）的平衡点的存在性， 令I_{t+\mathrm{1}}=I_t\triangleq I， 则式（7）是一个关于I的三次方程， 即

\begin{array}{c}f(I)\triangleq \frac{b{\beta }_{\mathrm{0}}}{N}{I}^{\mathrm{3}}-(a{\beta }_{\mathrm{0}}+b{\beta }_{\mathrm{0}}+\alpha \gamma )I^{\mathrm{2}}+\\ (a{\beta }_{\mathrm{0}}N-\gamma )I=\mathrm{0}\mathrm{。}\end{array}（8）

显然I=\mathrm{0}是方程（7）的一个解， 即系统（5）有无病平衡点E_{\mathrm{0}}=(S_{\mathrm{0}},I_{\mathrm{0}})=(N,\mathrm{0})。

下面利用下一代生成矩阵^［12］的方法来计算出系统（5）的基本再生数为R_{\mathrm{0}}=\frac{a{\beta }_{\mathrm{0}}N}{\gamma }\mathrm{。}

1.2 正平衡点

当I\ne \mathrm{0}时， 将式（8）化简整理得

\begin{array}{c}{f}_{\mathrm{1}}(I)\triangleq \frac{b{\beta }_{\mathrm{0}}}{N}{I}^{\mathrm{2}}-(a{\beta }_{\mathrm{0}}+b{\beta }_{\mathrm{0}}+\alpha \gamma )I+\\ \gamma (R_{\mathrm{0}}-\mathrm{1})=\mathrm{0}\mathrm{。}\end{array}（9）

记

{\Delta }_{\mathrm{1}}=(a{\beta }_{\mathrm{0}}+b{\beta }_{\mathrm{0}}{+\alpha \gamma )}^{\mathrm{2}}-\mathrm{4}\frac{b{\beta }_{\mathrm{0}}}{N}\gamma (R_{\mathrm{0}}-\mathrm{1})，

经验证，

{\Delta }_{\mathrm{1}}=[({\beta }_{\mathrm{0}}^{\mathrm{2}}{(a-b)}^{\mathrm{2}}+\mathrm{2}\alpha \gamma {\beta }_{\mathrm{0}}(a+b)+{\alpha }^{\mathrm{2}}{\gamma }^{\mathrm{2}})N+\mathrm{4}b{\beta }_{\mathrm{0}}\gamma ]N>\mathrm{0}\mathrm{。}

系统（5）正平衡点的存在性问题可以转化为分析方程（9）正根的存在性问题。下面分两种情形分别讨论。

情形 1　当R_{\mathrm{0}}\le \mathrm{1}时， 方程（9）存在一个正根

I_{\mathrm{*}}=\frac{(a+b){\beta }_{\mathrm{0}}+\alpha \gamma +\sqrt[]{{\Delta }_{\mathrm{1}}}}{\mathrm{2}b{\beta }_{\mathrm{0}}}N\mathrm{。}

容易判别当R_{\mathrm{0}}\le \mathrm{1}时， I_{\mathrm{*}}>N， 它不在可行域D内， 故当R_{\mathrm{0}}\le \mathrm{1}时， 系统（5）不存在正平衡点。

情形 2　当R_{\mathrm{0}}>\mathrm{1}时， 方程（9）存在两个正根I_{\mathrm{1}}^{\mathrm{*}},I_{\mathrm{2}}^{\mathrm{*}}。容易验证

故系统（5）存在正平衡E_{\mathrm{1}}^{\mathrm{*}}=(S_{\mathrm{1}}^{\mathrm{*}},I_{\mathrm{1}}^{\mathrm{*}})=(N-I_{\mathrm{1}}^{\mathrm{*}},I_{\mathrm{1}}^{\mathrm{*}})\mathrm{。}

综上所述， 有下面的定理：

定理 1　系统（5）始终存在一个无病平衡点E_{\mathrm{0}}； 当R_{\mathrm{0}}\le \mathrm{1}时， 系统（5）在可行域D内不存在正平衡点； 当R_{\mathrm{0}}>\mathrm{1}时， 系统（5）存在一个正平衡点E_{\mathrm{1}}^{\mathrm{*}}。

2 平衡点的稳定性

下面首先分析无病平衡点的稳定性^［13-14］， 计算f(I_t)在无病平衡点I_{\mathrm{0}}=\mathrm{0}处的导数

{\left|\frac{\mathrm{d}}{\mathrm{d}I}f(I_t)\right|}_{I_t=I_{\mathrm{0}}}=\left|\mathrm{1}+Na{\beta }_{\mathrm{0}}-\gamma \right|=\left|\mathrm{1}+\gamma (R_{\mathrm{0}}-\mathrm{1})\right|\mathrm{。}

无病平衡点I_{\mathrm{0}}局部渐近稳定的充要条件为通过简单计算得到下面定理：

定理 2　当时， 无病平衡点I_{\mathrm{0}}局部渐近稳定； 当R_{\mathrm{0}}>\mathrm{1}时， I_{\mathrm{0}}不稳定。

定理 3　当时， 无病平衡点I_{\mathrm{0}}全局渐近稳定。

证明　构造Lyapunov函数V_t=I_t， 对模型求差分， 得

\mathrm{\Delta }{V}_t=V_{t+\mathrm{1}}-V_t=I_{t+\mathrm{1}}-I_t=\frac{{\beta }_{\mathrm{0}}(a-b\frac{I_t}{N})S_t{I}_t}{\mathrm{1}+\alpha I_t}-

\gamma I_t\le [Na{\beta }_{\mathrm{0}}-\gamma ]I_t=\gamma (R_{\mathrm{0}}-\mathrm{1})I_t\mathrm{。}

若， ， 易知当且仅当系统（5）在I_{\mathrm{0}}处， \mathrm{\Delta }{V}_t=\mathrm{0}。根据Krasovkii-LaSalle定理^［15］可得， 无病平衡点I_{\mathrm{0}}是全局渐近稳定的。证毕。

下面分析当R_{\mathrm{0}}=\mathrm{1}时I_{\mathrm{0}}的稳定性情况。当R_{\mathrm{0}}=\mathrm{1}时， {\left|\frac{\mathrm{d}}{\mathrm{d}I}f(I_t)\right|}_{I_t=I_{\mathrm{0}}}=\mathrm{1}，故系统（5）在I_{\mathrm{0}}处可能发生分岔。先在{\beta }_{\mathrm{0}}=\mathrm{0}的某个小邻域内分析， 令I_t=x_t， 将I_t平移至原点， 则式（7）变为

\begin{array}{c}{x}_{t+\mathrm{1}}=\frac{b{\beta }_{\mathrm{0}}}{N}{x}_t^{\mathrm{3}}-(a{\beta }_{\mathrm{0}}+b{\beta }_{\mathrm{0}}+\alpha \gamma )x_t^{\mathrm{2}}+\\ (\mathrm{1}+a{\beta }_{\mathrm{0}}N-\gamma )x_t ，\end{array}（10）

将式（10）在(x_t,{\beta }_{\mathrm{0}})=(\mathrm{0,0})处泰勒展开得\begin{array}{c}x\to f(x,{\beta }_{\mathrm{0}})=(\mathrm{1}-\gamma )x+aNx{\beta }_{\mathrm{0}}-\gamma x^{\mathrm{2}} +\\ {\rm O}({\left|x+{\beta }_{\mathrm{0}}\right|}^{\mathrm{3}})，\end{array}

通过计算发现， f(x,{\beta }_{\mathrm{0}})在R_{\mathrm{0}}=\mathrm{1}处满足

\begin{array}{c}f(\mathrm{0,0})=\mathrm{0},\frac{\partial f}{\partial I}(\mathrm{0,0})=\mathrm{1},\frac{\partial f}{\partial {\beta }_{\mathrm{0}}}(\mathrm{0,0})=\mathrm{0},\\ \frac{{\partial }^{\mathrm{2}}f}{\partial I\partial {\beta }_{\mathrm{0}}}(\mathrm{0,0})\ne \mathrm{0},\frac{{\partial }^{\mathrm{2}}f}{\partial I^{\mathrm{2}}}(\mathrm{0,0})\ne \mathrm{0}\mathrm{。}\end{array}

因此， 当R_{\mathrm{0}}=\mathrm{1}时， 系统（5）在无病平衡点I_{\mathrm{0}}处发生跨临界分岔。

接下来分析正平衡点I_{\mathrm{1}}^{\mathrm{*}}的稳定性， 计算f(I_t)在正平衡点I_{\mathrm{1}}^{\mathrm{*}}处的导数

\begin{array}{c}{\left|\frac{\mathrm{d}}{\mathrm{d}I}f(I_t)\right|}_{I_t=I_{\mathrm{1}}^{\mathrm{*}}}=\mathrm{1}+\frac{\alpha b{\beta }_{\mathrm{0}}{I_{\mathrm{1}}^{\mathrm{*}}}^{\mathrm{3}}}{N(\mathrm{1}+\alpha I_{\mathrm{1}}^{\mathrm{*}}{)}^{\mathrm{2}}}+\\ \frac{\mathrm{2}b{\beta }_{\mathrm{0}}{I_{\mathrm{1}}^{\mathrm{*}}}^{\mathrm{2}}-N{\beta }_{\mathrm{0}}(b+a+a\alpha N)I_{\mathrm{1}}^{\mathrm{*}}}{N(\mathrm{1}+\alpha I_{\mathrm{1}}^{\mathrm{*}}{)}^{\mathrm{2}}}\mathrm{。}\end{array}

若正平衡点局部渐近稳定， 只要满足， 即。

下面分别从以及来寻找平衡点稳定的条件。

1）　 的情况。

下面分别分析{\left[\frac{\mathrm{d}}{\mathrm{d}I}f(I_t)\right]}_{I_t=I_{\mathrm{1}}^{\mathrm{*}}}>\mathrm{0}及时需满足的条件。

i）　 {\left[\frac{\mathrm{d}}{\mathrm{d}I}f(I_t)\right]}_{I_t=I_{\mathrm{1}}^{\mathrm{*}}}>\mathrm{0}等价于

\begin{array}{c}{G}_{\mathrm{1}}(I_{\mathrm{1}}^{\mathrm{*}})=\alpha b{\beta }_{\mathrm{0}}{I_{\mathrm{1}}^{\mathrm{*}}}^{\mathrm{3}}+(\mathrm{2}b{\beta }_{\mathrm{0}}+N{\alpha }^{\mathrm{2}}){I_{\mathrm{1}}^{\mathrm{*}}}^{\mathrm{2}}-\\ [{\beta }_{\mathrm{0}}(b{\beta }_{\mathrm{0}}+a{\beta }_{\mathrm{0}}+N^{\mathrm{2}}a\alpha )-\mathrm{2}N\alpha ]I_{\mathrm{1}}^{\mathrm{*}}+N>\mathrm{0}\mathrm{。}\end{array}

由于平衡点满足方程（9）， 因此， 利用f_{\mathrm{1}}(I_{\mathrm{1}}^{\mathrm{*}})=\mathrm{0}可得G_{\mathrm{1}}(I_{\mathrm{1}}^{\mathrm{*}})=n_{\mathrm{1}}{I}_{\mathrm{1}}^{\mathrm{*}}+n_{\mathrm{2}}， 其中

\begin{array}{c}{n}_{\mathrm{1}}=N [N\alpha (a^{\mathrm{2}}+b^{\mathrm{2}})+(a+b)b]{\beta }_{\mathrm{0}}^{\mathrm{2}}+\\ N\alpha [\mathrm{2}N(\gamma +\mathrm{1})\left(a+b\right)\alpha +b\left(\mathrm{3}\gamma +\mathrm{4}\right)]{\beta }_{\mathrm{0}}+\\ N^{\mathrm{2}}{\alpha }^{\mathrm{3}}\gamma \left(\gamma +\mathrm{2}\right)，\end{array}\begin{array}{c}{n}_{\mathrm{2}}=-N^{\mathrm{2}}a[\alpha N(a+b)+\mathrm{2}b]{\beta }_{\mathrm{0}}^{\mathrm{2}}-\\ N[N^{\mathrm{2}}a(\gamma +\mathrm{2}){\alpha }^{\mathrm{2}}-N\gamma \left(a+b\right)\alpha -\\ \mathrm{2}b\left(\gamma +\mathrm{1}\right)]{\beta }_{\mathrm{0}}+N^{\mathrm{2}}{\alpha }^{\mathrm{2}}\gamma \left(\gamma +\mathrm{2}\right)\mathrm{。}\end{array}

显然， 在I_{\mathrm{1}}^{\mathrm{*}}=-\frac{n_{\mathrm{2}}}{n_{\mathrm{1}}}处， G_{\mathrm{1}}(I_{\mathrm{1}}^{\mathrm{*}})=\mathrm{0}， 且G_{\mathrm{1}}(I_{\mathrm{1}}^{\mathrm{*}})关于I_{\mathrm{1}}^{\mathrm{*}}单调递增， 故当I_{\mathrm{1}}^{\mathrm{*}}=\frac{(a+b){\beta }_{\mathrm{0}}+\alpha \gamma -\sqrt[]{\Delta }}{\mathrm{2}b{\beta }_{\mathrm{0}}}>-\frac{n_{\mathrm{2}}}{n_{\mathrm{1}}}时， G_{\mathrm{2}}(I_{\mathrm{1}}^{\mathrm{*}})>\mathrm{0}， 即{\left[\frac{\mathrm{d}}{\mathrm{d}I}f(I_t)\right]}_{I_t=I_{\mathrm{1}}^{\mathrm{*}}}>\mathrm{0}。

ii） 等价于G_{\mathrm{2}}(I_{\mathrm{1}}^{\mathrm{*}})={\beta }_{\mathrm{0}}[\alpha b{I_{\mathrm{1}}^{\mathrm{*}}}^{\mathrm{2}}。下面只需分析的条件。很容易判别G_{\mathrm{2}}(I_{\mathrm{1}}^{\mathrm{*}})在[\mathrm{0},N]上关于I_{\mathrm{1}}^{\mathrm{*}}单调递增， 且G_{\mathrm{2}}(m)=\mathrm{0}， 其中， m=\frac{-b+\sqrt[]{b^{\mathrm{2}}+\alpha bN(b+a+Na\alpha )}}{\alpha b}。显然当I_{\mathrm{1}}^{\mathrm{*}}=即\mathrm{2}b{\beta }_{\mathrm{0}}m>(a+b){\beta }_{\mathrm{0}}+\alpha \gamma -\sqrt[]{\Delta }时， 。

因此， 当时，

。

2）　 的情况。

下面分别分析条件{\left[\frac{\mathrm{d}}{\mathrm{d}I}f(I_t)\right]}_{I_t=I_{\mathrm{1}}^{\mathrm{*}}}\le \mathrm{0}以及{\left[\frac{\mathrm{d}}{\mathrm{d}I}f(I_t)\right]}_{I_t=I_{\mathrm{1}}^{\mathrm{*}}}>-\mathrm{1}时需满足的条件。由上面的分析易知， 当I_{\mathrm{1}}^{\mathrm{*}}\le -\frac{n_{\mathrm{2}}}{n_{\mathrm{1}}}时， {\left[\frac{\mathrm{d}}{\mathrm{d}I}f(I_t)\right]}_{I_t=I_{\mathrm{1}}^{\mathrm{*}}}\le \mathrm{0}。 {\left[\frac{\mathrm{d}}{\mathrm{d}I}f(I_t)\right]}_{I_t=I_{\mathrm{1}}^{\mathrm{*}}}>-\mathrm{1}等价于

F_{\mathrm{1}}(I_{\mathrm{1}}^{\mathrm{*}})=\alpha b{\beta }_{\mathrm{0}}{I_{\mathrm{1}}^{\mathrm{*}}}^{\mathrm{3}}+\mathrm{2}(N{\alpha }^{\mathrm{2}}+b{\beta }_{\mathrm{0}}){I_{\mathrm{1}}^{\mathrm{*}}}^{\mathrm{2}}-(a\alpha N^{\mathrm{2}}{\beta }_{\mathrm{0}}+aN{\beta }_{\mathrm{0}}+bN{\beta }_{\mathrm{0}}-\mathrm{4}N\alpha )I_{\mathrm{1}}^{\mathrm{*}}+\mathrm{2}N>\mathrm{0}，

利用上面的分析方法F_{\mathrm{1}}(I_{\mathrm{1}}^{\mathrm{*}})=n_{\mathrm{1}\mathrm{1}}{I}_{\mathrm{1}}^{\mathrm{*}}+n_{\mathrm{22}}， 其中，

\begin{array}{c}{n}_{\mathrm{11}}=N(N{a}^{\mathrm{2}}\alpha +N\alpha b^{\mathrm{2}}+ab+b^{\mathrm{2}}){\beta }_{\mathrm{0}}^{\mathrm{2}}+\\ N\alpha [\mathrm{2}N(\gamma +\mathrm{1})(a+b)\alpha +b(\mathrm{3}\gamma +\mathrm{4})]{\beta }_{\mathrm{0}}+\\ N^{\mathrm{2}}{\alpha }^{\mathrm{3}}(\gamma +\mathrm{2})\gamma ，\end{array}

\begin{array}{c}{n}_{\mathrm{22}}=-N^{\mathrm{2}}a(a\alpha N+\alpha bN+\mathrm{2}b){\beta }_{\mathrm{0}}^{\mathrm{2}}+\\ N[-N^{\mathrm{2}}a(\gamma +\mathrm{2}){\alpha }^{\mathrm{2}}+N\gamma (a+b)\alpha +\\ \mathrm{2}b(\gamma +\mathrm{1})]{\beta }_{\mathrm{0}}+N^{\mathrm{2}}{\alpha }^{\mathrm{2}}\gamma (\gamma +\mathrm{2})\mathrm{。}\end{array}

在I_{\mathrm{1}}^{\mathrm{*}}=-\frac{n_{\mathrm{22}}}{n_{\mathrm{11}}}处， F_{\mathrm{1}}(I_{\mathrm{1}}^{\mathrm{*}})=\mathrm{0}， 且F_{\mathrm{1}}(I_{\mathrm{1}}^{\mathrm{*}})关于I_{\mathrm{1}}^{\mathrm{*}}单调递增， 故当

I_{\mathrm{1}}^{\mathrm{*}}=\frac{(a+b){\beta }_{\mathrm{0}}+\alpha \gamma -\sqrt[]{{\Delta }_{\mathrm{1}}}}{\mathrm{2}{\beta }_{\mathrm{0}}}>-\frac{n_{\mathrm{22}}}{n_{\mathrm{11}}}时，

F_{\mathrm{1}}(I_{\mathrm{1}}^{\mathrm{*}})>\mathrm{0}， 即{\left[\frac{\mathrm{d}}{\mathrm{d}I}f(I_t)\right]}_{I_t=I_{\mathrm{1}}^{\mathrm{*}}}>-\mathrm{1}。下面利用\mathrm{M}\mathrm{a}\mathrm{p}\mathrm{l}\mathrm{e}软件， 经过计算可以判断I_{\mathrm{1}}^{\mathrm{*}}恒大于-\frac{n_{\mathrm{22}}}{n_{\mathrm{11}}}。

因此， 当I_{\mathrm{1}}^{\mathrm{*}}\le -\frac{n_{\mathrm{2}}}{n_{\mathrm{1}}}时， 。

综合以上分析， 有下面的定理：

定理 4　当R_{\mathrm{0}}>\mathrm{1}时， 若\mathrm{2}b{\beta }_{\mathrm{0}}m>(a+b){\beta }_{\mathrm{0}}+\alpha \gamma -\sqrt[]{{\Delta }_{\mathrm{1}}}， 系统（5）的正平衡点I_{\mathrm{1}}^{\mathrm{*}}是局部渐近稳定的； 若， 正平衡点I_{\mathrm{1}}^{\mathrm{*}}不稳定。

3 数值模拟

本节将通过数值模拟分析时间延迟\Delta（如信息获取和行为反应的延迟）对系统性态的影响。依据文献［10］参数的选择， 固定N=\mathrm{10}\mathrm{000}， {\beta }_{\mathrm{0}}=\mathrm{0.05}， b=\mathrm{0.1}， \alpha =\mathrm{0.375}，a=\mathrm{0.001}\mathrm{5}， 分别模拟\Delta =\mathrm{0,1},\mathrm{2,3},\mathrm{4,5}时， 系统（6）的性态变化情况。由图 2（a） 可以看出： 当\Delta =\mathrm{0}时， 随着时间的增加， 染病者I单调趋于正平衡点I_{\mathrm{1}}^{\mathrm{*}}； 当\Delta =\mathrm{1}和\Delta =\mathrm{2}时， 系统（6）会振荡收敛到地方病平衡点I_{\mathrm{1}}^{\mathrm{*}}（见图 2（b） 和图 2（c））； 当\Delta =\mathrm{3}时， 系统（6）的轨线围绕正平衡点周期性振荡， 如图 2（d） 所示； 当\Delta =\mathrm{4}及\Delta =\mathrm{5}时， 如图 2（d） 和图 2（e）， 系统（6）失去稳定， 染病者数量在某个时刻急剧下降， 变为负值， 系统性态崩塌（见图 2（e） 和图 2（f））。

为了全面分析时间延迟对传染病传播的影响， 本文将所提出的模型与文献［10］中的模型进行了对比。尽管两种模型都考虑了社会行为的影响并引入了时间延迟， 但在模型结构上存在显著差异。首先， 本文模型采用了混合指数型损失厌恶函数， 这一函数设计能够更加准确地反映人们在面对传染病时的风险感知以及相应的行为决策。特别是在感染风险较高的情况下， 该函数能够有效捕捉人们采取预防措施的行为， 从而更有效地控制感染人数的增长。相比之下， 文献［10］中的模型采用了线性效用函数， 其对感染人数变化的描述相对平缓。其次， 随着时间延迟的增加， 本文模型展现出了更为稳定的动力学行为， 而文献［10］的模型则更容易出现振荡和周期性流行现象， 本文在时间延迟\Delta =\mathrm{3}时发生振荡， 文献［10］在\Delta =\mathrm{2}时发生振荡。这一特性使得本文模型在面对复杂多变的疫情发展时能够提供更为可靠的预测和决策支持。

4 结　论

本文通过时间延迟的效用函数得到平衡预期成本和感染风险的最优接触率， 进而建立了具有最优接触率的离散传染病模型。通过分析模型动力学性态， 给出了系统无病平衡点和正平衡点的存在和稳定的条件。利用Lyapunov函数的方法证明， 当时， 无病平衡点是全局渐近稳定的； 利用Jury判据得到， 当R_{\mathrm{0}}>\mathrm{1}时， 在正平衡点局部渐近稳定的条件。通过数值模拟发现， 时间延迟由小变大时， 系统性态出现了从单调稳定到振荡稳定以及出现周期解的变化， 甚至存在峰值和激增， 增加了系统性态的复杂性。因此， 在传染病流行期间， 政府和公共卫生部门应及时准确发布疫情信息， 避免由于信息获取不及时等增加疾病传播的复杂性。政府还应该考虑在传染病爆发早期， 严格实施防控措施， 提供明确的行动指南， 以确保人们能够快速采取有效防护措施。

在传染病爆发期间， 信息延迟对防控措施的实施非常重要。进一步可以考虑多群体模型， 如不同年龄段、 职业等群体间的差异； 还可以引入动态政策干预以及多政策组合， 如防控措施随时间的变化。

参考文献

[1] 王婉菁， 张兵， 朱红兵.经济活动与传染病的传播——基于流行性感冒的研究［J］.中国经济问题， 2023（6）： 176-191.

[2] WANG Wanjing， ZHANG Bing， ZHU Hongbing. Economic activities and the spread of infectious diseases： Research based on influenza［J］. China Economic Studies， 2023（6）： 176-191. （in Chinese）

[3] FUNK S， SALATHÉ M， JANSEN V A A. Modelling the influence of human behaviour on the spread of infectious diseases： A review［J］. Journal of the Royal Society Interface， 2010， 7（50）： 1247-1256.

[4] ARTHUR R F， GURLEY E S， SALJE H， et al. Contact structure， mobility， environmental impact and behaviour： The importance of social forces to infectious disease dynamics and disease ecology［J］. Philosophical Transactions of the Royal Society of London， Series B， Biological Sciences， 2017， 372（1719）： 20160454.

[5] NAIDOO M， SHEPHARD W， KAMBEWE I， et al.Incorporating social vulnerability in infectious disease mathematical modelling： A scoping review［J］.BMC Medicine， 2024， 22（1）： 125.

[6] SMITH R D， KEOGH-BROWN M R， BARNETT T， et al. The economy-wide impact of pandemic influenze on the UK： A computable general equilibrium modelling experiment［J］. BMJ， 2009， 339： b4571.

[7] KEOGH-BROWN M R， WREN LEWIS S， EDMUNDS W J， et al. The possible macroeconomic impact on the UK of an influenza pandemic［J］. Health Economics， 2010， 19（11）： 1345-1360.

[8] FENICHEL E P， CASTILL CHAVEZ C， CEDDIA M G， et al. Adaptive human behavior in epidemiological models［J］. Proceedings of the National Academy of Sciences of the United States of America， 2011， 108（15）： 6306-6311.

[9] PANICKER A， SASIDEVAN V. Social adaptive behavior and oscillatory prevalence in an epidemic model on evolving random geometric graphs［J］. Chaos， Solitons and Fractals， 2024， 178： 114407.

[10] 吴鑫育， 周海林. 定价核、 市场效用函数与投资者偏好［J］. 系统工程学报， 2017， 32（1）： 44-54.

[11] WU Xinyu， ZHOU Hailin. Pricing kernels， market utility functions and investor preferences［J］. Journal of Systems Engineering， 2017， 32（1）： 44-54. （in Chinese）

[12] RONAN F A， JAMES J H， MATTHEW H B， et al. Adaptive social contact rates induce complex dynamics during epidemics ［J］. PLoS Computational Biology， 2021， 17（2）： 1-17.

[13] 温利民， 冯会珍， 李俊雪， 等.基于混合指数型损失厌恶函数的投资组合模型［J］. 江西师范大学学报（自然科学版）， 2023， 47（1）： 1-7.

[14] WEN Limin， FENG Huizhen， LI Junxue， et al. The portfolio model based on mixed exponential loss aversion function ［J］. Journal of Jiangxi Normal University（Natural Science Edition）， 2023， 47（1）： 1-7.（in Chinese）

[15] VAN DEN DRIESSCHE P， WATMOUGH J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission［J］. Mathematical Biosciences， 2002， 180（1/2）： 29-48.

[16] BUONOMO B， ONOFRIO A， LACITIGNOLA D. Global stability of an SIR epidemic model with information dependent vaccination［J］. Mathematical Biosciences. 2008， 216（1）： 9-16.

[17] LI M Y， MULDOWNEY J S. Global stability for the SEIR model in epidemiology［J］. Mathematical Biosciences. 1995， 125（2）： 155-164.

[18] MARTCHEVA M. An introduction to mathematical epidemiology［M］. Boston： Springer， 2015.

基金资助

山西省自然科学基金资助项目(201901D111179)

山西省自然科学基金资助项目(202203021212130)