带有个体自发行为变化的传染病模型的研究

The Epidemic Model with Individual Spontaneous Behavior Change

摘要: 本文考虑了一类带有个体自发行为变化的SIR模型。在某些传染病流行期间，易感人群可以采取正常行为或谨慎行为(如戴口罩、保持社交距离等)，个体通过比较两种行为的回报(包括感染风险和经济成本)自发地选择其中之一。这种个体行为变化可以通过进化博弈论中的模仿动力学建模。在个体行为变化的时间尺度比传染病传播的时间尺度快得多的情况下，用几何奇异摄动理论分析我们建立的带有个体自发行为变化的SIR模型。借助快–慢结构和入–出积分得到模型的奇异轨道，并进行数值模拟。

Abstract: In this paper, we consider a SIR model with individual spontaneous behavior change. During some epidemics, susceptibles can adopt normal or altered behaviors (such as wearing masks, social dis-tancing, etc.). Individuals spontaneously choose one or the other by comparing the payoff (including the risk of infection and the economic costs) of the two behaviors. This individual behavior changes are modeled by imitation dynamics in evolutionary game theory. Under the condition where the time scale of behavior change is much faster than the time scale of infectious disease transmission, we use the geometric singular perturbation theory to analyze the SIR model with individual spontaneous behavior change. The singular orbit of the model is obtained by slow-fast structure and entry-exit integration, and simulated numerically.

1. 引言

传染病模型在控制传染病传播和治疗方面有重要意义，如SIR模型、SIS模型、SEIR模型等 [1] [2] [3] 。实际上，个体行为变化也会对传染病模型产生影响，为此很多学者开始研究带有个体行为变化的传染病模型 [4] [5] [6] [7] 。Reluga等学者使用进化博弈论将SIR模型与疫苗接种相结合，研究个体利益对疫苗摄取的影响 [8] ；Bauch等学者基于进化博弈论研究个体疫苗接种行为对疫苗恐慌的影响 [9] ；Poletti等学者基于进化博弈论研究风险认知对个体行为变化和传染病传播的影响 [10] 。个体行为变化体现在：某些传染病流行期间，易感者会随着人群中感染者数量的不断增加，而选择呆在家里、保持社交距离、戴口罩等谨慎行为，这些行为变化使得传染病的传播率降低。而当人群中的感染者数量不断减少，易感者又会自发的恢复到正常行为 [11] 。

2009年，Poletti等学者建立了带有个体自发行为变化的SIR模型 [12] ：在传染病流行期间，易感个体可以采用正常行为或谨慎行为。谨慎行为的个体能够降低单位时间内的接触人数，因此对两种行为的人群考虑两种不同的传播率。设正常行为的易感个体 $S_n$ 的传播率为 ${\beta }_n$ ，谨慎行为的易感个体 $S_a$ 的传播率为 ${\beta }_a$ 。假设易感个体可以根据成本或利益的考虑自发的改变自己的行为，这个现象可以用进化博弈论(进化博弈论指的是整个群体都可以采用某种策略，回报高的策略在群体中传播)中的模仿动力学(模仿动力学指的是策略可以通过模仿在群体中传播)描述，即行为对应于合适的博弈策略，有一定的预期回报 [12] [13] [14] 。模型使用双策略形式 $\frac{\text{d}x}{\text{d}t}\left(t\right)=x\left(t\right)\left(1-x\left(t\right)\right)\Delta p$ ，其中 $x\left(t\right)$ 为正常行为的易感人群比例， $1-x\left(t\right)$ 为谨慎行为的易感人群比例， $\Delta p$ 为两种行为的回报差。谨慎行为可以降低感染风险，但经济成本较高。正常行为的经济成本较低，但感染风险更高。Poletti等学者将个体行为变化引入传染病模型，建立带有个体自发行为变化的SIR模型：

$\begin{array}{l}\frac{\text{d}S}{\text{d}t}\left(t\right)=-\left({\beta }_{n}x\left(t\right)+{\beta }_a\left(1-x\left(t\right)\right)\right)S\left(t\right)I\left(t\right),\\ \frac{\text{d}I}{\text{d}t}\left(t\right)=\left({\beta }_{n}x\left(t\right)+{\beta }_a\left(1-x\left(t\right)\right)\right)S\left(t\right)I\left(t\right)-\gamma I\left(t\right),\\ \frac{\text{d}R}{\text{d}t}\left(t\right)=\gamma I\left(t\right),\\ \frac{\text{d}x}{\text{d}t}\left(t\right)=x\left(t\right)\left(1-x\left(t\right)\right)\left({\beta }_a-{\beta }_n\right)I\left(t\right)+\frac{\rho }{\alpha }x\left(t\right)\left(1-x\left(t\right)\right)\left(k-\left(m_n-m_a\right)I\left(t\right)\right).\end{array}$ (1)

2021年，Schecter用几何奇异摄动理论分析Poletti等建立的传染病模型，借助快–慢结构，构造出了模型的奇异轨道，并借助奇异轨道捕捉快时间尺度上的行为变化 [15] 。2021年，Mendonça等学者对Poletti等采用的回报函数进行了有意义的扩展，建立了新的回报函数，这种回报函数不仅与社会经济成本和感染者有关，而且与康复者和谨慎行为的易感者有关 [16] 。

在本文中，我们将SIR模型中的易感人群S分为正常行为的易感人群x和谨慎行为的易感人群 $1-x$ 两类，并将Mendonça等 [16] 提出的新回报函数引入模型，建立新的带有个体自发行为变化的SIR模型。我们建立的模型是多时间尺度(快时间尺度和慢时间尺度)模型，而几何奇异摄动理论是研究含有多时间尺度的复杂系统的工具，因此在Schecter工作 [15] 的启发下，我们用几何奇异摄动理论分析新的带有个体自发行为变化的SIR模型，并借助快–慢结构和入–出积分得到模型的奇异轨道。在第2节建立带有个体自发行为变化的SIR模型，第3节利用几何奇异摄动理论分析模型，第4节数值模拟并得出相关结果。

2. 模型建立

本文中我们建立如下带有个体自发行为变化的SIR模型：

$\begin{array}{l}\frac{\text{d}S}{\text{d}t}\left(t\right)=-\left({\beta }_{n}x\left(t\right)+{\beta }_a\left(1-x\left(t\right)\right)\right)S\left(t\right)I\left(t\right),\\ \frac{\text{d}I}{\text{d}t}\left(t\right)=\left({\beta }_{n}x\left(t\right)+{\beta }_a\left(1-x\left(t\right)\right)\right)S\left(t\right)I\left(t\right)-\gamma I\left(t\right),\\ \frac{\text{d}R}{\text{d}t}\left(t\right)=\gamma I\left(t\right),\\ \frac{\text{d}x}{\text{d}t}\left(t\right)=x\left(t\right)\left(1-x\left(t\right)\right)\left({\beta }_a-{\beta }_n\right)I\left(t\right)+\rho x\left(t\right)\left(1-x\left(t\right)\right)\left[k_0-k_{1}I\left(t\right)+k_{2}R\left(t\right)+k_3\left(1-x\left(t\right)\right)\right].\end{array}$ (2)

其中， ${\beta }_n$ 为正常行为的易感个体 $S_n\left(t\right)$ 的传播率， ${\beta }_a$ 为谨慎行为的易感个体 $S_a\left(t\right)$ 的传播率，且 ${\beta }_n>{\beta }_a$ 。 $S\left(t\right)=S_n\left(t\right)+S_a\left(t\right)$ 为易感人群比例， $I\left(t\right)$ 为感染人群比例， $R\left(t\right)$ 为康复人群比例，则 $S\left(t\right)+I\left(t\right)+R\left(t\right)=1$ 。 $x\left(t\right)=S_n\left(t\right)/\left(S_n\left(t\right)+S_a\left(t\right)\right)$ 为正常行为的易感人群比例， $1-x\left(t\right)=S_a\left(t\right)/\left(S_n\left(t\right)+S_a\left(t\right)\right)$ 为谨慎行为的易感人群比例。 $\gamma$ 为康复率， $\rho$ 为行为之间的变化率。

这里我们引入Mendonça等提出的新回报函数 [16] ，其中正常行为的回报函数为

$p_n\left(t\right)=-k_1^{n}I\left(t\right)+k_2^{n}R\left(t\right),$ (3)

谨慎行为的回报函数为

$p_a\left(t\right)=-k_0-k_1^{a}I\left(t\right)-k_2^{a}R\left(t\right)-k_3\left(1-x\left(t\right)\right),$ (4)

两种行为的回报差为

$\Delta p=p_n\left(t\right)-p_a\left(t\right)=k_0-k_{1}I\left(t\right)+k_{2}R\left(t\right)+k_3\left(1-x\left(t\right)\right).$ (5)

其中 $k_1=k_1^n-k_1^a>0$ ， $k_2=k_2^n+k_2^a>0$ 。

所有人都有感染的风险，而正常行为的感染风险比谨慎行为更高，(3)中 $-k_1^n$ 为正常行为人群的感染成本，(4)中 $-k_1^a$ 为谨慎行为人群的感染成本，且 $k_1^n>k_1^a>0$ 。如果有相当高比例的个体曾接触过传染病并恢复健康，那么感染的可能性就会降低，因此易感人群会随着康复人群比例的增加改变自己的行为，(3)中 $k_2^n$ 为正常行为人群的利益，(4)中 $-k_2^a$ 为谨慎行为人群的成本，其中 $k_2^n>0$ ， $k_2^a>0$ 。(4)中 $-k_0$ 为谨慎行为人群因个人生存问题产生的固定经济成本，其中 $k_0>0$ 。(4)中 $-k_3$ 为谨慎行为人群因粮食短缺和大规模裁员等产生的额外经济成本，这些成本会随着人群的增加而增加，其中 $k_3>0$ [16] 。模型(2)的仓室图如图1所示。

Figure 1. SIR model with individual spontaneous behavioral change diagram

图1. 带有个体自发行为变化的SIR模型图

为计算方便，对模型(2)进行简化。设存在参数 $0<\alpha <1$ ，使得 $1-x=\alpha S$ 且 $\alpha k_2=k_3$ 。因为 $R=1-I-S$ ，省略R的方程。则模型(2)简化为

$\begin{array}{l}\frac{\text{d}S}{\text{d}t}\left(t\right)=-\left({\beta }_{n}x\left(t\right)+{\beta }_a\left(1-x\left(t\right)\right)\right)S\left(t\right)I\left(t\right),\\ \frac{\text{d}I}{\text{d}t}\left(t\right)=\left({\beta }_{n}x\left(t\right)+{\beta }_a\left(1-x\left(t\right)\right)\right)S\left(t\right)I\left(t\right)-\gamma I\left(t\right),\\ \frac{\text{d}x}{\text{d}t}\left(t\right)=x\left(t\right)\left(1-x\left(t\right)\right)\left({\beta }_a-{\beta }_n\right)I\left(t\right)+\rho x\left(t\right)\left(1-x\left(t\right)\right)\left[k_0+k_2-\left(k_1+k_2\right)I\left(t\right)\right].\end{array}$ (6)

这里，模型(6)的状态空间是 $P=\left\{\left(S,I,x\right):S\ge 0,I\ge 0,S+I\le 1,0\le x\le 1\right\}$ 。

此时，模型(6)的基本再生数为 $R_0=\frac{{\beta }_n{x}^{\ast }+{\beta }_a\left(1-x^{\ast }\right)}{\gamma }=R_0^n{x}^{\ast }+R_0^a\left(1-x^{\ast }\right)$ ，其中 $R_0^n=\frac{{\beta }_n}{\gamma }$ ， $R_0^a=\frac{{\beta }_a}{\gamma }$ 。

对模型(6)做两个假设：

1) $R_0^a<1<R_0^n$ ；

2) $\frac{k_0+k_2}{k_1+k_2}<1$ ，虽然两种行为都不能保证总是带来更高的回报，但是这个假设保证，在 $0\le I\le 1$ 中，谨慎行为在I足够大时带来更高的回报，正常行为在I足够小时带来更高的回报。

3. 几何奇异摄动理论分析模型

在Schecter工作的启发下 [15] ，我们接下来用几何奇异摄动理论来分析模型(6)。由于传染病传播和行为变化可能以不同的速度发展，所以设传染病传播的时间尺度为t，行为变化的时间尺度为 $\tau =\frac{t}{\epsilon }$ 。在下文中假设行为变化的时间尺度比传染病传播的时间尺度快得多，即 $\epsilon >0$ 很小。假设变化率 $\rho$ 很大且与时间尺度有关，则令 $\rho =\frac{1}{\epsilon }$ 。

3.1. 快–慢时间

将 $\rho =\frac{1}{\epsilon }$ 代入模型(6)，得到时间尺度为慢时间t的慢系统

$\begin{array}{l}\stackrel{\dot{}}{S}=-\left({\beta }_{n}x+{\beta }_a\left(1-x\right)\right)SI,\\ \stackrel{\dot{}}{I}=\left({\beta }_{n}x+{\beta }_a\left(1-x\right)\right)SI-\gamma I,\\ \epsilon \stackrel{\dot{}}{x}=\epsilon x\left(1-x\right)\left({\beta }_a-{\beta }_n\right)I+x\left(1-x\right)\left[k_0+k_2-\left(k_1+k_2\right)I\right].\end{array}$ (7)

其中 ${}^{·}=\frac{\text{d}}{\text{d}t}$ 。

借助链式法则 $\frac{\text{d}}{\text{d}t}=\frac{\text{d}}{\text{d}\tau }\frac{\text{d}\tau }{\text{d}t}=\frac{1}{\epsilon }\frac{\text{d}}{\text{d}\tau }$ ，将时间尺度 $\tau =\frac{t}{\epsilon }$ 代入模型(6)，得到时间尺度为快时间 $\tau$ 的快系统

$\begin{array}{l}S\text{'}=-\epsilon \left({\beta }_{n}x+{\beta }_a\left(1-x\right)\right)SI,\\ I\text{'}=\epsilon \left({\beta }_{n}x+{\beta }_a\left(1-x\right)\right)SI-\epsilon \gamma I,\\ x\text{'}=\epsilon x\left(1-x\right)\left({\beta }_a-{\beta }_n\right)I+x\left(1-x\right)\left[k_0+k_2-\left(k_1+k_2\right)I\right].\end{array}$ (8)

其中 $\text{'}=\frac{\text{d}}{\text{d}\tau }$ 。

对于 $\epsilon >0$ ，慢系统(7)和快系统(8)具有相同的轨道，但沿轨道运动的速度不同。然而，在 $\epsilon =0$ 处的极限是不同的。因此，我们需要考虑慢系统(7)和快系统(8)的极限系统。

令慢系统(7)中 $\epsilon =0$ ，得到慢极限系统

$\begin{array}{l}\stackrel{\dot{}}{S}=-\left({\beta }_{n}x+{\beta }_a\left(1-x\right)\right)SI,\\ \stackrel{\dot{}}{I}=\left({\beta }_{n}x+{\beta }_a\left(1-x\right)\right)SI-\gamma I,\\ 0=x\left(1-x\right)\left[k_0+k_2-\left(k_1+k_2\right)I\right].\end{array}$ (9)

令快系统(8)中 $\epsilon =0$ ，得到快极限系统

$\begin{array}{l}S\text{'}=0,\\ I\text{'}=0,\\ x\text{'}=x\left(1-x\right)\left[k_0+k_2-\left(k_1+k_2\right)I\right].\end{array}$ (10)

3.2. 快极限系统

快极限系统(10)有三个平衡面： $x=0$ 、 $x=1$ 和 $I=\frac{k_0+k_2}{k_1+k_2}$ 。在 $x=0$ 或 $x=1$ 上一点的线性化矩阵有特征值0、0和 $\frac{\partial x\text{'}}{\partial x}=\left(1-2x\right)\left[k_0+k_2-\left(k_1+k_2\right)I\right]$ 。因此，只要 $I\ne \frac{k_0+k_2}{k_1+k_2}$ ，这两个平面上的平衡是双曲平衡。如果特征值是负的，平衡通常为吸引；如果特征值是正的，平衡通常为排斥 [17] 。因此，快极限系统(10)，在三角形 $x=0$ 上，当 $I<\frac{k_0+k_2}{k_1+k_2}$ 时，平衡通常是排斥的(正特征值)；当 $I>\frac{k_0+k_2}{k_1+k_2}$ 时，平衡通常是吸引的(负特征值)。在三角形 $x=1$ 上，当 $I>\frac{k_0+k_2}{k_1+k_2}$ 时，平衡通常是排斥的(正特征值)，当 $I<\frac{k_0+k_2}{k_1+k_2}$ 时，平衡通常是吸引的(负特征值)。

综上，我们可以得到 $0<x<1$ 的快极限系统(10)的如下结果：

1) S和I沿解不变。

2) 如果 $0<I<\frac{k_0+k_2}{k_1+k_2}$ ，则 $\stackrel{\dot{}}{x}>0$ ，且 $\underset{t\to -\infty }{\lim }x\left(t\right)=0$ ， $\underset{t\to +\infty }{\lim }x\left(t\right)=1$ 。

3) 如果 $\frac{k_0+k_2}{k_1+k_2}<I<1$ ，则 $\stackrel{\dot{}}{x}<0$ ，且 $\underset{t\to -\infty }{\lim }x\left(t\right)=1$ ， $\underset{t\to +\infty }{\lim }x\left(t\right)=0$ 。

上述结果说明，如果I很小，易感人群很快就会朝着正常行为的方向发展；如果I很大，易感人群很快就会朝着谨慎行为的方向发展。

3.3. 慢极限系统

慢极限系统(9)在状态空间的顶部和底部才有意义，即在三角形 $x=0$ 和 $x=1$ 上。

在三角形 $x=0$ (所有易感人群都是谨慎行为)上，慢极限系统(9)简化为

$\begin{array}{l}\stackrel{\dot{}}{S}=-{\beta }_{a}SI,\\ \stackrel{\dot{}}{I}={\beta }_{a}SI-\gamma I.\end{array}$ (11)

这是一个 $\beta ={\beta }_a$ ，基本再生数 $0<R_0^a<1$ 的SIR模型。

在三角形 $x=1$ (所有易感人群都是谨慎行为)上，慢极限系统(9)简化为

$\begin{array}{l}\stackrel{\dot{}}{S}=-{\beta }_{n}SI,\\ \stackrel{\dot{}}{I}={\beta }_{n}SI-\gamma I.\end{array}$ (12)

这是一个 $\beta ={\beta }_n$ ，基本再生数 $R_0^n>1$ 的SIR模型。

系统(11)或(12)的在SI平面上的相空间为三角形 $T=\left\{\left(S,I\right):S\ge 0,I\ge 0,S+I\le 1\right\}$ 。

设 $T_{+}=\left\{\left(S,I\right)\in T:S>0,I>0\right\}$ 。在 $T_{+}$ 中，(11)的轨道满足微分方程 $\frac{\text{d}I}{\text{d}S}=-1+\frac{\gamma }{{\beta }_{a}S}$ ，因此轨道是曲线

$I+S-\frac{\gamma }{{\beta }_a}\ln S=C.$ (13)

在 $T_{+}$ 中，(12)的轨道满足微分方程 $\frac{\text{d}I}{\text{d}S}=-1+\frac{\gamma }{{\beta }_{n}S}$ ，因此轨道是曲线

$I+S-\frac{\gamma }{{\beta }_n}\ln S=C.$ (14)

系统(11)和(12)有平衡线段 $I=0$ ， $0\le S\le 1$ 。每个解都趋于平衡 $\left(S_f,0\right)$ ，即当疫情结束时，没有人感染， $R=1-S_f$ 是在疫情期间感染的人口比例。

3.4. 入–出积分

3.4.1. 三角形 $x=0$ 的入–出积分

快极限系统(10)在点 $\left(S\left(t\right),I\left(t\right)\right)$ 处对 $x=0$ 的吸引或排斥由以下方程决定：

$\frac{\partial x\text{'}}{\partial x}\left(S\left(t\right),I\left(t\right),0\right)=k_0+k_2-\left(k_1+k_2\right)I\left(t\right).$ (15)

在三角形 $x=0$ 中，设 $\left(S_0,I_0\right)\in T_{+}$ ，其中 $I_0>\frac{k_0+k_2}{k_1+k_2}$ ，因此 $\left(S_0,I_0\right)$ 位于三角形的吸引部分。设 $\left(S\left(t\right),I\left(t\right)\right)$ 是系统(11)的解，令 $t_1>0$ ，且令 $\left(S_1,I_1\right)=\left(S\left(t_1\right),I\left(t_1\right)\right)$ 。解 $\left(S\left(t\right),I\left(t\right)\right)$ 描绘出曲线 $\Gamma$ ，由(13)得到方程 $I+S-\frac{\gamma }{{\beta }_a}\ln S=v_0$ ，则有

$v_0=I_0+S_0-\frac{\gamma }{{\beta }_a}\ln S_0=I_1+S_1-\frac{\gamma }{{\beta }_a}\ln S_1.$ (16)

三角形 $x=0$ 的入–出积分为：

$\begin{array}{c}{\mathcal{I}}_0\left(\left(S_0,I_0\right),\left(S_1,I_1\right)\right)={\displaystyle {\int }_0^{t_1}\frac{\partial x\text{'}}{\partial x}\left(S\left(t\right),I\left(t\right),0\right)\text{d}t}={\displaystyle {\int }_0^{t_1}{k}_0+k_2-\left(k_1+k_2\right)I\left(t\right)\text{d}t}\\ ={\displaystyle {\int }_{S_0}^{S_1}\frac{k_0+k_2-\left(k_1+k_2\right)\left(v_0-S+\frac{\gamma }{{\beta }_a}\ln S\right)}{{\beta }_{a}S\left(v_0-S+\frac{\gamma }{{\beta }_a}\ln S\right)}\text{d}S}.\end{array}$ (17)

将 $S=S\left(t\right)$ ， $\text{d}S=-{\beta }_{a}S\left(t\right)I\left(t\right)\text{d}t$ 和 $I=v_0-S+\frac{\gamma }{{\beta }_a}\ln S$ 代入第二个积分得到第三个积分。当 $I>\frac{k_0+k_2}{k_1+k_2}$ 时，第二个积分的被积函数为负，积分表示对平面 $x=0$ 的累积吸引力；当 $I<\frac{k_0+k_2}{k_1+k_2}$ 时，被积函数为正，积分表示对平面 $x=0$ 的累积排斥力。

在三角形 $x=0$ 中，对于 $T_{+}$ 中 $I_0>\frac{k_0+k_2}{k_1+k_2}$ 的每个点 $\left(S_0,I_0\right)$ ，都存在一个 $t_1>0$ ，使得 ${\mathcal{I}}_0\left(\left(S_0,I_0\right),\left(S_1,I_1\right)\right)=0$ 。则 $\left(S_1,I_1\right)$ 位于 $I<\frac{k_0+k_2}{k_1+k_2}$ 的区域，即在 $\left(S_1,I_1\right)$ 处，来自平面 $x=0$ 的累积排斥力平衡了对平面的累积吸引力。

综上可得，在三角形 $x=0$ 中，对于小的 $\epsilon >0$ ，快系统(8)的解在 $\left(S_0,I_0\right)$ 附近进入平面 $x=0$ 的邻域，然后遵循系统(11)，直到它离开 $\left(S_1,I_1\right)$ 附近的邻域。

3.4.2. 三角形 $x=1$ 的入–出积分

快极限系统(10)在点 $\left(S\left(t\right),I\left(t\right)\right)$ 处对 $x=1$ 的吸引或排斥由以下方程决定：

$\frac{\partial x\text{'}}{\partial x}\left(S\left(t\right),I\left(t\right),1\right)=-\left[k_0+k_2-\left(k_1+k_2\right)I\left(t\right)\right].$ (18)

在三角形 $x=1$ 中，设 $\left(S_0,I_0\right)\in T_{+}$ ，其中 $I_0<\frac{k_0+k_2}{k_1+k_2}$ ，因此 $\left(S_0,I_0\right)$ 位于三角形的吸引部分。设 $\left(S\left(t\right),I\left(t\right)\right)$ 是系统(12)的解，令 $t_1>0$ ，且令 $\left(S_1,I_1\right)=\left(S\left(t_1\right),I\left(t_1\right)\right)$ 。解 $\left(S\left(t\right),I\left(t\right)\right)$ 描绘出曲线 $\Gamma$ ，由(14)得到方程 $I+S-\frac{\gamma }{{\beta }_n}\ln S=v_0$ ，则有

$v_0=I_0+S_0-\frac{\gamma }{{\beta }_n}\ln S_0=I_1+S_1-\frac{\gamma }{{\beta }_n}\ln S_1.$ (19)

三角形 $x=1$ 的入–出积分为：

$\begin{array}{c}{\mathcal{I}}_1\left(\left(S_0,I_0\right),\left(S_1,I_1\right)\right)={\displaystyle {\int }_0^{t_1}\frac{\partial x\text{'}}{\partial x}\left(S\left(t\right),I\left(t\right),1\right)\text{d}t}=-{\displaystyle {\int }_0^{t_1}{k}_0+k_2-\left(k_1+k_2\right)I\left(t\right)\text{d}t}\\ ={\displaystyle {\int }_{S_0}^{S_1}\frac{k_0+k_2-\left(k_1+k_2\right)\left(v_0-S+\frac{\gamma }{{\beta }_n}\ln S\right)}{{\beta }_{n}S\left(v_0-S+\frac{\gamma }{{\beta }_n}\ln S\right)}\text{d}S}.\end{array}$ (20)

当 $I<\frac{k_0+k_2}{k_1+k_2}$ 时，第二个积分的被积函数为负，积分表示对平面 $x=1$ 的累积吸引力；当 $I>\frac{k_0+k_2}{k_1+k_2}$ 时，被积函数为正，积分表示对平面 $x=1$ 的累积排斥力。

在三角形 $x=1$ 中，与 $x=0$ 不同，不一定存在一个 $t_1>0$ ，使得 ${\mathcal{I}}_1\left(\left(S_0,I_0\right),\left(S_1,I_1\right)\right)=0$ 。系统(12)的解可能不会进入平面的排斥部分，或者只在排斥部分停留很短的时间随后就进入吸引部分。

综上可得，在三角形 $x=1$ 中，对于小的 $\epsilon >0$ ，快系统(8)的解在 $\left(S_0,I_0\right)$ 附近进入平面 $x=1$ 的邻域。

1、如果存在时间 $t_1>0$ ，使得 ${\mathcal{I}}_1\left(\left(S_0,I_0\right),\left(S_1,I_1\right)\right)=0$ ，那么解将离开点 $\left(S_1,I_1\right)$ 附近的邻域。

2、如果不存在这样的 $t_1$ ，那么解将永远不会离开邻域，并将继续遵循系统(12)。

3.5. 奇异轨道

基于前面的分析，借助快–慢结构和入–出积分，构建了快系统(8)(或慢系统(7))的奇异轨道。

考虑起点为 $P_{start}=\left(S_0,I_0,x_0\right)$ 的奇异轨道 $\mathcal{S}$ (见图2)，其中 $\left(S_0,I_0\right)\in T_{+}$ ， $I_0<\frac{k_0+k_2}{k_1+k_2}$ ，且 $0<x_0<1$ 。奇异轨道 $\mathcal{S}$ 可以用序列 $\left(P_{start},P_1,P_2,\cdots ,P_{2k},P_{end}\right)$ 表示，其中

-第一个快轨道从 $P_{start}=\left(S_0,I_0,x_0\right)$ 到 $\left(S_0,I_0,1\right)$ 。由于 $I_0<\frac{k_0+k_2}{k_1+k_2}$ ，所有个体迅速采用正常行为；

-第一个慢轨道从 $\left(S_0,I_0,1\right)$ 到 $P_1=\left(S_1,I_1,1\right)$ 。在平面 $x=1$ 中，存在 $t_1>0$ ，使得 ${\mathcal{I}}_1\left(\left(S_0,I_0\right),\left(S_1,I_1\right)\right)=0$ 且 $I_1>\frac{k_0+k_2}{k_1+k_2}$ ；

-第二个快轨道从 $P_1=\left(S_1,I_1,1\right)$ 到 $\left(S_1,I_1,0\right)$ 。由于 $I_1>\frac{k_0+k_2}{k_1+k_2}$ ，所有个体迅速采用谨慎行为；

-第二个慢轨道从 $\left(S_1,I_1,0\right)$ 到 $P_2=\left(S_2,I_2,0\right)$ 。在平面 $x=0$ 中，存在 $t_1>0$ ，使得 ${\mathcal{I}}_0\left(\left(S_1,I_1\right),\left(S_2,I_2\right)\right)=0$ 且 $I_2<\frac{k_0+k_2}{k_1+k_2}$ ；

-第三个快轨道从 $P_2=\left(S_2,I_2,0\right)$ 到 $\left(S_2,I_2,1\right)$ 。由于 $I_2<\frac{k_0+k_2}{k_1+k_2}$ ，所有个体迅速采用正常行为；

$⋮$

-最后一个快轨道从 $P_{2k}=\left(S_{2k},I_{2k},0\right)$ 到 $\left(S_{2k},I_{2k},1\right)$ 。由于 $I_{2k}<\frac{k_0+k_2}{k_1+k_2}$ ，所有个体迅速采用正常行为；

-最后一个慢轨道从 $\left(S_{2k},I_{2k},1\right)$ 到 $P_{end}=\left(S_f,0,1\right)$ 。这个轨道接近(12)的平衡，奇异轨道 $\mathcal{S}$ 的构造终止。

换句话说， $P_1,\cdots ,P_{2k}$ 是快速跳跃的起点；当i为奇数时， $P_i$ 在 $x=1$ 中；当i为偶数时， $P_i$ 在 $x=0$ 中。图2展示了一个奇异轨道的开始部分，其中蓝色部分表示快轨道，黑色部分表示慢轨道。

Figure 2. A singular orbit

图2. 奇异轨道

4. 数值模拟

模型假设： $\epsilon =0.005$ ， ${\beta }_n=0.5$ 、 ${\beta }_a=0.05$ 、 $\gamma =0.1$ ，则正常行为 $R_0^n=5$ ，谨慎行为 $R_0^a=0.5$ 。

数值模拟1： $P_{start}=\left(0.99,0.01,0.98\right)$ 。

1) $I=\frac{k_0+k_2}{k_1+k_2}=\frac{1}{20}$ ，其中 $k_0=0.2$ 、 $k_1=5.9$ 、 $k_2=0.1$ ，产生六次跳跃的奇异轨道。

$\begin{array}{l}{P}_{start}=\left(0.99,0.01,0.98\right)\\ P_1=\left(0.8304487197,0.1344035278,1\right)\\ P_2=\left(0.7499955134,0.0110588211,0\right)\\ P_3=\left(0.5872560692,0.1248770121,1\right)\end{array}$

$\begin{array}{l}{P}_4=\left(0.5432992953,0.0132325820,0\right)\\ P_5=\left(0.3741281508,0.1077913328,1\right)\\ P_6=\left(0.3540923350,0.0177457975,0\right)\\ P_{end}=\left(0.0839349074,0,1\right)\end{array}$

(a) 三角形 $x=1$ 中慢极限系统(12)的相图。其中垂直线 (b) 三角形 $x=0$ 中慢极限系统(11)的相图。其中水平线 $S=\frac{\gamma }{{\beta }_n}=\frac{1}{5}$ ，水平线 $I=\frac{k_0+k_2}{k_1+k_2}=\frac{1}{20}$ 。 $I=\frac{k_0+k_2}{k_1+k_2}=\frac{1}{20}$ 。

Figure 3. Phase portraits of Numerical simulation 1 (1)

图3. 数值模拟1 (1)相图

从图3数值模拟1 (1)的相图可以看到：从 $P_{start}$ 到 $P_1$ 、从 $P_2$ 到 $P_3$ 、从 $P_4$ 到 $P_5$ 和从 $P_6$ 到 $P_{end}$ 的慢轨道在三角形 $x=1$ 中慢极限系统(12)的相图中显示；从 $P_1$ 到 $P_2$ 、从 $P_3$ 到 $P_4$ 和从 $P_5$ 到 $P_6$ 的慢轨道在三角形 $x=0$ 中慢极限系统(11)的相图中显示。

2) $I=\frac{k_0+k_2}{k_1+k_2}=\frac{1}{15}$ ，其中 $k_0=0.3$ 、 $k_1=8.7$ 、 $k_2=0.3$ ，产生四次跳跃的奇异轨道。

$\begin{array}{l}{P}_{start}=\left(0.99,0.01,0.98\right)\\ P_1=\left(0.7519573308,0.1930375968,1\right)\end{array}$

$\begin{array}{l}{P}_2=\left(0.6542259343,0.0123153447,0\right)\\ P_3=\left(0.4116559272,0.1622323749,1\right)\\ P_4=\left(0.3764635621,0.0186915191,0\right)\\ P_{end}=\left(0.0765316159,0,1\right)\end{array}$

(a) 三角形 $x=1$ 中慢极限系统(12)的相图。其中垂直线 (b) 三角形 $x=0$ 中慢极限系统(11)的相图。其中水平线 $S=\frac{\gamma }{{\beta }_n}=\frac{1}{5}$ ，水平线 $I=\frac{k_0+k_2}{k_1+k_2}=\frac{1}{15}$ 。 $I=\frac{k_0+k_2}{k_1+k_2}=\frac{1}{15}$ 。

Figure 4. Phase portraits of Numerical simulation 1 (2)

图4. 数值模拟1 (2)相图

从图4数值模拟1 (2)的相图可以看到：从 $P_{start}$ 到 $P_1$ 、从 $P_2$ 到 $P_3$ 和从 $P_4$ 到 $P_{end}$ 的慢轨道在三角形 $x=1$ 中慢极限系统(12)的相图中显示；从 $P_1$ 到 $P_2$ 和从 $P_3$ 到 $P_4$ 的慢轨道在三角形 $x=0$ 中慢极限系统(11)的相图中显示。

3) $I=\frac{k_0+k_2}{k_1+k_2}=\frac{1}{10}$ ，其中 $k_0=0.4$ 、 $k_1=9.4$ 、 $k_2=0.6$ ，产生两次跳跃的奇异轨道。

$\begin{array}{l}{P}_{start}=\left(0.99,0.01,0.98\right)\\ P_1=\left(0.5966296974,0.3020886436,1\right)\\ P_2=\left(0.4905599561,0.0166601001,0\right)\\ P_{end}=\left(0.0498286735,0,1\right)\end{array}$

从图5数值模拟1 (3)的相图可以看到：从 $P_{start}$ 到 $P_1$ 和从 $P_2$ 到 $P_{end}$ 的慢轨道在三角形 $x=1$ 中慢极限系统(12)的相图中显示；从 $P_1$ 到 $P_2$ 的慢轨道在三角形 $x=0$ 中慢极限系统(11)的相图中显示。

结论1：由数值模拟1可知，在起始点不变的情况下，随着 $I=\frac{k_0+k_2}{k_1+k_2}$ 越来越大，跳跃点的数量减少。

数值模拟2： $I=\frac{k_0+k_2}{k_1+k_2}=\frac{1}{10}$ ，其中 $k_0=0.4$ 、 $k_1=9.4$ 、 $k_2=0.6$ 。

1) $P_{start}=\left(0.98,0.02,0.98\right)$ ，产生四次跳跃的奇异轨道。

(a) 三角形 $x=1$ 中慢极限系统(12)的相图。其中垂直线 (b) 三角形 $x=0$ 中慢极限系统(11)的相图。其中水平线 $S=\frac{\gamma }{{\beta }_n}=\frac{1}{5}$ ，水平线 $I=\frac{k_0+k_2}{k_1+k_2}=\frac{1}{10}$ 。 $I=\frac{k_0+k_2}{k_1+k_2}=\frac{1}{10}$ 。

Figure 5. Phase portraits of Numerical simulation 1 (3)

图5. 数值模拟1 (3)相图

$\begin{array}{l}{P}_{start}=\left(0.98,0.02,0.98\right)\\ P_1=\left(0.6743548111,0.2508859542,1\right)\\ P_2=\left(0.5726498790,0.0256270867,0\right)\\ P_3=\left(0.2533532428,0.1818257696,1\right)\\ P_4=\left(0.2346667203,0.0472754212,0\right)\\ P_{end}=\left(0.0897789174,0,1\right)\end{array}$

从图6数值模拟2 (1)的相图可以看到：从 $P_{start}$ 到 $P_1$ 、从 $P_2$ 到 $P_3$ 和从 $P_4$ 到 $P_{end}$ 的慢轨道在三角形 $x=1$ 中慢极限系统(12)的相图中显示；从 $P_1$ 到 $P_2$ 和从 $P_3$ 到 $P_4$ 的慢轨道在三角形 $x=0$ 中慢极限系统(11)的相图中显示。

2) 起始点为 $P_{start}=\left(0.96,0.04,0.98\right)$ ，产生六次跳跃的奇异轨道。

$\begin{array}{l}{P}_{start}=\left(0.96,0.04,0.98\right)\\ P_1=\left(0.7638923304,0.1904063827,1\right)\\ P_2=\left(0.6808061750,0.0431941322,0\right)\\ P_3=\left(0.4783899887,0.1750400441,1\right)\end{array}$

$\begin{array}{l}{P}_4=\left(0.4410805067,0.0499518019,0\right)\\ P_5=\left(0.2005958777,0.1328494132,1\right)\\ P_6=\left(0.1940499608,0.0730419904,0\right)\\ P_{end}=\left(0.0738347689,0,1\right)\end{array}$

(a) 三角形 $x=1$ 中慢极限系统(12)的相图。其中垂直线 (b) 三角形 $x=0$ 中慢极限系统(11)的相图。其中水平线 $S=\frac{\gamma }{{\beta }_n}=\frac{1}{5}$ ，水平线 $I=\frac{k_0+k_2}{k_1+k_2}=\frac{1}{10}$ 。 $I=\frac{k_0+k_2}{k_1+k_2}=\frac{1}{10}$ 。

Figure 6. Phase portraits of Numerical simulation 2 (1)

图6. 数值模拟2 (1)相图

(a) 三角形 $x=1$ 中慢极限系统(12)的相图。其中垂直线 (b) 三角形 $x=0$ 中慢极限系统(11)的相图。其中水平线 $S=\frac{\gamma }{{\beta }_n}=\frac{1}{5}$ ，水平线 $I=\frac{k_0+k_2}{k_1+k_2}=\frac{1}{10}$ 。 $I=\frac{k_0+k_2}{k_1+k_2}=\frac{1}{10}$ 。

Figure 7. Phase portraits of Numerical simulation 2 (2)

图7. 数值模拟2 (2)相图

从图7数值模拟2 (2)的相图可以看到：从 $P_{start}$ 到 $P_1$ 、从 $P_2$ 到 $P_3$ 、从 $P_4$ 到 $P_5$ 和从 $P_6$ 到 $P_{end}$ 的慢轨道在三角形 $x=1$ 中慢极限系统(12)的相图中显示；从 $P_1$ 到 $P_2$ 、从 $P_3$ 到 $P_4$ 和从 $P_5$ 到 $P_6$ 的慢轨道在三角形 $x=0$ 中慢极限系统(11)的相图中显示。

3) $P_{start}=\left(0.94,0.06,0.98\right)$ ，产生八次跳跃的奇异轨道。

$\begin{array}{l}{P}_{start}=\left(0.94,0.06,0.98\right)\\ P_1=\left(0.8218870678,0.1512575567,1\right)\\ P_2=\left(0.7629421725,0.0613609325,0\right)\\ P_3=\left(0.6421425442,0.1476861751,1\right)\\ P_4=\left(0.6041058394,0.0636010855,0\right)\end{array}$

$\begin{array}{l}{P}_5=\left(0.4794446750,0.1420380687,1\right)\\ P_6=\left(0.4565484184,0.0670668287,0\right)\\ P_7=\left(0.3205524806,0.1323330165,1\right)\\ P_8=\left(0.3095287465,0.0733666691,0\right)\\ P_{end}=\left(0.0623101777,0,1\right)\end{array}$

(a) 三角形 $x=1$ 中慢极限系统(12)的相图。其中垂直线 (b) 三角形 $x=0$ 中慢极限系统(11)的相图。其中水平线 $S=\frac{\gamma }{{\beta }_n}=\frac{1}{5}$ ，水平线 $I=\frac{k_0+k_2}{k_1+k_2}=\frac{1}{10}$ 。 $I=\frac{k_0+k_2}{k_1+k_2}=\frac{1}{10}$ 。

Figure 8. Phase portraits of Numerical simulation 2 (3)

图8. 数值模拟2 (3)相图

从图8数值模拟2 (3)的相图可以看到：从 $P_{start}$ 到 $P_1$ 、从 $P_2$ 到 $P_3$ 、从 $P_4$ 到 $P_5$ 、从 $P_6$ 到 $P_7$ 和从 $P_8$ 到 $P_{end}$ 的慢轨道在三角形 $x=1$ 中慢极限系统(12)的相图中显示；从 $P_1$ 到 $P_2$ 、从 $P_3$ 到 $P_4$ 、从 $P_5$ 到 $P_6$ 和从 $P_7$ 到 $P_8$ 的慢轨道在三角形 $x=0$ 中慢极限系统(11)的相图中显示。

结论2：由数值模拟2可知，在 $I=\frac{k_0+k_2}{k_1+k_2}$ 一定的情况下， $\left(S_0,I_0\right)$ 距离相空间T中的点 $\left(1,0\right)$ 越近，跳跃点的数量越少。

5. 结论

本文建立的带有个体自发行为变化的SIR模型可以产生有反复波的传染病。在某些传染病流行期间，易感个体可以选择经济成本低但感染风险高的正常行为，或者感染风险低但经济成本高的谨慎行为。当 $I>\frac{k_0+k_2}{k_1+k_2}$ 时，谨慎行为回报更高，易感人群倾向于采取谨慎行为，使得传染病传播减慢；当 $I<\frac{k_0+k_2}{k_1+k_2}$ 时，正常行为回报更高，易感人群会自发的恢复到正常行为，传染病可能就会产生反复波。在个体行为变化的时间尺度比传染病传播的时间尺度快得多的情况下，我们用几何奇异摄动理论分析带有个体自发行为变化的SIR模型，并借助快–慢结构和入–出积分得到模型的奇异轨道。对模型的奇异轨道数值模拟得到： $I=\frac{k_0+k_2}{k_1+k_2}$ 越大，跳跃点越少，疫情反复波越少；感染人群比例越大，跳跃点越多，疫情反复波越多。我们建立的模型还可以扩展到有症状和无症状的感染人群，如将无症状感染人群看作与个体行为变化有关的易感人群。还可以把潜伏期的个体引入一个仓室，从而延缓传染病传播并影响个体行为变化。

基金项目

国家自然科学基金项目(批准号11901052)。