基于数值与分数阶方法的巴基斯坦结核病战略最优控制 SVITR 模型

《Computational Biology and Chemistry》：A strategic optimal control SVITR Model for tuberculosis (TB) disease in Pakistan by using numerical and fractional schemes

　　

• 针对巴基斯坦开发了确定性和分数阶SVITR-TB模型。• 稳定性和敏感性分析验证了结核病的传播动力学。• NSFD和分数阶RK4方案提高了数值模型精度。• 最优控制识别出提高认识和疫苗接种是关键干预措施。• 策略1是最具成本效益的结核病控制方法。引言分枝杆菌结核复合群（MTB）是导致结核病（TB）的病原体，结核病是一种广泛存在且可能致命的传染病。1882年，罗伯特·科赫发现分枝杆菌结核复合群是结核病的主要致病原（Cambau和Drancourt, 2014）。结核分枝杆菌主要侵犯肺部及其他器官。空气是结核病的传播媒介。当活动性结核病个体通过咳嗽、打喷嚏、说话或吐痰将含有结核分枝杆菌的飞沫释放到周围环境中时，传播便发生。吸入此类空气可能使易感者受到感染。与结核病患者经常接触的人风险更高。肺部是结核分枝杆菌的寄生场所，该菌可扩散至身体其他部位。患者的免疫系统不足以将其清除。活动性结核病患者的传染概率为10%。在数月乃至数年时间跨度内，大多数感染者保持潜伏感染状态，即不具有传染性。然而，这一时间框架可能因个体所合并感染的其他疾病类型而发生变化。合并感染会因削弱宿主免疫系统而增加感染风险并降低免疫力，其中HIV感染显著加速结核病向活动期的进展。据世界卫生组织2024年10月报告，2023年全球结核病病例超过1080万例，其中男性600万例，女性360万例，儿童130万例。同年，125万人死于该病。据世卫组织全球结核病规划部报告，2023年全球结核病估计死亡人数为：15岁及以上成年男性568,000人，15岁及以上成年女性352,000人，儿童及青少年166,000人。近期研究显示，结核病对健康负担显著，尤其是合并HIV感染的孕妇和产后妇女。据WHO 2024年报告，全球将有1070万例结核感染和123万例死亡。巴基斯坦在2022年位列全球30个高负担国家第5位，与印度（280万例）、印度尼西亚（110万例）、中国（748,000例）、菲律宾（737,000例）和巴基斯坦（608,000例）共同承担全球超过一半的结核病病例，2022年全球结核病病例总计1060万例。此外，据世卫组织2026年3月24日最新报告，巴基斯坦有669,000例结核病病例，年死亡51,000人。该国在全球受结核病影响最严重的国家中排名第5，占东地中海地区结核病病例的73%。每天巴基斯坦有超过1800例新发结核病病例和140例结核病死亡。世界银行集团报告指出，2021年至2023年间，每10万人结核病发病率为273、278和277。数学模型采用定量和分析方法来理解传染病传播、评估关键流行病学决定因素以及预测未来趋势。通过数值模拟和基于微分方程的分析，它们还可评估免疫接种和治疗等控制手段的有效性。Waaler等人（1962）基于结核病传播特征建立了一个将人群分为三个类别的模型系统，对结核病流行病学做出了重要贡献。Kalizhanova等人（2024）利用2014至2019年的数据，采用SARIMA和SIR模型研究了哈萨克斯坦结核病传播模式。SIR模型显示流行病处于下降趋势，基本再生数低于1，而SARIMA模型准确刻画了季节性变化。两种模型均表现出较强的预测性能，有助于未来的结核病控制工作。Side等人（2018）采用四阶Runge-Kutta方法对结核病传播的SIR模型进行数值求解，使用了2015年南苏拉威西省马卡萨尔的结核病病例数据。研究结果表明该方法可支持结核病预防工作，并准确模拟发病率的上升和下降。Ullah等人（2019）利用2002至2017年数据，阐明了巴基斯坦开伯尔-普什图省结核病传播动力学，并建立了数学框架。研究表明当基本再生数小于1时，无病平衡状态得以维持。模拟和敏感性分析揭示了影响疾病消除的关键因素。Castillo-Chavez和Feng（1997a, 1997b）研究了疾病动力学，提供了多个模型用于深入观察并生成了相关结果。Castillo-Chavez和Song（2004）对结核病动力学的文献进行了全面综述，收集了多个时间序列系统并提供了严谨的框架。Tan和Isa（2025）考察了含有人口学因素和不考虑人口学因素两种情形的SIR模型以预测马来西亚结核病传播。他们发现人口学SIR模型产生了更准确的预测，且传播率和恢复率对结核病动力学具有重大影响。Adi（2024）利用印度尼西亚数据建立并检验了一个包含心理影响和假性恢复的SEIR结核病模型，分析了稳定性、敏感性、分岔以及影响R0的关键因素。1994年，Mickens提出使用数值技术是确定微分方程精确解的唯一途径。Algehyne和ud Din（2021）研究了基本再生数，即一次主感染在感染期内引发的预期继发感染数。Gurski（2013）开发了一个SQIR模型以考察COVID-19传播的平衡点、基本再生数和稳定性特征，并采用非标准有限差分（NSFD）方案进行数值模拟。Memarbashi等人（2020）提出采用Mickens离散化方案引入NSFD技术，证明了所提出的NSFD方案能够保持相关连续流行病学模型的平衡点。作者还考察了所提出模型的分类特征，发现其保持了平衡点稳定性、正性和Neimark-Sacker分岔。Ahmed等人（2018）为饱和传播率SEIQV疫情模型开发了一种NSFD方案，证明其具有优于RK4方法的稳定性和正性保持能力。Puntani Pongsumpun（2021）通过求解数值解和平衡点来研究一个数学模型，计算了基本再生数以检验各因素的影响，并计算了单个感染者在感染期内可能引发的继发感染数。Fatima等人（2018）首先将结核病的SEIQR模型转化为微分方程，随后确定了传播潜力、无病平衡状态、 endemic 平衡点及局部稳定性评估，最后对阈值数的数值模拟进行了快速考察。Ahmed等人（2019c）为反应-扩散SEIR疫情模型的数值解提出了两个算子的NSFD方案。Ahmed等人（2019a）于2019年提出了Brusselator反应-扩散模型并使用Neumann稳定性准则评估其稳定性，其中包括实施通用有限差分方案、前向欧拉显式有限差分法和半隐式Crank-Nicolson有限差分法。Ahmed等人（2019b）还研究了自催化糖酵解模型及其数值解，该模型得出了阳性结果，因为未知条件揭示了化学组分的浓度。Kanwal等人（2022）创建并对比了前向欧拉、RK4和NSFD方案用于结核病传播模型，证明NSFD在保持正性和稳定性方面优于传统数值方法。Jeelani等人（2024）针对独特疾病开发了SARS-CoV-2双变异株传播的计算模型，针对变异株一进行疫苗接种，并采用随机方法确定平衡点、再生数和随机系统的稳定性，最后使用NSFD方法获得数值模拟结果。潜伏期长、进展延迟和治疗相关副作用是结核病的特征，表明疾病当前状态可能既取决于当前条件也取决于既往流行病学历史。这种记忆和遗传效应可以通过非局部算子自然地融入分数阶模型中。因此，与传统整数阶模型相比，分数阶模型通过记忆效应为表征疾病动力学提供了更大的灵活性。此外，最优控制理论能够确定既有效降低疾病流行率又最小化实施成本的干预策略。因此，最优控制为评估有效控制策略和理解结核病传播提供了有用的框架。Ullah等人（2020）使用Caputo-Fabrizio导数证明了开伯尔-普什图省结核病数据中复杂TB行为解的存在性和唯一性。Ghoreishi等人（2023）利用Caputo广义Taylor级数提出了一类新的分数阶Runge-Kutta（FORK）方法以求解分数微分方程，发展了显式和隐式方案，并对分数导数的非局部特性进行了调整，收敛性分析和数值实验证明了该方法的准确性和鲁棒性。Al-Shimmary（2021）推导了一种新的4阶分数阶Runge-Kutta方法（4sFRKM）用于求解时间分数阶初值问题，获得了良好的数值解，使该方法具有实用价值。作者研究了所提方法的稳定性并与精确解进行对比，描绘了相应的稳定区域，还通过示例数值算例验证了该方法的准确性和效率，用表格和图形证实了结果。Jawias等人（2010）针对线性常微分方程（LODEs）创建了一种新型四阶对角隐式Runge-Kutta（DIRK）方法，通过利用Butcher误差方程并利用LODE特性省略其中一阶，简化了推导过程。与类似的显式Runge-Kutta（ERK）方法相比，稳定性分析和数值结果表明所提方法具有更优的稳定性和更低的整体误差。Flaig（2023）探讨了用于演化方程最优控制问题时间离散化的隐式Runge-Kutta方案，确保离散化与优化可交换，并为特定配置方法（如Gauss、Radau和Lobatto方案，最高六阶）建立了高阶要求自然满足的条件。尽管已有大量关于数学建模的研究，包括聚焦于开伯尔-普什图省结核病动力学的研究（Ullah et al., 2019; Khan et al., 2019），但在结核病传播的理解上仍存在显著空白。此前区域建模研究主要关注确定性传播动力学，既未考虑疫苗接种、治疗和恢复路径，也未将参数估计、分数阶动力学、最优控制和全局敏感性分析（如偏序相关系数PRCC）充分整合到单一建模框架中。具体而言，在巴基斯坦数据校准的结核病模型中，记忆依赖性疾病动力学的综合效应和理论技术尚未得到充分重视。因此，需要一个更全面的框架，不仅刻画结核病传播并确立其数学性质，还需使用实证流行病学数据校准模型、识别关键参数、考察分数阶效应并评估可行的控制措施。本研究开发了一个SVITR-TB框架，结合分数阶建模、最优控制分析、归一化敏感性和PRCC分析、确定性建模以及真实数据参数估计，以填补这些空白。这种综合方法为巴基斯坦的结核病传播模式、参数不确定性、记忆效应以及可能的预防和控制措施提供了互补性洞见。

本研究的主要贡献在于开发并评估了一种基于经典微分方程和分数阶形式建立的SVITR-TB结核病模型，以及一个理想控制框架，用于评估巴基斯坦可能的干预策略：

• 开发一个确定性的SVITR-TB模型，并确定其基本数学性质，如适定性、正性和有界性，以解释巴基斯坦结核病传播的动力学。
• 采用基于优化的参数估计过程估算未知模型参数，并使用巴基斯坦实际的结核病发病率数据（2000-2020年）对建议的模型进行校准。
• 通过归一化敏感性分析和偏秩相关系数（PRCC）确定影响结核病传播的最重要因素。
• 通过将确定性模型扩展到分数阶框架，探索记忆依赖动力学对结核病传播的可能影响。
• 使用合适的数值方案（如非标准有限差分（NSFD）方案和分数阶RK4方案）研究建议模型的数值动力学。
• 建立一个包含意识努力（u?）和加强疫苗接种（u?）的理想控制模型，以评估其在降低结核病传播方面的潜在效果。
• 通过将理论分析、真实数据校准、敏感性分析、分数阶动力学和最优控制整合到一个综合建模框架中，全面理解巴基斯坦结核病的传播与控制。

本研究的目标兼具方法论和流行病学性质。在方法论方面，我们开发并评估了一个确定性的SVITR-TB模型，将其扩展到分数阶框架，从实际的结核病数据中估计参数，进行PRCC和归一化敏感性分析，并创建理想控制模型。在流行病学方面，我们的目标是了解巴基斯坦结核病传播的动力学，识别关键因素，并评估意识宣传活动（u?）和加强疫苗接种（u?）对减少结核病传播的潜在影响。

**SVITR-TB模型构建**

为了评估疾病传播的动力学并评估不同干预策略的效果，我们利用数学模型。我们研究一种通常由结核分枝杆菌（MTB）引起的传染性疫情，并设计了一个微分比较模型，以建立SVITR-TB数学模型形式（1）。建议的SVITR-TB传播动力学模型在任何时间t将总人口N(t)分为五个互斥的类别：易感者（S(t)）、已接种疫苗者（V(t)）、潜伏感染者（I(t)）、活跃感染者（T(t)）和康复者（R(t)）。

**SVITR-TB模型数学分析**

在本节中，我们分析系统（1）在t≥0时的基本动力学性质，建立其解的正性、有界性和适定性。

设x:[0→∞)?R?由x=(S+V+I+T+R)?给出，五个状态变量包含在R?中的列向量里。设g(x)=(g?(x)+g?(x)+g?(x)+g?(x)+g?(x))?由g:R??R?给出，一个函数，其中R?定义向量值为，g?(x)=β?kSIN?(ρ+d?)S，g?(x)=ρS?(φ+d?)V，g?(x)=kSIN?(λ+μ+d?)I，g?(x)=μI?(η+d?)T，g?(x)=φV+ηT?d?R。

系统（2）可表达为...

**R?的表达式**

本节严格分析了结核病流行模型的三个基本数学特征：阈值参数、系统平衡点及其稳定性。基本再生数R?是衡量单个初始病例在完全易感人群中引起的平均二次感染数的关键参数，它控制着传染病传播的动力学。关键在于，这个阈值决定了疾病是否会持续...

**SVITR-TB模型稳定性分析**

本节研究系统平衡点的局部和全局稳定性特征。提供了这些稳定性性质的全面分析，详细的数学表述在定理5.1-5.4中给出。

**模型拟合**

为了确定未知模型参数并评估建议的SVITR-TB模型对观测结核病发病率数据的再现程度。使用有界非线性优化方法，通过最小化观测数据与模型预测之间的差异来进行参数估计。

为了估算未知参数k、ρ、μ、φ和η，最小化了观测结核病发病率数据与模型预测结核病发病率数据之间的平方差之和。有限记忆Broyden Fletcher方法...

**SVITR-TB模型参数敏感性分析**

本研究旨在评估重要因素对特定人群中疾病传播发展或管理的影响。我们将SVITR-TB模型再生数关于任何特定参数的敏感性系数记为ε。可表示为如下形式（Aalam and Pongsumpun, 2026b, Ahmad et al., 2025）：R?∩ε=(?R?/?ε)·(ε/R?)。其中，R?=k d?/((ρ+d?)(λ+μ+d?))，对R?中每个相关参数进行敏感性确定，并将结果展示在图3中。

**数值分析**

Ronald Mickens发明了非标准有限差分（NSFD）方案，旨在离散化连续的生物和流行病学模型，在这些模型中，非物理行为（如负种群密度或数值混沌）经常导致标准数值求解器（如向前Euler或经典Runge-Kutta方法）失败。NSFD框架的主要优势是它能够保持与基础连续微分方程系统的动力学一致性。

**分数阶模型**

经典SVITR-TB模型被扩展到Caputo分数阶框架，并使用分数阶RK4技术进行数值求解，以更好地研究结核病传播的时间动力学。引入分数阶导数主要是因为分数阶形式为表达记忆依赖的疾病动力学提供了更灵活的框架。相比之下，经典模型假设当前的变化率仅取决于当前状态。先前的感染...

**SVITR-TB模型的最优控制分析**

为了评估其对模型（1）中SVITR-TB动力学的影响，建立了一个最优控制框架。我们纳入了两个时间依赖的控制函数：u?(t)代表针对易感人群S的公众意识宣传活动（平面媒体、数字媒体、社交媒体），u?(t)代表针对已接种疫苗人群V的加强疫苗接种策略。通过定义目标泛函来最小化疾病流行率...来建立最优控制问题。

**结论**

本研究引入了一种确定性的SVITR框架来分析巴基斯坦的结核病传播动力学和缓解策略。通过确认系统的适定性、正性和有界性，我们确保数学解具有生物学意义和物理可靠性。此外，利用下一代矩阵技术导出了基本再生数（R?）。同时，采用了Routh-Hurwitz准则和Lyapunov函数...

**局限性与未来工作**

建议的研究提供了关于结核病传播和控制的见解；然而，应注意几个局限性。首先，该模型没有明确考虑社会经济、地区或人口统计异质性，尽管它假设了同质的人口混合。由于巴基斯坦各地生活条件、医疗保健获取、人口特征和结核病传播模式存在巨大差异，这些变量可能影响疾病动力学和控制结果。

**CRediT作者贡献声明**

Balal Aalam：写作-审阅与编辑，写作-初稿撰写，形式分析，数据管理，概念化。Omar Kihel：写作-审阅与编辑，写作-初稿撰写。Puntani Pongsumpun：写作-初稿撰写，可视化，验证，监督，方法学，资金获取，形式分析，数据管理，概念化。

**资金声明**

作者未收到本研究报告出版的任何资金。

**未引用的参考文献**

Aalam and Pongsumpun (2026a)

**利益冲突声明**

作者声明没有已知的可能影响本研究报告工作的财务利益或个人关系。

**致谢**

本研究得到了泰国玛格达理工学院科学学院的资金支持。Balal Aalam感谢玛格达理工学院科学学院提供的研究生奖学金。额外资金支持来自玛格达理工学院研究助理（R.A.）补助金（KREF026724）。

Balal Aalam | Omar Kihel | Puntani Pongsumpun