细菌感染与噬菌体防御的协同动力学模型：序参量视角下的多物种网络构建与爆发机制分析

本文从协同学视角构建并分析了一个包含易感细菌(S)、感染细菌(I)、噬菌体(P)和抗性细菌(R)的四物种动力学模型，通过推导序参量方程揭示感染爆发的自组织机制。研究表明，序参量幅值决定感染动态的时间特性，而抗性细菌通过CRISPR系统产生的防御效应（γR项）能有效降低有效接触率β=β0/(1+γR)。模型首次建立了线性预测方程（如I=(F0/(λ2+k1))P），为多物种网络调控提供理论框架。

　　

1 引言

流行病学系统通常由不同物种相互作用的复杂网络构成。虽然孤立种群通常表现出相对简单的动力学行为，但网络生理学领域的挑战在于理解不同类型物种间的相互作用如何塑造整体网络动力学。在此背景下，重要的第一步是考虑基于常微分方程(ODE)和耦合微分方程模型的平均场近似，这些模型由于其相对简单性常常允许解析求解方法。

对于病毒感染，ODE三物种TIV模型捕捉了靶细胞、感染细胞和病毒颗粒相互作用的基本动力学。同样，在噬菌体感染背景下，我们处理的是靶细菌、感染细菌和噬菌体——后者充当病毒角色。研究细菌感染和噬菌体的作用是一项重要任务，也是尝试在现代医学中使用噬菌体治疗人类某些疾病不可或缺的一步。

为此，文献中研究了简化和广义的三物种模型。特别是，作为其综合研究的一部分，Skanata和Kussell (2021)考虑了细菌可以存在对给定入侵噬菌体的抗性和非抗性表型。增加抗性表型的浓度是对抗噬菌体攻击的一种防御机制，因为这种机制降低了噬菌体与非抗性细菌之间的有效接触率，使得在适当条件下感染会消失。

虽然动力学系统视角通常是分析噬菌体感染模型不可或缺的工具，但很少有人关注利用更具体的协同学动力学系统概念在这方面。然而，在COVID-19大流行之后，已经表明协同学可以补充现有的动态系统方法来理解流行病学和病毒动力学模型。

在本简报中，将考虑一个包含如Skanata和Kussell (2021)中抗性表型的简约四物种噬菌体感染和细菌防御模型。模型的相对简单性将允许分析方法。研究的目的是明确推导模型的序参量，并展示它如何决定噬菌体攻击的初始组织和相应的防御反应。此外，目的是推导确定系统沿序参量和剩余方向演化的振幅方程。在此背景下，还将推导两个系统动力学的近似模型。示例模拟将说明一些分析结果。

2 材料与方法

如上所述，我们的起点是包含易感(S)和感染(I)细菌以及噬菌体负载(P)的三物种模型。与病毒动力学的TIV模型一致，模型方程如下：

dS/dt = -k₀(P)S + μS(1-S/K)

dI/dt = k₀(P)S - k₁I

dP/dt = qI - k₂P

其中k₁, k₂, q > 0，其中t表示时间，k₀描述了P依赖的S→I转变速率，k₁和k₂分别表示感染细胞和噬菌体的衰减率，q描述了每个感染细菌的噬菌体产生率。S的演化方程涉及一个具有生长速率μ和容量K > 0的逻辑增长项。

为了考虑引言中提到的主动防御机制，模型(1)以两种方式进行了修改。首先，假设当感染细菌I出现在细菌种群中时，抗性细菌突变(R)会像这样生长：

dR/dt = αI(1-R/R_m)

其中α ≥ 0和R_m > 0分别表示抗性细菌的生长速率和最大浓度。从机制角度来看，方程(2)通过所谓的CRISPR系统捕捉了细菌的适应性防御机制。CRISPR系统允许细菌记忆攻击的噬菌体，从而使它们对未来攻击产生免疫。这样，在入侵噬菌体存在的情况下，出现了噬菌体抗性细菌。

其次，通常有几种机制可以使抗性细菌R减缓或停止噬菌体感染。为简洁起见，仅考虑了R对S→I转变速率的影响。这样做，转变速率k₀成为R和P的函数，并表示为：

k₀(P,R) = (β₀/(1+γR))P

其中γ ≥ 0测量主动防御机制的有效性。基本上，方程(3)指出抗性细胞的存在降低了噬菌体与易感细菌之间有效接触的机会，因此R依赖的有效接触率读作β = β₀/(1+γR)。

从协同学视角来看，对于模型(1-3)，推导了序参量和振幅方程。为此，使用状态向量X = (S,I,P,R)，考虑了从初始固定点X_st,0 = (S_st,0, I_st,0, P_st,0, R_st,0)附近开始的噬菌体感染。随后，借助从线性稳定性分析获得的特征向量v_j，振幅A_j被隐式定义为：

X = X_st,0 + Σ_j=1⁴ A_jv_j

通过构建由满足w_iv_k = δ_ik（Kronecker符号）的向量w_j张成的双正交基，振幅被明确表达为：

A_j = w_ju = w_j(X - X_st,0)

其中u表示差向量u = X - X_st,0。从模型方程(1-3)和方程(5)，最终推导出形式如下的模型振幅方程：

dA_j/dt = λ_jA_j + N_j(A)

其中j=1,2,3,4，A构成振幅向量A = (A₁, A₂, A₃, A₄)。在方程(6)中，λ_j表示系统的特征值，N_j是振幅中的非线性函数。模型的序参量及其序参量幅值被确定为对应于模型潜在正特征值λ_j的特征向量v_j及其幅值A_j。

3 结果与讨论

3.1 振幅方程视角

固定点分析表明，模型(1-3)表现出无噬菌体固定点，定义为：

S_st ≥ 0, I_st = 0, P_st = 0, R_st ∈ [0, R_m]

对于μ = 0。对于μ > 0，方程(7)成立且S_st = K。如方法部分所述，假设在时间t < 0，即感染发生之前，系统保持在固定点(7)。该固定点被称为初始固定点，记为X_st,0。在时间t = 0时，细菌种群被浓度为P(0) > 0的噬菌体感染，使得状态被移出其固定点。因此，该模型描述了从时间t = 0初始噬菌体感染P(0) > 0开始并以方程(7)定义的无噬菌体状态或如果存在的地方性状态结束的感染爆发。

在X_st,0处的线性稳定性分析表明，对于μ ≥ 0的模型表现出特征值λ₁ = -μ，λ₄ = 0，以及：

λ_2,3 = -(k₁+k₂)/2 ± √[(k₁+k₂)²/4 - k₁k₂ + qβ₀f_st,0S_st,0]

其中f_st,0 = 1/[1+γR_st,0]，上（下）号对λ₂（λ₃）成立。对于μ > 0，在方程(8)及后续内容中我们必须代入S_st,0 = K。可以证明对于任意模型参数都有λ₃ < 0。相反，λ₂可以取正值或负值。

案例I定义为qβ₀f_st,0S_st,0 < k₁k₂ ? λ₂ < 0，这等价于R₀ < 1，因此固定点X_st是对于μ = 0的中性稳定/对于μ > 0的渐近稳定固定点。没有感染爆发。

相反，案例II characterized by qβ₀f_st,0S_st,0 > k₁k₂ ? λ₂ > 0，这等价于说R₀ > 1成立。固定点不稳定。感染动力学描述了感染爆发。

3.2 影响

3.2.1 序参量：基本构建块和初始组织

模型最多表现一个正特征值。因此，在案例II情景下，具有λ₂ > 0 (R₀ > 1)，系统表现出由v₂给出的序参量及其序参量幅值A₂。让我们将状态动力学分为两部分X = X_out + X_s，其中X_s = A₃(t)v₃（对于μ = 0）和X_s = A₁(t)v₁ + A₃(t)v₃（对于μ > 0）描述了沿稳定方向朝向序参量的动力学，而X_out捕获了剩余动力学。

最初，即对于t ≈ 0，我们有：

X_out ≈ K + v₂A₂(0)exp(λ₂t)

其中K是常数。方程(14)描述了沿不稳定方向远离初始固定点的动力学（即向外动力学）。相反，X_s描述了朝向不稳定方向的动力学，即朝向序参量。因此，序参量v₂描述了多物种生理网络的新兴组织，其幅值A₂描述了这种组织如何随时间演化。

由于X_s最初在量级上随时间衰减，当考虑初始感染动力学时，我们可以忽略其对状态动力学的贡献。如果这样，那么任何由ΔX = X(t) - X_st,0定义的状态变化近似由下式给出：

ΔX ≈ v₂ΔA₂ ≈ v₂(exp(λ₂t)-1)

方程(15)再次说明序参量描述了多物种网络感染爆发组织的基本构建块，包括防御反应（参见分量v_2,R）。

3.2.2 停止机制

沿v₂的指数增长如方程(15)所述，在某个时间点减慢并最终停止。与稳定性分析一致，让我们假设δ, I, P, ω是ε阶的小量。那么N_j当将N_j解释为差变量的函数时可以展开，使得A₂的振幅方程读作：

dA₂/dt = λ₂A₂ - [v_3,Pβ₀f_st,0][γS_st,0f_st,0Pω - Pδ] + O(ε³)

注意，对于任何t > 0，有δ < 0。因此，网络生理学产生了两种减缓A₂指数增长的机制：噬菌体存在下易感细菌的衰减，由相互作用项-Pδ > 0测量；以及噬菌体存在下抗性细菌数量的增加，由相互作用项Pω > 0测量。前一种机制是一种被动机制，仅表示由于资源（即易感细菌）的衰减，指数感染爆发减慢。后一种机制是一种主动机制，表明将抗性细菌突变引入细菌菌落具有减缓噬菌体入侵的期望效果。

3.2.3 线性预测方程

方程(15)暗示所有物种最初满足形式如下的线性回归方程：

X_i ≈ a_i,j + r_i,jX_j, r_ij = v_2,i/v_2,j

因此，网络的任何物种都可以用于预测任何其他网络物种（假设所有i,j的r_ij ≠ 0）。网络组件通过线性序参量链接耦合。例如，噬菌体种群大小P可用于构建回归模型，如：

S = a_S,P - (F₀/λ₂)P

I = (F₀/(λ₂+k₁))P

R = a_R,P + (αF₀/((λ₂+k₁)λ₂))P

其中a_X(i),X(j)是取决于X_st,0的截距参数。如方程(18)所示，由于I_st,0 = P_st,0 = 0，有a_I,P = 0。

3.2.4 细菌生长项的有限影响

显然，细菌生长项可能影响S动力学。然而，它不影响序参量特征值λ₂，也不影响序参量v₂在(I,P,R)空间中的方向，这是主要关注点。因此，由序参量动力学确定的在(I,P,R)空间中的初始爆发动力学完全不受细菌生长项的影响。

3.3 二维感染/传染物种动力学和双指数动力学近似

感染和传染物种I和P的动力学完全由振幅A₂和A₃描述。这样做的原因是特征向量v₁和v₄在I-P子空间中不表现分量。从A₂和A₃到I和P的映射读作：

[I; P] = v′₂A₂ + v′₃A₃

其中v′₂和v′₃表示v₂和v₃到I-P子空间的投影。从方程(19)可以看出，I和P的初始演化满足双指数函数，如方程(20)所示：

[I; P] ≈ v′₂A₂(0)exp(λ₂t) + v′₃A₃(0)exp(λ₃t)

方程(20)说明最初，即对于较小的t，I和P的演化满足双指数函数。

3.4 缩放模型

使用变量变换s = S/S_st,0, i = I/S_st,0, p = k₁P/(qS_st,0), r = R/R_m，模型(1-3)变为：

ds/dt = -k′₀s + μs(1-s)

di/dt = k′₀s - k₁i

dp/dt = k₁i - k₂p

dr/dt = α′I(1-r)

其中k′₀(r,p) = (β′₀/(1+γ′r))p

且β′₀ = β₀S_st,0q/k₁, γ′ = γR_m, 和α′ = αS_st,0/R_m。缩放模型表现出以下两个特性。首先，变量s, i, p, r是无量纲的。其次，变量变换P → p将噬菌体变量P转变为类细菌变量。也就是说，正如模型描述易感细菌转变为感染细菌一样，缩放模型描述当以变量p（而不是P）表示噬菌体时，感染细菌以1:1方式转变为噬菌体单位。从数学上讲，从方程(21)可以看出，p由于项k₁I而增加的速率与感染细菌数量由于项-k₁I而衰减的速率相同，这意味着模型描述了上述1:1转变。

3.5 模拟

使用前向欧拉模拟方案求解方程(21)和(22)。在第一次模拟中，仅考虑了被动防御机制，其中α = 0和γ = 0。在第二次模拟中，考虑了主动防御机制，其中α = 1和γ = 5。其余参数和初始条件为：k₁ = 0.2/天，k₂ = 1/天，β′₀ = 5/天，s_0,st = 1, i_st,0 = 0, p_st,0 = 0, r_st,0 = 0.01, 和p(t=0) = 0.001。为简洁起见，仅呈现最相关的相曲线，并且仅考虑初始爆发阶段的第一周。

如面板(a)和(b)所示，相曲线迅速收敛到序参量v₂，随后沿v₂演化。根据定义，序参量v₂不捕获沿稳定方向v₃的动力学。如面板(c)和(d)所示，双指数近似可以描述朝向序参量的瞬态初始动力学（即沿v₃的动力学）以及随后沿序参量v₂的动力学。

比较可见，由于噬菌体抗性细菌的影响，实际动力学与序参量动力学的差异更大。同样，实际动力学更早地偏离双指数近似。这些观察结果并不令人惊讶，因为主动防御机制减缓并最终停止感染爆发。因此，实际动力学将更早地偏离序参量动力学，一方面，和双指数动力学，另一方面。粗略地说，主动机制削弱了网络组件之间的线性序参量链接。

4 结论与局限性

我们得出结论，在方法和模拟部分规定的适当条件下，序参量及其幅值表征了噬菌体感染和相应细菌防御的（自）组织。在物理学中，已经进行了各种实验来专门测试上述自组织理论（和协同学）的预测。因此，正如在物理学中一样，上述结果可能作为进行噬菌体感染实验室测试以检验序参量假设的基础。此外，我们得出结论，如上推导的线性回归模型可用于基于更便于测量的物种来估计难以观察的物种。

为简洁起见，在当前研究中，未详细研究地方性固定点的性质。同样，当前研究仅限于考虑一种可能的防御机制，而忽略了替代机制。更全面的研究（超出本简报范围）可以通过推广上述结果来克服此类局限性。