对一类猴痘模型的最优控制

Optimal Control of One Kind of Mpox Model

摘要: 作为一类以动物宿主接触及性传播为主要途径的病毒性传染病，猴痘在2022~2023年全球范围内呈现出突发性与快速传播特征。本研究聚焦于无动物宿主输入的国家中，女性与男性患者感染特征的异质性，进而构建性别分组的SEIR动力学模型以解析病毒传播动态。通过引入庞特里亚金极大值原理(Pontryagin’s maximum principle)，系统推导了模型的最优控制策略。基于美国疫情数据的参数化数值模拟表明，干预措施需优先针对高风险人群——尤其是男性群体，且疫苗早期覆盖率对疫情控制具有显著敏感性。研究结果为政府防控策略提供了理论依据，强调在疫情初期即需优化男性群体疫苗接种优先级，以有效抑制传播链。

Abstract: As a viral infectious disease primarily transmitted through animal-to-human contact and sexual transmission, monkeypox exhibited sudden and rapid global spread characteristics during 2022~2023. This study focuses on the heterogeneity of infection characteristics between female and male patients in countries without animal host introduction and subsequently constructs a gender-stratified SEIR dynamic model to analyze the transmission dynamics of the virus. By applying Pontryagin’s maximum principle, the study systematically derives the model’s optimal control strategies. Parameterized numerical simulations based on U.S. outbreak data indicate that intervention measures should prioritize high-risk populations—particularly the male group—and that early vaccine coverage rates are highly sensitive in controlling the epidemic. The findings provide theoretical support for government prevention and control policies, emphasizing the need to optimize vaccine prioritization for the male population in the early stage of the epidemic to effectively curb transmission chains.

1. 引言

2022年，多国暴发的猴痘疫情使这一沉寂多年的病毒重新引发关注[1]。作为痘病毒科正痘病毒属的双链DNA人畜共患病病原体，猴痘病毒可分为西非分支与刚果盆地分支[2]。该病毒最早于1958年在实验猴群痘样疾病暴发研究中被发现，故得名monkeypox (简称Mpox) [3]。首例人类感染病例发现于刚果民主共和国，患者曾被误诊为天花。2021年前，猴痘主要在中西非地区呈地方性流行，通过接触感染动物传播，人际传播链较短，以儿童及青少年的散发病例和聚集性病例为主，偶尔通过家庭或旅行方式扩散至其他国家和地区[4]。自2022年5月起，全球多国出现以男男性行为者(MSM)间性接触为主要传播途径的疫情，病例多通过大型集会及后续社区MSM网络传播，并扩散至其他国家。2022年后各国猴痘患者症状普遍较轻，死亡病例多见于未接受抗逆转录病毒治疗的HIV感染者等免疫缺陷人群[5]。

目前认为猴痘病毒宿主主要为非洲啮齿类动物，感染啮齿类、猴类及猿类等灵长类动物均可成为传染源[6]。该病毒通过黏膜及破损皮肤侵入人体，主要传播途径包括：直接接触病例病变皮肤或黏膜、接触病毒污染物品、长时间近距离吸入病例呼吸道飞沫，以及接触感染动物呼吸道分泌物、病灶渗出液、血液等体液或被其抓咬伤[7]。此外，猴痘病毒对干燥与低温具有较强抗性，可在痂皮、土壤及衣物寝具等物体表面存活数月，存在通过污染寝具感染的案例。T. T. Pattiyakumbura等研究表明，猴痘存在垂直传播风险，需重视孕妇健康管理[8]。人群普遍易感，但天花疫苗接种者可能获得一定交叉保护。

猴痘潜伏期为5~21天(多为6~13天) [9]。患者在出现症状至皮疹自然脱痂形成新皮肤期间均具有传染性。早期症状以发热、头痛、背痛、肌痛为主，可伴淋巴结肿大，皮肤黏膜皮疹多出现于热退后，少数可先于全身症状[10]。皮疹通常经历斑疹、丘疹、水疱、脓疱、结痂至脱落等阶段，不同形态皮疹可同时存在，常伴明显瘙痒疼痛。脱痂后可遗留红斑或色素沉着甚至瘢痕，瘢痕可持续数年。症状持续约2~4周，免疫缺陷患者病程可能延长[10]。作为自限性疾病，多数病例可自愈，但存在重症及死亡病例，主要集中于儿童、孕妇及免疫缺陷人群[8]。现有JYNNEOS® (Imvamune/Imvanex)与ACAM2000两种疫苗可降低感染风险，但尚未实现全球普及[11] [12]。

在数学模型研究方面，C. P. Bhunu与S. Mushayabasa [13]于2011年建立包含人–人及啮齿类–人传播的SIR模型，揭示人群免疫状态与感染后恢复方式的关系。2012年，C. P. Bhunu等[14]通过HIV/AIDS与猴痘共感染研究，发现两种感染的相互促进效应。2017年，S. Usman与I. I. Adamu [15]构建考虑潜伏期的SEIR模型，将疫苗接种与外科干预作为控制策略。2018年，P. C. Emeka等[16]通过SVIR模型探讨免疫系统差异、感染效应及接种率对传播的影响。2022年非流行国突发疫情后，相关研究显著增加：O. J. Peter与F. A. Oguntolu等既考虑病例隔离[17]又采用分数阶模型分析传播动力学[17]；2023年Yang等[18]研究高风险人群(主要为MSM)与普通人群间的传播差异，Batiha等[19]同时考虑休假模式与分数阶模型，Al-Shomrani等[20]建立包含继发感染的模型，O. C. Collins与K. J. Duffy [21]则分析传播过程中的随机性。

本文结构如下：第二节构建包含垂直传播的无动物宿主国家猴痘传播模型；第三节引入最优控制理论探讨疫苗干预效应；第四节基于美国2022年6月13日至2024年1月29日疫情数据进行数值模拟；最后总结研究结论。

2. 具有垂直传播的猴痘模型

本文构建的猴痘传播模型基于O. J. Peter等[22]与F. B. Agusto等[23]的研究框架。根据世界卫生组织(WHO)数据[1]，尽管女性患者数量显著少于男性，但其临床症状更为严重。为此，本研究将人群划分为两个独立群体：女性(标记为f)与男性(标记为m)。模型中定义以下状态变量：$S_f$ ($S_m$ )表示女性(男性)易感人群，$E_f$ ($E_m$ )人群，$I_f$ ($I_m$ )表示女性(男性)易感人群，$R_f$ ($R_m$ )表示女性(男性)易感人群。基于WHO数据[1]，模型设定以下关键假设：

(1) 新生儿性别比恒定为1:1；

(2) 现有证据未显示恢复率存在性别差异，故模型中两群体恢复率参数设为一致。

Figure 1. Schematic of Mpox transmission

图1. 猴痘传播的流程图

图1中，${\gamma }_f$ 和${\gamma }_m$ 分别代表女性和男性的感染力度，有

$\begin{array}{l}{\gamma }_f\left(t\right)=\frac{{\beta }_{fm}{I}_m\left(t\right)+{\beta }_{ff}{I}_f\left(t\right)}{N_f\left(t\right)},\\ {\gamma }_m\left(t\right)=\frac{{\beta }_{mm}{I}_m\left(t\right)+{\beta }_{mf}{I}_f\left(t\right)}{N_m\left(t\right)}.\end{array}$ (1)

且$N_f,N_m$ 代表女性和男性的总人口。

$N_f\left(t\right)=S_f\left(t\right)+E_f\left(t\right)+I_f\left(t\right)+R_f\left(t\right),$

$N_m\left(t\right)=S_m\left(t\right)+E_m\left(t\right)+I_m\left(t\right)+R_m\left(t\right).$

图1中的模型的微分方程形式由下面的系统(2)给出。

$\frac{\text{d}{S}_f}{\text{d}t}=\left({\pi }_B-q_E{E}_f-q_I{I}_f-q_R{R}_f\right)-{\gamma }_f\left(t\right)S_f-\mu S_f,$ (2)

$\frac{\text{d}{E}_f}{\text{d}t}=q_E{E}_f+{\gamma }_f\left(t\right)S_f-\left(\mu +\alpha \right)E_f,$

$\frac{\text{d}{I}_f}{\text{d}t}=q_I{I}_f+\alpha E_f-\left(\mu +\sigma \right)I_f,$

$\frac{\text{d}{R}_f}{\text{d}t}=q_R{R}_f+\sigma I_f-\mu R_f,$

$\frac{\text{d}{S}_m}{\text{d}t}=\left({\pi }_B-q_E{E}_f-q_I{I}_f-q_R{R}_f\right)-{\gamma }_m\left(t\right)S_m-\mu S_m,$

$\frac{\text{d}{E}_m}{\text{d}t}=q_E{E}_f+{\gamma }_m\left(t\right)S_m-\left(\mu +\alpha \right)E_m,$

$\frac{\text{d}{I}_m}{\text{d}t}=q_I{I}_f+\alpha E_m-\left(\mu +\sigma \right)I_m,$

$\frac{\text{d}{R}_m}{\text{d}t}=q_R{R}_f+\sigma I_m-\mu R_m\text{.}$

系统(2)中的所有参数值由表1给出。

Table 1. Definition of parameters

表1. 参数的定义

参数	定义
${\pi }_B$	新生女性、男性婴儿的数量
$q_E$	一个潜伏期女性生出一个潜伏期女性、男性婴儿的概率
$q_I$	一个感染期女性生出一个感染期女性、男性婴儿的概率
$q_R$	一个治愈期女性生出一个有抵抗力的女性、男性婴儿的概率
${\beta }_{ff}$	感染期女性对易感期女性的疾病传播率
${\beta }_{fm}$	感染期男性对易感期女性的疾病传播率
${\beta }_{mm}$	感染期男性对易感期男性的疾病传播率
${\beta }_{mf}$	感染期女性对易感期男性的疾病传播率
$\sigma$	恢复率
$\mu$	自然死亡率
$\alpha$	潜伏期个体到疾病下一阶段的进展率

系统(2)的最优控制

在本节中，将原系统(2)扩展引入两个控制变量：$v_1$ 和$v_2$ 。其中$v_1$ 表示女性群体的人均疫苗接种率，$v_2$ 表示男性群体的人均疫苗接种率。通过这一扩展，可进一步探讨已实施疫苗接种的国家现状，以及政府应推行的政策建议。扩展后的模型如系统(3)所示。

$\frac{\text{d}{S}_f}{\text{d}t}=\left({\pi }_B-q_E{E}_f-q_I{I}_f-q_R{R}_f\right)-{\gamma }_f{S}_f-v_1\left(t\right)S_f-\mu S_f,$ (3)

$\frac{\text{d}{E}_f}{\text{d}t}=q_E{E}_f+{\gamma }_f{S}_f-\left(\mu +\alpha \right)E_f,$

$\frac{\text{d}{I}_f}{\text{d}t}=q_I{I}_f+\alpha E_f-\left(\mu +\sigma \right)I_f,$

$\frac{\text{d}{R}_f}{\text{d}t}=q_R{R}_f+\sigma I_f+v_1\left(t\right)S_f-\mu R_f\text{,}$

$\frac{\text{d}{S}_m}{\text{d}t}=\left({\pi }_B-q_E{E}_f-q_I{I}_f-q_R{R}_f\right)-{\gamma }_m{S}_m-v_2\left(t\right)S_m-\mu S_m,$

$\frac{\text{d}{E}_m}{\text{d}t}=q_E{E}_f+{\gamma }_m{S}_m-\left(\mu +\alpha \right)E_m,$

$\frac{\text{d}{I}_m}{\text{d}t}=q_I{I}_f+\alpha E_m-\left(\mu +\sigma \right)I_m,$

$\frac{\text{d}{R}_m}{\text{d}t}=q_R{R}_f+\sigma I_m+v_2\left(t\right)S_m-\mu R_m\text{.}$

为界定控制策略的边界，设定$v_{1,\max }$ 与$v_{2,\max }$ 分别为$v_1\left(t\right)$ 和$v_2\left(t\right)$ 的最大值。基于此可构建策略$v_1\left(t\right)$ 与$v_2\left(t\right)$ 的控制集$U$ ，其数学定义为：

$U=\left\{\left(v_1,v_2\right)\in {\left[L^{\inf }\left(0,T\right)\right]}^2|0\le v_i\left(t\right)\le v_{i,\max },i=1,2\right\}.$ (4)

在控制集$U$ 中，$T$ 表示策略的实施周期。通过将系统(3)的解定义为状态变量，可构造目标泛函$J\left(v_1,v_2\right)$ 如下：

$\begin{array}{c}J\left(v_1,v_2\right)=J\left(v_1\right)+J\left(v_2\right)\\ ={\displaystyle {\int }_0^T\left[b_1\left(q_E{E}_f+q_I{I}_f+{\gamma }_f{S}_f\right)+A_1{v}_1{S}_f+{ϵ}_1{v}_1^2\right]\text{d}t}\\ +{\displaystyle {\int }_0^T\left[b_2\left(q_E{E}_f+q_I{I}_f+{\gamma }_m{S}_m\right)+A_2{v}_2{S}_m+{ϵ}_2{v}_2^2\right]\text{d}t}.\end{array}$ (5)

其中，参数$b_i$ 表示女性或男性群体中每新增病例的成本；$A_i$ 表示女性或男性群体中单次疫苗接种的成本；${ϵ}_i$ 表征女性或男性群体疫苗接种的非线性成本系数。根据庞特里亚金极大值原理[22]，为获得满足条件的最优控制向量$\left(v_1^{*},v_2^{*}\right)$ ，

$J\left(v_1^{*},v_2^{*}\right)={\min }_{U}J\left(v_1,v_2\right),$ (6)

需构建相应的哈密顿函数：

$\begin{array}{c}H\left[x\left(t\right),v\left(t\right),\lambda \left(t\right),t\right]\triangleq L\left[x\left(t\right),v\left(t\right),t\right]+{\lambda }^{\text{T}}\left(t\right)f\left[x\left(t\right),v\left(t\right),t\right]\\ =b_1\left(q_E{E}_f+q_I{I}_f+{\gamma }_f{S}_f\right)+A_1{v}_1{S}_f+{ϵ}_1{v}_1^2\\ +b_2\left(q_E{E}_f+q_I{I}_f+{\gamma }_m{S}_m\right)+A_2{v}_2{S}_m+{ϵ}_2{v}_2^2\\ +{\lambda }_1\left[\left({\pi }_B-q_E{E}_f-q_I{I}_f-q_R{R}_f\right)-{\gamma }_f{S}_f-v_1\left(t\right)S_f-\mu S_f\right]\\ +{\lambda }_2\left[q_E{E}_f+{\gamma }_f{S}_f-\left(\mu +\alpha \right)E_f\right]\\ +{\lambda }_3\left[q_I{I}_f+\alpha E_f-\left(\mu +\sigma \right)I_f\right]\\ +{\lambda }_4\left[q_R{R}_f+\sigma I_f+v_1\left(t\right)S_f-\mu R_f\right]\\ +{\lambda }_5\left[\left({\pi }_B-q_E{E}_f-q_I{I}_f-q_R{R}_f\right)-{\gamma }_m{S}_m-v_2\left(t\right)S_m-\mu S_m\right]\\ +{\lambda }_6\left[q_E{E}_f+{\gamma }_m{S}_m-\left(\mu +\alpha \right)E_m\right]\\ +{\lambda }_7\left[q_I{I}_f+\alpha E_m-\left(\mu +\sigma \right)I_m\right]\\ +{\lambda }_8\left[q_R{R}_f+\sigma I_m+v_2\left(t\right)S_m-\mu R_m\right].\end{array}$ (7)

该哈密顿函数$H$ 需满足伴随方程与状态方程：

$\dot{\lambda }\left(t\right)=-\frac{\partial H}{\partial x\left(t\right)},\dot{x}\left(t\right)=\frac{\partial H}{\partial \lambda \left(t\right)}$ (8)

其中

$x\left(t\right)=\left(S_f\left(t\right),E_f\left(t\right),I_f\left(t\right),R_f\left(t\right),S_m\left(t\right),E_m\left(t\right),I_m\left(t\right),R_m\left(t\right)\right),$ (9)

$v\left(t\right)=\left(v_1\left(t\right),v_2\left(t\right)\right),$

$\lambda \left(t\right)=\left[{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_8\right],$

$f\left[x\left(t\right),v\left(t\right),t\right]=\dot{x}\left(t\right)$ 由系统(3)得到。

依据庞特里亚金极大值原理，若上述条件成立，则使目标泛函(5)最小化的$\left(v_1^{*},v_2^{*}\right)$ 在每一时刻均使哈密顿函数关于控制变量取得极小值。因此，在控制集内域上求解$\left(v_1^{*},v_2^{*}\right)$ 需满足以下最优性条件：

$\frac{\partial H}{\partial v_i}=0,i=1,2.$ (10)

在计算方程(10)前，首先通过求解下式检验$H$ 关于$v_i$ 的凸性：

$\frac{{\partial }^{2}H}{\partial v_i^2}=2{ϵ}_i>0,i=1,2.$ (11)

方程(11)保证最优控制$\left(v_1^{*},v_2^{*}\right)$ 确为极小值。继而通过求解方程(10)可获得最优控制向量$\left(v_1^{*},v_2^{*}\right)$ 。

$A_1{S}_f+2{ϵ}_1{v}_1-{\lambda }_1{S}_f+{\lambda }_4{S}_f=0,$

$A_2{S}_m+2{ϵ}_2{v}_2-{\lambda }_5{S}_m+{\lambda }_8{S}_m=0.$

由此可得

$v_1=\frac{\left({\lambda }_1-{\lambda }_4-A_1\right)S_f}{2{ϵ}_1},$

$v_2=\frac{\left({\lambda }_5-{\lambda }_8-A_2\right)S_m}{2{ϵ}_2}.$

最终获得最优控制向量表达式为

$\begin{array}{l}{v}_1^{*}=\min \left\{v_{1,\max },\max \left\{0,\frac{\left({\lambda }_1-{\lambda }_4-A_1\right)S_f}{2{ϵ}_1}\right\}\right\},\\ v_2^{*}=\min \left\{v_{2,\max },\max \left\{0,\frac{\left({\lambda }_5-{\lambda }_8-A_2\right)S_m}{2{ϵ}_2}\right\}\right\}.\end{array}$ (12)

为构建系统(12)，需满足方程(8)中给定的伴随方程与状态方程，其对应的伴随微分方程组为

$\begin{array}{l}{{\lambda }^{\prime }}_1=-[b_1\frac{{\beta }_{fm}{I}_m+{\beta }_{ff}{I}_f}{N_f}+A_1{v}_1+{\lambda }_1\left(-\frac{{\beta }_{fm}{I}_m+{\beta }_{ff}{I}_f}{N_f}-v_1-\mu \right)\\ +{\lambda }_2\frac{{\beta }_{fm}{I}_m+{\beta }_{ff}{I}_f}{N_f}+{\lambda }_4{v}_1],\end{array}$ (13)

${{\lambda }^{\prime }}_2=-\left[\left(b_1+b_2-{\lambda }_1+{\lambda }_2-{\lambda }_5+{\lambda }_6\right)q_E-{\lambda }_2\left(\mu +\alpha \right)+{\lambda }_3\alpha \right],$

$\begin{array}{l}{{\lambda }^{\prime }}_3=-[\left(b_1+b_2-{\lambda }_1+{\lambda }_3-{\lambda }_5+{\lambda }_7\right)q_I+\left(b_1-{\lambda }_1+{\lambda }_2\right)\frac{{\beta }_{ff}{S}_f}{N_f}\\ +\left(b_2-{\lambda }_5+{\lambda }_6\right)\frac{{\beta }_{mf}{S}_m}{N_m}-{\lambda }_3\left(\mu +\sigma \right)+{\lambda }_4\sigma ],\end{array}$

${{\lambda }^{\prime }}_4=-\left[\left(-{\lambda }_1+{\lambda }_4-{\lambda }_5+{\lambda }_8\right)q_R-{\lambda }_4\mu \right],$

$\begin{array}{l}{{\lambda }^{\prime }}_5=-[b_2\frac{{\beta }_{mm}{I}_m+{\beta }_{mf}{I}_f}{N_m}+A_2{v}_2+{\lambda }_5\left(-\frac{{\beta }_{mm}{I}_m+{\beta }_{mf}{I}_f}{N_m}-v_2-\mu \right)\\ +{\lambda }_6\frac{{\beta }_{mm}{I}_m+{\beta }_{mf}{I}_f}{N_m}+{\lambda }_8{v}_2],\end{array}$

${{\lambda }^{\prime }}_6=-\left[-{\lambda }_6\left(\mu +\alpha \right)+{\lambda }_7\alpha \right],$

${{\lambda }^{\prime }}_7=-\left[\left(b_1-{\lambda }_1+{\lambda }_2\right)\frac{{\beta }_{fm}{S}_f}{N_f}+\left(b_2-{\lambda }_5+{\lambda }_6\right)\frac{{\beta }_{mm}{S}_m}{N_m}-{\lambda }_7\left(\mu +\sigma \right)+{\lambda }_8\sigma \right],$

${{\lambda }^{\prime }}_8={\lambda }_8\mu ,$

${\lambda }_j\left(T\right)=0,j=1,2,\cdots ,8.$

3. 对系统(2)的数值模拟

本节采用MCMC算法估计系统(2)的未知参数，并利用四阶龙格–库塔方法求解最优控制值。数据选取方面采用美国2022年6月至2024年1月期间的每周新增猴痘病例数据。鉴于美国政府自2022年8月16日起开始实施Mpox疫苗接种计划(使用JYNNEOS® (Imvamune或Imvanex)和ACAM2000疫苗) [12]，选取2022年6月至8月数据用于系统(2)的参数估计，并以2022年8月至2024年1月数据验证控制策略$v_1$ 和$v_2$ 的有效性。

3.1. 数据拟合和参数估计

在应用MCMC算法前，先基于现实数据对部分初始值和参数进行先验估计：

(i) 初始值$S_{f0}$ 和$S_{m0}$ 采用美国截至2021年末的统计女性及男性总人口数据；

(ii) 参数${\pi }_B$ 取2022年美国周平均新生儿数量，由当年总出生人口除以52周获得；

(iii) 参数$\mu ,$ 取2022年美国周平均死亡率，由当年总死亡人数除以52周计算；

(iv) 由于女性平均妊娠周期约为40周，超过美国Mpox疫情的持续时间，在模拟过程中我们将垂直传播参数$q_E,q_I,q_R$ 设为0。

下面使用MCMC方法来估计其余的初始值和参数。参数的先验区间根据其生物学意义进行设置，时间步长为1周。对于参数${\beta }_{ff}$ ，${\beta }_{fm}$ ，${\beta }_{mf}$ 和${\beta }_{mm}$ ，它们代表了猴痘通过性接触传播的概率，以${\beta }_{fm}$ 为例，其范围选择$0\le {\beta }_{fm}\le 10$ 以捕捉广泛的传播情景，从无传播(${\beta }_{fm}=0$ )到高效传播(${\beta }_{fm}=10$ )。这个范围考虑了病原体传播性和接触模式的变化，同时确保了模型在不同条件下动态的稳健性。参数$\sigma$ 代表恢复率，根据世界卫生组织的报告，猴痘病毒的恢复期从2周到4周不等，因此将$\sigma$ 的范围设为$\frac{1}{2}\le \sigma \le \frac{1}{4}$ 。此外，根据美国疾病控制和预防中心(CDC)的报告，美国猴痘病毒的潜伏期从3天到17天不等，因此取$\alpha$ 的范围为$\frac{7}{17}\le \alpha \le \frac{7}{3}$ 。考虑到美国直到2022年6月13日的总计病例为240，为了确保模型具有足够的灵活

度，模拟时将其他初值的上限设为1000，最终模拟得到的值如表2所示。

Table 2. Values of initial values and parameters estimated from MCMC and the reality

表2. 从MCMC中获得的以及根据真实数据估算的初始值和参数值

参数	数值	标准差	范围	来源
${\pi }_B$	70,534	-	-	(ii)
$q_E$	0	-	-	(iv)
$q_I$	0	-	-	(iv)
$q_R$	0	-	-	(iv)
${\beta }_{ff}$	0.120188866	0.057988	$\left[0,10\right]$	MCMC
${\beta }_{fm}$	0.208074914	0.09142	$\left[0,10\right]$	MCMC
${\beta }_{mm}$	1.379101637	0.22812	$\left[0,10\right]$	MCMC
${\beta }_{mf}$	0.031029074	0.02861	$\left[0,10\right]$	MCMC
$\sigma$	0.278073917	0.12085	$\left[\frac{1}{4},\frac{1}{2}\right]$	MCMC
$\mu$	0.0048575	-	$\left[0,1\right]$	(iii)
$\alpha$	0.823686743	0.15905	$\left[\frac{7}{17},\frac{3}{7}\right]$	MCMC
$S_{f0}$	$168\times {10}^8$	-	-	(i)
$E_{f0}$	2	176.34	$\left[0,1000\right]$	MCMC

续表

$I_{f0}$	4	30.86	$\left[0,1000\right]$	MCMC
$R_{f0}$	0	3.16	$\left[0,1000\right]$	MCMC
$S_{m0}$	$165.37\times {10}^8$	-	-	(i)
$E_{m0}$	120	191.98	$\left[0,1000\right]$	MCMC
$I_{m0}$	97	24.787	$\left[0,1000\right]$	MCMC
$R_{m0}$	4	253.38	$\left[0,1000\right]$	MCMC

图2展示了MCMC模拟后的拟合结果，可以看到模型得到的数值解和数据拟合得较好。

Figure 2. MCMC fitting result

图2. MCMC拟合结果

3.2. 最优控制

本节为契合美国实际情况，将控制变量($v_1,v_2$ )于2022年6月13日8周后投入应用。初始值与参数均取自表2。根据E. Howerton等[3]的研究，设定各成本系数为：单例病例成本$b_1=b_2=1$ ，单次疫苗接种成本$A_1=A_2=0.125$ ，疫苗接种非线性成本${ϵ}_1={ϵ}_2=5\times {10}^8$ ，女性与男性最大疫苗接种率分别为$v_{1,\max }=0.35$ 和$v_{2,\max }=0.55$ 。图3中曲线表示应用不同控制策略后模型的输出结果，红点代表美国2022年6月13日至2024年1月29日实际新增病例数据。图中展示了美国2022年6月至2024年1月猴痘实际新增病例数据(红点)与系统(3)自2022年8月起应用不同控制策略$v_1,v_2$ 的模拟结果(曲线)。控制策略选取了四种形式：$v_1=\frac{1}{2}{v}_{1,\max }$ 且$v_2=\frac{1}{2}{v}_{2,\max }$ ；$v_1=v_{1,\max }$ 且$v_2=\frac{1}{4}{v}_{2,\max }$ ；$v_1=0$ 且$v_2=\frac{1}{2}{v}_{2,\max }$ ；$v_1=v_{1,\max }$ 且$v_2=v_{2,\max }$ 。图3(a)、图3(c)、图3(e)显示固定控制策略的结果，图3(b)、图3(d)、图3(f)展示第3节计算所得最优控制策略的结果。

(a) (b)

(c) (d)

(e) (f)

Figure 3. Different levels of control

图3. 不同水平下的控制

图3(a)显示，当固定女性周接种率$v_1$ 时，调整男性周接种率$v_2$ 对女性群体中猴痘病毒传播具有显著影响。而图3(c)中深蓝色虚线(对应$v_1=\frac{1}{2}{v}_{1,\max },v_2=\frac{1}{2}{v}_{2,\max }$ )与紫色实线(对应$v_1=0,v_2=\frac{1}{2}{v}_{2,\max }$ )基本重合，表明调整女性接种率$v_2$ 对男性群体中病毒传播影响甚微。值得注意的是，采用$v_1=v_{1,\max },v_2=v_{2,\max }$ 策略的模拟结果与实际数据最为接近，而最优控制结果在2周后的下降速率较实际数据更为平缓。图4(a)显示最优控制强度的时间演化过程，其中红色曲线表征女性接种率$v_1$ ，绿色曲线表征男性接种率$v_2$ ；图4(b)对比不同控制策略的实施成本。其中J表示最优控制总成本，Jb对应$v_1=\frac{1}{2}{v}_{1,\max }$ 且$v_2=\frac{1}{2}{v}_{2,\max }$ 策略，Jc对应$v_1=v_{1,\max }$ 且$v_2=\frac{1}{4}{v}_{2,\max }$ 策略，Jd对应$v_1=0$ 且$v_2=\frac{1}{2}{v}_{2,\max }$ 策略，Jmax对应$v_1=v_{1,\max }$ 且$v_2=v_{2,\max }$ 策略。图4(a)表明，在最优控制策略下，$v_1$ 仅需维持最大值至2.6周后，$v_2$ 则需维持至3.3周后，随后均衰减至零。图4(b)显示，虽然最优控制对疫情抑制效果弱于实际情况，但其成本显著降低。这一现象提示政府仅需在疫情初期集中资源推进疫苗接种。

(a) (b)

Figure 4. Level of optimal control and cost

图4. 最优控制强度和成本

4. 讨论

本文构建了包含女性与男性群体的SEIR模型，用以研究无动物宿主国家猴痘病毒的传播动力学。基于庞特里亚金极大值原理，建立以新增病例成本与疫苗接种成本为指标的目标泛函$J$ ，推导出政府可采取的最优控制策略。选择美国2022年6月至2024年1月的数据进行验证(猴痘疫苗自2022年8月起可获取)，马尔可夫链蒙特卡洛模拟结果表明：若无干预措施，猴痘病毒将在美国持续流行；控制策略分析建议政府应重点针对男性群体，且在疫情初期即需全力推进疫苗接种。

作为曾局限于非洲的地方性传染病，猴痘在2022~2023年间于无动物宿主国家的突发性流行值得警惕。本文虽探讨了性别因素与疫苗接种的影响，但年龄结构、人口流动、自我隔离等其他因素仍需深入研究。唯有纳入最关键的影响因子，方能准确刻画猴痘病毒的真实传播动态。