模拟具有再感染机制的双菌株竞争：数学分析与流行病学意义

《Infectious Disease Modelling》：Modeling Two-Strain Competition with Reinfection: Mathematical Analyses and Epidemiological Implications

　　杰白|金王
辽宁大学数学与统计学院，沈阳110036，中国

摘要
我们提出了一种新的流行病模型，用于描述由相同病原体引起的两种不同菌株之间的相互作用。我们关注的是多种菌株之间的竞争以及再感染对疾病传播、扩散和进展的影响。我们分析了该数学模型的动态行为，特别是各种平衡点的存在及其稳定性。此外，我们利用公开可用的COVID数据进行了详细的数值研究以校准模型。我们的模拟结果揭示了新出现菌株与原有菌株之间相互作用的复杂模式，并强调了有效控制原发感染和再感染途径的重要性。

1. 引言
由病原体（细菌、病毒、真菌、寄生虫等）引起的人类感染通常涉及同一病原体的多个菌株。菌株代表病原微生物的遗传变体或亚型，这些变体在免疫原性、传染性、毒力、生长率和耐药性等特性上可能存在显著差异。尽管多菌株感染非常普遍，但其重要性和临床相关性仍不清楚[3]，[8]。数学建模作为流行病学研究的理论工具，可以提供有关不同菌株竞争、互利共生和进化的宝贵见解，并有助于我们更好地理解疾病动态和宿主-病原体相互作用[17]，[18]，[23]，[26]，[37]，[38]，[39]，[41]，[43]。
已经有许多关于多菌株传染病的数学模型被发表[5]，[13]，[20]，[22]，[25]，[26]，[28]，[29]，[32]，[33]，[34]，[35]，[39]，[40]，[42]。例如，Minayev和Ferguson开发了一套相对简单的确定性模型，用于描述多菌株病原体的传播动态，从而提高了生物学上的真实性[28]。Martcheva提出了一种具有周期系数的非自主多菌株流行病模型，并进行了详细分析，包括解的一般适定性和稳定性，特别是对于双菌株版本的模型[25]。Rong等人基于双菌株数学模型研究了丙型肝炎病毒感染，该模型纳入了野生型和耐药变体，并考察了反向突变和肝细胞增殖对突变病毒存在及变体之间竞争的影响[34]。Murall等人对人类乳头瘤病毒（HPV）感染进行了建模，并研究了哪些HPV类型的相互作用及其免疫反应产生了在患者中频繁观察到的相似模式[29]。
在过去几年中，许多建模研究致力于COVID-19。自2019年12月首次报告SARS-CoV-2病毒株以来，该病毒经历了大量突变。几种变体，如Alpha（B.1.1.7）、Beta（B.1.351）、Gamma（P.1）、Delta（B.1.617.2）和Omicron（B.1.1.529），导致了大规模感染，并被世界卫生组织列为“关注变异株”[49]。利用双菌株模型，Li等人评估了SARS-CoV-2变体对大流行进展和社会恢复的影响，并估计了相关疫苗的有效性[21]。Massard等人基于多菌株建模研究显示，Alpha、Beta和Gamma变体比原始病毒株更具传染性[27]。Thompson等人发现，如果之前的SARS-CoV-2感染导致了完全的交叉反应性免疫，那么新的变体需要比之前传播的变体更具传染性才能侵入人群[36]。Tang等人提出了一个包含多种病毒株的COVID-19模型，这些病毒株既有症状性传播途径也有无症状传播途径，并研究了Omicron变体在美国的出现[14]。其他多菌株COVID-19建模研究包括但不限于[2]，[6]，[9]，[19]，[24]。
尽管做出了这些努力，我们对涉及多个菌株的感染动态的定量理解仍然不足。实证研究表明，与原始菌株相比，新变体可能更具传染性、对抗疫苗的抵抗力更强，并且更快产生耐药性[8]，[10]。然而，这些差异对病原体传播和流行病/大流行进展（特别是疾病发病率、患病率和攻击率）的影响尚未得到充分研究。同时，从一个菌株的感染和康复通常会为同一菌株提供免疫力，但并不一定保护宿主免受其他菌株（尤其是新变体）的感染[30]。因此，量化由于病原体突变和进化而导致的再感染动态对于更好地理解其流行病学和地方性以及设计更有效的控制措施至关重要。
在本文中，我们提出了一个基于微分方程的数学模型，用于描述与传染病相关的两种不同菌株之间的相互作用。这两种菌株在传播率、康复率和疾病引起的死亡率方面有所不同。特别强调了再感染的影响，即从一种菌株感染中康复的个体会对另一种菌株易感。此外，我们纳入了可能包括迁移和自然出生与死亡在内的关键动态因素。结合这些因素，该模型使我们能够研究多种传播模式之间的相互作用、不同菌株的竞争、再感染的影响以及人口统计效应。
本文的其余部分组织如下：第2节描述了模型，第3节总结了主要数学结果，附录中提供了详细证明。第4节展示了数据拟合和模拟结果，第5节给出了结论和讨论。

2. 模型构建
我们提出了一个数学模型，该模型考虑了两种菌株之间的竞争以及与我们关注的病原体相关的再感染的影响。我们将人类群体划分为以下几类：易感个体（用S表示）、被菌株j感染的个体（用Ij表示）以及从菌株j感染中康复的个体（用Rj表示），其中j = 1, 2。我们假设来自Rj类别的康复个体对同一菌株j的再感染具有抵抗力，但仍有可能被另一种菌株感染。
模型描述如下：
(2.1) 易感个体以流入率Λ进入群体，并分别以初级感染率β1和β2感染菌株1和菌株2。来自菌株1和菌株2的感染个体分别具有康复率γ1和γ2以及疾病引起的死亡率w1和w2。从菌株1感染中康复的个体可能以再感染率再次感染菌株2，而从菌株2感染中康复的个体可能以再感染率再次感染菌株1。标准双线性发病率用于表示初次感染和再感染。此外，个体还具有自然死亡率μ。所有这些模型参数都是正常数。
许多实证研究[4]，[16]，[31]发现，由于自然免疫的作用，再感染率通常低于同一病原体的初次感染率（特别是SARS-CoV-2），因此我们在模型中假设（j = 1, 2）。
我们还强调，为了简化，我们的模型没有考虑免疫力减弱的影响。因此，免疫力主要被视为一种静态分层，即从感染中康复的个体将对同一病原体菌株具有永久免疫力，但他们仍有可能被另一种菌株感染。

3. 数学结果
我们数学模型的主要动态特性可以通过以下定理总结，所有证明都提供在附录A中。定理3.1涉及无病平衡状态（DFE）的存在性和基本再生数的推导。定理3.2解决了边界平衡状态的存在问题。定理3.3和定理3.4分别分析了DFE和边界平衡状态的稳定性。此外，定理3.5建立了内部平衡状态的存在性和唯一性的结果。

定理3.1
系统(2.1)在正不变域内具有唯一的无病平衡状态（DFE）(3.1)
模型的基本再生数由下式给出：
(3.3)
其中 和 分别是菌株1和菌株2的个体再生数。

定理3.2
当 时，DFE x0是系统(2.1)在Ω中的唯一平衡状态。当 时，除了DFE之外还存在边界平衡状态。具体来说，如果 ，系统(2.1)具有菌株-1边界平衡状态(3.4)；如果 ，系统(2.1)具有菌株-2边界平衡状态(3.5)

定理3.3
当 时，DFE x0在Ω中全局渐近稳定。

定理3.4
当 和 时，菌株-1边界平衡状态 在Ω中局部渐近稳定。当 和 时，菌株-2边界平衡状态 在Ω中局部渐近稳定。
系统和 的动态尚未完全解决。下面的定理提出了一个充分条件，以确保存在唯一的正内部平衡状态。从生物学角度来看，这个正平衡状态代表一种稳态，在这种状态下两种菌株共存，竞争排斥原理不适用。

定理3.5
设 和 。在以下条件下，
(3.6)
其中a = μ + γ1 + w1且b = μ + γ2 + w2，系统(2.1)具有唯一的正内部平衡状态。
我们对该条件(3.6)的流行病学含义进行了一些评论。该条件的左侧表示跨菌株传播（包括初次感染和再感染）的综合效应，按两个康复率的乘积进行了标准化。右侧的第一项表征了包括初次感染和再感染在内的菌株-1传播效应，按其康复率和菌株1导致的疾病暴发风险进行了标准化。类似的解释也适用于右侧的第二项。需要注意的是，在实际情况下，康复率通常远大于自然死亡率和疾病引起的死亡率，因此a ≈ γ1且b ≈ γ2。总体而言，条件(3.6)表明，当标准化的跨菌株传播强度小于两种标准化菌株特异性传播强度之和时，系统(2.1)具有唯一的正平衡状态，表示两种菌株共存。另一方面，如果跨菌株传播效应足够强，以至于不满足条件(3.6)，优势菌株可能会完全排除并取代另一种菌株，从而消除共存状态。
我们在附录A中提供了定理3.5的证明以及两种特殊情景的更详细分析。接下来，我们将使用数值模拟来说明这些数学结果，并通过一个现实案例研究来深入了解该模型的复杂动态。

4. 模拟结果
4.1. 平衡状态和稳定性
我们首先通过数值方法验证了第3节中提出的关于各种平衡状态及其稳定性预测。为此，我们对系统(2.1)进行了多次不同初始条件的数值模拟，每个初始条件确定了一个代表相图中轨道的解。然后我们绘制这些解轨迹随时间的变化，以检查它们的行为，从而确定相关平衡状态的稳定性。
图1展示了四种典型情景下I1与I2的相图：(a) ，所有轨道都趋向于DFE，其中I1 = I2 = 0；(b) ，解收敛于菌株-1边界平衡状态，其中I1 > 0且I2 = 0；(c) ，解收敛于菌株-2边界平衡状态，其中I1 = 0且I2 > 0；(d) ，存在一个内部平衡状态，其中I1 > 0且I2 > 0，并且至少吸引了一些轨道。请注意，这些图中没有显示所有与给定初始条件相关的解轨迹。图1(a)中的观察结果证实了定理3.3的结论。同时，图1(b)和图1(c)不仅证明了定理3.4，还表明在这些条件下这些边界平衡状态确实是全局渐近稳定的。此外，图1(d)中的参数设置满足条件(3.6)，因此系统(2.1)具有唯一的正内部平衡状态，这与定理3.5一致。

下载：下载高分辨率图像（717KB）
下载：下载全尺寸图像
图1. 四种不同情景下系统(2.1)的I1与I2的相图。在每个面板中，黑点代表关注的平衡点，红点代表初始点。4.2. 应用于COVID数据在过去的几年中，几种SARS-CoV-2变体导致了大规模的感染浪潮。作为我们研究的实际应用，我们将把我们的模型（2.1）应用于SARS-CoV-2变体的数据。我们考虑了2021年末到2022年初期间的美国数据，这段时间标志着Delta毒株的主导地位向Omicron毒株的出现和主导地位的过渡。具体来说，我们关注美国田纳西州的COVID数据。我们使用了两个主要的数据集来分析2021年11月13日至2022年4月16日这段时间的数据：（1）CDC提供的变体数据[45]，其中包含了田纳西州每周的变体比例报告；（2）田纳西州卫生部门提供的COVID病例数据[46]，其中包含了田纳西州每天新增病例和活跃病例的数量报告。在我们的模型中，我们将Omicron毒株标记为毒株1，它包括所有其亚变体，如在变体数据中报告的B.1.1.529、BA.1.1、BA.2和BA.5。模型中的毒株2主要包含Delta（B.1.617.2）。此外，在这段时间内还有其他一些毒株，但它们的比例非常低。为了简化，我们将所有这些变体与Delta合并为毒株2。结合上述两个数据集，我们能够确定与毒株1和毒株2相关的活跃病例数量，这将用于拟合我们模型中的变量I1和I2。对于参数校准，我们将使用已发布的来源来确定自然死亡率、人口流入率、恢复率和疾病引起的死亡率的值（详见第4.4节）。然而，传播率β1、β2的值在文献中并不可得，因为它们通常会因地点、人口、具体毒株和时间段而异。因此，我们将使用数据拟合来估计这四个传播率。在展示拟合结果之前，我们首先对模型参数进行了可识别性分析。4.3. 可识别性分析我们采用尺度不变局部结构可识别性方法[7]来分析模型参数的可识别性。在我们的模型和数据设置中，需要识别的四个传播率是β1、β2、...（此处数据缺失），而其他所有参数以及两个变量I1和I2都是可以观察到的。根据这种方法，我们将未观察到的变量S、R1、R2以及参数β1、β2、...（此处数据缺失）重新缩放如下（4.1），其中uS、...（此处数据缺失）是需要确定的缩放因子。对于系统（2.1）中的每个变量X（X = S、I1、I2、R1和R2），让fX表示相应方程的右侧。我们将每个fX表示为两个线性独立函数fX1和fX2的总和，如下（4.2）。系统（2.1）的不变性方程变为（4.3）系统4.3第一行的第二个方程立即得出uS = 1。同时，第四行的两个方程得出...（此处数据缺失），最后两行的两个方程得出...（此处数据缺失）。将这些结果代入第二行和第三行后得到...（此处数据缺失）。因此，系统4.3只有一个解。根据尺度不变局部结构可识别性方法[7]、[44]的理论，四个参数β1、β2、...（此处数据缺失）都是可识别的。4.4. 数据拟合我们现在继续通过将我们的模型（2.1）同时拟合到第4.2节描述的两个COVID数据集（变体比例和病例数据）来进行模型校准和参数估计。先前的实证和建模研究表明，不同的SARS-CoV-2变体，特别是Delta和Omicron，主要在初次感染和再感染的传播率上存在差异[4]、[6]、[16]、[27]。因此，我们将重点关注四个传播率β1、β2、...（此处数据缺失），并通过数据拟合来估计它们的值。为了简化，我们假设在这项研究中，毒株1和毒株2的恢复率相同，即γ1 = γ2 = 1/14每天，基于COVID-19的平均恢复期为2周[48]，以及疾病引起的死亡率相同，即ω1 = ω2 = 7.5×10^-4每天，基于美国三个州的平均估计[4]。2021年美国的平均预期寿命约为76岁[1]。因此，我们在模型中设置自然出生率和死亡率为μ = 1/(76年) = 1/(76×365)每天。此外，我们将人口流入率设置为Λ = μ×N，其中N = 7.05×10^6是2021-2022年田纳西州的人口规模[47]。所有这些参数的值，包括四个传播率的拟合估计值和置信区间，列在表1中。表1. 模型（2.1）的参数值参数描述值95%置信区间单位来源μ自然死亡率3.605×10^-5每天[1]Λ人口流入率254.1525人/天^-1[47]γ1毒株1的恢复率0.071428每天[48]γ2毒株2的恢复率0.071428每天[48]w1毒株1的死亡率7.5×10^-4每天[4]w2毒株2的死亡率7.5×10^-4每天[4]β1毒株1的初次感染率3.6607×10^-8[3.5329e-08, 3.7885e-08]（人/天）^-1拟合值β2毒株1的再感染率1.7069×10^-8[1.4727e-08, 1.9411e-08]（人/天）^-1拟合值β2毒株2的初次感染率1.1693×10^-8[1.0267e-08, 1.3118e-08]（人/天）^-1拟合值β2毒株2的再感染率1.0273×10^-8[8.8602e-09, 1.1686e-08]（人/天）^-1拟合值从拟合值可以看出，...（此处数据缺失）与我们的假设一致，即每个初次感染率都高于相应的再感染率。此外，β1 > β2，表明在这段时间内，毒株1（Omicron）的传播性和侵袭性更强。另外，...（此处数据缺失）的值也高于β2和...（此处数据缺失）。这表明Omicron的再感染在塑造COVID疫情中可能发挥了重要作用，而Delta的再感染作用则相对较小。这些结果与先前的实证发现[4]、[16]一致。图2显示了2021年11月13日至2022年4月16日活跃病例数量和累计病例数量的拟合曲线。活跃病例和累计病例的模拟曲线与报告的数据非常吻合。活跃病例数量的峰值（约1.9×10^5）出现在2022年1月。之后，活跃病例数量减少，在时间期末降至非常低的水平，而累计病例数量稳定在略高于2×10^6的水平。下载：下载高分辨率图像（295KB）下载：下载全尺寸图像图2. 2021年11月13日至2022年4月16日活跃病例数量（左）和累计病例数量（右）的拟合结果。蓝色圆圈代表报告数据，红色实线代表数值结果。使用拟合后的参数值评估方程（3.3），我们得到了毒株1和毒株2的再生数...（此处数据缺失）。这清楚地表明...（此处数据缺失），量化了COVID疫情的持续性。同时，...（此处数据缺失）再次表明，在我们关注的时期，Omicron的作用比Delta（和其他变体）更为重要。基本再生数可用于估计群体免疫阈值（HIT），即为了控制疾病所需具备的免疫宿主人口的比例。这个阈值可以通过一个简单的公式[15]来评估：...（此处数据缺失）。根据我们在田纳西州的数据拟合结果，我们得到HIT ≈ 72%。截至2021年11月，田纳西州大约55%的人口至少接种了一剂COVID疫苗，约49%的人口完全接种了疫苗[46]。这些比率远低于估计的72%的HIT。这一事实，加上Omicron的强大免疫逃逸能力，可能解释了这一时期的感染浪潮。然而，鉴于存在具有不同感染性的多种毒株以及COVID疫苗接种和感染后提供的保护期限有限，实际的人群免疫水平可能难以评估。使用拟合后的参数值，我们可以很容易地验证在这种情况下确实满足共存条件（3.6），表明存在一个唯一的正平衡。数值上，我们发现这个平衡大约位于...（此处数据缺失）。注意，在这个平衡状态下，毒株1的活跃病例数量是毒株2的10倍以上，毒株1的恢复病例数量是毒株2的20倍以上。接下来，我们对条件（3.6）进行更深入的数值研究，并探索其在参数空间中的变化。为此，我们引入了表达式（4.4）显然，条件（3.6）成立当且仅当D < 0。然后我们使用从数据中拟合得到的参数值来进行以下数值实验来表征D。1.我们将D象征性地视为...（此处数据缺失）和...（此处数据缺失）的函数，同时保持所有其他参数不变。当...（此处数据缺失）和...（此处数据缺失）在1到+∞之间变化时，D的值和符号也会改变。图3(a)在第一个象限绘制了曲线D = 0，其中...（此处数据缺失）和...（此处数据缺失）分别是横轴和纵轴。曲线以下的区域是D < 0，表明如果两个个体再生数中的任何一个足够高，则可能由于较弱毒株的完全排除和替代而不满足共存条件。下载：下载高分辨率图像（121KB）下载：下载全尺寸图像图3. 对于（4.4）中定义的量D的数值模拟：(a) D作为再生数...（此处数据缺失）和...（此处数据缺失）的函数；(b) D作为比率...（此处数据缺失）和...（此处数据缺失）的函数；(c) D作为比率...（此处数据缺失）和...（此处数据缺失）的函数。每个面板中的蓝色曲线代表D = 0.2。我们将D象征性地视为比率...（此处数据缺失）和...（此处数据缺失）的函数，同时保持所有其他参数不变。这两个比率比较了毒株1和2在初次感染和恢复方面的优势。当这些比率在0到+∞之间变化时，D的值和符号也会改变。我们再次在第一象限绘制曲线D = 0，其中两个比率分别代表横轴和纵轴。图3(b)显示了结果，其中D < 0的区域位于曲线D = 0之上。显然，对于每个感染比率值，恢复比率...（此处数据缺失）需要高于某个阈值才能确保两种毒株的共存。否则，非常低的比率...（此处数据缺失）只会促进毒株1的竞争力，导致毒株2的完全排除。3.同样地，我们还可以将D视为比率...（此处数据缺失）和...（此处数据缺失）的函数，比较毒株1和2在再感染和恢复方面的优势。结果在图3(c)中呈现，解释与情景2类似。我们注意到，上述每一个数值实验都是以某种简化的方式进行的，其中各种量（如再生数和各种比率）和参数实际上彼此依赖。尽管如此，通过关注每个情景中的特定量，我们获得了关于参数空间中条件（3.6）的有用见解。本质上，每个情景中的曲线D = 0区分了共存状态（D < 0）和竞争排斥状态（D > 0），代表了两种毒株在各种特征方面的平衡：对疾病爆发风险的贡献（情景1）、初次感染和恢复的综合效应（情景2），以及再感染和恢复的综合效应（情景3）。此外，我们还可以数值上检验正内平衡的稳定性。系统（2.1）在此平衡下的雅可比矩阵由...（此处数据缺失）给出。通过代入从数据中拟合得到的参数值和通过数值模拟找到的平衡，我们得到一个常数矩阵，其五个特征值分别为?6.6848×10^-5 ± 2.6065×10^-3i、?4.2874×10^-5 ± 4.3653×10^-4i和?3.6840×10^-5。所有这些特征值的实部都是负的，表明正平衡至少是局部渐近稳定的。4.5. 敏感性分析我们对系统（2.1）中的参数进行了敏感性分析，以量化它们对我们模型结果的影响。我们采用了[12]中提出的方法来计算每个参数相对于状态变量的相对敏感性。我们将系统（2.1）写成向量形式（4.5），其中表示未知向量，α是包含所有模型参数的向量，F是表示系统（2.1）右侧的向量值函数。让Xi和Fi分别表示X和F的第i个分量，αj表示α的第j个分量。对系统（4.5）中的第i个方程关于αj求偏导数，我们得到（4.6）。结合方程（4.5）和（4.6），我们可以数值求解每个敏感性函数。然后，参数αj相对于状态变量Xi的相对敏感性被定义为整个时间段内标准化敏感性函数的最大值：（4.7）这种敏感性分析的结果绘制在图4中，其中正值和负值分别表示状态变量和参数之间的正相关和负相关。我们清楚地看到，β1（菌株1的初次感染率）是对所有状态变量最敏感的参数。参数β2（菌株2的初次感染率）对与菌株2相关的感染和康复病例敏感。同时，菌株1的再感染率显示出显著的敏感性，而菌株2的再感染率的敏感性在所有变量中保持非常低。这些结果表明，除了之外，其他三个拟合参数β1、β2和都可能以显著的方式影响模型结果（详见第4.7节的进一步讨论）。下载：下载高分辨率图像（556KB）下载：下载完整尺寸图像图4. 基于方程（4.7）的系统（2.1）的敏感性分析结果。4.6. 菌株比较为了量化菌株1和菌株2的贡献，以及初次感染和再感染的贡献，我们现在比较了我们模型中这些不同组分的发病率和攻击率。图5显示了发病率比较的结果。这里的发病率是指在所关注的时间段内，由特定菌株和/或感染模式引起的新的病例数量。图5(a)比较了我们模型中由β1SI1表示的菌株1的初次感染发病率，以及由表示的菌株1的再感染发病率。初次感染发病率从几乎0开始，由于活跃感染（I1）的增加而不断上升，在时间间隔的中期达到峰值。之后，由于易感个体和活跃病例的减少（参见图2），它开始下降，在接近结束时降至非常低的水平。再感染发病率落后于初次感染发病率，并且大部分时间都低于它。然而，在达到峰值后，再感染发病率超过了初次感染发病率，并在其后一直保持较高水平。这一结果与先前的实证观察相符，即Omicron菌株的再感染显著增加[4]，[16]。下载：下载高分辨率图像（386KB）下载：下载完整尺寸图像图5. 发病率比较：(a) 菌株1的初次感染发病率β1SI1（绿色）与再感染发病率（红色）；(b) 菌株2的初次感染发病率β2SI2（绿色）与再感染发病率（红色）；(c) 菌株1的综合发病率（绿色）与菌株2的综合发病率（红色）；(d) 菌株1和菌株2的综合初次感染发病率（绿色）与综合再感染发病率（红色）。图5(b)比较了由β2SI2表示的菌株2的初次感染发病率，以及由表示的菌株2的再感染发病率。初次感染发病率从峰值开始下降，直到降至非常低的水平，这反映了菌株2（即Delta）在整个时间段内的下降。图5(c)比较了由表示的菌株1的初次感染和再感染的综合发病率，以及由表示的菌株2的初次感染和再感染的综合发病率。菌株2的发病率在第一个月占主导地位，但在整个时间段内持续下降。菌株1的发病率虽然开始时落后，但迅速达到高峰并在最后几个月保持主导地位。这些结果再次证明了从Delta到Omicron的主导地位转变。图5(d)比较了由β1SI1 + β2SI2表示的两种菌株的综合初次感染发病率，以及由表示的两种菌株的综合再感染发病率。模式与图5(a)中的类似。接下来，我们在图6中绘制了每1000人的攻击率（4.8）。在我们的模型设定中，易感个体面临两种菌株的初次感染风险，而从一种菌株康复的个体面临另一种菌株的再感染风险。下载：下载高分辨率图像（435KB）下载：下载完整尺寸图像图6. 每1000人的攻击率比较：(a) 菌株1的初次感染攻击率β1I1?1000（绿色）与再感染攻击率（红色）；(b) 菌株2的初次感染攻击率β2I2?1000（绿色）与再感染攻击率（红色）；(c) 菌株1的综合攻击率（绿色）与菌株2的综合攻击率（红色）；(d) 菌株1和菌株2的综合初次感染攻击率（β1I1 + β2I2）?1000（绿色）与综合再感染攻击率（红色）。图6(a)比较了由（β1SI1/S）?1000 = β1I1?1000表示的菌株1的初次感染攻击率，以及由表示的菌株1的再感染攻击率。初次感染攻击率的峰值约为每1000人72次。同时，再感染攻击率的峰值约为每1000人33次，显示出菌株1再感染的显著威胁。与图5(a)相比，攻击率的峰值出现在菌株1初次感染发病率峰值之后，而攻击率和发病率的峰值几乎同时出现在菌株1再感染期间。图6(b)比较了由β2I2?1000表示的菌株2的初次感染攻击率，以及由表示的菌株2的再感染攻击率。两种攻击率的峰值分别出现在大约每1000人3次和2.6次，远低于菌株1的对应值。图6(c)比较了由表示的菌株1的初次感染和再感染的综合攻击率，以及由表示的菌株2的初次感染和再感染的综合攻击率。对于每种菌株，综合攻击率可以解释为初次感染和再感染攻击率之间的一种平均值。菌株1的综合攻击率峰值约为每1000人53次，虽然落后于菌株2的综合攻击率峰值（约为每1000人4次），但要高得多。这种行为再次显示了从Delta到Omicron的主导地位转变。图6(d)比较了由（β1I1 + β2I2）?1000表示的两种菌株的综合初次感染攻击率，以及由表示的两种菌株的综合再感染攻击率。可以看出，菌株1在初次感染和再感染攻击率中都占主导地位。结合图5和图6，我们观察到以下情况：(1) 对于每种菌株，再感染发病率的峰值都落后于初次感染发病率的峰值，而初次感染和再感染的攻击率峰值几乎同时出现。(2) 菌株2（Delta）的发病率和攻击率峰值分别出现在菌株1（Omicron）的发病率和攻击率峰值之前，并且低于它们，证实了菌株1在整个时间段内的出现和主导地位。(3) 再感染发病率和攻击率主要由菌株1贡献，它在塑造流行病波中发挥了重要作用。(4) 菌株1的发病率和攻击率峰值显著高于菌株2的对应值。4.7. 参数变化我们从数据拟合得到的四个传播率β1、β2、和的估计值被视为这些参数的基线值（见表1）。第4.5节中的敏感性分析表明，这些参数中的大多数（菌株2的再感染率除外）都是敏感的。为了进一步展示各种传播率对模型结果的影响，我们对β1、β2、和进行了参数变化测试。对于每次测试，我们将一个传播率分别扰动?50%、?25%、+25%和+50%，同时保持所有其他模型参数在其基线值。然后我们运行模型，并为2021年11月13日至2022年4月16日的相同时间间隔内的每个状态变量生成假设模拟曲线。图7、图8、图9和图10展示了这些结果（见附录B）。从图7中我们观察到菌株1的初次感染率β1与直接相关的两个状态变量；即活跃病例I1和康复病例R1之间存在强烈的正相关。β1的增加或减少会导致I1和R1同向的显著变化。相比之下，直接与菌株2相关的两个状态变量；即活跃病例I2和康复病例R2似乎与β1呈负相关。当β1降低时，由于菌株1的贡献减弱（因此菌株2的贡献增强），I2和R2增加，从而影响宿主群体内的整体感染动态。类似地，图8显示菌株2的初次感染率β2与两个菌株2变量I2和R2之间存在强烈的正相关，而β2与两个菌株1变量I1和R1之间存在负相关。图9中菌株1的再感染率的变化在定量变化程度上与图7中的β1相似，尽管程度较低。最后，图10中菌株2的再感染率的变化结果加强了我们之前的观察，即表现出非常低的敏感性，对模型结果的影响微乎其微。5. 讨论通过突变和进化，许多病原体会不断产生新的菌株，这些菌株可能对公共卫生构成威胁。在本文中，我们提出了一个新的数学模型来研究同时循环的两种菌株的传染病群体动态。这些菌株在宿主群体内的传播和扩散能力方面相互竞争。我们的模型包含了不同菌株之间的相互作用、再感染的影响以及多种传播模式的存在。我们分析了基于微分方程和动力系统理论的该模型的基本动态行为，包括各种类型平衡解的存在及其稳定性属性。我们特别关注导致两种菌株共存的独特正平衡条件，并对此共存条件进行了详细的数值研究。在我们的模拟研究中，我们使用真实数据校准了模型并进行了一系列数值模拟。我们的工作强调了多菌株病原体传播所涉及的复杂动态的全面理解。利用公开可用的数据，我们将模型应用于2021年11月至2022年4月美国田纳西州的COVID疫情。这一时期出现了两种主要的SARS-CoV-2变体，其中Delta菌株正在传播，而Omicron菌株正在出现。我们仔细比较和量化了每种菌株的贡献以及每种传播途径的重要性。我们的数据拟合和数值模拟清楚地表明，与德尔塔变异株相比，奥密克戎变异株的传染性更强、毒性更大，发病率和攻击率显著更高，再次感染的风险也大幅增加。同时，我们发现奥密克戎的初次感染和再次感染在这一时期对疫情波动的形成起到了重要作用。另一方面，德尔塔变异株的再次感染的影响似乎可以忽略不计，尽管其初次感染途径仍对整体疾病流行率有相当大的贡献。这些结果与先前的实证研究发现一致。一种新出现的变异株（尤其是奥密克戎）通常由于其较高的传染率、更强的感染性（体现在初次感染和再次感染的发病率和攻击率上）以及更强的免疫逃逸能力而具有竞争优势。这项工作的研究结果强调了针对这类新兴病原菌的有效预防和干预措施的重要性。有策略地部署疫苗，尤其是针对当时最流行变异株的更新疫苗和加强针，将是控制疾病最有效的措施之一。此外，非药物干预措施如卫生和清洁、社交距离、教育活动及提高意识的计划对于疾病暴发的管理也非常有用。所有这些方法都可以有效降低循环病原菌的传播率，尤其是新兴变异株的传播率。如图7和图9所示的假设模拟结果所展示的，如果将Omicron的初次感染率和再次感染率降低25%，疾病流行率将分别减少约33%和16%。此外，这些数据还表明，降低传播率可以显著推迟感染高峰的出现，为管理和控制新兴变异株争取更多时间。

本研究存在一些局限性。我们的模型没有设定潜伏期/暴露期状态来表示已感染但尚未具有传染性的个体，也没有区分无症状和有症状的传染者。模型没有考虑免疫力减弱的影响；相反，它假设从感染中恢复后会获得对该变异株的永久性免疫，这可能并不现实。这样设计的目的是为了使模型（2.1）尽可能简洁，以专注于本研究的主要关注点——不同变异株之间的竞争以及再次感染的影响。此外，模型中没有明确考虑接种疫苗的个体；相反，疫苗的影响是通过从真实数据中拟合得到的初次感染率和再次感染率来间接体现的。未来的一个有意义的方向是扩展当前模型，明确纳入感染潜伏期、无症状感染、免疫力减弱和疫苗接种等额外因素。这些因素的加入可能会改变模型的动态特性，特别是不同变异株的共存情况。此外，还可以将这个建模框架扩展到能够表示三种或更多同时传播的变异株之间的相互作用。

为了简化问题，我们假设所有模型参数都是常数，但实际上这些参数，尤其是传播率，很可能会随着疫情的发展和不同变异株之间的相互作用而发生变化。引入时间依赖参数可以使模型更加真实，但也会增加模型校准和数据拟合的难度。将机理建模与机器学习相结合可能会帮助克服这些挑战。

**作者贡献声明：**
白杰：撰写原始稿件、可视化展示、验证、软件开发、形式分析、数据整理。
王进：审稿与编辑、监督、方法论设计、数据分析、概念构建。

**未引用参考文献：[11]**