两地区人口流动肺结核动力学模型研究

Dynamic Model of Tuberculosis in Population Migration between Two Regions

摘要: 为研究人口迁移与疫苗接种对肺结核传播的影响，建立人口密度不同的两地区之间流动的SEIR流行病动力学模型，并考虑对各地区实施疫苗接种措施，计算其有效再生数，通过数值模拟作有效再生数关于迁移率与疫苗接种率的敏感性分析。结果表明，高人口密度地区向低密度地区迁移可以降低有效再生数，缩小染病规模，低密度地区向高密度地区迁移则结果相反，但均不能通过迁移达到消除肺结核的目的，而适当的对高密度人口地区施加疫苗可以消除肺结核。

Abstract: To study the effects of population migration and vaccination on the spread of tuberculosis, a dy-namic SEIR epidemic model of the two regions with different population densities is established. The effective reproduction number is calculated, and the vaccination strategy for each region is considered in model. We numerically simulate the sensitivity of the effective reproduction number with respect to migration rate and vaccination rate. The results show that migration from areas with high population density to areas with low population density can reduce the effective re-production number and disease scale. The migration from low-density areas to high-density areas has the opposite result, but they cannot achieve the purpose of eliminating tuberculosis through migration. Appropriate application of vaccines to people in the high-density area can eliminate tuberculosis.

1. 引言

传染病是致命性微生物在人与人之间，人与动物之间，或者动物与动物之间相互传播的疾病，给人类的健康带来极大的灾难，给社会经济发展带来极大的损失 [1] 。随着我国经济发展水平与医疗卫生水平得到提高，降低了大部分传染病的病死率，部分疾病也得到控制 [2] ，但是有些疾病仍然难以控制，例如肺结核。由于环境变化、流动人口的增加等原因使得肺结核依然在部分地区流行 [3] 。进城务工人员多数集体居住，劳动强度大，生活条件相对差，更容易诱发肺结核暴发，进城工人的流动对肺结核的传播影响较大 [4] 。近年来，国内外有很多研究，通过建立数学模型来研究流动人口对肺结核传播的影响，并探究治疗，疫苗等防控措施对其传播的影响。针对意大利的人口迁移问题，周义仓 [5] 等人建立了具有移民的肺结核SEIRT模型，考虑了移民以及不同地区治疗率对肺结核传播的影响。许 [4] 等人通过建立SEIR模型来研究广东省本地人口和外来移入人口之间肺结核的传播，探究了疫苗措施对其控制的作用，得出控制广东肺结核传播的最优疫苗分配策略。靳帧 [6] 等人通过建立具有出生与死亡的人口流动SIR流感模型，研究两地区的移出率对基本再生数的影响，得出人口高密度地区向低密度地区迁移可以有效控制流感暴发的结论。

2. 模型建立及再生数的计算

为研究人口迁移与疫苗接种对肺结核传播的影响，基于靳帧 [6] 等人模型，建立以下带有疫苗项的两城市之间人口流动的肺结核SEIR动力学模型。

$\begin{array}{l}\frac{\text{d}{S}_{11}}{\text{d}t}=A_{11}-{\sigma }_{12}{S}_{11}+{\rho }_{12}{S}_{12}-S_{11}{\displaystyle \underset{j=1}{\overset{2}{\sum }}{k}_1{\beta }_{11j}}\frac{I_{j1}}{N_1^t}-{\chi }_{11}{S}_{11}-d{S}_{11}\\ \frac{\text{d}{S}_{12}}{\text{d}t}=A_{12}+{\sigma }_{12}{S}_{11}-{\rho }_{12}{S}_{12}-S_{12}{\displaystyle \underset{j=1}{\overset{2}{\sum }}{k}_2{\beta }_{12j}}\frac{I_{j2}}{N_2^t}-{\chi }_{12}{S}_{12}-d{S}_{12}\\ \frac{\text{d}{S}_{22}}{\text{d}t}=A_{22}-{\sigma }_{21}{S}_{22}+{\rho }_{21}{S}_{21}-S_{22}{\displaystyle \underset{j=1}{\overset{2}{\sum }}{k}_1{\beta }_{22j}}\frac{I_{j2}}{N_2^t}-{\chi }_{22}{S}_{22}-d{S}_{22}\\ \frac{\text{d}{S}_{21}}{\text{d}t}=A_{21}+{\sigma }_{21}{S}_{22}-{\rho }_{21}{S}_{21}-S_{21}{\displaystyle \underset{j=1}{\overset{2}{\sum }}{k}_1{\beta }_{21j}}\frac{I_{j1}}{N_1^t}-{\chi }_{21}{S}_{21}-d{S}_{21}\end{array}$

$\begin{array}{l}\frac{\text{d}{E}_{11}}{\text{d}t}=S_{11}{\displaystyle \underset{j=1}{\overset{2}{\sum }}{k}_1{\beta }_{11j}}\frac{I_{j1}}{N_1^t}-{\sigma }_{12}{E}_{11}+{\rho }_{12}{E}_{12}-\left(\alpha +d\right)E_{11}\\ \frac{\text{d}{E}_{12}}{\text{d}t}=S_{12}{\displaystyle \underset{j=1}{\overset{2}{\sum }}{k}_2{\beta }_{12j}}\frac{I_{j2}}{N_2^t}+{\sigma }_{12}{E}_{11}-{\rho }_{12}{E}_{12}-\left(\alpha +d\right)E_{12}\\ \frac{\text{d}{E}_{22}}{\text{d}t}=S_{22}{\displaystyle \underset{j=1}{\overset{2}{\sum }}{k}_1{\beta }_{22j}}\frac{I_{j2}}{N_2^t}-{\sigma }_{21}{E}_{22}+{\rho }_{21}{E}_{21}-\left(\alpha +d\right)E_{22}\\ \frac{\text{d}{E}_{21}}{\text{d}t}=S_{21}{\displaystyle \underset{j=1}{\overset{2}{\sum }}{k}_1{\beta }_{21j}}\frac{I_{j1}}{N_1^t}+{\sigma }_{21}{E}_{22}-{\rho }_{21}{E}_{21}-\left(\alpha +d\right)E_{21}\end{array}$

$\begin{array}{l}\frac{\text{d}{I}_{11}}{\text{d}t}=\alpha E_{11}-{\sigma }_{12}{I}_{11}+{\rho }_{12}{I}_{12}-\left(\gamma +d\right)I_{11}\\ \frac{\text{d}{I}_{12}}{\text{d}t}=\alpha E_{12}+{\sigma }_{12}{I}_{11}-{\rho }_{12}{I}_{12}-\left(\gamma +d\right)I_{12}\\ \frac{\text{d}{I}_{22}}{\text{d}t}=\alpha E_{22}-{\sigma }_{21}{I}_{22}+{\rho }_{21}{I}_{21}-\left(\gamma +d\right)I_{22}\\ \frac{\text{d}{I}_{21}}{\text{d}t}=\alpha E_{21}+{\sigma }_{21}{I}_{22}-{\rho }_{21}{I}_{21}-\left(\gamma +d\right)I_{21}\end{array}$

$\begin{array}{l}\frac{\text{d}{R}_{11}}{\text{d}t}=\gamma I_{11}-{\sigma }_{12}{R}_{11}+{\rho }_{12}{R}_{12}+{\chi }_{11}{S}_{11}-d{R}_{11}\\ \frac{\text{d}{R}_{12}}{\text{d}t}=\gamma I_{12}+{\sigma }_{12}{R}_{11}-{\rho }_{12}{R}_{12}+{\chi }_{12}{S}_{12}-d{R}_{12}\\ \frac{\text{d}{R}_{22}}{\text{d}t}=\gamma I_{22}-{\sigma }_{21}{R}_{22}+{\rho }_{21}{R}_{21}+{\chi }_{22}{S}_{22}-d{R}_{22}\\ \frac{\text{d}{R}_{21}}{\text{d}t}=\gamma I_{21}+{\sigma }_{21}{R}_{22}-{\rho }_{21}{R}_{21}+{\chi }_{21}{S}_{21}-d{R}_{21}\end{array}$ (1)

其中， $S_{ij},E_{ij},I_{ij},R_{ij}$ 表示分别表示t时刻i户籍迁移到j地区的易感者，潜伏者，染病者，恢复者。 $N_i^t$ 表示t时刻在i地区的总人口，t时刻两地的总人数分别为

$\begin{array}{l}{N}_1^t=S_{11}\left(t\right)+S_{21}\left(t\right)+E_{11}\left(t\right)+E_{21}\left(t\right)+I_{11}\left(t\right)+I_{21}\left(t\right)+R_{11}\left(t\right)+R_{21}\left(t\right)\\ N_2^t=S_{12}\left(t\right)+S_{22}\left(t\right)+E_{12}\left(t\right)+E_{22}\left(t\right)+I_{12}\left(t\right)+I_{22}\left(t\right)+R_{12}\left(t\right)+R_{22}(t)\end{array}$

表1对模型(1)中参数进行说明。

Table 1. The meaning of parameter in model

表1. 模型参数含义表

由模型(1)得 $\{\begin{array}{l}\frac{\text{d}{N}^1\left(t\right)}{\text{d}t}=A_{11}+A_{21}-d{N}^1\\ \frac{\text{d}{N}^2\left(t\right)}{\text{d}t}=A_{12}+A_{22}-d{N}^2\end{array}$ ，

可以推出：当时间 $t\to \infty$ 时，1，2两地区的人口总数趋于稳定，分别趋于 $\frac{A_{11}+A_{21}}{d}$ ， $\frac{A_{12}+A_{22}}{d}$ 。

需要注意的是流动的人口(即户籍地与所在地不同的人群)流动的原因分是短期旅行和因工作长期居住两种情况，当短期旅行时，流动人口就不存在出生人口，假设不存在旅行者在异地生育的情况，此时有异地生育率 $A_{12}=A_{21}=0$ ，当长期居住时，则存在异地生育，则 $A_{12},A_{21}$ 不为0。

当没有染病者且不考虑注射疫苗措施时：模型(1)转化为模型(2)

$\begin{array}{l}\frac{\text{d}{S}_{11}}{\text{d}t}=A_{11}-{\sigma }_{12}{S}_{11}+{\rho }_{12}{S}_{12}-d{S}_{11}\\ \frac{\text{d}{S}_{12}}{\text{d}t}=A_{12}+{\sigma }_{12}{S}_{11}-{\rho }_{12}{S}_{12}-d{S}_{12}\\ \frac{\text{d}{S}_{22}}{\text{d}t}=A_{22}-{\sigma }_{21}{S}_{22}+{\rho }_{21}{S}_{21}-d{S}_{22}\\ \frac{\text{d}{S}_{21}}{\text{d}t}=A_{21}+{\sigma }_{21}{S}_{22}-{\rho }_{21}{S}_{21}-d{S}_{21}\end{array}$ (2)

求得模型(2)的平衡点为： $\left(S_{011}^{*},S_{012}^{*},S_{021}^{*},S_{022}^{*},0,0,0,0,0,0,0,0,0,0,0,0\right)$ ，其中有：

$S_{011}^{*}=\frac{{\rho }_{12}\left(A_{11}+A_{12}\right)+A_{11}d}{\left({\rho }_{12}+d\right)\left({\sigma }_{12}+d\right)-{\sigma }_{12}{\rho }_{12}},S_{012}^{*}=\frac{A_{12}d+{\sigma }_{12}\left(A_{11}+A_{12}\right)}{\left({\rho }_{12}+d\right)\left({\sigma }_{12}+d\right)-{\sigma }_{12}{\rho }_{12}}$

$S_{022}^{*}=\frac{{\rho }_{21}\left(A_{21}+A_{22}\right)+A_{22}d}{\left({\rho }_{21}+d\right)\left({\sigma }_{21}+d\right)-{\sigma }_{21}{\rho }_{21}},S_{021}^{*}=\frac{A_{21}d+{\sigma }_{21}\left(A_{21}+A_{22}\right)}{\left({\rho }_{21}+d\right)\left({\sigma }_{21}+d\right)-{\sigma }_{21}{\rho }_{21}}$

此时达到平衡点时的两地区总人口分别为 $N_1^{t*}=S_{011}^{*}+S_{021}^{*}$ ， $N_2^{t*}=S_{022}^{*}+S_{012}^{*}$

在平衡点处模型(2)的系数矩阵为

$\left(\begin{array}{cccc}-{\sigma }_{12}-d& {\rho }_{12}& 0& 0\\ {\sigma }_{12}& -{\rho }_{12}-d& 0& 0\\ 0& 0& -{\sigma }_{21}-d& {\rho }_{21}\\ 0& 0& {\sigma }_{21}& -{\rho }_{21}-d\end{array}\right)$

其对应的特征方程的特征根具有负实部 [6] ，则平衡点为全局渐近稳定的。

当没有染病者且考虑注射疫苗措施时：则有模型(3)

$\begin{array}{l}\frac{\text{d}{S}_{11}}{\text{d}t}=A_{11}-{\sigma }_{12}{S}_{11}+{\rho }_{12}{S}_{12}-{\chi }_{11}{S}_{11}-d{S}_{11}\\ \frac{\text{d}{S}_{12}}{\text{d}t}=A_{12}+{\sigma }_{12}{S}_{11}-{\rho }_{12}{S}_{12}-{\chi }_{12}{S}_{12}-d{S}_{12}\\ \frac{\text{d}{S}_{22}}{\text{d}t}=A_{22}-{\sigma }_{21}{S}_{22}+{\rho }_{21}{S}_{21}-{\chi }_{22}{S}_{22}-d{S}_{22}\\ \frac{\text{d}{S}_{21}}{\text{d}t}=A_{21}+{\sigma }_{21}{S}_{22}-{\rho }_{21}{S}_{21}-{\chi }_{21}{S}_{21}-d{S}_{21}\end{array}$

$\begin{array}{l}\frac{\text{d}{R}_{11}}{\text{d}t}={\sigma }_{12}{R}_{11}+{\rho }_{12}{R}_{12}+{\chi }_{11}{S}_{11}-d{R}_{11}\\ \frac{\text{d}{R}_{12}}{\text{d}t}={\sigma }_{12}{R}_{11}-{\rho }_{12}{R}_{12}+{\chi }_{12}{S}_{12}-d{R}_{12}\\ \frac{\text{d}{R}_{22}}{\text{d}t}={\sigma }_{21}{R}_{22}+{\rho }_{21}{R}_{21}+{\chi }_{22}{S}_{22}-d{R}_{22}\\ \frac{\text{d}{R}_{21}}{\text{d}t}={\sigma }_{21}{R}_{22}-{\rho }_{21}{R}_{21}+{\chi }_{21}{S}_{21}-d{R}_{21}\end{array}$ (3)

求得系统(3)的平衡点为： $\left(S_{111}^{*},S_{112}^{*},S_{121}^{*},S_{122}^{*},0,0,0,0,0,0,0,0,R_{111}^{*},R_{112}^{*},R_{121}^{*},R_{122}^{*}\right)$ ，其中有

$S_{111}^{*}=\frac{A_{11}\left({\chi }_{12}+d\right)+{\rho }_{12}\left(A_{11}+A_{12}\right)}{\left({\rho }_{12}+{\chi }_{12}+d\right)\left({\sigma }_{12}+{\chi }_{11}+d\right)-{\sigma }_{12}{\rho }_{12}},S_{112}^{*}=\frac{A_{12}\left({\chi }_{11}+d\right)+{\sigma }_{12}\left(A_{11}+A_{12}\right)}{\left({\rho }_{12}+{\chi }_{12}+d\right)\left({\sigma }_{12}+{\chi }_{11}+d\right)-{\sigma }_{12}{\rho }_{12}}$

$S_{122}^{*}=\frac{A_{22}\left({\chi }_{21}+d\right)+{\rho }_{21}\left(A_{22}+A_{21}\right)}{\left({\rho }_{21}+{\chi }_{21}+d\right)\left({\sigma }_{21}+{\chi }_{22}+d\right)-{\sigma }_{21}{\rho }_{21}},S_{121}^{*}=\frac{A_{21}\left({\chi }_{22}+d\right)+{\sigma }_{21}\left(A_{22}+A_{21}\right)}{\left({\rho }_{21}+{\chi }_{21}+d\right)\left({\sigma }_{21}+{\chi }_{22}+d\right)-{\sigma }_{21}{\rho }_{21}}$

$R_{111}^{*}=\frac{{\rho }_{12}{\chi }_{12}{S}_{112}^{*}+\left({\rho }_{12}+d\right){\chi }_{11}{S}_{111}^{*}}{d\left({\sigma }_{12}+{\rho }_{12}+d\right)},R_{112}^{*}=\frac{{\sigma }_{12}{\chi }_{11}{S}_{111}^{*}+\left({\sigma }_{12}+d\right){\chi }_{12}{S}_{112}^{*}}{d\left({\sigma }_{12}+{\rho }_{12}+d\right)}$

$R_{122}^{*}=\frac{{\rho }_{21}{\chi }_{21}{S}_{121}^{*}+\left({\rho }_{21}+d\right){\chi }_{22}{S}_{22}^{*}}{d\left({\sigma }_{21}+{\rho }_{21}+d\right)},R_{121}^{*}=\frac{{\sigma }_{21}{\chi }_{22}{S}_{122}^{*}+\left({\sigma }_{21}+d\right){\chi }_{21}{S}_{121}^{*}}{d\left({\sigma }_{21}+{\rho }_{21}+d\right)}$

此时达到平衡点时的两地区总人口分别为

$N_1^{t*}=S_{111}^{*}+S_{121}^{*}\text{+}{R}_{111}^{*}+R_{121}^{*},N_2^{t*}=S_{222}^{*}+S_{212}^{*}+R_{222}^{*}+R_{212}^{*}$

利用文献 [8] 中再生矩阵法求得有效再生数 $R_v$ ，由模型(1)得

$\mathcal{F}=\left(\begin{array}{c}{S}_{11}{\displaystyle \underset{j=1}{\overset{2}{\sum }}{k}_1{\beta }_{11j}}\frac{I_{j1}}{N_1^t}\\ S_{12}{\displaystyle \underset{j=1}{\overset{2}{\sum }}{k}_2{\beta }_{12j}}\frac{I_{j2}}{N_2^t}\\ S_{22}{\displaystyle \underset{j=1}{\overset{2}{\sum }}{k}_1{\beta }_{22j}}\frac{I_{j2}}{N_2^t}\\ S_{21}{\displaystyle \underset{j=1}{\overset{2}{\sum }}{k}_1{\beta }_{21j}}\frac{I_{j1}}{N_1^t}\\ 0\\ 0\\ 0\\ 0\end{array}\right)$ , $\mathcal{V}=\left(\begin{array}{c}{\sigma }_{12}{E}_{11}-{\rho }_{12}{E}_{12}+\left(\alpha +d\right)E_{11}\\ -{\sigma }_{12}{E}_{11}+{\rho }_{12}{E}_{12}+\left(\alpha +d\right)E_{12}\\ {\sigma }_{21}{E}_{22}-{\rho }_{21}{E}_{21}+\left(\alpha +d\right)E_{22}\\ -{\sigma }_{21}{E}_{22}+{\rho }_{21}{E}_{21}+\left(\alpha +d\right)E_{21}\\ -\alpha E_{11}+{\sigma }_{12}{I}_{11}-{\rho }_{12}{I}_{12}+\left(\gamma +d\right)I_{11}\\ -\alpha E_{12}-{\sigma }_{12}{I}_{11}+{\rho }_{12}{I}_{12}+\left(\gamma +d\right)I_{12}\\ -\alpha E_{22}+{\sigma }_{21}{I}_{22}-{\rho }_{21}{I}_{21}+\left(\gamma +d\right)I_{22}\\ -\alpha E_{21}-{\sigma }_{21}{I}_{22}+{\rho }_{21}{I}_{21}+\left(\gamma +d\right)I_{21}\end{array}\right)$

令

$\begin{array}{l}{m}_{12}={\sigma }_{12}+\alpha +d,n_{12}={\rho }_{12}+\alpha +d,m_{21}={\sigma }_{21}+\alpha +d,n_{21}={\rho }_{21}+\alpha +d,\\ p_{12}={\sigma }_{12}+r+d,q_{12}={\rho }_{12}+r+d,p_{21}={\sigma }_{21}+r+d,q_{21}={\rho }_{21}+r+d\end{array}$

$\mathcal{F},\mathcal{V}$ 分别关于潜伏项 $E_{ij}$ 与染病项 $I_{ij}$ 求导得 $F,V$

$F=\left(\begin{array}{cc}0& F_{12}\\ 0& 0\end{array}\right)$ ，其中 $F_{12}=\left(\begin{array}{cccc}\frac{S_{11}^{*}{k}_1{\beta }_{111}}{N_1^{t*}}& 0& 0& \frac{S_{11}^{*}{k}_1{\beta }_{112}}{N_1^{t*}}\\ 0& \frac{S_{12}^{*}{k}_2{\beta }_{121}}{N_2^{t*}}& \frac{S_{12}^{*}{k}_2{\beta }_{122}}{N_2^{t*}}& 0\\ 0& \frac{S_{22}^{*}{k}_2{\beta }_{221}}{N_2^{t*}}& \frac{S_{22}^{*}{k}_2{\beta }_{222}}{N_2^{t*}}& 0\\ \frac{S_{21}^{*}{k}_1{\beta }_{211}}{N_1^{t*}}& 0& 0& \frac{S_{21}^{*}{k}_1{\beta }_{212}}{N_1^{t*}}\end{array}\right)$

$V=\left(\begin{array}{cc}{V}_{11}& 0\\ V_{21}& V_{22}\end{array}\right)$ ，其中 $$

$V_{21}=\left(\begin{array}{cccc}-\alpha & 0& 0& 0\\ 0& -\alpha & 0& 0\\ 0& 0& -\alpha & 0\\ 0& 0& 0& -\alpha \end{array}\right)$ , $V_{22}=\left(\begin{array}{cccc}{p}_{12}& -{\rho }_{12}& 0& 0\\ -{\sigma }_{12}& q_{12}& 0& 0\\ 0& 0& p_{21}& -{\rho }_{21}\\ 0& 0& -{\sigma }_{21}& q_{21}\end{array}\right)$

求逆得 $V^{-1}=\left(\begin{array}{cc}{V}_{11}^{-1}& 0\\ -V_{22}^{-1}{V}_{21}{V}_{11}^{-1}& V_{22}^{-1}\end{array}\right)$ 其中：

$V_{11}^{-1}=\left(\begin{array}{cccc}\frac{n_{12}}{m_{12}{n}_{12}-{\rho }_{12}{\sigma }_{12}}& \frac{{\rho }_{12}}{m_{12}{n}_{12}-{\rho }_{12}{\sigma }_{12}}& 0& 0\\ \frac{{\sigma }_{12}}{m_{12}{n}_{12}-{\rho }_{12}{\sigma }_{12}}& \frac{m_{12}}{m_{12}{n}_{12}-{\rho }_{12}{\sigma }_{12}}& 0& 0\\ 0& 0& \frac{n_{21}}{m_{21}{n}_{21}-{\rho }_{21}{\sigma }_{21}}& \frac{{\rho }_{21}}{m_{21}{n}_{21}-{\rho }_{21}{\sigma }_{21}}\\ 0& 0& \frac{{\sigma }_{21}}{m_{21}{n}_{21}-{\rho }_{21}{\sigma }_{21}}& \frac{m_{12}}{m_{21}{n}_{21}-{\rho }_{21}{\sigma }_{21}}\end{array}\right)$

$V_{22}^{-1}=\left(\begin{array}{cccc}\frac{q_{12}}{p_{12}{q}_{12}-{\rho }_{12}{\sigma }_{12}}& \frac{{\rho }_{12}}{p_{12}{q}_{12}-{\rho }_{12}{\sigma }_{12}}& 0& 0\\ \frac{{\sigma }_{12}}{p_{12}{q}_{12}-{\rho }_{12}{\sigma }_{12}}& \frac{p_{12}}{p_{12}{q}_{12}-{\rho }_{12}{\sigma }_{12}}& 0& 0\\ 0& 0& \frac{q_{21}}{p_{21}{q}_{21}-{\rho }_{21}{\sigma }_{21}}& \frac{{\rho }_{21}}{p_{21}{n}_{21}-{\rho }_{21}{\sigma }_{21}}\\ 0& 0& \frac{{\sigma }_{21}}{p_{21}{q}_{21}-{\rho }_{21}{\sigma }_{21}}& \frac{p_{21}}{p_{21}{q}_{21}-{\rho }_{21}{\sigma }_{21}}\end{array}\right)$

令 $U_{12}=m_{12}{n}_{12}-{\rho }_{12}{\sigma }_{12},U_{21}=m_{21}{n}_{21}-{\rho }_{21}{\sigma }_{21},W_{12}=p_{12}{q}_{12}-{\rho }_{12}{\sigma }_{12},W_{21}=p_{21}{q}_{21}-{\rho }_{21}{\sigma }_{21}$ ，有

$-V_{22}^{-1}{V}_{21}{V}_{11}^{-1}=\left(\begin{array}{cccc}\frac{\alpha \left(n_{12}{q}_{12}+{\sigma }_{12}{\rho }_{12}\right)}{W_{12}{U}_{12}}& \frac{\alpha \left({\rho }_{12}{q}_{12}+m_{12}{\rho }_{12}\right)}{W_{12}{U}_{12}}& 0& 0\\ \frac{\alpha \left({\sigma }_{12}{\eta }_{12}+{\sigma }_{12}{p}_{12}\right)}{W_{12}{U}_{12}}& \frac{\alpha \left({\sigma }_{12}{\rho }_{12}+p_{12}{m}_{12}\right)}{W_{12}{U}_{12}}& 0& 0\\ 0& 0& \frac{\alpha \left(n_{21}{q}_{21}+{\sigma }_{21}{\rho }_{21}\right)}{W_{21}{U}_{21}}& \frac{\alpha \left({\rho }_{21}{q}_{21}+m_{21}{\rho }_{21}\right)}{W_{21}{U}_{21}}\\ 0& 0& \frac{\alpha \left({\sigma }_{21}{\eta }_{21}+{\sigma }_{21}{p}_{21}\right)}{W_{21}{U}_{21}}& \frac{\alpha \left({\sigma }_{21}{\rho }_{21}+p_{21}{m}_{21}\right)}{W_{21}{U}_{12}}\end{array}\right)$

令 $K_{11}=-F_{12}{V}_{22}^{-1}{V}_{21}{V}_{11}^{-1}$ ，则有

再生矩阵为 $K=F{V}^{-1}$ ，得 $K=\left(\begin{array}{cc}{K}_{11}& \ast \\ 0& 0\end{array}\right)$ ，由于*不影响K的特征值 [9] ，则有 $\rho \left(K\right)=\rho \left(K_{11}\right)$ ，因此有效再生数 $R_v=\rho \left(K_{11}\right)$ 。

3. 数值模拟

在本节中，我们研究两个地区之间人口流动与疫苗注射对肺结核传播的影响，对参数(迁移率与疫苗率)进行敏感性分析。我们选取人口密度较高的北京市和人口密度较低的呼和浩特市作为要研究的地区，即城市1为北京，城市2为呼和浩特。由统计局调查数据显示2017年北京人口出生率为473.45人/天，呼和浩特出生率为68.27人/天。模拟时参数见表2。

Table 2. The value of parameter in simulation

表2. 模拟参数取值表

注：模拟中取单位时间为月。

首先，探究迁移率对染病规模与基本再生数的影响。

Figure 1. Trend graph of total infected I with time at different migration rates ${\sigma }_{12}$

图1. 不同迁移率 ${\sigma }_{12}$ 下总染病量I时间序列图

染病规模用总染病量来表示。总染病量为 $I=I_{11}+I_{22}+I_{12}+I_{21}$ ，假设北京市每天人均接触的人数 $k_1$ 为12，呼和浩特市人均接触人数 $k_2$ 为6，从呼和浩特迁移到北京的迁移率 ${\sigma }_{21}=0.001$ ，另外，假设模型为短期旅行，则不存在外地户籍的新生儿，即 $A_{12}=A_{21}=0$ ，另外，不考虑两地感染率的差异，即 ${\beta }_{ijk}=0.05$ ，此时不考虑采取注射疫苗策略，即 ${\chi }_{11}={\chi }_{22}={\chi }_{12}={\chi }_{21}=0$ 。观察图1得到随着北京到呼和浩特的迁移率 ${\sigma }_{12}$ 增加，染病规模和染病峰值逐渐降低，并且染病峰值出现的时间逐渐延后，但是不能完全控制出染病量，持续在染病量大于零的状态，即不能消除肺结核。进一步说明，高密度地区人口向低密度地区的迁移会减少肺结核的染病规模，但难以消除肺结核。

Figure 2. Trend graph of total infected I with time at different migration rates ${\sigma }_{21}$

图2. 不同迁移率 ${\sigma }_{21}$ 下总染病量I时间序列图

与图1不同的是，假设 ${\sigma }_{12}=0.001$ ，其他条件相同。由图2可得出，随着由呼和浩特向北京迁移的迁移率 ${\sigma }_{21}$ 增加，染病规模与染病峰值逐渐升高，染病峰值出现的时间逐渐所提前，由此得出低人口密度地区人口向高密度地区迁移，会导致扩大肺结核的染病规模，使得其更加难以控制。

Figure 3. Graph of relation between ${\sigma }_{12}$ and $R_0$

图3. $R_0$ 与 ${\sigma }_{12}$ 关系图

Figure 4. Graph of relation between ${\sigma }_{21}$ and $R_0$

图4. $R_0$ 与 ${\sigma }_{21}$ 关系图

当不考虑任何控制措施得出的再生数为基本再生数 $R_0$ ，图3与图4描述的为不同迁移率下 $R_0$ ，由图3得到随着 ${\sigma }_{12}$ 增加 $R_0$ 先快速下降后趋于平稳，使得 $R_0$ 保持在4.5左右，由图3得到随着 ${\sigma }_{21}$ 增加 $R_0$ 逐渐增加后趋于平稳，使得 $R_0$ 保持在8.42左右。进一步得出高密度人口向低密度人口迁移会让肺结核基本再生数下降，低密度人口向高密度迁移会导致相反的结果，但是两种迁移均不能使得 $R_0<1$ ，即均不能控制肺结核的传播。

接下来考虑注射疫苗措施对染病规模与有效再生数的影响。假设两地区之间存在一定的人口流动，取各地各户籍的迁移率与返回率为。考虑只对人口密度较高的城市1的户籍1 的人口注射疫苗的情况，得到图5与图6。

当考虑控制措施时求得的基本再生数则被称为有效再生数 $R_v$ ，由图5得出对人口密度较高的城市1户籍1的人采取注射疫苗策略可以大幅度降低染病规模，肺结核趋于消亡，使得总染病量趋于0。由图6得出对城市1户籍1的人采取注射疫苗策略时，由6.85下降到1以下，同样得出可以通过控制密度较高的地区的本户籍人口注射疫苗能够达到使肺结核消亡的效果。

Figure 5. Graph of of relation between ${\chi }_{11}$ and total infected I

图5. ${\chi }_{11}$ 与总染病量I关系图

Figure 6. Graph of relation between ${\chi }_{11}$ and $R_v$

图6. $R_v$ 与 ${\chi }_{11}$ 关系图

Figure 7. Graph of relationship between vaccine rate in each group and $R_v$

图7. 各组疫苗率与 $R_v$ 关系图

图7所描述的是分别对所在两地区的两户籍人分别注射一定剂量的疫苗率时有效再生数的变化，在所用参数和假设条件与前面相同的条件下，城市1的1户籍疫苗接种率对再生数影响最大，可以通过对单独对该组的人注射一定剂量的疫苗，使得 $R_v$ 小于1，从而控制肺结核的传播。其次是在城市2户籍1的人群疫苗接种率对再生数影响较大，但始终使得 $R_v$ 保持大于1，不同通过单独对城市2户籍1的人群疫苗接种的方式控制肺结核的传播。另外对城市2户籍2与城市1户籍1的人群疫苗接种率对 $R_v$ 的减少影响不大，进一步说明可以通过对高密度地区的人群施加疫苗措施达到消除肺结核的目的。

4. 结论

本文研究了考虑疫苗策略的高密度城市与低密度城市之间人口流动的SEIR肺结核模型，得出该模型的有效再生数，通过调查数据与假设条件对模型进行数值模拟分析，通过做染病规模与有效再生数对迁移率与疫苗率的敏感性分析得出由高人口密度地区向低人口密度地区能够大幅度降低染病规模，推迟染病高峰期，而由低人口密度地区向高人口密度地区迁移会扩大染病规模，使染病高峰期提前，但两种迁移方式均不能使得有效再生数小于1，所以均不能使肺结核消亡。另外得出对高密度地区的人口采取疫苗措施可以使得有效再生数小于1，可达到使肺结核逐渐消亡的目的。

致 谢

感谢许传青老师与崔景安老师的耐心指导！

基金项目

国家自然科学基金项目(11871093, 11871093)，北京市教委科技面上项目(KM2016100116018)，北京建筑大学教学实践类项目(J1703)，市属高校基本科研业务费专项基金(X18017, X18080)以及研究生创新项目(PG2018095)。