基于精确贝叶斯推断的普通燕鸥高致病性禽流感流行病学参数估计

本文采用可逆跳跃马尔可夫链蒙特卡洛(RJ MCMC)方法，对德国Banter See普通燕鸥 colony在2022-2023年高致病性禽流感(HPAI)爆发期间的每日死亡数据建立随机个体传播模型，首次实现了野生鸟类种群中基本繁殖数(R0)和死亡率(pd)等关键参数的精确贝叶斯推断，为海鸟保护提供了重要流行病学依据。

　　

1 引言

高致病性禽流感(HPAI)是一种传染性病毒性疾病，已导致大量家禽被扑杀以及无数野生鸟类死亡。2021年开始的最新全球大流行因其对全球海鸟种群的毁灭性影响而尤为引人注目。虽然HPAI在家禽中的传播已得到相对充分的研究，但当前大流行的新近性加上数据收集的挑战，意味着对野生种群中关键流行病学参数（如基本繁殖数R₀）的了解要少得多。

数学和统计建模方法越来越多地用于支持控制和预防禽流感的有效决策。先进的统计推断方法为我们提供了估计关键流行病学参数的工具，例如R₀值，即使在数据稀疏的情况下也能告诉我们，在时间t将单个感染个体引入易感种群中预期会产生多少二次感染。这些统计工具非常灵活，能够使建模假设适应当前的基础问题，并已应用于应对一系列疾病挑战。更具体地说，多种统计方法已被应用于估计家禽养殖场中HPAI的群内和群间传播参数，通常是为了参数化模拟模型以帮助决策者管理疫情。类似的估计野生鸟类传播的努力更为有限，可能是由于数据可用性和资金有限。

我们调整了先前应用于估计非洲猪瘟病毒在猪群内传播参数的精确贝叶斯推断方法，并将其应用于估计德国Banter See一个经过充分研究的普通燕鸥 colony在2021年和2023年连续爆发的HPAI传播情况。我们的模型假设 colony内的疾病传播遵循随机个体基础的流行病模型。然后，我们使用数据增强的可逆跳跃马尔可夫链蒙特卡洛从每日死亡记录中推断关键流行病学参数。最后，我们强调了该方法超出此处给出的具体示例的广泛适用性，并结合传播模式和免疫性来解释结果，这对于预测对种群的长期影响非常重要。

2 材料与方法

2.1 Banter See 普通燕鸥 HPAI 爆发数据

数据收集于2022年和2023年高致病性禽流感爆发期间，地点在德国北部Banter See的一个普通燕鸥繁殖 colony。该 colony自1992年以来，对成年燕鸥及其后代进行了系统的环志和应答器标记，通过天线系统远程记录个体存在和繁殖性能。在2022年和2023年的HPAI爆发期间，对 colony及周边湖泊进行了每日监测，所有发现的死亡鸟类均从环境中移除并冷冻保存。 colony在HPAI致死前的总个体数估计2022年为2071只，2023年为1061只。两年内记录的累计死亡个体数分别为511只和111只。考虑到可能存在死亡个体未被完全检测到的情况，推断方法中对此进行了考虑。

2.2 模型构建与假设

我们使用了Ewing等人提出的推断框架，该框架旨在从疫情早期监测期间的每日死亡动物记录数来推断流行病学参数。与养殖场不同，海鸟 colony是一个基本未受管理的环境，这意味着在疾病诊断后扑杀动物所导致的观察窗口截断在这里不适用，我们可以观察到爆发的整个过程。然而，对野生种群的观察也有其自身的挑战：首先，我们无法控制不同海鸟物种之间潜在的混合，而且我们几乎不可能完全检测到所有死亡个体。此外，与某些可以合理假定在家猪中致死率为100%的ASFV毒株不同，一部分海鸟会康复并对禽流感产生免疫力。这种康复意味着，特别是在2023年，初始种群可能包含未知比例的易感和免疫个体。

我们专注于模拟禽流感在 colony内的传播，并简化假设与其他物种的混合对这种传播没有影响。虽然这显然是该方法应用于其他野生鸟类物种或存在大量混合的其他 colony的一个局限性，但对于像本研究这样不与其他 colony紧密相邻繁殖的种群来说，这可能是一个合理的假设。

通过扩展易感-潜伏-感染-移除模型框架来解决死亡个体检测不完善和包含康复免疫个体的问题，模型包括以下类别：历史免疫、易感、潜伏、感染、康复、已观察死亡和未观察死亡。初始种群以概率p_h和1-p_h分为历史免疫个体和易感个体。历史免疫个体被认为因在越冬地或前一年经历过疫情而具有免疫力。一旦易感个体被感染，它们进入潜伏状态，此时已感染但不具传染性。经过这段潜伏期后，个体进入传染期，之后要么死亡，要么康复。那些死亡的个体随后以概率1-p_u被观察到，或以概率p_u未被观察到。从先前疫情获得免疫力的历史免疫个体和从当前疫情获得免疫力的康复个体被分开建模。这使我们能够推断特定年份内的疫情参数，而不需要明确模拟年际种群动态。

我们专注于估计R₀，不过在假设种群中部分个体具有预先免疫力的案例中，更有意义的指标或许是有效繁殖数R_e，其计算公式为R_e = R₀(1 - p_h)。这考虑了由于易感个体比例减少而导致的后续传播减弱。在假设存在预先免疫力的案例中，我们同时给出R₀和R_e值；然而，我们专注于比较年份间和模型间的R₀值，既是为了确保同类比较，也是由于历史免疫水平的不确定性。

假设 colony内的疾病传播遵循随机个体基础的流行病模型。我们假设存在一个封闭且均匀混合的种群，其规模为N(t)，其中在时间t有I(t)个个体具有传染性。在每年疫情开始时，我们假设N(t)是该年估计的总种群大小，随后N(t)随着因HPAI导致的观察到的和推断出的死亡率而随时间减少。这假设病毒不是由早期抵达者引入的，并且忽略了个体在繁殖后或应对感染时离开 colony的速率。感染的发生要么通过直接的、频率依赖的 colony内传播β₂，要么通过来自环境储库的传播β₁，因此作用于易感个体的感染力由λ(t) = (β₂I(t)/N(t) + β₁)给出。因此，新感染的发生遵循一个非齐次泊松过程，其速率为λ(t)。通过纳入外部感染力，我们允许存在多个个体通过 colony内传播以外的过程被感染的可能性，从而隐含地允许多次传入事件进入 colony的可能性。这与病原体序列数据一致，表明在该 colony内存在多个起源点。我们还通过考虑去掉分母中N(t)项的情况，来检验我们选择频率依赖传播（相对于密度依赖传播）的影响。为简单起见，我们未尝试明确模拟环境感染储库，而是假设其引起恒定的背景感染力。我们假设潜伏期和传染期的持续时间服从伽马分布。