提示:因为课题需要,对以下文章进行总结提炼,方便对文章的理解。
论文标题:Estimating epidemiologic dynamics from cross-sectional viral load distributions(从横断面病毒载量分布估计流行病学动态)Science 2021
DOI: 10.1126/science.abh0635
《Estimating epidemiologic dynamics from cross-sectional viral load distributions》论文提炼
摘要
提示:这里可以添加本文要记录的大概内容:
估计流行病的轨迹对于制定传染病公共卫生对策至关重要,但用于此类估计的病例数据受到各种检测方法的干扰。我们表明,在随机或基于症状的监测下观察到的病毒载量的人群分布,以逆转录定量聚合酶链反应(RT-PCR)检测获得的循环阈值(Ct)值的形式,在流行期间发生变化。
因此,即使是有限数量的随机样本的Ct值也可以提供对流行病轨迹的更好估计。组合来自多个这样的样本的数据提高了这样的估计的精度和鲁棒性。
我们将我们的方法应用于SARS-CoV-2大流行期间在各种环境下进行的监测的Ct值,并为疫情管理和应对提供了实时估计流行轨迹的替代方法。
介绍
实时跟踪流行轨迹和感染发生率是大流行期间公共卫生规划和干预的基础。在SARS-CoV-2大流行中,主要流行病学参数(如Rt)通常使用观察到的病例计数、住院或死亡的时间序列进行估计,通常基于逆转录定量聚合酶链反应(RT-qPCR)检测。然而,有限的检测能力、检测可用性随时间的变化以及报告延迟都对结果产生一定影响。
RT-qPCR检测以周期阈值(Ct)的形式提供半定量结果,其与log10病毒载量呈负相关,但它们通常仅作为二进制“阳性”或“阴性”报告。对于SARS-CoV-2,Ct值可用于临床确定是否需要隔离、识别个体感染的阶段和预测疾病严重程度。然而,基于Ct值的个体层面决策尚未广泛实施,因为测试平台和样本之间的测量差异以及对无症状和症状前感染中SARS-CoV-2病毒动力学的了解有限。尽管单个高Ct值可能不能保证一个样本中的病毒载量低,例如,因为样本采集可变,但在许多样本中测量高Ct数值将表明病毒载量主要较低的人群。因此,Ct值的横截面分布应代表随着时间推移潜在人群中的病毒载量,这可能与流行轨迹的变化一致。例如,量化Ct值分布的系统性增加伴随着流行病的下降。
在这里,我们证明了来自RT-qPCR数据的单个或连续横截面样本的Ct值可以用于估计流行轨迹,而不需要来自测试阳性率或连续病例计数的额外信息。我们证明,观察到的Ct值分布中的人群水平变化可能作为一种流行病学现象出现,这对根据新出现的SARS-CoV-2变体解释RT-qPCR数据具有一定的意义。此外,我们还展示了如何将多个横断面样本结合起来,以提高对人口发病率的估计,如果没有血清学监测研究,这一指标往往难以实现。总体而言,我们使用收集到但目前通常被丢弃的Ct数据提供实时监测疫情的指标,我们的方法推动了旨在监测疫情的测试程序的开发。
Ct值与流行病动态的关系
模拟SEIR模型
首先,我们表明宿主内病毒动力学和流行病动力学的相互作用可以驱动Ct值随时间的变化,而不改变潜在病原体动力学。也就是说,Ct值分布的群体水平变化可以在个体水平上发生,而基本的感染后病毒载量轨迹没有系统变化。为了证明传播率和测量的病毒载量或Ct值之间的流行病学联系,我们首先模拟了SEIR模型下产生的感染(图1A(使用SEIR模型模拟的流行病中新增感染的人均日发病率(incidence)(直方图)和日增长率(蓝线)?
)和材料和方法:流行病传播模型),所用参数见表S1。在疫情期间的选定测试日,使用图S1和S2所示的材料和方法:Ct值模型中描述的Ct分布模型,从人群的随机横截面样本中观察模拟Ct值。通过在特定时间点从人群中抽取模拟样本进行测试,这些模拟再现了整个流行病过程中可检测病毒载量的真实横截面分布。在整个过程中,我们假设每个人最多感染一次,忽略了再感染,因为到目前为止,再感染似乎是疫情中微不足道的一部分。
SEIR传播模型:
模型具有易感(S)、暴露-非感染(E)、感染(I)和恢复(R)四个状态
对于Ct模型拟合,每日感染概率𝜋𝑡由累积暴露(E)状态的每日人均增加量给出:
Ct值模型
第a天模态(modal)Ct值:
随机选择的个体在其感染后a天的Ct值C(a)服从
其中
表S1:病毒载量和Ct值分布、SEIR传播模型、SEEIRR传播模型和高斯过程(GP)模型参数表
通过模拟得到感染率和群体Ct间关系
在疫情早期,感染率(infection incidence)迅速增长,因此大多数感染都是由最近的接触引起的。然而,随着疫情的减弱,新感染率的降低,受感染个体的平均暴露时间增加(图1,B(感染后中位天数与新增感染的日增长率。此处和(E-G)中的标记点显示了模拟流行病中的五个时间点
)和E(流行病期间随机选择的个体自感染以来的天数分布(小提琴图和直方图)。中位数以及第一四分位数和第三四分位数显示为绿色线和点。
))。这类似于增长人口与下降人口的平均年龄较低。虽然感染通常是未观察到的事件,但我们可以依赖一个可观察到的数量,例如病毒载量,作为感染后时间的一个指标。由于Ct值在感染宿主体内随时间不对称地变化(图1C(每天观察500例随机抽样感染个体的Ct值
)),病毒载量峰值出现在感染早期,因此在流行病生长期间对个体进行随机抽样更有可能对最近感染的个体进行抽样,因此病毒RNA的量更高。相反,在疫情下降期间随机抽样的感染个体更有可能处于感染的后期,通常抽样的病毒RNA量较低,尽管在所有时间点都存在大量抽样和病毒载量变异(图1D(病毒动力学模型(峰值后Ct值增加,随后达到检测限附近的平台期),显示Ct值(x轴,线显示平均值,条带显示95%分位数范围)相对于感染后天数(y轴)的时间进程。注意,y轴与(E)对齐。
))。因此,随机监测测试下观察到的Ct值的总体分布会随着时间的推移而变化,如中值、四分位数和偏度所示(图1G(按流行日列出的采样感染个体中观察到的Ct值分布(小提琴图和直方图)。中位数以及第一四分位数和第三四分位数显示为紫色线和点
))。虽然基于单个Ct值对个体感染后的时间估计将高度不确定,但观察到的Ct值的人口水平分布将随着新感染的增长率和Rt而变化(图1,F(观察到的Ct值分布偏度与感染每日增长率
)和H(根据观察到的Ct值分布的中位数和偏度绘制由SEIR模拟得出的Rt
))。
图1:Ct值分布反映了疫情暴发过程中的流行病学动态
为了在经典结果的背景下总结这一关键观察结果,我们发现快速增长的流行病学人群(Rt>1,增长率r>0)将以新感染为主,从而具有高病毒载量,而萎缩的流行病(Rt<1,r<0)将具有更多的老年感染,从而在给定的横截面上具有低病毒载量。并且r由生成间隔的分布来调制。相似的原则已应用于血清学数据,以推断未观察到的个体水平感染事件和传染病传播的人群水平参数。
我们发现,在基于症状的监测下获得的Ct值中,这种现象也可能存在,尽管不太明显,在症状出现后,对个体进行识别和测试。与随机监测测试的情况类似,测试最近有症状的个体获得的Ct值在疫情增长期间比在疫情下降期间更低(即,病毒载量更高)(图S3和S4)。然而,定义这种关系的确切性质和强度将取决于满足的许多条件(参见图S4标题)。
通过对个体水平的病毒生长/清除动力学和采样误差引起的观察到的Ct值的变化进行建模,随机样本中的观察到Ct值分布成为自感染以来的时间的可估计函数,然后可以从流行病增长率(epidemic growth rate)预测给定时间点的预期Ct值中值和偏差。然后,该函数可用于根据一组观察到的Ct值估计流行病增长率。如上所述,在大多数采样策略下,Ct值和流行病增长率之间存在关系,但校准精确映射是实现推断所必需的(例如,使用不同的RT-qPCR;见图S5)。这种映射可能会被检测偏差所混淆,例如,当检测能力有限时,感染和样本采集日期之间的延迟,或对病毒载量较高的样本(如来自严重疾病患者的样本)的系统性偏差。在这里,我们重点关注随机监测测试的情况,即在感染过程中的随机点对个体进行采样
用单个截面推断传染病的传播轨迹
拟合两个传播模型:
①指数增长模型
发病率(incidence)以r的速度指数增长(适合感染初期):
π
t
−
a
=
π
0
e
x
p
[
r
(
t
−
a
)
]
\pi_{t-a}=\pi_0exp[r(t-a)]
πt−a=π0exp[r(t−a)]
则似然函数为:
根据贝叶斯框架:
P
(
r
∣
X
1
,
.
.
.
,
X
n
)
=
P
(
X
1
,
.
.
.
,
X
n
∣
r
)
P
(
r
)
P
(
X
1
,
.
.
.
,
X
n
)
P(r|X_1,...,X_n)={P(X_1,...,X_n|r)P(r)\over P(X_1,...,X_n)}
P(r∣X1,...,Xn)=P(X1,...,Xn)P(X1,...,Xn∣r)P(r),其中P(r)服从均值为0,方差为0.25的正态分布,利用MCMC求得后验分布
P
(
r
∣
X
1
,
.
.
.
,
X
n
)
P(r|X_1,...,X_n)
P(r∣X1,...,Xn)
②SEIR模型
增长率(growth rate)根据剩余的易感个体数量每天变化,日感染概率
π
t
=
β
S
I
N
,
β
=
R
0
γ
\pi_t={\beta SI\over N},\beta={R_0\over \gamma}
πt=NβSI,β=γR0
从这些关系中,我们提出了一种方法,以推断给定随机抽样的RT-qPCR测试结果的单个截面的流行病增长率(epidemic growth rate)。该方法结合了两个模型:(1)观察到的Ct值的概率分布(以及阴性结果的概率),取决于感染和采样之间的天数;以及(2)在采样日期之前的给定日期被感染的可能性。对于(1),我们使用贝叶斯模型,并基于现有文献(材料和方法:Ct值模型和单截面模型)定义感染后Ct值的模式和范围的先验。对于(2),我们最初开发了两个模型来描述感染随时间的概率:(a)感染率(infection incidence)的持续指数增长;和(b)在SEIR模型下产生的感染。这两个模型都提供了流行病增长率的估计,但对爆发轨迹的可能形状作出了不同的假设:指数增长模型假设在前五周内的增长率是恒定的,几乎不需要事先的假设,而SEIR模型则假设增长率根据剩余的易感个体数量每天变化,但需要更多的先验信息。
Ct观测值与每日感染概率之间的关系之单截面模型:
对于一个测试日t,令 π t − A m a x , . . . , π t − 1 \pi_{t-A_{max}},...,\pi_{t-1} πt−Amax,...,πt−1分别为t天之前 A m a x A_{max} Amax天和前1天内整个人群日感染概率。其中 t − A m a x t-A_{max} t−Amax为检测日t时可检测出病毒的最早感染日。所以 π t − a \pi_{t-a} πt−a为是随机选择的个体在t − a日被感染的概率。假设数值可检测,则 p a ( x ) p_a(x) pa(x)为感染后第a天测试Ct值为x的概率(归一化到可观察值的Gumbel概率密度函数),即 p a ( x ) = P [ C ( a ) = x ] / P [ 0 ≤ C ( a ) < C L O D ] p_a(x)=P[C(a)=x]/P[0\leq C(a) < C_{LOD}] pa(x)=P[C(a)=x]/P[0≤C(a)<CLOD],其中 P [ C ( a ) = x ] P[C(a)=x] P[C(a)=x]为具有位置参数 C m o d e ( a ) C_{mode}(a) Cmode(a)和尺度参数 σ ( a ) \sigma(a) σ(a)的Gumbel概率密度函数。令 ϕ a \phi_a ϕa为感染a天后Ct值可检测的概率,取决于 C ( a ) C(a) C(a)和检测能力。
将n个个体样本的PCR检测结果记录为 X 1 , . . . , X n X_1,...,X_n X1,...,Xn,然后对于 x i < C L O D x_i < C_{LOD} xi<CLOD(可检测的Ct值),个体i具有Ct值 x i x_i xi的概率为: P ( X i = x i ∣ π t − A m a x , . . . , π t − 1 ) = ∑ a = 1 A m a x p a ( x i ) ϕ a π t − a P(X_i=x_i|\pi_{t-A_{max}},...,\pi_{t-1})=\sum_{a=1}^{A_{max}}p_a(x_i)\phi_a\pi_{t-a} P(Xi=xi∣πt−Amax,...,πt−1)=∑a=1Amaxpa(xi)ϕaπt−a.随机选择的个体在检测日t可被PCR检出的概率为: P ( X i < C L O D ∣ π t − A m a x , . . . , π t − 1 ) = ∑ a = 1 A m a x ϕ a π t − a P(X_i<C_{LOD}|\pi_{t-A_{max}},...,\pi_{t-1})=\sum_{a=1}^{A_{max}}\phi_a\pi_{t-a} P(Xi<CLOD∣πt−Amax,...,πt−1)=∑a=1Amaxϕaπt−a.所以n个PCR值的似然函数为: L ( X 1 , . . . , X n ∣ π t − A m a x , . . . , π t − 1 ) = ∏ i = 1 n [ ( ∑ a = 1 A m a x p a ( X i ) ϕ a π t − a ) I ( X i < C L O D ) ( 1 − ∑ a = 1 A m a x ϕ a π t − a ) I ( X i ≥ C L O D ) ] L(X_1,...,X_n|\pi_{t-A_{max}},...,\pi_{t-1})=\prod_{i=1}^{n}[(\sum_{a=1}^{A_{max}}p_a(X_i)\phi_a\pi_{t-a})^{I(X_i<C_{LOD})}(1-\sum_{a=1}^{A_{max}}\phi_a\pi_{t-a})^{I(X_i\geq C_{LOD})}] L(X1,...,Xn∣πt−Amax,...,πt−1)=∏i=1n[(∑a=1Amaxpa(Xi)ϕaπt−a)I(Xi<CLOD)(1−∑a=1Amaxϕaπt−a)I(Xi≥CLOD)]
如果仅将可检测的Ct值记录为 X 1 , . . . , X n X_1,...,X_n X1,...,Xn,则似然函数为
这些似然函数中的任何一个都可以被最大化,以获得每日感染概率的非参数估计,其约束条件为: ∑ a = 1 A m a x π t − a ≤ 1 \sum_{a=1}^{A_{max}}\pi_{t-a} \leq 1 ∑a=1Amaxπt−a≤1.然而,为了提高估计值的可解释性,我们考虑了基于上述流行病传播模型的两个参数模型:(i)在采样日之前的限定时期内感染发生率(incidence)呈指数增长的模型;和(ii)封闭有限群体中的SEIR模型,其中基本再生数R0是由模型估计的参数但不随时间变化(即,不存在降低可传播性的干预)。参见补充资料:拟合横断面病毒载量数据的参数模型.
以SEIR模型为例:
生成所有暴露、感染和康复个体的Ct值。SEIR模型参数的先验分布如表S1所示。对于任何一组参数值,通过发现新暴露个体的发病率(incidence)和使用上文给出的适当非参数似然方程(包括所有样本或仅具有可检测Ct值的样本)发现的似然性来确定每日感染概率 π t \pi_t πt。使用MCMC拟合获得SEIR模型参数的后验分布以及Ct动力学模型参数的后验分布。根据SEIR模型参数的后验估计值,还可以发现相关模型特征的后验分布,例如截至测试日t的Rt。在该模型下,参数 A m a x A_{max} Amax简单地取为模型周期的测试日t和第0天之间的天数。例如,在长期护理机构分析中,假设第0天为2020年2月1日,并估计了流行种子日期t0。
指数增长模型
假设在测试第t天之前的 A m a x A_{max} Amax天内,日发病率(incidence)以r的速度呈指数增长,因此 π t − a = π 0 e x p [ r ( t − a ) ] \pi_{t-a}=\pi_0exp[r(t-a)] πt−a=π0exp[r(t−a)]。因此,指数增长率β是 t − A m a x t-A_{max} t−Amax至t − 1天内发病率日增长率的对数.β值越大表示新感染的增长越快,而值为0表示每天新感染的数量没有增加或减少。对于正值𝑟, 新感染的倍增时间(以天为单位)由 l o g 2 r log2\over r rlog2给出。当易感个体的数量与感染数量相比很大时,该模型可能是疫情暴发初期的合理近似。然而,对于较长时间内的流行过程,恒定指数增长是一个糟糕的假设。因此,我们设置 A m a x A_{max} Amax, 并且等效地 t m a x t_{max} tmax 根据流行病传播模型一节中的命名法,选择35天,以在大多数个体在感染后保持PCR阳性的持续时间之间折中,而不假设在太长的时间内持续指数增长。
拟合该似然性,并忽略多余参数 π 0 \pi_0 π0,给出:
是否记录可检出和不可检出Ct值。如果仅记录可检测Ct值,则参数似然性由下式给出:
为了将不确定性纳入感染后每天病毒载量的分布中,我们构建了一个贝叶斯框架用于估计和推断。在本文中,我们使用r的均值为0且标准差为0.25的正态分布先验,使用马尔可夫链蒙特卡罗(MCMC)算法拟合以获得后验分布。注意,r = 0.1对应于新感染的倍增时间约为1周。Ct动力学模型参数的先验分布见表S1,如上所述。
该方法还产生Ct动力学模型的参数的后验分布。这些是估计流行轨迹的多余参数,但可能有助于改善该模型未来使用的先验。
使用SEEIRR模型拟合马萨诸塞州长期护理机构真实数据,提供近似ground truth
为了证明这种方法在封闭人群中具有单一横断面的潜力,我们首先调查了在2020年3月和4月经历SARS-CoV-2疫情的四个马萨诸塞州长期护理机构中,Ct值的分布和PCR阳性率随时间的变化。在每个机构中,我们都有疫情开始后三个时间点对居民和工作人员进行的近乎全员的PCR检测结果,包括阳性样本数量、阳性样本的Ct值和阴性样本数量(材料和方法:长期护理设施数据)。为了对我们基于Ct值的流行病轨迹估计进行基准化,我们首先使用标准房室模型来估计轨迹,该方法适合每个机构中随时间变化的测量点流行率(图2A(SEEIRR模型拟合3个采样时间点患病率的估计患病率(prevalence)(淡橙子线显示后验样本,橙色实线显示后验中位数,橙色丝带显示95%置信区间(CrI))和发生率(incidence)(红线显示后验中位数,红色丝带显示95% CrI)
))。具体而言,我们将一个简单的扩展SEIR(SEEIRR)模型与每个机构的三个观察点流行率(prevalence)值进行拟合,其中包括描述PCR阳性持续时间的额外暴露和恢复分区(材料和方法:流行病传播模型)。由于测试几乎是通用的,因此该方法提供了流行病轨迹的近似ground truth,我们可以根据该ground truth评估基于Ct值的方法的准确性。我们称之为基线估计。图2显示了其中一个长期护理机构的结果和数据,而图S6和S7显示了其他三个机构的结果。
马萨诸塞州长期护理机构的数据是从麻省理工学院布罗德研究所和哈佛CRSP CLIA实验室的工作人员和住院医生采集的鼻咽标本,采用FDA紧急使用授权实验室开发的检测试剂盒进行处理。提供N1和N2基因靶标的Ct值沿着样本采集日期、随机试管ID和唯一匿名机构ID,以反映标本来自不同机构。此处使用的标本源于2020年初,当时马萨诸塞州州的公共卫生工作导致对高级护理机构进行了全面的系列检测。这些公共卫生工作的拭子经过处理,用于临床诊断。样本采集日期范围为2020年4月6日至2020年5月5日,每家机构进行三轮采样。完成每轮研究的中位时间为2天(范围1-6天)。提供匿名Ct数据,并将N2 Ct值用于这些分析。对于此处列出的所有分析,将样本采集日期分组为采样轮次,并基于该轮次的平均采集日期进行分析(即,日期如图2和图S6和S7)。
SEEIRR模型
模型具有状态:易感(S)、暴露、无传染性、不可检出(E1)、暴露、无传染性、可检出(E2)、传染性且可检出(I)、恢复后仍可检出(R1)、恢复后不可检出(R2).
这些额外的暴露和恢复状态说明了个体检测PCR阳性但尚未或不再具有传染性的时期。我们使用马尔可夫链蒙特卡罗(MCMC)将该模型拟合到马萨诸塞州4家长期护理机构的数据中,以生成图2A中的轨迹,信息先验如表S1所示,将一致先验置于机构特定的基本生殖数R0和有效种子时间t0上。
随着时间的推移,长期护理机构中每个时间点观察到的Ct值的分布(图2B(模型预测的Ct分布(蓝色)与三个横截面样本中每个样本的实测Ct值(灰色条)拟合。所示为预期Ct分布的后中位数(黑线)和95% CrI(深蓝色条带)以及基于模拟观察结果的95%预测区间(浅蓝色条带)。请注意,预测区间比可信区间要宽得多,因为它们是由小样本量的模拟观测结果得出的
))变得更高(较低的病毒载量),变得更左偏。我们观察到,这些变化与机构中流行率(prevalence)的变化(即下降)有关。为了评估Ct值分布的这些变化是否确实反映了流行病增长率的潜在变化,我们使用Ct似然性对Ct值的每个单独横截面拟合指数增长和简单SEIR模型,以获得截至该时间点的流行病轨迹的后验分布(图2C(每个图显示了将基于Ct的SEIR模型分别拟合到病毒学数据的三个横截面的结果。显示的是发生率曲线的随机后验样本(红线)和最大后验概率轨迹(黑线
))。注意,每种拟合的唯一机构特定数据是Ct值和每个横截面样本的阴性试验次数。其他辅助信息包括流行种子时间(3月1日之后)和宿主内病毒动力学的先前分布。为了评估拟合,我们将每个拟合的预测Ct分布(图2B)和点流行率(图2D(基于Ct模型的检测阳性样本比例与观察到的检测阳性样本比例(灰色交叉)的拟合中位数(蓝色点)和95% CrI(蓝色误差线)
))与数据进行比较,并将这些拟合的增长率与基线估计值进行比较。所有Ct值模型参数的后验分布如图S8所示。
图2:观察到的Ct值的单横截面分布可用于重建马萨诸塞州州长期护理机构的流行轨迹
虽然这两组结果都是拟合模型,因此两者都不能被视为事实,但我们发现,拟合到一个数据截面的Ct方法提供了与基线估计相似的后验中值轨迹,这需要在每个时间点进行三个独立的点流行率和近乎普遍的测试。特别是,基于Ct的模型似乎可以准确地辨别样本是在感染率峰值后不久还是很久采集的。这两种方法在过去平均值和最近每日增长率的方向上一致(即,疫情当前是在增长还是在下降,以及增长率相对于过去平均值是否有所下降)。在大多数时间点,仅流行率模型和Ct值模型之间的平均增长率估计值非常相似,尽管仅流行率房室模型中的每日增长率似乎更早下降。然而,这些估计值具有很大的可变性,应该在这种背景下进行解释。这在图S7中尤其明显,其中其他机构在两种方法的估计值之间表现出更大的可变性。总体而言,这些结果表明,当流行轨迹受到限制时,如在封闭人群中,Ct值的单个横截面可以提供类似于三个不同采样轮的点流行率估计的信息。
为了确保我们的方法能够准确估计整个流行病曲线,我们使用经历随机SEIR流行病的合成封闭人群进行了广泛的模拟恢复实验。图S9显示了一个此类模拟的结果,展示了在试图估计爆发期间不同点的真实感染率曲线时,使用单个病毒学测试数据截面获得的信息。我们使用各种模拟来评估性能,包括仅使用正Ct值而不使用正分数信息的方法版本,不同人口规模的模拟,以及使用强度降低的先验分布的估计;详细信息见补充材料:模拟长期护理设施暴发,结果见图S10至S12。
模拟数据
所有模拟数据都是在相同的框架下生成的,但对潜在的流行病轨迹使用了不同的模型和假设。对于每种模拟,分四步生成数据:1)使用确定性SEIR、随机SEIR模型或高斯过程模型计算每日传播概率{ π t \pi_t πt};2)在模拟的每一天,在模型 I t B i n o m i a l ( N , π t ) I_t~Binomial(N,\pi_t) It Binomial(N,πt)下模拟新感染,其中N为模拟的群体大小,It为第t天新感染的数量(假设所有其他个体均已逃脱感染);3)在由下述测试方案和补充材料确定的模拟的特定日期对个体子集进行抽样:分析方法的比较;(4)对于第u天采集的每个样本,Ct值 X i G u m b e l ( C m o d e ( u − t i n f ) , σ ( u − t i n f ) ) , t i n f X_i~Gumbel(C_{mode}(u-t_{inf}),\sigma(u-t_{inf})),t_{inf} Xi Gumbel(Cmode(u−tinf),σ(u−tinf)),tinf为个体i的感染时间。
尽管在疫情暴发的早期阶段,没有真实的长期护理机构数据可用于评估该方法的准确性,但模拟实验表明,该方法可用于疫情的所有阶段。此外,尽管增长率估计存在很大的不确定性,但这些分析表明,可以使用单一的数据截面来确定疫情最近是在增加还是在减少。增长与下降的后验概率可用于此评估,当可信区间排除零时,其作用类似于假设检验,如果不排除零,则以更广泛的推断方式进行。尽管这对于SARS-CoV-2发病率来说是一个微不足道的结果,在许多情况下,病例、住院或死亡已经提供了疫情增长或下降的清晰图景,但对于检测能力受限的地点和未来疫情,我们的结果表明,一个由数百名受试个体组成的单一横断面随机样本,结合合理的先验知识(例如,将疫情种子时间限制在一到两个月的窗口内),可以用来立即估计疫情爆发的阶段。此外,该推断方法为组合多个测试日的横截面提供了基础。
基于多截面的流行病轨迹推断
虽然Ct值的单个横截面可以合理地估计由房室模型表示的简单暴发的轨迹,但更复杂的流行轨迹将需要更多的横截面来进行适当的估计。在这里,我们扩展了我们的方法来组合来自多个横截面的数据,使我们能够更可靠地估计整个流行病轨迹(材料和方法:多横截面模型和马尔可夫链蒙特卡罗框架)。在许多情况下,使用报告的病例数来监测疫情轨迹。将报告的病例限制为检测结果为阳性的病例,每日新增阳性病例数可用于计算Rt。然而,当病例的定义在流行过程中发生变化时,这种方法可能会变得模糊。此外,此类数据通常代表阳性检测的增长率,这种增长率可能会根据检测能力的变化而不是感染的发生率发生显著变化,需要仔细监测和调整,以说明检测能力的改变、感染与检测报告日期之间的延迟以及从流行率(prevalence)到发病率(incidence)的转换。死亡人数也被用来估计疫情轨迹,但这些数字被大大推迟,病例与死亡之间的关系并不稳定。相反,当监测采样的Ct值可用时,我们的方法可以通过提供Ct值分布和感染率之间的直接映射来克服这些限制。虽然病例计数方法由于测试率的变化而显示出偏差,但我们的方法提供了一种仅使用一个或几个监测样本来估计Rt的方法,并且这种方法可以适应随时间增加或减少的随机抽样方案,随测试可用性而增加或减少。
Ct观测值与每日感染概率之间的关系之多截面模型:
我们考虑有多个测试日 t 1 , . . . , t T t_1,...,t_T t1,...,tT,再次使用 π t \pi_t πt表示第t天感染概率,使用 X i t j X_{i}^{t_j} Xitj表示第i个个体在测试日 t j t_j tj采样的Ct值,其中 i ∈ 1 , . . . , n j i\in 1,...,n_j i∈1,...,nj,测试日 t j ∈ 1 , . . . , T t_j \in 1,...,T tj∈1,...,T。请注意,个体i可能指不同检测日的不同个体。设{ π t \pi_t πt}为任何第t天的每日感染概率,其中t天的感染可以在其中一个测试日使用PCR测试检测到。通过直接扩展单横截面模型的似然,当包括具有和不具有可检测Ct值的样本时,感染概率集{ π t \pi_t πt}的非参数似然由下式给出:
其中 n j − n_{j}^{-} nj−为检测日 t j t_j tj未检出的样本数量。
仅考虑具有可检测Ct值的样本,则似然函数:
这些似然函数中的任何一个都可以使用上述指数增长率模型来参数化。然而,指数增长率模型不太可能是较长时间段内真实发生概率的良好近似值,因此对于涵盖较长时间段的多个检测日,该模型可能不是良好模型。
多横截面似然主要用于拟合高斯过程模型,以观察到的Ct值集为条件估计每日感染概率{ π t \pi_t πt}。(See补充材料:用于拟合横截面病毒载量数据的参数模型)。SEIR模型也可用于多个测试日。如单横截面模型所述进行拟合,但使用其中一个可能性代替单横截面模型可能性,通过MCMC拟合获得后验分布估计值。
马尔可夫链蒙特卡罗框架
所有模型,包括使用Ct值的模型(SEIR模型、指数增长模型和高斯过程模型)和仅使用患病率的模型(SEIR模型、SEEIRR模型),均使用马尔可夫链蒙特卡罗框架进行拟合。我们使用Metropolis-Hastings算法来生成多元高斯或单变量样本。对于所有单横截面分析(图2和3),我们使用具有并行回火的该框架的修改版本:算法的扩展,其使用多个并行链来改进多峰后验分布的采样(47)。对于包括图4中的多个横截面分析,我们使用未修改的Metropolis-Hastings算法,因为并行回火算法的计算时间要长得多,并且这些分析得到更多数据的支持,受多模态的影响较小。在所有分析中,三个链运行80,000次迭代(高斯过程模型为500,000次迭代)。基于有效样本量大于200且潜在比例缩减因子小于1.1的所有估计参数评估收敛性,并使用coda R包进行评价(55)。所有假设的先验分布见表S1。
高斯模型
我们根据 π t \pi_t πt(每日发病率值)的后验分布给出了高斯过程模型的结果。对于可辨识性,将t的所有可能值上的 π t \pi_t πt的和设置为等于1;因此,所得到的估计值应被认为是每天相对于可能的天数的相对感染概率。当推断中仅包括阳性PCR检测结果时,我们直接估计 π t \pi_t πt。当阴性PCR检测也包括在内时,我们将 π t \pi_t πt乘以0 - 1之间的估计比例因子,这是整个发病率曲线中感染的绝对概率。Xu et al.描述并举例说明了GP作为先验分布用于各种传染病环境的发病率非参数推断的用途。在该模型下, π t \pi_t πt被考虑用于模拟的每一天,并且不像在SEIR模型中那样存在单独的流行病种子时间参数。考虑到第t天采样的Ct值将由首次采样时间前Amax天的感染产生,我们将Amax设定为35天。这包括了在第一个采样日期之前感染并且仍然可以检测到的大多数个体。因此,在对多轮测试数据进行高斯过程模型拟合时,在Tmin-Amax至Tmax的时间段内,每天估计一个日感染概率参数 π t \pi_t πt,其中Tmin和Tmax分别为最早和最晚测试日。最后,我们注意到高斯过程超参数v和p在模型拟合中是固定的。虽然理论上可以联合估计v和p(具有指数分布的先验),但我们发现,在拟合模型时,不约束这些参数会导致非常差的收敛链。
为了证明这些基于Ct的方法的性能,我们使用基于SEIR的模拟和来自暴发的样本Ct值(材料和方法:模拟测试方案)在各种测试方案下模拟暴发。我们通过R包EpiNow2使用报告的病例数(基于测试方案)对Rt估计的性能进行了比较,其中报告取决于检测能力和感染个体的症状状态,当一个、两个或三个监测样本具有观察到的Ct值时,总共约0.3%的样本(3000次测试在样本中分布)。
模拟测试方案:
估计倍增时间、生长速率或Rt的标准方法容易因试验政策的变化而被错误估计。为了评估这些变化对我们方法的影响,我们模拟了检测率的变化,并评估了对几种Rt估计方法的影响:使用EpiNow 2和报告的病例计数,使用基于Ct的方法和随机监测样本,以及使用PCR检测阳性单独和监测样本。我们在疾病爆发的两个阶段测试这些方法:一次是在疫情上升时,一次是在疫情下降时。对于每个分析时间点的随机样本,我们在1 - 3天的采样期内进行病毒学检测,检测日的样本量不同。结果见图3和图S13;详情请参阅补充资料:分析方法比较。
模拟长期护理机构疫情
为了确保我们的方法仅使用可检测的Ct值或所有PCR检测结果提供流行轨迹的准确估计,我们使用经历随机SEIR流行的合成封闭群体进行了广泛的模拟恢复实验(图S9、S10A)。我们评价了SEIR和指数生长模型拟合流行病增长、下降和峰值附近观察到的横截面Ct值后生长速率估计值的准确度和精密度(图S10B-D)。当使用所有样本结果时,SEIR模型在所有3个时间点均一致提供了无偏、约束的日生长率估计值。当仅使用可检测Ct值的分布时,生长和下降阶段的估计值准确,但显示出较宽的可信区间,而峰阶段的估计值略微向上偏倚。无论是使用全部阳性样本还是仅使用阳性样本,指数模型得出的平均增长率估计值都是一致的–略高于SEIR模型在高峰期和下降期得出的日增长率,反映出随着疫情开始下降,日增长率相对于过去平均值有所下降。
我们还发现,将模拟群体大小从300增加到5000不会改变我们估计的准确性,对95%可信区间宽度仅产生适度影响(图S11)。同样,对病毒动力学参数使用信息量逐渐减少的先验不会改变推断生长速率的准确度,但会增加生长速率估计值的不确定性(图S12)。
基于症状的报告
报告的病例计数来自感染后症状发作和报告延迟的模型。我们假设对数正态潜伏期平均值为log(5)天,标准差为0.418,如Lauer等人所估计。在模拟中,我们假设35%的个体有症状,潜伏期来自对数正态分布。然后,每个有症状的个体都有一定的概率在症状发作和检测报告日期之间延迟进行检测,其中该概率可能随疾病爆发的日期而变化。考虑了三种情况:平坦测试(固定测试概率为10%);增加检测率(检测概率从分析日前36天的10%线性增加到分析日前1天的20%);降低检测率(检测概率从分析日前36天的10%线性降低至分析日前1天的1%)。
使用R包EpiNow2(32,33)从这些模拟数据估计Rt,该R包可在https://github.com/epiforecasts/EpiNow2中获取,输入数据如下:
1.每日新增确诊病例数时间序列数据
2.指定的潜伏期分布,给出了感染和症状发作之间的延迟分布
3.指定的报告延迟分布,给出症状发作和病例确认之间的延迟分布
4.世代间隔分布的先验,指定感染者-受感染者对中感染之间时间的平均值和标准差
对于报告延迟分布,我们假设离散gamma分布的形状和尺度参数分别为5和2(平均值为2.5天,标准差为1.12天)。对于世代间隔,在我们使用SEIR模型的模拟中,平均世代间隔由下式给出: T c = 1 / σ + 1 / γ T_c=1/\sigma+1/\gamma Tc=1/σ+1/γ, σ \sigma σ和 γ \gamma γ分别为平均潜伏期和感染期(天)的倒数.世代间隔分布的方差为 V a r = 2 T c 2 2 Var=2{T_c \over 2}^2 Var=22Tc2,因此,𝑇𝑐=8天,标准差为5.66天。在EpiNow2中,将正态先验置于这些量上,两者的标准差均为3。
随机监测
对于监测Ct样本分析,我们使用随机抽样,其中每个个体在暴发的某个时间点接受检测的概率为0.3%。我们考虑了一个检测日(抽样0.3%的群体)、间隔一周的两个检测日(每天抽样0.15%的群体)和间隔一周的三个检测日。对于三个检测日,我们考虑了在其中任何一天采样的概率在三天内持平于约0.1%、从0.05%上升至0.10%再上升至0.15%或从0.15%下降至0.10%再下降至0.05%的情况。在所有这些设置中,总共进行大约3000次测试。图3B绘制了使用这些方法在流行病的两个时间点进行100次模拟的估计值:一个在峰值入射之前(真实Rt〉1),一个在峰值入射之后(真实Rt〈1)。
图S13比较了使用Ct值的监测样本分析结果与仅使用阳性率的结果,总监测样本量范围为100 - 3000。样本量分为1 - 3个测试日,多个测试日间隔一周,3个测试日进一步包括3天内的上升测试、下降测试或平坦测试。图S13A显示了中位数和四分位距图S13B显示了100个样本的中位后验Rt估计值的均方误差,图S13C显示了100个样本的Rt估计值的95%置信区间的平均宽度,图S13D显示了100个样本中95%可信区间完全大于1(流行性增长)或小于1(流行性下降)的比例。
与仅使用阳性百分比检验的模型拟合的比较
我们还考虑仅使用阳性率估计SEIR模型(即,Ct值低于40的检测比例)。在这种情况下,我们将SEIR模型拟合到流行率估计值,假设当且仅当个体处于感染状态时PCR阳性发生。与模型拟合的其余部分一样,我们使用相同的MCMC框架和先验,但对SEIR模型轨迹条件下观察到的阳性和阴性检测结果的数量使用二项式似然。
图3描绘了当疫情增长(第60天)和下降(第88天)时,每种方法的100个模拟中的后验中值Rt。除非只使用一个样本,否则即使测试数量在样本日内发生了实质性变化,拟合SEIR模型的基于Ct的方法也显示出最小的偏差。对于生长阶段的单样本估计,后验中值估计值高于真实值,因为R0值的范围与数据一致,R0的先验密度在1和10之间一致,中值为5.5,这使后验中值权重高于真实值。另一方面,基于报告病例数的方法在观察期内测试率增加时始终表现出明显的向上偏差,而测试率下降时则表现出显著的向下偏差。然而,基于Ct的方法显示出更高的可变性。这由贝叶斯推断模型捕获,因为所有基于Ct的方法在这100个模拟中至少实现了95%可信区间的标称覆盖(图S13)。
使用病例计数来估计Rt的另一种方法是将标准房室模型拟合到随机样本中阳性检测的观察比例。为了证明合并Ct值的价值,而不是简单地使用监测样本的阳性率,我们还将结果与在相同样本时间观察到的点流行率的SEIR模型进行比较,假设PCR阳性代表疾病的感染阶段。在这种替代方法中,SEIR模型的这种不规范导致在模拟的下降阶段不准确的Rt估计(图3B(每个流行病采样时间、检验方案和估计方法的100次模拟的Rt估计值。每个点是来自单个模拟的后验中位数。报告病例计数的Rt估计值使用EpiNow2估计值,监测Ct样本使用拟合至SEIR模型的一个或多个横截面的基于Ct的似然性。图右侧的半透明点是仅使用检测的二元结果拟合SEIR模型的监测样本,假设PCR阳性反映感染区室。采样日基于模型的真实Rt用黑星星和水平虚线表示,而Rt为1表示爆发平稳,用水平实线表示。
))。虽然一个更准确的模型可能会区分PCR阳性的感染阶段和持续时间,如SEEIRR模型,但这个简单的模型代表了一种方法,可以在感染状态和PCR阳性之间缺乏量化关系的情况下,从流行率数据推断发病率变化。
图3:根据观察到的Ct值从横断面监测样本推断流行轨迹,可以得到随时间变化的有效再生数Rt的几乎无偏估计值,而改变检测率会导致使用报告病例计数进行有偏估计。
A:模拟SEIR流行病中,按流行病采样时间和报告病例计数(顶行)和监测Ct采样(底行)的检测方案列出的每日阳性检测数量。对应于图(B)的分析时间由虚垂直线示出
我们还使用较小的样本量和给定样本量的测试日之间不同的测试部署来评估我们的估计精度。这些比较如图S13所示,图S13还比较了基于Ct的方法和基于正性的估计。基于Ct的方法在许多情况下表现良好,样本量低至200–500个测试。当测试稳定时,报告的病例计数提供了更精确的轨迹估计。然而,当报告的病例计数可能有偏差时,用于两个或三个监测样本的少量测试(例如,与常规病例检测一天所用的测试数量相同)可以提供无偏估计。
利用Ct值重建复杂的发生(Incidence)曲线
当数据稀少或在相对封闭的人群中,流行病开始时间大致已知时,简单的流行病模型有助于了解最近的发病趋势(补充材料:流行病种子时间优先级)。然而,在现实中,流行病通常遵循更复杂的轨迹,难以进行参数化建模。例如,SEIR模型不考虑影响病原体传播的非药物干预和行为改变的实施/放松,除非模型中明确规定。为了更灵活地从多个横截面估计疫情轨迹,我们开发了第三个感染发生率模型,使用高斯过程(GP)预测潜在的每日感染概率。GP方法提供了估计的每日感染概率,而没有对流行轨迹做出强有力的假设,仅假设同时期的感染概率是相关的,随着时间距离的增加,相关性降低(补充材料:流行传播模型)。电影S1展示了随着时间的推移,随着新的样本变得可用,如何使用该模型依次更新对整个流行病轨迹的估计(代表非药物干预的实施和随后放松的模拟)。电影S2显示了在较小样本量下估计的流行曲线的精度如何降低,其中每周200个样本足以可靠地跟踪流行曲线。电影S3显示了如果采样仅在疫情期间开始,估计如何保持准确。
流行病种子时间优先级
关于流行病学背景的外部信息可用于进一步限制估计的流行病轨迹。当将SEIR模型拟合到马萨诸塞州BWH的Ct值的单个横截面时,我们估计了观察时间之前的单个流行高峰的动态。因此,我们对流行病种子时间t0设置了统一先验,以反映两个流行病增长阶段开始的先验知识。具体而言,对于2020年6月1日之前采集的样本,我们假设接种时间为2020年2月1日至2020年4月1日。对于2020年6月1日至2020年8月1日之间采集的样本,我们假设2020年2月1日至采样时间前两周之间的接种时间未知。这就抓住了一个假设,即我们不知道这段时间的感染是由第一波的下降阶段还是第二波的增长阶段主导的。如果采样时间在2020年8月1日之后,我们假设第二波的种子时间在2020年6月1日和采样时间前2周之间未知。这确保了我们基于这些后续样本的第二波来估计发生率。
作为敏感性分析,我们假设2020年2月1日至所有采样时间的采样日期之间的流行种子时间未知(图S14D、E)。估计的增长率非常相似,尽管当流行种子时间限制较少时,前四周的后验密度更宽。一些后验估计值为双峰(图S15),其中相同的Ct分布可解释为是由非常近期和快速的流行增长(大多数高Ct值来自近期感染的上升)或流行下降的下降(大多数高Ct值来自清除期)引起的。虽然我们的方法能够准确地估计这些双峰后验分布,但重要的是要结合流行病学背景来解释它们。没有合适的先验,数学上正确的估计不一定在流行病学上可信。
高斯模型
对于我们方法的高度通用版本,我们开发了一个灵活的每日感染概率模型,表示为 π t \pi_t πt,任何一天𝑡 当天的感染可能导致在测试日(或其中一天)的PCR测试呈阳性。具体地说,我们使用高斯过程先验对每一天t的向量值 k t k_t kt;则第t天的发生率(incidence)为 π t = ( 1 + e − k t ) − 1 \pi_t=(1+e^{-k_t})^{-1} πt=(1+e−kt)−1.GP先验 k M V N ( ( 0 , . . . , 0 ) , K ) k~MVN((0,...,0),K) k MVN((0,...,0),K)的协方差矩阵为 K i j = η 2 e x p ( − ρ 2 D i j 2 ) K_{ij}=\eta^2exp(-\rho^2D_{ij}^2) Kij=η2exp(−ρ2Dij2),其中 D i j D_{ij} Dij为是i和j之间的差值,v和p是固定值为1.5和0.03的超参数。(见表S1)。参数p决定了协方差随天数增加而下降的速率,因此p值越高,相关性越小。更多详情请参见McElreath。
这种方法比基于SEIR模型的方法更灵活,使每日发病率能够反映传播率和接触模式的变化(尽管模型没有单独确定这些变化),以及易感个体库的消耗。因此,它不需要SEIR模型的严格参数假设。我们用它来模拟马萨诸塞州的疫情过程,因为各种政策和行为变化会随着时间的推移影响轨迹。
为了使用常规收集的RT-qPCR数据重建整个发病率曲线,我们使用了4月15日至11月10日在马萨诸塞州波士顿布里格姆妇女医院(BWH)对所有入院和未入院的急诊室患者进行的近全员检测中测得的阳性样本的匿名Ct值,2020(材料和方法:Brigham和妇女医院数据)。我们将这些与基于马萨诸塞州病例数的Rt估计值相一致(图4,A至C)。可检测的Ct分布的中值和偏度与Rt相关(图4B(按每周采样时间列出的病例计数得出的Rt与观察到的Ct值分布的中位数和偏度
)),与我们的理论预测一致(如图1所示)。Ct值中值上升(对应于病毒载量中值下降)。Ct分布偏度在春季末和夏季初下降,随着就地庇护所订单和其他非药物干预措施的推出(图4C(按采样周列出的Ct观测值的分布(小提琴图和点)和平滑中位数(蓝线)。红框突出显示了用于为子图(D)中的估计提供信息的数据。
))。但在夏末和初秋,随着这些措施的放松,中值下降,偏度上升,同时该州观察到的病例数也增加(图4A(马萨诸塞州每日新增确诊病例(灰色条)和估计的随时间变化的有效再生数Rt
))。
Brigham和妇女医院数据
马萨诸塞州波士顿布里格姆妇女医院的数据是用Hologic Panther Fusion SARS-CoV-2检测试剂处理的患者鼻咽标本。ORF1ab基因的Ct值与样本采集日期一同提供,采集日期范围为2020年4月3日至2020年11月10日。对于这些分析,我们按照流行病学日历上的采集周对样本进行分组,并使用每周的中点进行分析,如图4所示。2020年4月前两周的检测仅限于症状与COVID-19一致且需要住院的患者。4月15日之后,该平台的检测标准扩大到包括所有无症状住院患者、急诊室中未入院的有症状患者以及未临产需要检测的住院患者。入院的症状性ER患者在不同的PCR平台上进行检测,此处不予考虑。在本文的分析中,我们仅使用4月15日之后采集的样本。虽然这不是一个完全具有代表性的监测样本,但对没有寻求新冠肺炎治疗的住院患者进行常规检测所产生的队列比基于症状的检测更少偏倚,并代表了医院覆盖区域内病例的总体上升和下降。每日数据按周汇总。马萨诸塞州每日确诊病例数来自《纽约时报》,基于州和地方卫生机构的报告。
使用观察到的Ct值,我们使用单个横截面上的SEIR模型估计了感染的日增长率(图4D(从2020年6月14日开始的一周,SEIR模型拟合观察到的Ct值数据的单个横截面,估计的新发感染日增长率的后中位数(黄色箭头)和分布(蓝色阴影区域)。阴影密度与后验密度成正比。图S14显示了所有单周横截面的拟合
)和图S14和S15)和使用GP模型的完整流行轨迹(图4E(来自高斯过程(GP)模型的感染相对概率的后验分布(按日期)拟合所有观察到的Ct值(条带显示95%和50%置信区间(CrI),线显示后验中位数)。请注意,y轴显示的是相对而非绝对的感染概率,因为基础发病率曲线之和必须为1:只有阳性样品包括在估计中,因此假定所有样品都来自感染。
)和图S16)。在两种模型下推断出相似的时间趋势(图S17),GP模型提供的生长率估计值遵循使用观察到的病例计数估计的值(图4F(按日期比较GP模型(蓝线和阴影带显示后中位数和95% CrI)与使用观察到的病例计数(红线和绿色线和阴影带显示后中位数和95% CrI)进行Rt估计的新发感染估计日增长率。请注意,感染发生率(incidence)估计值是在首次观察到的采样日期2020年4月15日之前的日期进行的,最早可追溯至2020年3月1日,但x轴在2020年4月1日截断(图S19)。
))。虽然这些数据并非严格意义上的社区随机样本,并且观察到的病例计数不一定提供Rt值的真实数据,但这些结果证明了该方法仅使用常规检测收集的阳性Ct值重建流行轨迹和估计病例增长或下降的能力。我们通过在数据集中子采样不同数量的Ct值后重新拟合模型,评估了估计GP轨迹对较小样本量的稳健性(图S18)。有趣的是,我们仅使用一家医院常规生成的Ct值估计的流行轨迹与从废水数据中获得的社区水平病毒载量变化显著相关(图S19)。
图4:观察到的Ct值的横断面分布可以从马萨诸塞州布里格姆妇女医院基于医院的监测中估计复杂的全州流行轨迹
讨论
Ct值对公共卫生决策的有用性目前是许多讨论的主题。在许多地方持续观察到的一个无法解释的结果是,观察到的Ct值的分布在当前SARS-CoV 2大流行的过程中发生了变化,这导致了病毒适应性是否发生了变化的问题。我们的结果表明,这可以解释为流行病学现象,而不引起个体水平病毒动力学或检测实践的任何变化。然而,单独使用该方法无法证明任何特定设置的情况,因为病毒特性的变化或检测可用性的变化也可能导致Ct值分布发生此类偏移。我们发现,群体水平Ct分布的性质与现实环境中有效生殖数量或生长率的估计值密切相关,这与我们的理论预测一致。
使用在一次横断面调查中进行的多个不同测试的定量诊断测试结果,以前从病毒学数据推断出流行趋势。我们在这里描述的方法使用在当前大流行中观察到的现象以及发病率、自感染以来的时间和病毒学测试结果之间的关系,根据使用单一病毒学测试的一个或多个横断面调查的数据,在各种流行轨迹模型下,估计社区在流行曲线中的位置。模拟Ct值和观察到的Ct值与生长速率和Rt估计值的比较验证了这种通用方法。尽管存在采样变异性、病毒动力学、发病率、感染后时间和病毒学测试结果之间的个体水平差异等挑战,但根据一次或多次使用单一病毒学测试的横断面调查数据,在不同的流行轨迹模型下,评估社区在流行曲线中的位置。模拟Ct值和观察到的Ct值与生长速率和Rt估计值的比较验证了这种通用方法。尽管存在采样变异性、病毒动力学的个体水平差异以及比较不同实验室或仪器结果的局限性等挑战,但我们的结果表明,RT-qPCR Ct值及其对个体的所有变异性可以对种群水平动态提供高度信息。当测量结果被简化为二进制阳性或阴性分类时,这一信息就丢失了,SARS-CoV-2大类中的大多数都是如此。我们的结果表明,RT-qPCR Ct值及其对个体的所有变异性,可以高度了解种群水平动态。如SARS-CoV-2大流行中的大多数情况一样,当测量结果减少为二进制阳性或阴性分类时,这些信息就会丢失。
在这里,我们关注的是从总体中随机抽取个体的情况。因此,这种方法在可以独立于COVID-19症状获得代表性监测样本的情况下最有用,例如英国的REACT研究。即使是相对较小的横断面调查,例如在特定城市进行的调查,也可能非常有助于了解疫情的发展方向。跨区域的标准化数据收集和管理,沿着更广泛地使用随机抽样,将进一步提高这些方法的有用性,这证明了此类监测的另一种使用情况。这些方法使市政当局能够实时评价和监测各种流行病缓解干预措施的作用,例如,作为监测的一部分,进行哪怕是一次或少量的随机病毒学检测样本,而不是简单地依赖例行检测结果。
将这些结果外推至通过人口普查或大部分随机样本以外的策略获得的Ct值需要额外考虑。当检测主要基于症状的存在或接触者追踪工作时,感染个体更有可能在感染后的特定时间采样,这将影响Ct测量值的分布。当感染或症状出现与样本采集之间的延迟在流行病期间发生变化时,例如由于检测能力紧张,就会出现进一步的并发症。尽管如此,我们的模拟结果表明,流行轨迹仍然可以影响基于症状监测下测量的Ct值,尽管这种关联的强度将取决于图S4所述的许多额外考虑。还需要做更多的工作来扩展这里提出的推断方法,以使用非随机监测样本。
流行病增长率和Ct分布之间联系的总体发现对于根据新出现的SARS-CoV-2变异体解释病毒学数据非常重要。当通过全人群检测获得样本时,较低Ct值与新出现变异体之间的相关性可部分解释为与既存下降变异体相比,这些变异体具有较高的生长速率,且近期感染占优势。例如,最近对巴西马瑙斯P.1和非P.1变异样本的Ct值进行的分析最初发现P.1样本的Ct值显著较低。然而,在考虑症状发作与样本采集日期之间的时间(较短的延迟应导致较低的Ct值)后,该差异的显著性丧失。我们警告说,这一发现并不排除新的变异引起更高病毒载量感染的可能性;相反,它强调需要除监测测试数据之外的证据线。
这些结果对感染后每天观察到的病毒载量的真实分布敏感。不同的拭子类型、样本类型、仪器或Ct阈值可能改变Ct分布的变异性,导致特定Ct分布与流行轨迹之间的不同关系。在可能的情况下,设置特定校准(例如基于Ct值的参考范围)将有助于生成精确估计值。对于检测中使用的仪器和样品,该方法在可通过直接验证或与参比标准品比较估计群体水平病毒载量动力学的情况下最有用。在本文中,我们基于文献中测量的病毒载量的观察特性(症状发作后随时间推移的可检出比例、阳性标本的Ct值分布)生成了病毒动力学模型,并在估计生长速率时使用这些结果为关键参数的先验提供信息。因此,可以通过选择与模型拟合期间使用的观测相关的更精确、准确的先验来改进生长速率估计。如果结果来自多个检测平台,则应根据每个平台的特性指定不同的分布,调整模型以考虑这一点,或者,如果可能,应将Ct值转换为通用尺度,如log病毒拷贝。如果试验的这些特征随时间发生显著变化,则包含多个横截面的结果可能会出现偏倚,并且不可靠。
如果数据中包含可能影响病毒载量的个体水平特征,如症状状态、年龄和抗病毒治疗,并将其纳入Ct值模型,则结果也会得到改善。类似的方法也可能使用血清学调查,作为将感染后时间与其他传染病抗体滴度联系起来的工作的延伸。如果多种类型的测试(例如,抗原和PCR)同时进行,则组合信息可显著降低这些估计中的不确定性。如果变异株与不同的病毒载量动力学相关并变得常见,则也应将其纳入模型中。病原体的其他特征,如感染者和被感染者病毒载量之间的关系,也可能影响人群水平随时间的变异性。使用病毒学数据作为监测信息来源将需要投资,以更好地了解Ct值分布,因为新的仪器和技术上线,以及变异出现,并迅速表征这些分布,以应对未来出现的传染病。剩余的不确定性可以合并到贝叶斯先验分布中。
这种方法有几个局限性。虽然贝叶斯框架将病毒载量分布的不确定性纳入生长率的推断中,但参数假设和这些分布的合理强先验有助于识别。如果违反了这些参数假设,例如,当在干预可能影响传播率的时间段内使用SEIR模型时,推断可能不可靠。此外,当阳性病例非常少时,本文描述的方法以及发生率和Ct分布偏度之间的关系变得不太可靠,因此应谨慎解释结果,并在发生率低的时期增加样本量。在某些情况下,使用一个或少量横截面,观察到的Ct分布可能是由所有个体在快速流行增长开始时的感染早期、流行下降期间的感染恢复期或两者混合产生的(图4E和图S15)。因此,我们使用并行回火马尔可夫链蒙特卡罗(MCMC)算法进行单横截面估计,该算法可以准确估计这些多峰后验分布(47)。应在适当的流行病学背景下解释估计的中位增长率和可信区间:与其他数据严重不一致的估计增长率可以安全地排除。
如果同时使用不同机器或方案的结果来通知先验,则该方法也可能夸大病毒载量分布的不确定性。更精确地理解病毒载量动力学-特别是,以解释测量的流行病学和技术设置的方式对这些动力学建模-将有助于改进这种方法,并确定不同设置的Ct分布参数是否具有可比性。因此,对于SARS-CoV-2病例,应定期报告RTqPCR的半定量测量结果,对于未来出现的病原体,应优先早期评估病原体负荷动力学。使用对照测量,如使用检测到的病毒RNA与检测到的人RNA的比值,也可以提高Ct测量的可靠性和可比性。
Ct值是一种量级测量值,可提供潜在病毒动力学信息。尽管依赖单个Ct值进行个体水平决策存在挑战,但来自群体的许多此类测量值的聚集包含大量信息。这些结果表明,一次或少量随机病毒学调查如何能够最好地用于流行病监测。总的来说,Ct值的人群水平分布和一般的定量病毒学数据可以提供重要流行病学问题的信息,即使是来自单一的横断面调查。然后,可以在这种调查的基础上实施更好的流行病规划和更有针对性的流行病学措施,或者可以合并重复样本的Ct值,以最大限度地利用现有证据。