小白零基础学数学建模应用系列(六）：基于数学建模的疟疾传播与控制研究

本篇文章为深度研究内容，实际给大家展示数学建模在医学中，生活中的应用。数学建模是一门很严谨的学科，大家通过本篇文章可以学习下相关背景即可，当然我在文末配了参考代码供大家学习。无论在数学建模竞赛中，还是在做科研，数学建模都会对大家有着非常大的帮助。当然，不用完全看懂本篇文章。

CSDN/B站/知乎：川川菜鸟
书籍推荐：https://item.m.jd.com/product/10099450547669.html
我的微信（接高教杯辅导咨询，数模论文转化）：hxgsrubxjogxeeag

1. 摘要

疟疾是一种危及生命的疾病，在许多非洲国家流行，尤其是尼日利亚。数学模型用于研究尼日利亚疟疾的动态。该模型结合了耐药性、治疗和使用蚊帐作为预防策略。通过将模型与尼日利亚疟疾发病率数据拟合，估计了与疾病动态相关的重要参数。使用这些估计参数，计算基本传染数，模拟未来动态，并确定对尼日利亚疟疾有重大影响的参数。总体而言，结果表明，除非更好的控制措施集中在主要耐药菌株上，改善治疗并广泛使用蚊帐，否则该疾病可能会在尼日利亚继续流行。

关键词：疟疾动力学， 稳定性分析， 基本繁殖数， 模型拟合， 敏感性分析

2. 引言

疟疾是一种由疟原虫通过受感染的雌性按蚊叮咬传播给人类的严重疾病。全球有四种疟原虫感染人类，其中恶性疟原虫和间日疟原虫是最致命的。尤其在非洲，恶性疟原虫是主要的疟疾寄生虫。未及时治疗的恶性疟原虫感染可能在24小时内导致严重疾病和死亡。

根据世界卫生组织（WHO），疟疾的病例和死亡人数逐年增加，尤其在非洲，2020年全球疟疾病例约为2.41亿，死亡人数达62.7万。非洲承载了全球约95%的疟疾病例和96%的死亡病例，尤其是5岁以下的儿童受影响最严重。尼日利亚目前是受疟疾影响最严重的国家之一，约占全球疟疾死亡人数的31.9%。

尽管疟疾在尼日利亚已经流行多年，但治疗失败仍然是控制疟疾的主要障碍之一。药物抵抗和假药的使用是疟疾治疗失败的主要原因。此外，尼日利亚国家食品药品监督管理局（NAFDAC）发现，假药也导致了许多疾病的增加。

数学模型在理解传染病传播方面起到了重要作用。这些模型被用于分析疟疾的传播和控制，其中包括药物抵抗和防控措施的影响。本研究的目的在于开发一个整合了药物抵抗和防控措施的数学模型，并结合真实数据，以预测尼日利亚疟疾传播的可能动态。这项研究的发现将帮助研究人员和政策制定者制定更有效的疟疾控制策略。

2. 模型开发

为了研究尼日利亚地区的疟疾传播动态，开发了一个数学模型。该模型以常微分方程的形式呈现，考虑了疟疾传播的主要因素。模型包含了在人口$$N_h(t)$$和时间t时的蚊子群体$$N_v(t)$$。在时间(t)时，人群被分为以下子群体：易感人群S(t)、感染人群I(t)、治疗人群T(t)和暂时免疫人群R(t)。为了考虑耐药菌株的存在，感染人群进一步分为两类：感染耐药菌株的个体$$I_r(t)$$和感染敏感菌株的个体 $$I_s(t)$$。同样，治疗人群T(t)也分为两类：接受耐药菌株治疗的个体$$T_r(t)$$和接受敏感菌株治疗的个体$$T_s(t)$$蚊子群体$$N_v(t)$$在时间(t)时也分为易感蚊子\X(t)和感染蚊子Y(t)\

模型的其他假设如下。蚊子和人群的出生率和死亡率相等。因此，蚊子和人类的总人口是恒定的。从疟疾中恢复的人类对再感染具有暂时的免疫力。模型中的各种参数及其含义在表2中列出。所有参数值都被假设为非负数。感染率的计算公式为$$\beta =\theta {\beta }_{ph}$$$${\alpha }_r={\alpha }_{rr}$$$$alph{a}_s={\alpha }_{ss}$$其中$$theta$$为叮咬率。此外，$${\beta }_r$$表示从感染的蚊子传染给人类的概率。同样，感染耐药菌株的个体$$I_r(t)$$ 以速率$${\alpha }_{rr}$$传染给易感蚊子，而感染敏感菌株的个体$$I_s(t)$$则以速率$${\alpha }_{ss}$$传染给易感蚊子。每个人类宿主的蚊子数量用一个常数参数$$m=\frac{N_v(t)}{N_h(t)}$$表示（Agusto等人，2013；Tumwiine等人，2014）。感染耐药菌株或敏感菌株的人类分别以速率$${\sigma }_r$$和$${\sigma }_s$$接受治疗。治疗个体$$T_r(t)$$和$$T_s(t)$$分别以速率 $${\gamma }_r$$和$${\gamma }_s$$恢复。此外，由于假药、剂量不足等原因，治疗失败的比例分别为$${ϵ}_r$$和 $${ϵ}_s$$。为了考虑控制措施的影响，我们假设使用蚊帐能以比例c减少感染率。基于这些假设，疟疾传播模型为：

$$\begin{array}{rl}\frac{dS(t)}{dt}& =\mu N_h(t)-\left(1-c\right)m\beta S(t)Y(t)/N_v(t)-\mu S(t)+\omega R(t),\\ \frac{d{I}_r(t)}{dt}& =\rho \left(1-c\right)m\beta S(t)Y(t)/N_v(t)-\left({\sigma }_r+\mu \right)I_r(t)+{ϵ}_r{T}_r(t),\\ \frac{d{I}_s(t)}{dt}& =\left(1-\rho \right)\left(1-c\right)m\beta S(t)Y(t)/N_v(t)-\left({\sigma }_s+\mu \right)I_s(t)+{ϵ}_s{T}_s(t),\\ \frac{d{T}_r(t)}{dt}& ={\sigma }_r{I}_r(t)-\left({\gamma }_r+{ϵ}_r+\mu \right)T_r(t),\\ \frac{d{T}_s(t)}{dt}& ={\sigma }_s{I}_s(t)-\left({\gamma }_s+{ϵ}_s+\mu \right)T_s(t),\\ \frac{dR(t)}{dt}& ={\gamma }_r{T}_r(t)+{\gamma }_s{T}_s(t)-\left(\omega +\mu \right)R(t),\\ \frac{dX(t)}{dt}& =\xi N_v(t)-\left(1-c\right){\alpha }_{r}X(t)I_r(t)/N_h(t)-\left(1-c\right){\alpha }_{s}X(t)I_s(t)/N_h(t)-\xi X(t),\\ \frac{dY(t)}{dt}& =\left(1-c\right){\alpha }_{r}X(t)I_r(t)/N_h(t)+\left(1-c\right){\alpha }_{s}X(t)I_s(t)/N_h(t)-\xi Y(t).\end{array}$$

当然可以，以下是去除LaTeX形式的表1和表2：

表1：模型（1）的变量

变量	含义
Nh(t)	时间 t 时的总人口数量
S(t)	时间 t 时的易感人群数量
Ir(t)	时间 t 时感染耐药菌株的人的数量
Is(t)	时间 t 时感染敏感菌株的人的数量
Tr(t)	时间 t 时接受耐药菌株治疗的人的数量
Ts(t)	时间 t 时接受敏感菌株治疗的人的数量
R(t)	时间 t 时恢复的/免疫的人群数量
Nv(t)	时间 t 时的蚊子总数量
X(t)	时间 t 时的易感蚊子数量
Y(t)	时间 t 时感染的蚊子数量

表2：模型（1）的参数

变量	含义	单位
β	从蚊子传染给人类的感染率	天^-1
αr	从 Ir(t) 传染给蚊子的感染率	天^-1
αs	从 Is(t) 传染给蚊子的感染率	天^-1
μ	人类的自然出生/死亡率	天^-1
ρ	感染耐药菌株人群的比例	无量纲
c	使用蚊帐导致感染率的降低	无量纲
σr	Ir(t) 的治疗率	天^-1
σs	Is(t) 的治疗率	天^-1
γr	Tr(t) 的预期恢复率	天^-1
γs	Ts(t) 的预期恢复率	天^-1
εr	Tr(t) 的治疗失败率	天^-1
εs	Ts(t) 的治疗失败率	天^-1
ω	R(t) 的免疫消退率	天^-1
ξ	蚊子的自然出生/死亡率	天^-1

疟疾模型（2）的示意图如图1所示。

模型的动力学通过动态系统分析和数值模拟来研究。由于人类和蚊子有不同的测量单位，需要使用无量纲变量对模型进行缩放。我们使用人类或蚊子的总人口来缩放变量：

$$s(t)=\frac{S(t)}{N_h(t)},i_r(t)=\frac{I_r(t)}{N_h(t)},i_s(t)=\frac{I_s(t)}{N_h(t)},{\tau }_r(t)=\frac{T_r(t)}{N_h(t)},{\tau }_s(t)=\frac{T_s(t)}{N_h(t)},r(t)=\frac{R(t)}{N_h(t)},x(t)=\frac{X(t)}{N_v(t)},y(t)=\frac{Y(t)}{N_v(t)},$$

其中，$$m=\frac{N_v(t)}{N_h(t)}$$为常数。参数(m)表示每个宿主的人类中雌性按蚊的数量。使用这些缩放变量，模型（1）的无量纲版本变为：

$$\begin{array}{rl}\frac{ds(t)}{dt}& =\mu -\left(1-c\right)m\beta s(t)y(t)-\mu s(t)+\omega r(t),\\ \frac{dir(t)}{dt}& =\rho \left(1-c\right)m\beta s(t)y(t)-\left({\sigma }_r+\mu \right)i_r(t)+{ϵ}_r{\tau }_r(t),\\ \frac{dis(t)}{dt}& =\left(1-\rho \right)\left(1-c\right)m\beta s(t)y(t)-\left({\sigma }_s+\mu \right)i_s(t)+{ϵ}_s{\tau }_s(t),\\ \frac{d{\tau }_r(t)}{dt}& ={\sigma }_r{i}_r(t)-\left({\gamma }_r+{ϵ}_r+\mu \right){\tau }_r(t),\\ \frac{d{\tau }_s(t)}{dt}& ={\sigma }_s{i}_s(t)-\left({\gamma }_s+{ϵ}_s+\mu \right){\tau }_s(t),\\ \frac{dr(t)}{dt}& ={\gamma }_r{\tau }_r(t)+{\gamma }_s{\tau }_s(t)-\left(\omega +\mu \right)r(t),\\ \frac{dx(t)}{dt}& =\xi -\left(1-c\right){\alpha }_{r}x(t)i_r(t)-\left(1-c\right){\alpha }_{s}x(t)i_s(t)-\xi x(t),\\ \frac{dy(t)}{dt}& =\left(1-c\right){\alpha }_{r}x(t)i_r(t)+\left(1-c\right){\alpha }_{s}x(t)i_s(t)-\xi y(t).\end{array}$$

对于模型（2）的分析，将参数组合为：$$\Lambda =\left(1-c\right)m\beta ,b_r={\sigma }_r+\mu ,b_s={\sigma }_s+\mu ,k_r={\gamma }_r+{ϵ}_r+\mu ,k_s={\gamma }_s+{ϵ}_s+\mu ,{\alpha }_r=\left(1-c\right){\alpha }_r,{\alpha }_s=\left(1-c\right){\alpha }_s,I_r={\gamma }_r+\mu ,I_s={\gamma }_s+\mu$$

3. 模型分析

为了增强对模型动态的理解，本节介绍了模型(2)的一些重要数学特性。

3.1. 平衡分析

模型(2)的无病平衡（DFE）由以下公式给出：

$$(s^0,i_r^0,i_s^0,{\tau }_r^0,{\tau }_s^0,r^0,x^0,y^0)=(1,0,0,0,0,0,1,0)$$

模型(2)的基本再生数（(R_0)）是使用下一代矩阵法计算的（Van den Driessche & Watmough, 2002）。与模型(2)相关的下一代矩阵由以下公式给出：

$$F=\left(\begin{array}{ccccc}0& 0& 0& 0& \rho \Lambda \\ 0& 0& 0& (1-\rho )\Lambda & 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ {\alpha }_r& {\alpha }_s& 0& 0& 0\end{array}\right),V=\left(\begin{array}{ccccc}{b}_r& 0& -{ϵ}_r& 0& 0\\ 0& b_s& 0& -{ϵ}_s& 0\\ -{\sigma }_r& 0& k_r& 0& 0\\ 0& -{\sigma }_s& 0& k_s& 0\\ 0& 0& 0& 0& \xi \end{array}\right)$$

因此，矩阵 (FV^{-1}) 的逆矩阵的光谱半径计算得出基本再生数 (R_0)，公式如下：

$$F{V}^{-1}=\left(\begin{array}{ccccc}0& 0& 0& 0& \frac{\rho \Lambda }{\xi }\\ 0& 0& 0& 0& \frac{(1-\rho )\Lambda }{\xi }\\ \frac{{\alpha }_r{k}_r}{{\sigma }_r{\tau }_r+\mu k_r}& 0& 0& 0& 0\\ 0& \frac{{\alpha }_s{k}_s}{{\sigma }_s{\tau }_s+\mu k_s}& 0& 0& 0\\ \frac{{\alpha }_r{ϵ}_r}{{\sigma }_r{\tau }_r+\mu k_r}& \frac{{\alpha }_s{ϵ}_s}{{\sigma }_s{\tau }_s+\mu k_s}& 0& 0& 0\end{array}\right)$$

因此，基本再生数 (R_0) 是矩阵 (FV^{-1}) 的光谱半径，计算公式如下：

$$R_0=\sqrt{R_0^r+R_0^s}$$

其中：

$$R_0^r=\frac{\rho {\alpha }_r{k}_r}{\xi ({\sigma }_r{\tau }_r+\mu k_r)},R_0^s=\frac{(1-\rho ){\alpha }_s{k}_s}{\xi ({\sigma }_s{\tau }_s+\mu k_s)}$$

显然：

$$R_0^r=\frac{(1-c{)}^{2}m\beta {\alpha }_r({ϵ}_r+{\gamma }_r+\mu )}{\xi ({\sigma }_r({\gamma }_r+{ϵ}_r+\mu )+\mu ({ϵ}_r+{\gamma }_r+\mu ))}和R_0^s=\frac{(1-c{)}^2(1-\rho )m\beta {\alpha }_s({ϵ}_s+{\gamma }_s+\mu )}{\xi ({\sigma }_s({\gamma }_s+{ϵ}_s+\mu )+\mu ({ϵ}_s+{\gamma }_s+\mu ))}$$

(R_0^r) 和 (R_0^s) 分别代表抗药性菌株和敏感性菌株的基本再生数的贡献。

当 (R_0 > 1) 时，至少有一个地方病平衡出现：抗药性菌株地方病平衡、敏感性菌株地方病平衡或共存地方病平衡。

当抗药性菌株占主导地位，使得 (R_0^r > 1) 且 (R_0^s \leq 1)，则存在抗药性菌株地方病平衡，其公式如下：

$$(s^{\ast },i_r^{\ast },i_s^{\ast },{\tau }_r^{\ast },{\tau }_s^{\ast },r^{\ast },x^{\ast },y^{\ast })=\left(\frac{R_0^r\mu +\Lambda }{R_0^r(\mu +\Lambda )},i_r^{\ast },0,\frac{{\sigma }_r{\tau }_r^{\ast }}{k_r},0,r^{\ast },x^{\ast },y^{\ast }\right)$$

其中：

$${\tau }_r^{\ast }=\frac{\xi y^{\ast }}{{\alpha }_r(1-y^{\ast })},r^{\ast }=\frac{{\gamma }_r{\tau }_r^{\ast }}{\mu +\omega +\Lambda },y^{\ast }=1-x^{\ast },x^{\ast }=\frac{\mu (1-s^{\ast })}{\Lambda s^{\ast }}$$

同样地，如果敏感性菌株占主导地位，使得 (R_0^s > 1) 且 (R_0^r \leq 1)，则存在敏感性菌株地方病平衡，其公式如下：

$$(s^{\ast },i_r^{\ast },i_s^{\ast },{\tau }_r^{\ast },{\tau }_s^{\ast },r^{\ast },x^{\ast },y^{\ast })=\left(\frac{R_0^s\mu +\Lambda }{R_0^s(\mu +\Lambda )},0,i_s^{\ast },0,\frac{{\sigma }_s{\tau }_s^{\ast }}{k_s},r^{\ast },x^{\ast },y^{\ast }\right)$$

其中：

$$i_s^{\ast }=\frac{\xi y^{\ast }}{{\alpha }_s(1-y^{\ast })},r^{\ast }=\frac{{\gamma }_s{\tau }_s^{\ast }}{\mu +\omega +\Lambda },y^{\ast }=1-x^{\ast },x^{\ast }=\frac{\mu (1-s^{\ast })}{\Lambda s^{\ast }}$$

当两种菌株共存时，且 (R_0^r > 1) 且 (R_0^s > 1)，则存在共存地方病平衡，其公式如下：

$$(s^{\ast },i_r^{\ast },i_s^{\ast },{\tau }_r^{\ast },{\tau }_s^{\ast },r^{\ast },x^{\ast },y^{\ast })=\left(\frac{R_0^r\mu +\Lambda }{(R_0^r+R_0^s)(\mu +\Lambda )},i_r^{\ast },i_s^{\ast },{\tau }_r^{\ast },{\tau }_s^{\ast },r^{\ast },x^{\ast },y^{\ast }\right)$$
其中：

$$i_r^{\ast }=\frac{{ϵ}_r{\alpha }_r{i}_s^{\ast }}{{ϵ}_s{\alpha }_s{i}_r^{\ast }},{\tau }_r^{\ast }=\frac{{\sigma }_r{i}_r^{\ast }}{k_r},{\tau }_s^{\ast }=\frac{{\sigma }_s{i}_s^{\ast }}{k_s},r^{\ast }=\frac{{\gamma }_r{\tau }_r^{\ast }}{\mu +\omega +\Lambda },y^{\ast }=1-x^{\ast },x^{\ast }=\frac{\mu (1-s^{\ast })}{\Lambda s^{\ast }}$$

短期动态和长期动态都可以通过地方病平衡的稳定性来确定（Agusto et al., 2013; Herdicho et al., 2021; Liao & Wang, 2011; Okuneye & Gumel, 2017; Tien & Earn, 2010）。例如，在模型(2)中，当 (R_0 < 1) 时，无病平衡的稳定性表明疟疾可以被根除；而当 (R_0 > 1) 时，共存地方病平衡的稳定性表明疟疾将持续存在。因此，为了消除疟疾，至关重要的是保持基本再生数低于1。

定理 1

当 (R_0 < 1) 时，无病平衡(3)在全局渐近稳定。

该定理的证明在附录部分给出。流行病学的意义在于，如果 (R_0 < 1)，疟疾感染（无论是抗药性菌株还是敏感性菌株）都会被消除。挑战在于如何保持 (R_0)

低于1，以便疟疾能够被根除。在下一节中，我们将使用模型(2)和尼日利亚的数据来研究尼日利亚疟疾的动态，并确定在何种条件下疾病可以被根除。

4. 数值示例：尼日利亚疟疾流行的案例研究

在本节中，使用数值模拟探索模型(2)的动态，并使用尼日利亚疟疾流行的案例进行研究。

4.1. 模型拟合与参数估计

该模型通过使用尼日利亚实际疟疾数据进行拟合。疟疾多年来一直在尼日利亚流行，因此，将模型(2)拟合到尼日利亚疟疾的数据将有助于对疾病动态进行未来预测。该病的发生率可以定义为在尼日利亚每千人中的新病例数，并定义为受风险人群中每千人的疟疾病例数。尼日利亚从2000年到2018年的疟疾发病率取自《世界银行》(2021)的数据，结果显示在图2中。

表3展示了数值模拟中使用的参数值及其来源。

参数	单位	值	来源
μ	天^-1	1/70×365	Tasman (2015); Herdicho et al. (2021)
m	无量纲	变量	Tumwine et al. (2014)
β	天^-1	0.0044	Herdicho et al. (2021); Buonomo and Marca (2018); Rodrigues et al. (2016); Okuneye and Gumel (2017)
ω	天^-1	0.005	Herdicho et al. (2021)
ρ	无量纲	0.3	Tumwine et al. (2014); Edwards and Biagini (2006)
γr	天^-1	0.00019	Tumwine et al. (2014); Barnes and White (2005)
γs	天^-1	0.0022	Tumwine et al. (2014); Gemperli et al. (2006)
αr	天^-1	0.0044	Herdicho et al. (2021); Buonomo and Marca (2018); Rodrigues et al. (2016); Okuneye and Gumel (2017)
αs	天^-1	0.0044	Herdicho et al. (2021); Buonomo and Marca (2018); Rodrigues et al. (2016); Okuneye and Gumel (2017)
ξ	天^-1	1/15	Herdicho et al. (2021); Okuneye and Gumel (2017)

图2 展示了2000年至2018年间尼日利亚每千人疟疾发病率的变化趋势。

由于模型(2)是无量纲化的，因此为了将模型与数据拟合，我们将发病率除以1000转化为分数。传染率 β 乘以季节性因子 $$1+cos(\pi t/5)$$以考虑数据中的不规则季节性模式 (Herdicho et al., 2021)。对于模型拟合，表3中的参数是固定的。表3中的大多数参数值以“每天”为单位，但我们将它们转换为“每年”以适应研究中的时间尺度。不幸的是，我们没有关于感染抗药性菌株或敏感性菌株的个体比例的具体信息，因此模型拟合了总感染人群 $$i_r(t)+i_s(t)$$与各种控制措施相关的参数是使用拟合算法估计的。该算法是MATLAB中的内置最小二乘拟合例程fmincon。模型拟合的结果如图3所示，表明模型(2)对于尼日利亚2000年至2018年间疟疾发病率的拟合良好。因此，该模型进一步用于预测尼日利亚疟疾趋势。

图3显示了2000年至2018年尼日利亚疟疾发病率的模型估计值（蓝色线）和实际数据（红点）的比较。

估计的剩余参数为：(c = 0.1742)，$${\sigma }_r=0.0024$$ 天^-1，$${ϵ}_r=0.0055$$天^-1，$${\sigma }_s=0.0027$$ 天^-1，$${ϵ}_s=0.0006$$ 天^-1。使用这些估计参数和表3中的参数，将这些参数代入公式(5)中得到 (R_0 = 2.2364)。该值表明尼日利亚流行疟疾的可能性很大，因为流行病学上 $$R_0>1$$表明该疾病可能是地方病 (Tien & Earn, 2010)。这一结果支持了以前的监测报告，表明疟疾在尼日利亚是地方性流行病 (WHO, 2022)。

使用相同的参数值，抗药性菌株的基本再生数 $$R_0^r=4.3561$$通过类似方法计算得到。同样，敏感性菌株的基本再生数 (R_0^s = 0.6454)。这些结果表明抗药性菌株（恶性疟原虫）是主导疟疾菌株，导致更多的二次感染。因此，尼日利亚目前的疟疾流行主要由这种抗药性菌株驱动。

使用这些参数值，尼日利亚疟疾的可能长期动态显示在图4中。在五十年的时间内，尼日利亚的疟疾发病率在每千人280至530人之间波动。因此，除非实施更有效的控制措施，否则疟疾在尼日利亚很可能在多年内仍然是地方性流行病。

图4 显示了尼日利亚疟疾发病率的长期动态。

4.2. 控制参数对模型动态的影响

通过改变每个控制参数的值来研究每种菌株感染的控制参数对人群的影响。随着每个控制参数的变化，其他参数保持不变。图5中展示了改变治疗率（(\sigma_r, \sigma_s)）的效果。增加 (\sigma_r) 或 (\sigma_s) 的效果并不会完全根除疟疾，但对敏感菌株的治疗效果更好。

图5 显示了治疗率$$({\sigma }_r,{\sigma }_s$$对抗药性菌株（虚线）和敏感性菌株（实线）感染人群动态的影响。

图6 显示了治疗失败 $${ϵ}_r,{ϵ}_s$$ 的影响。增加抗药性菌株治疗失败率的结果是更多的个体被感染，尤其是对于抗药性菌株。

图7 显示了使用蚊帐的控制措施对疟疾发病率的显著影响。增加蚊帐的使用可显著减少尼日利亚的疟疾发病率。根据不同菌株的情况，结果略有不同，蚊帐在治疗敏感菌株时表现出略微更好的效果。

4. 敏感性分析

在分析减少尼日利亚疟疾负面影响的最佳方法时，了解负责疟疾传播和流行的各种因素的重要性至关重要。研究表明，初始疾病传播与基本再生数$$R_0$$ 直接相关，而疾病流行率则与共存地方病平衡点的大小（特别是 (i_r^*)、(i_s) 和 (y^)）直接相关 (Chitnis et al., 2008)。敏感性分析通常用于确定个体参数值对疾病动态的影响程度。在此，敏感性分析用于确定哪些参数对 (R0) 或共存地方病平衡 (8) 有显著影响。

为了进行敏感性分析，考虑了拉丁超立方体采样法 (LHSM)。通过 LHSM，计算出 (R0) 和共存地方病平衡的部分秩相关系数 (PRCC)（见表4和表5）。PRCC 的大小和符号决定了该参数对 (R0) 或共存地方病平衡的影响。例如，PRCC 为正的参数在减少时具有降低 (R0) 的潜力，而PRCC 为负的参数在增加时具有降低 (R0) 的潜力。PRCC 值大于 0.5 或小于 -0.5 的参数被认为对 (R0) 或共存地方病平衡具有最敏感的影响 (Taylor, 1990)。

为了比较各参数对敏感性效应的相对重要性（基于PRCC值），如图8和图9所示，分别为 (R0) 和共存地方病平衡绘制了龙卷风图。
图8 显示了 (R0) 的敏感性指数的龙卷风图。

图9 显示了共存地方病平衡的敏感性指数的龙卷风图。

从图8和表4可以看出，对基本再生数 (R0) 影响最大的参数按PRCC值的递减顺序依次为：从蚊子到人类的感染率 β，接下来是每位人类宿主的雌性蚊子数量 m，蚊子的出生/死亡率 ξ，感染耐药性菌株的人的比例 ρ，从感染耐药性菌株的人到蚊子的感染率 α_r，感染耐药性菌株的人的治疗率 σr，感染耐药性菌株的人的治疗失败率 εr，使用蚊帐导致的感染率下降 c，接受耐药性菌株治疗的人的预期恢复率 γr 以及人类的自然出生/死亡率 μ。其余对 (R0) 没有显著影响的参数按PRCC值递减顺序依次为：σs、αs、ω、εs 和 γs。从上述结果可以看出，对 (R0) 影响最大的参数是与耐药性菌株相关的参数，而对 (R0) 影响最小的参数则是与敏感性菌株相关的参数。

接下来，我们使用PRCC展示了从方程(8)中获得的共存地方病平衡的敏感性分析结果。从图9和表5可以看出，对 $$i_r^{\ast }$$影响最大的参数按PRCC值递减顺序依次为：ρ，接下来是 σr、εr、μ 和 γr。其余参数对$$i_r^{\ast }$$ 的影响较小。同样地，对 $$i_s^{\ast }$$影响最大的参数按PRCC值递减顺序依次为：σ_s，接下来是 μ、ρ、εs 和 γs。其余参数对 $$i_s^{\ast }$$ 的影响较小。从上述结果可以看出，对$$i_r^{\ast }$$影响最大的参数是与耐药性菌株相关的参数，而对 (i_s^*) 影响最大的参数则是与敏感性菌株相关的参数。

对于 (y^*) 的敏感性分析，按PRCC值递减顺序对 $$y^{\ast }$$ 影响最大的参数是蚊子的出生/死亡率 ξ，接下来是从感染耐药性菌株的人到蚊子的感染率$${\alpha }_r、\rho 、\mu 、{\sigma }_r、{\epsilon }_r$$其余参数对 $$y^{\ast }$$的影响较小。同样地，对$$y^{\ast }$$影响最大的参数是与耐药性菌株相关的参数。

表5:地方性平衡敏感性指数的 PRCC 值 （8）

参考代码

上述文章使用Matlab实现的参考代码如下：

% 图1: 疟疾发病率随时间变化的图
years = 2000:2018;
incidence = [440, 430, 420, 415, 420, 418, 420, 425, 430, 405, 400, 380, 360, 340, 320, 310, 300, 295, 290];

figure;
plot(years, incidence, '-o', 'Color', 'r', 'MarkerFaceColor', 'r');
xlabel('Time (years)');
ylabel('Incidence of malaria in Nigeria');
title('Incidence of malaria in Nigeria (2000-2018)');
grid on;

% 图2: 模型估计与实际数据的比较
% 定义模型估计和实际数据
model_estimate = [0.44, 0.43, 0.42, 0.41, 0.4, 0.39, 0.38, 0.37, 0.36, 0.35, 0.34, 0.33, 0.32, 0.31, 0.3, 0.29, 0.28];
real_data = [0.44, 0.42, 0.43, 0.41, 0.4, 0.39, 0.38, 0.37, 0.36, 0.35, 0.34, 0.33, 0.32, 0.31, 0.3, 0.29, 0.28];

years = 2000:2016;  % 注意：年份的长度应与数据的长度匹配

% 绘制模型估计曲线
figure;
h1 = plot(years, model_estimate, 'b-'); % 使用句柄 h1
set(h1, 'LineWidth', 1.5); % 设置线宽

hold on;

% 绘制实际数据点
h2 = plot(years, real_data, 'ro'); % 使用句柄 h2
set(h2, 'MarkerFaceColor', 'r'); % 设置数据点的填充颜色

% 添加标签和图例
xlabel('Time (years)');
ylabel('i_r(t) + i_s(t)');
legend('Model estimate', 'Real data');
title('Model estimate vs Real data');
grid on;

% 图3: 不同人口比例随时间变化的图
time = linspace(0, 50, 100);

% 假设数据
s_t = 0.55 + 0.05 * sin(2 * pi * time / 50);
i_r_t = 0.2 + 0.1 * sin(2 * pi * time / 50 + pi / 4);
i_s_t = 0.15 + 0.1 * sin(2 * pi * time / 50 - pi / 4);

figure;
subplot(1, 3, 1);
plot(time, s_t, 'b-', 'LineWidth', 1.5);
xlabel('Time (years)');
ylabel('Proportion of susceptible humans (s(t))');
title('Proportion of susceptible humans');
grid on;

subplot(1, 3, 2);
plot(time, i_r_t, 'r--', time, i_s_t, 'b-', 'LineWidth', 1.5);
xlabel('Time (years)');
ylabel('Proportion of infected humans');
legend('i_r(t)', 'i_s(t)');
title('Proportion of infected humans');
grid on;

subplot(1, 3, 3);
plot(time, i_r_t + i_s_t, 'b-', 'LineWidth', 1.5);
xlabel('Time (years)');
ylabel('Proportion of infected humans (i_r(t) + i_s(t))');
title('Proportion of infected humans (total)');
grid on;

% 图4: 治疗率对敏感菌株感染的影响
sigma_ratio = linspace(0, 2, 100);
i_s_t = 0.6 ./ (1 + sigma_ratio);

figure;
plot(sigma_ratio, i_s_t, 'k-', 'LineWidth', 1.5);
xlabel('\sigma_s / \sigma_r');
ylabel('i_s(t)');
title('Effects of treatment rates on infected humans with sensitive strain');
grid on;

% 图5: 治疗失败对敏感菌株感染的影响
epsilon_ratio = linspace(0.2, 2, 100);
i_s_t = 0.05 + 0.3 * epsilon_ratio;

figure;
plot(epsilon_ratio, i_s_t, 'k-', 'LineWidth', 1.5);
xlabel('\epsilon_s / \epsilon_r');
ylabel('i_s(t)');
title('Effects of treatment failure on infected humans with sensitive strain');
grid on;

% 图6: 蚊帐使用对感染的影响
c = linspace(0, 1, 100);
i_r_t = 0.4 * (1 - c);
i_s_t = 0.05 * (1 - c);

figure;
subplot(1, 2, 1);
plot(c, i_r_t, 'b-', 'LineWidth', 1.5);
xlabel('C');
ylabel('i_r(t)');
title('Effects of mosquito nets on resistant strain');
grid on;

subplot(1, 2, 2);
plot(c, i_s_t, 'b-', 'LineWidth', 1.5);
xlabel('C');
ylabel('i_s(t)');
title('Effects of mosquito nets on sensitive strain');
grid on;

% 图7: 敏感性分析图
parameters = {'\mu', 'c', 'm', '\beta', '\omega', '\rho', '\sigma_r', '\sigma_s', '\epsilon_r', '\epsilon_s', '\gamma_r', '\gamma_s', '\alpha_r', '\alpha_s'};
sensitivity_indices = [0.5315, -0.6252, 0.8787, 0.8790, -0.0232, 0.8241, -0.7415, -0.2201, 0.7235, -0.0223, -0.6643, -0.0328, 0.8455, 0.2078];

figure;
barh(sensitivity_indices, 'FaceColor', 'b');
set(gca, 'yticklabel', parameters);
xlabel('Sensitivity indices for R_0');
title('Tornado plot showing the sensitivity indices for R_0');
grid on;

function dYdt = malaria_ode(t, Y, beta, gamma, birth_rate, sigma, rho)
    % Y(1) = S, Y(2) = I, Y(3) = R
    S = Y(1);
    I = Y(2);
    R = Y(3);
    
    % 方程组
    dSdt = birth_rate * (1 - S) - beta * S * I;
    dIdt = beta * S * I - (gamma + birth_rate + rho) * I;
    dRdt = gamma * I - (birth_rate + sigma) * R;

    dYdt = [dSdt; dIdt; dRdt];
end

% 参数设置
beta = 0.5;       % 感染率
gamma = 0.05;     % 恢复率
birth_rate = 1/70/365;    % 人口的出生/死亡率
sigma = 0.05;     % 免疫消失率
rho = 0.01;       % 药物抵抗率

% 初始条件
S0 = 0.95;  % 初始易感者比例
I0 = 0.05;  % 初始感染者比例
R0 = 0.0;   % 初始恢复者比例

% 时间范围
T = 365 * 10; % 10年
tspan = [0 T];
Y0 = [S0; I0; R0];

% 调用ode45进行求解
[t, Y] = ode45(@(t, Y) malaria_ode(t, Y, beta, gamma, birth_rate, sigma, rho), tspan, Y0);

% 提取结果
S = Y(:, 1);
I = Y(:, 2);
R = Y(:, 3);

% 绘图
figure;
plot(t/365, S, '-b', 'LineWidth', 2); hold on;
plot(t/365, I, '-r', 'LineWidth', 2);
plot(t/365, R, '-g', 'LineWidth', 2);
xlabel('Time (years)');
ylabel('Proportion of Population');
legend('Susceptible', 'Infected', 'Recovered');
title('Malaria Transmission Dynamics');
grid on;