摘要
细菌群落对于维持生态功能和保持生物多样性至关重要。环境测序技术（如宏基因组学）的快速发展极大地提升了我们探测这种多样性的能力。然而，尽管取得了这些进展，将实证数据与生态系统建模联系起来的理论理解，特别是在类似于自旋玻璃的无序系统框架内，仍处于起步阶段。在这里，我们提出了一种使用无序系统理论解码微生物组数据的综合框架，该框架提供了关于塑造宏观生态状态的生态力量的见解。通过应用淬火无序广义Lotka-Volterra模型，我们分析了健康和患病人类肠道微生物组中的物种丰度数据。结果揭示了两种不同的物种相互作用网络模式，阐明了菌群失调可能导致微生物组不稳定的途径。相互作用模式因此提供了从健康到疾病转变过程中系统变化的窗口，为微生物群落动态提供了新的视角。我们的研究结果表明，无序系统理论有可能通过捕捉生态相互作用的本质及其对稳定性和功能的影响，来表征微生物组，从而加深我们对菌群失调与微生物动态之间联系的理解。
引言
细菌群落是生态功能和生物多样性的基本储库，它们在任何环境和宿主相关生态系统中都具有重要意义，从土壤到人类肠道。特别是，人类肠道微生物组通过与宿主代谢和免疫系统的相互作用，成为人类健康的关键调节器。菌群失调被定义为与胃肠道相关的健康微生物组组成的改变。最近的一些研究表明，肠道菌群失调与微生物物种相互作用模式的变化有关，而关注单个物种属性（如物种丰度分布）的新兴模式未能表征微生物组的异常状态。微生物组包含各种相互作用，这些相互作用已在体外群体中测量，并预计会调控自然群体。
近年来，随机非交互中性模型和逻辑斯谛模型已被应用于描述不同类型的社区（从人类到土壤微生物组）中的经验模式。然而，这些模型完全忽略了物种间的相互作用，因此似乎没有合适的建模工具来从菌群失调模型中获得新的见解。
因此，近年来提出了几种方法来推断微生物相互作用。然而，这种推断仍然非常具有挑战性和问题性：a) 典型微生物组数据集的非常高维性和缺乏纵向（长期）实验；b) 微生物物种相互作用的时间依赖性；c) 诸如资源、温度和pH值等非生物因素可以引入环境过滤效应并诱导有效相互作用；d) 实验和技术挑战引入了采样效应、假阳性物种、虚假相关性和数据偏差。因此，即使在元基因组样本无噪声和偏差的情况下，重建底层微生物动力学的物种相互作用网络仍然是一个艰巨的任务。
一种理论上更为优雅的方法是由May提出的，建议使用随机矩阵来模拟物种相互作用网络，即邻接矩阵的每个条目都是从给定分布中随机抽取的。由于无法经验测量相互作用强度，这种方法的优势在于参数数量从~ S 2（如果S是物种数量）减少到仅几个参数（例如2个），用于参数化分布。近年来，许多研究已将这一框架应用于通过广义Lotka-Volterra (gLV) 方程研究种群动力学，采用来自统计物理学的淬火无序系统技术（也称为“玻璃状”）。主要思想是专注于gLV的集合性质，这将只依赖于少数相关控制参数（例如，随机相互作用强度的均值 μ 和方差 _σ_2）。然而，尽管取得了一系列非常有趣的结果，这种方法主要局限于纯理论领域，直到最近才被用于解释受控实验中的微生物群落。
在这项工作中，我们提供了一个概念验证，证明了该框架适用于人类微生物组。具体来说，我们将从健康和不健康的肠道微生物组队列中推断淬火dgLV模型的参数。我们将展示人类肠道微生物组的不同生理状态对应于dgLV模型的不同相互作用模式。然后，我们将引入定量指标来量化稳定性和生态力对动态的贡献。
结果
2.1 无序广义Lotka-Volterra模型
dgLV模型描述了局部池中 S 个相互作用物种的浓度丰度随时间的演化。其中 N i 是第 i 个物种的数量， K i 是其承载能力， r i 是生长率， ρi 设置了物种动态的时间尺度，我们假设它不依赖于物种，即 ρi = ρ。系数 αij 是独立同分布的随机变量，且 <αij> = 0 和 <αij²> = σ²。我们还包括一个物种无关的移民率 λ，将其视为数学上规范问题的反射壁，并引入一个人口噪声项 ηi，其中 δ 表示狄拉克δ函数， T 设置噪声强度的尺度。为了方便表示，我们还引入其倒数 β = T -1。最后，我们要求相互作用是对称的，即 αij = αji。这样，我们可以将这个问题映射到一个精确可解的平衡热力学问题。如补充信息（SI）A部分所示，可以通过考虑MacArthur消费者-资源模型的准静态近似来证明这种对称假设和生长率与承载能力之间的线性关系 r i = ρKi。特别地，我们设置 ρ = 1。噪声项的存在使我们能够写出系统的Fokker-Planck方程并研究其稳态解（详见Altieri et al. [2021]，Altieri and Biroli [2022] 中类似哈密顿形式主义的详细推导）。在这里，我们采用Itô规定来处理随机动力学，以防止噪声引起的物种复苏。
副本形式主义是一种在无序系统中广泛应用的技术，允许我们推导出对应于dgLV对称相互作用 αij 的非交互哈密顿量（见方法和SI B部分的完整推导）。简而言之，与其尝试解决一个随机设置的问题，该形式主义考虑了许多系统的副本并平均它们的行为。这种方法有助于平滑随机性并揭示系统的典型行为。我们还假设了一个单一均衡情景，称为副本对称（RS）体制。尽管经验微生物群落的物种丰度分布（SADs）显示出胖尾特征，这种特征更符合不对称1RSB Mallmin et al. [2024] 情况或具有退火无序的GLV Suweis et al. [2023]，但我们考虑这种简化是因为它是可以得到模型参数与数据之间明确解析关系并且模型反演可行的体制。
另一方面，通过使用腔体论证 Mezard and Montanari [2009], Altieri and Baity-Jesi [2024]，可以解析地推导出模型的SAD（见方法和SI D部分）：
其中辅助变量 h 考虑了无序相互作用； z 是标准化高斯变量， q d 和 q 0 分别是复制自由能的二阶和零阶矩， m = 1 – _βσ_2(q d – q 0) 应大于零以保证理论的一致性 Biroli et al. [2018]。
RS假设还可以通过其对外部扰动（如外部温度、移民率或相互作用异质性）的稳定性来表征。为了研究RS相的稳定性，我们通过对复制自由能进行谐波展开来考虑理论的海森矩阵。当海森矩阵在所谓的“复制子部门”适当对角化后的主导特征值低于零时，会出现不稳定的平衡点：要么系统进入具有多个局部稳定极小值的副本对称破缺（RSB）相，要么发展成边际全RSB相，导致状态的分层结构和非常复杂的景观（也称为玻璃相）。因此，在后一种情况下，RS假设不再代表热力学稳定的相。基于无序系统中的标准计算，这个海森矩阵的主导特征值（称为“复制子” R）可以解析地评估为：
其中 χ 是物种对外部扰动的响应，仅跟踪未灭绝的物种。更多细节见SI C部分。
2.2 数据透视
我们现在考虑两个不同队列的横截面采样的肠道微生物组：一个是健康个体，另一个是患病个体。这里我们关注慢性炎症综合征，但原则上，我们所呈现的内容适用于表型不同的群体。我们可以将每个组的样本理想化为由dgLV Eq.(1) 的平稳分布生成，具有 α 的不同实现，但共享 μ 和 ρ（见方法）。我们提出的方法因此认为同一组的所有样本来自相同的统计系综，但通过不同的随机物种相互作用 α 的统计数据（例如，不同的 μ 和 σ）区分表征不同表型的生态制度。我们旨在测试每个队列将由一组不同的生态参数（θ = (μ, σ, T, λ））描述的假设。
为了从数据中计算顺序参数 (h, q d, q 0, K)，我们需要指定如何从数据中执行系综平均和无序平均。鉴于我们通常缺乏时间序列数据，并利用有效的平均场描述，我们认为系综平均可以估计为物种的平均值（即，不同物种是相同潜在随机过程的实现 Azaele et al. [2016]，称为中性假设），即 。另一方面，我们认为无序平均可以计算为样本平均值。换句话说，在给定表型（健康/不健康）内，每个测量的微生物组配置是从具有特定无序实现的dgLV模型的平稳分布中抽取的样本，即 ，其中 R 是样本数量且通常 R ≫ 1。此外，由于我们对每个微生物组中相互作用的精细细节了解有限，假设所有具有给定宏观状态（即 (h, q 0, q d, K) 空间中的一个点）的社区经历相同的人口噪声 T、移民率 λ 和无序参数 μ 和 σ 是合理的。最终，我们可以从数据中计算顺序参数如下：
其中 N a,j 给出了样本 a 中物种 j 的种群密度。
为了从数据中评估承载能力 K i，如同大多数使用dGLV的研究 Bunin [2017], Biroli et al. [2018], Altieri et al. [2021], Suweis et al. [2023], Mallmin et al. [2024] 所做的那样，我们假设在每个队列中，所有物种具有相同的承载能力 K i = K。在这里，我们将 K 值分配为每个队列中可用样本中物种的最大相对丰度的平均值，即 。
还要注意的是，从宏基因组学中，我们只能获得物种丰度的组成数据 Pasqualini et al. [2023]。这对于正确处理不同样本并最终得到一致的分析至关重要。此外，这也允许我们对一些顺序参数进行生态解释。特别是，由于数据的组成性，我们有 和 （其中 称为平均局部多样性），而 q d 提供了每个副本中物种丰度之间成对乘积的信息。图2显示了健康和不健康队列之间的这些顺序参数值以及两个随机化队列的值。在数据中，我们观察到健康和不健康队列之间存在系统差异，指出了健康样本中较高的平均局部多样性（面板A、B），正如Pasqualini等人 [2023] 所观察到的那样。面板C突出了不健康患者中较高的 q d，这是这些样本中物种相互作用减弱的标志。在接下来的部分中，我们将通过推断两个队列中的物种相互作用进一步研究这一方面。
从建模的角度来看，可以严格证明将丰度归一化为1相当于在自由能表达式中引入额外的拉格朗日乘数。然而，这仅作用于线性项，贡献于移动随机相互作用的均值 μ，但对异质性 σ 和相图没有影响。更多细节请参阅SI（B1节）以及Altieri et al. [2021]。在Altieri et al. [2021]中，作者还展示了物种丰度的全局归一化约束如何反映在随机复制品模型 Biscari and Parisi [1995] 中的一对一映射。（Diederich and Opper [1989] 最初引入的复制子方程并在副本形式主义中重新表述 Biscari and Parisi [1995], Altieri et al. [2021] 描述了通过随机耦合进化的复制品集合的动力学。）
2.3 物种相互作用模式表征微生物组的状态
因此，我们在向量 ρ = (h, q d, q 0, K) 中收集了从数据中估计的所有参数。正如我们将在方法部分详细说明的那样，我们开发了一种矩匹配推断算法来推断模型参数 θ。因此，我们生成了尽可能接近数据的模型的顺序参数。该方法的思想是引入一个成本函数 C(θ|ρ)，表示某些自洽方程的总相对误差。如果参数 θ 使得自洽方程的右侧等于左侧，则可以认为问题得到了解决。由于该成本函数的景观呈现多个最小值，我们进行了多次优化程序以收集一组可能的解决方案，并从中保留了前30个。
首先，我们发现不同解决方案 θ_* 的优化问题提供了对底层微生物组种群的生态洞察。如图3A所示，健康和患病微生物组的推断人口噪声强度（T）和相互作用的异质性（σ）在 T - σ 平面上聚类。特别是，健康队列的SAD在推断过程的不同解决方案中是稳健的，如图3A插图中的不同曲线叠加所示。另一方面，从不健康患者推断出的SAD对不同解决方案具有高度敏感性。特别是，其中一些显示高丰度物种的模式（图3A中的浅红色线条），这是肠道中优势菌株的标志。一致地，通过推断参数 μ 和 σ 生成的相互作用 P(αi,j) 的分布在健康和患病队列之间不同，给出了不同的相互作用模式（见图3B），这一结果与Bashan等人 Bashan et al. [2016] 发现的结果一致。特别是，我们发现菌群失调降低了相互作用强度的异质性，这一结果也在将相关性作为相互作用代理时观察到 Seppi et al. [2023]。
我们随后评估了推断的 σ 和 β = 1/T（统计物理方法中的逆温度）与dgLV的临界副本对称破缺（RSB）线（R = 0）的距离，其他参数保持不变。我们再次发现，我们的优化协议每个解决方案对应的复制子值 R 对于所调查的两种微生物组表型显著不同（见图4A）。特别是，患病微生物组在副本对称假设下更接近边缘稳定性 Altieri et al. [2021], MÃ(C)zard et al. [1987], de Almeida and Thouless [1978]。此外，通过研究由公式（2）给出的SAD形状，我们可以估计生境（由物种相互作用表示）和中性（由出生/死亡和移民表示）生态力量之间的比率，这可以通过量 ψ 来捕获 Wu et al. [2021]。它检测SAD中内部模式的出现作为生境过程的标志（见SI，E节）。
受场论论证的启发（见方法部分 自由能景观探索：副本公式），我们可以将理论中的 m 参数称为“质量”，如上文公式（2）中定义的那样。当 m → 0 时，理论的顺序参数发散，系统进入无限增长的病态区域。因此，“质量”项可以被视为补充稳定性度量，能够捕捉过渡到无限增长区域。
在模型中，有两种效应竞争塑造群落结构。一方面，我们有生境效应，编码在无序相互作用中，因此由参数 μ, σ, 和 K 跟踪。它们的整体效应是选择性的，倾向于将SAD集中在典型的丰度值周围。另一方面，我们有中性效应，编码在随机动力学和移民中，控制SAD的低丰度区域。当人口噪声幅度强于移民（ν < 1，如我们的情况）时，SAD表现出低丰度的可积发散。相反，在 ν > 1 的情况下，没有发散，SAD是模态的。由于相互作用是随机的，观察到内部模式的概率可以估计为具有非平凡稳态解的SAD实现的分数。这种量，称为生境-中性比率，可以解析地评估为：
当 ψ ≈ 1 时，生境和中性力量对动力学的贡献相当，因为SAD同时存在低丰度发散和有限丰度模式。最后，如果典型丰度发散，我们进入了无限增长阶段，这意味着质量和生境-中性比率 ψ 不是独立的，如 ψ 的解析表达式所示。对于这一结果的详尽推导，请参见SI（E节）。通过获得的模型参数，我们能够评估健康和患病微生物组的 m 和 ψ。同样在这种情况下，健康和患病微生物组明显聚类，如图4所示。不健康的微生物组更接近无限增长阶段，生境-中性比率比健康情况高出五个数量级 ψU ≈ 10^5 ψH。这使我们得出结论，选择压力在疾病状态下要大得多，而在健康状态下，出生和死亡效应是动力学的关键驱动因素。这些结果也由图3（面板a）插图中的SAD形状证实。
讨论和结论
在通过无序系统和随机矩阵理论的视角探索肠道微生物组的过程中，我们开启了一项雄心勃勃的旅程，旨在将复杂的理论概念与实际的环境微生物组数据分析联系起来。特别是，我们使用无序广义Lotka-Volterra模型（dgLV）来表征健康和不健康的肠道微生物组。这种将无序系统领域的先进理论方法应用于现实世界生物学数据的尝试，代表了理解微生物群落复杂相互作用的重要一步。
我们研究的一个微妙之处在于解决固有的局限性，并在未来的研究中寻求改进的空间。一个显著的挑战是经验物种丰度分布（SAD）与淬火dgLV模型预测的理论SAD之间的差异。虽然经验数据显示了广泛的物种丰度，遵循幂律分布，但模型预测往往在SAD尾部表现出指数衰减。这种不匹配强调了需要进一步完善模型以准确捕捉现实数据中观察到的细微模式。例如，最近有人提出引入退火无序可以生成更接近经验数据的SAD Suweis et al. [2023]。另一种可能性是设置dgLV模型参数以重现多吸引子相。事实上，在这一区域，SAD表现出更异质的形状 Arnoulx de Pirey and Bunin [2024]。我们已经打算通过结合一阶副本对称破缺（1RSB）计算和不对称相互作用情况下的动力学平均场理论分析来探索这一后续方向。
另一个相关限制是生成从dgLV模型中提取的物种丰度样本，这些样本反映了观测到的经验数据的统计特性。在我们当前的框架中，每个微生物组样本都可以从 p(N |ξ) 中提取，其中 ξ 是无序的一个实现。然而，p(N |ξ) 在 N → 0 附近呈现出幂律指数 ν - 1，其中 νH,D ≈ 10^-3。这会产生数值不稳定并主导抽样过程，从而难以生成代表性的合成样本。此外，虽然微生物组数据本质上是组合的 Pasqualini et al. [2023]，但dgLV模型中的物种种群 N i 是正实数。然而，正如已经指出的那样，可以在这些方程中引入归一化 ，不会改变所提出解决方案的结构 Altieri et al. [2021]，因此不应影响我们工作的结论。该研究还揭示了dgLV模型背景下人口噪声的作用。实际上，在 T → 0 的极限下，SAD转变为截断高斯分布，峰值位于零处，反映了灭绝物种的比例。与零噪声情况下观察到物种灭绝的情况相反，dgLV模型中包括人口噪声表明不存在物种灭绝，支持“一切都在各地”的假设 Grilli [2020], Pigani et al. [2024]。换句话说，在这个框架中，数据中的零值是由于采样效应而不是本地物种灭绝造成的 Grilli [2020], Pasqualini et al. [2023]。
总之，我们的工作为未来的研究开辟了新的途径，特别是在改进理论模型以更好地与经验观察一致以及探索微生物组背景下的物种丰度分布细微差别方面。此外，整合其他形式的环境变异性及物种特异性特征可以提供生态系统动态的更全面视图，正如“One Health-One Microbiome”框架 Tomasulo et al. [2024] 所提议的那样。本研究获得的见解不仅增强了我们对微生物生态系统复杂动态的理解，还为生态建模提供了更先进的分析工具，大大提高了我们解释和管理这些重要生物系统的能力。
(全文结束)
