根据世界卫生组织的最新研究, 传染病依旧是人类死亡的第一杀手.由于对传染病的研究不宜采用实验的形式, 因此理论分析与数值模拟常被用于传染病的机理研究上.此时选择合适的传染病的动力学模型显得十分重要, 常见的传染病模型有SI, SIR, SEIR等等.

在传染病的动力学研究中, 许多学者做出有意义的结果, 特别是传染病模型的空间斑图动力学^[1-4].孙桂全, 靳祯等人研究了一类带有时滞的SIR空间传染病模型的Turing失稳, 通过数值模拟得出了条状与点状共存的斑图^[1].王玮明等人研究了一类SI传染病模型, 通过推导模型的振幅方程做出斑图选择, 得到了不同类型的Turing斑图^[2].马知恩和周义仓在研究中提出了这样的一个比例依赖型的模型^[5]

$ \left\{ \begin{array}{l} \frac{dS}{dt}=\gamma (1-\frac{S}{K})-\frac{\beta SI}{S+I}+cI, \\ \frac{dI}{dt}=\frac{\beta SI}{S+I}-(\mu+d+c)I(t), \\ \frac{dR}{dt}=\mu I-dR, \end{array} \right. $

(1.1)

由于$R$与$S$, $I$无关, 所以作者利用稳定性理论讨论了下面一个微分方程组的平衡点的存在性与稳定性

$ \left\{ \begin{array}{l} \frac{dS}{dt}=\gamma (1-\frac{S}{K})-\frac{\beta SI}{S+I}+cI, \\ \frac{dI}{dt}=\frac{\beta SI}{S+I}-(\mu+d+c)I(t). \end{array} \right. $

(1.2)

考虑到疾病对易感者与感染者心理上的影响, 本文在此基础上研究了一类带有负交叉扩散项的SI传染病模型

$ \left\{ \begin{array}{l} \frac{\partial S(t)}{\partial t}=\gamma S(t)(1-\frac{S(t)}{K})-\frac{\beta S(t)I(t)}{\alpha_{1}S(t)+\alpha_{2}I(t)}+cI(t)-\mu S(t)+d_{11}\Delta S(t)+d_{12}\Delta I(t), \\ \frac{\partial I(t)}{\partial t}=\frac{\beta S(t)I(t)}{\alpha_{1}S(t)+\alpha_{2}I(t)}-(\mu+d+c)I(t)+d_{21}\Delta S(t)+d_{22}\Delta I(t), \\ \frac{\partial S}{\partial \textbf{n}}=\frac{\partial I}{\partial \textbf{n}}=0, ~~ \mbox{在} \ \partial\Omega\times(0, \infty), \end{array} \right. $

(1.3)

在这里$S(t)$和$I(t)$分别代表时间$t$时易感者和感染者的密度, $\gamma$和$K$分别代表了内禀增长率和环境的承载力, $\mu$代表人口的自然死亡率, $d$代表了因病死亡率, $c$代表了染病者的恢复率, $\beta$是传染病系数, $\frac{\beta}{\alpha_{1}S(t)+\alpha_{2}I(t)}$是传染率, $\alpha_{1}$和$\alpha_{2}$心理影响系数, 并且$\alpha_{1}$和$\alpha_{2}$都是常数, $\Delta$代表空间$\Omega$中的Laplace算子, Neumann边界条件表明模型(1.3)是自我封闭的, 并且是零流量的, $\textbf{n}$代表着光滑边界$\partial\Omega$上的单位外法向量.

这里需要指出的是, 在之前的文献[1-4]中, 研究的传染病模型所带的交叉扩散系数都是正数, 其生物学意义是人群总是从另一人群的高密度区域向低密度区域移动^[6-7].而在实际生活中, 人群从低密度区域向高密度区域移动的现象也是存在的.一方面, 考虑到在生活中易感者有辨别染病者的能力并且会远离染病者, 同时染病者也会远离易感者^[6].另一方面, 在疾病爆发初期, 由于人们的心理因素, 觉得人多的地方就是安全的地方, 染病者反而会尽量接近易感者.所以在某种特定的情况下, 这种染病者向易感者移动的现象在模型中则表现为交叉扩散系数为负数.据我们所知, 带有负交叉扩散系数的传染病模型的Turing斑图在生物模型中很少被研究.因此本文将研究带有负交叉扩散系数的二维模型中的Turing斑图的生成问题.

首先考虑模型(1.3)

$ \left\{ \begin{array}{l} \frac{\partial S(t)}{\partial t}=\gamma S(t)(1-\frac{S(t)}{K})-\frac{\beta S(t)I(t)}{\alpha_{1}S(t)+\alpha_{2}I(t)}+cI(t)-\mu S(t)+d_{11}\Delta S(t)+d_{12}\Delta I(t), \\ \frac{\partial I(t)}{\partial t}=\frac{\beta S(t)I(t)}{\alpha_{1}S(t)+\alpha_{2}I(t)}-(\mu+d+c)I(t)+d_{21}\Delta S(t)+d_{22}\Delta I(t), \\ \frac{\partial S}{\partial \textbf{n}}=\frac{\partial I}{\partial \textbf{n}}=0, ~~\mbox{在} \ \partial\Omega\times(0, \infty). \end{array} \right.$

易知系统(1.3)有很多平衡点, 包括$E_0=(0, 0)$, 稳定节点$E_1=(K(1-\frac{\mu}{\gamma}), 0)$, 和正平衡点$E^{\ast}=(S^{\ast}, I^{\ast})$, 从生物学上考虑, 正平衡点更加有讨论的意义, 其中

$ \begin{eqnarray} && S^{\ast}=\frac{K(-\beta d-\beta \mu+\mu^{2} \alpha_1+d^{2}\alpha_1-\mu^{2}\alpha_2+2\mu d\alpha_1+c\mu\alpha_1+cd\alpha_1+\gamma \alpha_2\mu+\gamma\alpha_2 d+\gamma\alpha_2c-\mu\alpha_2d-\mu\alpha_2 c)}{\gamma \alpha_2(u+d+c)}, \\ && I^{\ast}=-K(-\beta d-\beta \mu+\mu^{2} \alpha_1+d^{2} \alpha_1-\mu^{2} \alpha_2+2 \mu d \alpha_1+c \mu \alpha_1+c d \alpha_1+\gamma \alpha_2 \mu+\gamma \alpha_2 d+\gamma \alpha_2 c-\mu \alpha_2 d-\mu \alpha_2 c)\\ &&~~~~~~~~~~ (-\beta+\mu \alpha_1+d \alpha_1+c \alpha_1)/ {\gamma \alpha_2^{2} (u+d+c)^{2}}, \\ && (\mu+d+c) \alpha_{1}<\beta<\frac{[(\mu+d) \alpha_{1}+(\gamma-\mu) \alpha_{2}](\mu+d+c)}{\mu+d}. \end{eqnarray} $

定理 2.1  不带有扩散项的系统(1.3)的平衡点$(S^{\ast}, I^{\ast})$是局部渐近稳定的.

证  为了讨论方便, 记

$ \left\{ \begin{array}{l} f(S, I)=\gamma S(t)(1-\frac{S(t)}{K})-\frac{\beta S(t)I(t)}{\alpha_{1}S(t)+\alpha_{2}I(t)}+cI(t)-\mu S(t), \\ g(S, I)=\frac{\beta S(t)I(t)}{\alpha_{1}S(t)+\alpha_{2}I(t)}-(\mu+d+c)I(t).\\ \end{array} \right.$

设

$ J=\left( \begin{matrix} a_{11}&a_{12}\\ a_{21}&a_{22} \end{matrix} \right), $

(2.1)

其中

$ \begin{align*} a_{11}=\frac{\partial f}{\partial S}|_{E^{\ast}}= &-(d^{3} \alpha_{1}^{2}+\mu^{3} \alpha_{1}^{2}-\beta \mu^{2} \alpha_{2}+c^{3} \alpha_{1} ^{2}-\beta \mu \alpha_{2} d-2 \beta c \mu \alpha_{1}-2 \beta c d\alpha_{1}-\beta \mu \alpha_{2} c \\ &+\beta \gamma \alpha_{2} \mu+\beta \gamma \alpha_{2} d+\beta \gamma \alpha_{2} c+6 \mu \alpha_{1}^{2}dc-2 \beta c^{2} \alpha_{1}+3 \mu^{2} \alpha_{1}^{2} d+3 \mu^{2} \alpha_{1}^{2} c+ 3 \mu \alpha_{1}^{2} d^{2}\\ &+3 \mu \alpha_{1}^{2} c^{2}+3 d^{2} \alpha_{1}^{2} c+3 d \alpha_{1}^{2} c^{2}-\beta^{2} \mu-\beta^{2} d+\beta^{2} c)\frac{1}{\alpha_{2} \beta (\mu+d+c)}, \\ a_{12}=\frac{\partial f}{\partial I}|_{E^{\ast}}=&-(-c \beta+\mu^{2} \alpha_{1}+2 \mu d \alpha_{1}+2 c \mu \alpha_{1}+d^{2} \alpha_{1}+2 c d \alpha_{1}+c^{2} \alpha_{1})\frac{1}{\beta}, \\ a_{21}=\frac{\partial g}{\partial S}|_{E^{\ast}}=&(\mu \alpha_{1}+d \alpha_{1}-\beta+c \alpha_{1})^{2}\frac{1}{\beta \alpha_{2}}, \\ a_{22}=\frac{\partial g}{\partial I}|_{E^{\ast}}=&(\mu^{2} \alpha_{1}+2 \mu d \alpha_{1}+2 c \mu \alpha_{1}+d^{2} \alpha_{1}+2 c d \alpha_{1}-\beta \mu-\beta d-c \beta+c^{2} \alpha_{1})\frac{1}{\beta}. \end{align*} $

$J$的特征方程是

$ \mid \lambda E-J\mid= \left|\begin{array}{cccc} \lambda-a_{11}& &-a_{12} \\ -a_{21}& &\lambda-a_{22} \end{array}\right|=\lambda^{2}-(a_{11}+a_{22})\lambda^{2}-\lambda \hbox{tr}_{0}+\hbox{det}J=0, $

其中${\hbox{tr}}_{0}=a_{11}+a_{22}, $ ${\hbox{det}}J=a_{11}a_{22}-a_{12}a_{21}.$易知${\hbox{tr}}_{0}<0$和${\hbox{det}}J>0$当且仅当$\max\{t_{1}, t_{3}\}<\beta<\min\{t_{2}, t_{4}\}$, 其中

$ \begin{align*} t_{1}= &\frac{(\mu+d+c)}{2(\mu-c+d)} \{\alpha_{2}(d+c+\gamma)-2\alpha_{1}c-[2\alpha_{2}^{2}(dc+dr+cr)+\alpha_{2}^{2}(d^{2}+c^{2}+\gamma^{2})\\ &+8d\mu(\alpha_{1}^{2}-\alpha_{1}\alpha_{2})+4\alpha_{1}^{2}(\mu^{2}+d^{2})-4\alpha_{1}\alpha_{2}(dc+rc+\mu^{2}+d^{2})]^{\frac{1}{2}}\}, \\ t_{2}= &\frac{(\mu+d+c)}{2(\mu-c+d)} \{\alpha_{2}(d+c+\gamma)-2\alpha_{1}c+[2\alpha_{2}^{2}(dc+dr+cr)+\alpha_{2}^{2}(d^{2}+c^{2}+\gamma^{2})\\ &+8d\mu(\alpha_{1}^{2}-\alpha_{1}\alpha_{2})+4\alpha_{1}^{2}(\mu^{2}+d^{2})-4\alpha_{1}\alpha_{2}(dc+rc+\mu^{2}+d^{2})]^{\frac{1}{2}}\}, \\ t_{3}=&\alpha_{1}(\mu+d+c), \\ t_{4}=&\frac{[(\mu+d)\alpha_{1}+(\gamma-\mu)\alpha_{2}](\mu+d+c)}{\mu+d}. \end{align*} $

通过Routh-Hurwitz, 可知$(S^{\ast}, I^{\ast})$是局部渐进稳定的.

现在考虑系统(1.3), 并且对平衡点$(S^{\ast}, I^{\ast})$进行线性化分析.如下面所示, 在平衡点$(S^{\ast}, I^{\ast})$处作微扰:

$ \left\{ \begin{array}{l} S(x, y, t)=S^{\ast}+\varepsilon S(x, y, t), |\varepsilon S(x, y, t)|\ll S^{\ast}, \\ I(x, y, t)=I^{\ast}+\varepsilon I(x, y, t), |\varepsilon I(x, y, t)|\ll I^{\ast}, \end{array} \right. $

(2.2)

其中

$\varepsilon S(x, y, t)=S_{0}e^{\lambda t}e^{i(k_{x}x+k_{y}y)}, \varepsilon I(x, y, t)=I_{0}e^{\lambda t}e^{i(k_{x}x+k_{y}y)}, $

这里$\lambda$是在时间$t$上的扰动增长率; $k_{x}$和$k_{y}$是相应的振幅; $i$是虚数单位并且有$i^{2}=-1$; $k=\sqrt{k_{x}^{2}+k_{y}^{2}}$是波数; $S_{0}$和$I_{0}$是两个正常数.把(2.2)式带入系统(1.3), 并且省略所有的非线性项, 可以得到特征方程

$|J-k^{2}D-\lambda I|=0, $

(2.3)

其中

$I=\left( \begin{matrix} 1&0\\ 0&1 \end{matrix} \right), D=\left( \begin{matrix} d_{11}&d_{12}\\ d_{21}&d_{22} \end{matrix} \right). $

特征方程(2.3)的解为如下形式

$ \lambda_{k}=\frac{{\hbox{tr}}_{k}\pm\sqrt{{\hbox{tr}}_{k}^{2}-4\delta_{k}}}{2}, $

(2.4)

其中

$ \begin{eqnarray} {\hbox{tr}}_{k}=a_{11}+a_{22}-k^{2}(d_{11}+d_{22}), \end{eqnarray} $

(2.5)

$ \delta_{k}=(d_{11}d_{22}-d_{12}d_{21})k^{4}-(d_{22}a_{11}+d_{11}a_{22}-d_{12}a_{21}-d_{21}a_{12})k^{2}+{\hbox{det}}J. $

(2.6)

选择$\beta$作为分支参数.当Im$(\lambda_{k})\neq0$和Re$(\lambda_{k})=0$在$k=0$时成立, 系统出现Hopf分支, 这样能得到Hopf分支曲线

$ a_{11}+a_{22}=0. $

(2.7)

利用稳定性定理^[9-11], 可知道当Im$(\lambda_{k})=0$和Re$(\lambda_{k})=0$在$k=k_{T}\neq0$时成立, Turing分支出现, 且波数$k_{_{T}}$满足$k_{T}^{2}=\sqrt{\frac{a_{11}a_{22}-a_{12}a_{21}}{d_{11}d_{22}-d_{12}d_{21}}}.$因此分支参数$\beta_{T}$满足如下Turing分支曲线

$ a_{11}d_{22}+a_{22}d_{11}-a_{12}d_{21}-a_{21}d_{12}-2\sqrt{(d_{11}d_{22}-d_{12}d_{21})\hbox{det}J}=0. $

(2.8)

根据Hopf和Turing分支曲线^{[12, 13]}, 能得到Hopf分支区域和Turing不稳定区域.

在图 1中, 可以看到系统(1.3)的分支图包含了Turing分支线和Hopf分支线, 并且它们把$\gamma-\beta$参数空间分成了四个区域, 区域$D_{11}$被称为Turing空间, 在这里发生Turing失稳, 区域$D_{12}$被称为Hopf空间, 在这里发生Hopf失稳.

为了更好地理解参数对系统稳定性的影响作用, 在图 2, 给出了随参数$d_{21}$变化的色散关系图.线(3)对应着Turing临界值$d_{21}=-1.13$, 当$d_{21}=-0.9>-1.13$时, Turing失稳发生; 当$d_{21}=-1.2<-1.13$时, Turing失稳消失, 也就是说, 此刻系统(1.3)的平衡解为稳定态.

在这一部分, 我们将通过Matlab对系统(1.3)进行一系列的数值模拟.所有的数值模拟均运用齐次Neumann边界条件.将空间区域离散为$200\times200$个格子.对空间的离散采用有限差分法, 设定空间步长为$\Delta h=0.25$, 对时间的离散采用欧拉方法, 取定时间步长为$\Delta t=0.01$.

首先设$d_{11}=0.02$, $d_{22}=5$, $d_{12}=0.1$, $\gamma=0.2$, $\mu=0.12$, $\beta=0.1380$, $d=0.08$, $c=0.04, $ $\alpha_{1}=0.4$, $\alpha_{2}=0.5$.现在研究参数$d_{21}$的不同值所产生的斑图.

在图 3中, $d_{21}=-0.05$, 这时可见: (A)中颜色条数值基本不变, 初值选取为平衡解加上一个随机扰动; (B)中出现类条状斑图; (C)中出现条状斑图; (D)条状斑图几乎占据了整个区域, 且系统的动力学行为不再发生变化.

图 4, 图 5分别是$d_{21}=-0.7$和$d_{21}=-0.9$时, 染病者的时间演化图.由图 4和图 5可见:随着时间的演化, 最终点状斑图和条状斑图共存.但图 4中条状斑图占优; 而当$d_{21}$达到$-0.9$时, 点状斑图会占优(图 5(D)).由图 6可见, 当$d_{21}$增至$-1.1$时, 最终点状斑图几乎占满整个空间.

本文研究了在Neumann边界条件下, 负交叉扩散对带有非线性传染率的传染病模型的影响.具体表现为负交叉扩散可引起系统(1.3)在平衡点$E^{\ast}$处的Turing失稳, 并由此得到了不同类型的斑图, 包括点状斑图、条形斑图和点条混合斑图.