化学反应网络的热力学空间

《Physical Review X》：Thermodynamic Space of Chemical Reaction Networks

　　生命系统在非平衡状态下运作，持续消耗能量以维持有序、功能性的状态。它们的涌现行为通常依赖于一系列相互连接的化学反应网络（CRNs），这些网络由外部通量驱动，使某些物种的浓度保持固定。因此，揭示控制这些CRNs运作的原理对于理解生命系统如何产生和调节复杂性至关重要。虽然动力学在塑造详细动态现象方面起着关键作用，但CRNs的操作范围从根本上受到热力学的限制。在这里，我们引入并分析推导了CRN的“热力学空间”，即在给定能量预算下可以实现的稳定浓度范围。我们为反应亲和力建立了类似的界限，揭示了全局热力学性质（如总的非平衡驱动力）如何限制局部非平衡量，并进一步表明这些界限可以用来测试候选网络拓扑结构并揭示隐藏的多分子路径。我们在各种典型示例中展示了我们的结果，证明了复杂行为的出现与非平衡条件的存在紧密相关。通过提供一个分析CRNs的通用工具，所提出的框架为加深我们预测复杂非平衡现象和设计人工化学系统的能力奠定了基础，这一切都基于对潜在热力学性质的单独了解。

I. 引言
生命系统远离平衡状态运作，消耗能量以维持有序状态，对抗熵增加的趋势[1]。这种生命的内在非平衡性质是由基于化学反应的生物过程维持的[2]。化学反应网络（CRNs）涵盖了众多这些生化过程，包括基因调控、新陈代谢和细胞信息处理[3–6]。化学反应网络理论为表征CRNs提供了有力的框架，并揭示了复杂行为是如何从潜在的分子相互作用中产生的[4,7]。这种复杂性表现为各种形式，如非平衡稳态、多稳态、混沌和振荡[8–12]。这种广泛的可能性对于建立稳健的生物功能至关重要，同时与不断从环境中吸收和耗散能量的需求紧密相关。实际上，非平衡条件在决定生物系统的运作中起着主导作用；例如，通过能量消耗可以放大所需的生化状态的丰度[13–19]，自由能耗散触发复杂动态的出现[10,20]，化学驱动调节大规模结构的形成[21–23]。由于其重要性，已经发展出一种热力学一致的化学网络理论来描述CRNs在平衡状态内外的热力学特性——通过量化通量的不可逆性并将其与自由能变化相关联[24–28]。最近的研究集中在CRNs的拓扑结构对其非平衡行为的影响上，例如，研究了它如何影响由于外部化学计量器（即固定选定化学物种浓度的装置）引起的耗散通量的性质[29,30]，以及它在确定扰动响应的热力学界限中的作用[31]，并表征了化学系统如何调整其效率[32]。尽管理解非平衡现象的重要性毋庸置疑，但热力学驱动力对生化系统所能施加的一般限制仍然难以捉摸。特别是，揭示全局非平衡性质如何限制特定的物理可观测量并塑造CRNs的操作范围是一个巨大的挑战。最近的研究对非平衡放大、信息处理的准确性以及反应-扩散模式的对比提供了一些启示[22,33–36]，但这些分析仅限于线性（或最多是催化性的）CRNs，这些CRNs总是可以用简单图表示并用拟平衡量来理解[37,38]。然而，尽管许多生物功能是由多分子反应实现的，但这些原理和限制如何适用于由超图描述的通用CRNs[30,39]仍是一个未解决的问题。

在这项工作中，我们解决了这个问题，并推导出了适用于任何CRN的反应亲和力和物种浓度的基本热力学界限。我们的方法超越了以前的研究，通过解决考虑任意化学反应网络拓扑结构的挑战性任务。这确实是一个关键步骤，使得我们的结果可以应用于其中复杂的多分子反应起关键作用的系统，如代谢途径、自组装过程和循环神经化学网络[40–42]。我们引入了热力学空间的关键概念，即在稳态下生化系统必须操作的浓度范围。通过提供一个理解通用CRNs的热力学限制的通用框架，我们的工作为预测和分析广泛的生物和合成化学系统中的复杂行为开辟了新的途径，并为的反应网络拓扑结构提供了推断工具[43]。作为教育性但非平凡的例子，我们在各种情况下展示了我们发现的应用性，即双稳态的Schl?gl模型、一个涉及化学复合物形成的简单自组装方案、一个展示手性对称性破缺的典型模型，以及反应-扩散模式的出现。最后，为了加强我们结果的广泛适用性，我们开发了TACOS（热力学可访问的非平衡稳态化学反应网络浓度），这是一个开源的Python包，用于确定任何CRN的热力学空间（见附录C）。

本文的结构如下：第二节详细介绍了这项工作的主要结果和核心信息。第三节介绍了对通用CRNs进行热力学一致建模的基本概念，以表征它们的性质和守恒定律。掌握了这些理论工具后，第四节我们展示了主要结果的推导，特别强调了“化学探针”的概念——一种在电路中起到类似电压表作用的人工反应。这种探针允许我们评估与守恒定律兼容的任何可能反应中涉及的物种浓度的热力学可访问空间。第五节我们将我们的结果应用于几个代表性的CRNs，以展示它们的含义和广泛适用性。最后，第六节提供了两个例子，其中我们的界限可以从实验数据中提供有用信息。

II. 主要结果
CRN描述了一组相互连接的反应，这些反应将化学物种转化为其他物种。这些反应可以通过外部化学计量器保持非平衡状态，这些装置提供特定化学物种的连续通量，使其浓度保持固定。任何非平衡CRN都由其内部物种的特性、其浓度动态演变、外部（化学计量控制的）物种以及相关的反应速率和一组守恒定律所表征。热力学通过将给定反应的前向和反向通量之比（）与涉及的物种的化学势（即沿的 free-energy change）联系起来，被纳入CRNs的建模中。这种关系的确切形式被称为局部详细平衡（LDB）[第三节A，方程（8）]。因此，每个反应都受到相应的热力学力的驱动，即“亲和力”。这种亲和力来自内部（）和外部（）物种的贡献，其中编码了非平衡驱动力[第三节A，方程（10）]。

这项工作提出了两个主要的热力学结果，这些结果独立于CRNs的结构特性，只要它们达到稳态。因此，我们的框架不适用于持续振荡的情况，但它适用于表现出多稳态的CRNs。这些结果将全局网络拓扑和能量驱动与局部反应性质以及系统的整体可访问浓度空间联系起来。

A. 反应亲和力的热力学界限
第一个主要结果确定，任何反应的稳态亲和力，记为，受到网络全局热力学性质的限制。我们首先将CRN分解为最小的一组基本循环——这些是基本通量模式（EFMs），其确切定义在第三节B中介绍。这些最小循环中的每一个，记为，可能由化学计量控制的物种驱动，并与循环亲和力相关联，代表这些外部物种引起的净自由能变化[第三节B，方程（14）]。我们证明，亲和力受到包含的所有与基本循环相关的亲和力的极值的限制（这个子集由表示）：（1）其中每个亲和力必须通过进入相应循环的有效次数进行归一化。在第四节，我们推导出了这个界限的符号解析形式，区分了是在前向还是后向方向上考虑，并包括了零亲和力的情况[方程（20）]。用更直观的语言来说，任何局部反应步骤的热力学驱动力从根本上受到全局非平衡驱动的限制——在整个CRN的水平上评估——决定了沿基本循环的极端可用亲和力。这些界限同样适用于与守恒定律兼容但不明确存在于CRN中的反应的有效亲和力[第四节，方程（24）。

B. 热力学空间：可访问浓度的界限
第二个主要结果将反应亲和力的界限转化为对化学物种稳态浓度的限制，从而确定了CRN的热力学空间。这是在给定非平衡驱动水平下生化系统必须操作的浓度范围。

考虑内部化学物种的净转化——无论是通过现有的反应还是与守恒定律兼容的有效反应的路径——例如，图1(a)中图解表示的四种不同组分的组装。这种转化可以通过不同的路径发生[图1(a)中的红色和蓝色路径]，其中每一个都由表示，其中表示从反应的左侧到右侧（或从右侧到左侧的特定转化涉及的内部物种集）。与每条路径相关的自由能变化包括了与内部物种转化相关的标准吉布斯自由能变化以及沿该路径的外部化学计量控制物种引起的累积吉布斯自由能变化，这些变化编码了非平衡驱动的效果。这些量的详细信息在第四节，方程（26）和（28）中提供。通过将反应路径的有效平衡常数定义为（减去）其总自由能变化的指数，并用表示，我们可以识别出导致极端贡献的路径（这些分别是图1(a)草图中的红色和蓝色路径）。我们证明，涉及的物种的稳态浓度必须满足以下界限：（2）其中是定义这种特定转化的内部物种变化的净化学计量系数（第三节）。这些限制规定了给定转化的可访问浓度范围由决定极端有效平衡常数的全局非平衡驱动控制（见图1(b)中的绿色区域，由相应路径的极端红色和蓝色贡献界定）。

图1. CRN的热力学空间。 (a) 开放CRN的示意图，其中内部物种通过多个相互连接的路径转化和组装，由外部化学储库驱动的非平衡燃料到废物的转化。突出显示了与最高（红色）和最低（蓝色）自由能变化相关的极端路径。(b) 对数轴上浓度热力学空间的概念性说明。它显示，可访问的稳态浓度受到来自极端路径的贡献的上限和下限限制，这些在(a)中以红色和蓝色显示。这个空间的宽度（垂直箭头）通过非平衡驱动力量化了能量预算。

这些界限可以用于任何与守恒定律兼容的内部物种的转化。为了评估那些不是通过现有反应直接连接的物种的界限，我们引入了“化学探针”的概念——一个辅助的、无限缓慢的反应，用于测量任意两组物种之间的有效亲和力（第四节A）。注意，当动态主要由相应的反应路径支配时（图1(b)中的红色和蓝色曲线），这些界限（上下）可能会在一般非平衡条件下饱和，即当这条路径比连接相同化学物种的所有其他路径都要快得多时。因此，方程（2）包含了优化特定浓度的相关路径的信息。

写下方程（2）可以确定给定全局能量预算下CRN的整个热力学空间。在热力学平衡状态下，所有来自外部物种的驱动力都消失[例如，方程（1）中的]。因此，所有路径上来自外部物种的能量贡献均为零，导致某些值与另一值重合，并且其值仅由内部物种的标准吉布斯自由能变化决定。方程（2）因此成为一个等式，系统稳定在一个由经典热力学和守恒定律决定的独特浓度集合中。相反，热力学空间的非零宽度是非平衡条件的标志，也是复杂行为出现的先决条件。

III. 化学反应网络的热力学与几何学

A. 动力学与LDB
首先，我们总结了用于描述通过化学计量器维持非平衡状态的开放CRN（Chemical Reaction Networks）的动态和热力学的框架的相关方面。我们将化学物种分为两种类型：内部物种，其浓度受网络内反应动态的控制；以及外部（受化学计量器控制）物种，其浓度通过化学计量器的作用保持恒定。用\(i\)表示CRN中某个可能的反应的索引。我们可以用以下形式表示一个通用反应：
\[ \text{Reaction} = a_iX_i + b_jY_j + \ldots \]
其中\(a_i\)和\(b_j\)分别表示内部物种在正向和反向反应中的化学计量系数。对于外部物种，\(a_i\)和\(b_j\)也表示相同的量。正向和反向反应的速率常数分别由\(k_i\)和\(k_j\)给出。如后文所述，为了确保热力学一致性，我们认为所有化学反应都是双向的，即可逆的CRN。网络中所有反应引起的化学计量变化由化学计量矩阵\(\mathbf{M}\)捕捉：
\[ \mathbf{M} = \begin{bmatrix} a_{ij} & b_{ij} \\ a_{ij} & b_{ij} \\ \vdots & \vdots \end{bmatrix} \]
其中\(a_{ij}\)和\(b_{ij}\)分别表示由于第\(i\)个反应导致内部和外部物种的分子数净变化。矩阵\(\mathbf{M}\)的维度为\(m \times n\)。为了通用性，用\(S_i\)表示一个通用化学物种（内部或外部），用\(m_i\)表示由于第\(i\)个反应导致其分子数变化对应的矩阵中的条目。矩阵\(\mathbf{M}\)的第\(i\)列，表示第\(i\)个反应中所有物种的分子数变化。在图2(a)中，我们展示了一个由四个反应组成的简单CRN的例子，涉及五种化学物种，其中三种是内部的，\(X_1, X_2, X_3\)，两种是外部的，\(Y_1, Y_2\)。其化学计量矩阵也在图中给出。

图2. 一个CRN的例子。(a) 一组带有相应速率和化学计量矩阵的化学反应。这些量定义了CRN。\(S_1\)和\(S_2\)表示两种不同形式的单体（例如，非活性和活性），而\(S_3\)是一个二聚体。化学势差\(\Delta\Phi\)驱动\(S_1\)转化为\(S_2\)，使系统脱离平衡状态。(b) 反应网络的超图表示。这里，指向每个化学物种的边数表示其在反应中的化学计量系数。箭头指示每个反应边的方向，由化学计量矩阵中相应条目的符号决定，从具有负条目的反应物指向具有正条目的产物。(c) 从化学计量矩阵得出的守恒定律垂直于内部物种的化学计量子空间。如正文中所述，浓度用小写字母表示。(d) 涉及反应4的两个受驱动的EFMs（Efficient Flux Modes）。EFM代表一个循环反应路径，不会导致内部物种的净变化。

对于稀薄的CRN，浓度演变可以通过质量作用定律来表达。特别地，我们可以引入正向和反向通量向量\(\vec{f}_i\)和\(\vec{f}_j\)，其中第\(i\)个分量量化了给定反应在正向或反向的倾向。这些通量的显式形式遵循质量作用动力学：
\[ \vec{f}_i = k_iS_i - k_jS_j \]
\(\vec{f}_j = k_jS_j - k_iS_i \]
通过进一步引入净通量向量\(\mathbf{F}\)，其元素为\(\vec{f}_i + \vec{f}_j\)，我们可以通过化学计量矩阵表达网络中物种浓度的动态：
\[ \mathbf{C} = \mathbf{M}\vec{F} \]
这里我们使用了简写符号\(\vec{C}_i\)来表示内部物种的浓度向量，其元素\(\vec{c}_i\)，类似地用于\(\vec{F}\)。我们用\(\mathbf{C}\)表示所有物种浓度的完整集合。外部通量\(\vec{F}\)由化学计量器引起，通过补偿内部反应的通量来维持受化学计量器控制物种的恒定浓度。如果系统达到稳态，它将满足条件：
\[ \mathbf{C} = \mathbf{F}_\text{steady} \]
其中\(\mathbf{F}_\text{steady}\)表示稳态通量向量，即所有浓度都被替换为其静止值的通量。然而，重要的是要注意，并非所有CRN都必然达到稳态。一些系统可能会表现出极限环、混沌或其他类型的吸引子。在这项工作中，我们将讨论限制在能够达到稳态解的CRN上。

热力学通过LDB条件[27]纳入，该条件限制了正向通量和反向通量的比率：
\[ \frac{\Delta\Phi_i}{\Delta\Phi_j} = \lambda \]
其中\(\Delta\Phi_i\)是沿反应\(\text{反应}_i\)的自由能变化，\(\lambda\)是反应亲和力。亲和力与净反应通量\(\vec{F}_i\)的方向一致；即，正的\(\lambda\)与正的\(\Delta\Phi_i\)相关联，因为它表明正向通量大于反向通量。一个物种的化学势由下式给出：
\[ \Phi_i = \Phi_0 + \delta_iS_i \]
其中\(\Phi_0\)是标准化学势，\(\delta_i\)是熵贡献。然后我们引入化学势向量\(\mathbf{\Phi}\)，其元素为\(\delta_i\)。通过使用这种简化的符号，可以计算网络中所有反应的亲和力向量：
\[ \mathbf{\Phi} = (\mathbf{M}\mathbf{\delta})^T \]
通过将化学势向量分解为来自内部和外部物种的贡献，亲和力可以表示为：
\[ \mathbf{\Phi} = (\mathbf{M}\mathbf{\delta}_i + \mathbf{\delta}_j) \]
亲和力\(\alpha\)表示与内部物种变化相关的热力学驱动力，而\(\beta\)表示与外部物种变化相关的热力学驱动力。

利用方程(5)中的质量作用动力学，LDB条件可以用反应速率常数来重写：
\[ \frac{\Delta\Phi_i}{\Delta\Phi_j} = k_ik_j \]
这种形式对应于随机热力学[28]中通常提到的LDB，即使方程(8)中的表达式通常更普遍适用。我们还指出，我们的结果仅基于方程(8)，而不依赖于动力学的具体形式。

B. 循环与EFMs
一个通用的CRN可以被视为一个超图，其中每个反应都是连接多个反应物种的超边[30,44]。化学计量矩阵定义了超图的几何结构，我们可以从中提取超循环，即系统保持不变的系列反应。从代数角度来看，如果向量\(\mathbf{v}\)是一个CRN的循环——即它包含不导致分子数净变化的反应集合——那么\(\mathbf{v}\)属于化学计量矩阵的右零空间，也称为矩阵的核\(\mathcal{K}\)。在图2(b)中，我们展示了图2(a)中示例的化学计量矩阵相关的超图。

在一个开放的CRN中，外部物种受到化学计量器的控制，而内部物种受到反应动态的影响。因此，我们对内部反应的化学计量矩阵的循环感兴趣。这些循环，称为“通量模式”，是大小为\(m \times n\)的向量，并满足：
\[ \mathbf{v} \cdot \mathbf{Q} = 0 \]
这个方程由所有仅保持内部物种分子数不变的反应集合解决。所有通量模式构成的空间是内部化学计量矩阵的核\(\mathcal{K}\)。我们可以立即看到，根据方程(7)，CRN的稳态解是一个通量模式。这意味着我们可以将通量模式解释为通量的所有可能稳态解，而实际解取决于速率的具体值。

理解CRN可以操作的所有可能路径对于引入EFMs的概念至关重要。EFMs是支持最小的非零通量模式，在这个意义上，所有可能的CRN中的通量模式都可以表示为EFMs的线性组合。更准确地说，一个EFM是一个大小为\(m \times n\)的向量，其支持，即非零条目的索引集合\(\mathcal{S}\)，在非零向量中是最小的。因此，EFM不能进一步分解为支持更小的通量模式。换句话说，EFMs是保持内部物种分子数的最小循环。所有EFMs的集合由\(\mathcal{EFM}\)表示。我们用\(\mathbf{EFM}_i\)表示EFM的第\(i\)个元素。向量\(\mathbf{EFM}\)的模表示反应\(\text{反应}_i\)发生的次数。然而，\(\mathbf{EFM}\)的符号表示反应是在EFM中的正向（正）还是反向（负）方向进行的。请注意，这个定义本质上依赖于我们选择给定反应的正方向的任意性，我们将其视为正向并与速率\(k_i\)相关联。不同的选择将导致相同的结果。

EFMs与通量模式之间的关系可以通过共形性概念来表征。我们说一个EFM\(\mathbf{EFM}\)与通量模式\(\mathbf{v}\)是共形的，如果对于每个反应\(r_j\)，要么\(\mathbf{EFM}_i \cdot \mathbf{v} \geq 0\)，要么\(\mathbf{EFM}_j \cdot \mathbf{v} \leq 0\)。这个性质确保了EFM和通量模式在共享反应时的对齐。EFMs的一个基本性质是任何通量模式都允许共形分解；即，它可以写成一个不相互抵消的EFMs之和[32,45]，每个EFM都与它共形。这样的分解通常不是唯一的，不同的通量模式可能涉及完整集合\(\mathcal{EFM}\)的不同子集。因此，由于稳态通量\(\mathbf{F}_\text{steady}\)是一个通量模式，它也允许共形分解。在第IV节中，我们利用这个性质来证明我们的主要结果。

在这里，我们还简要评论了EFMs的数值枚举。对于完全不可逆的网络，EFMs与由内部化学计量矩阵定义的非负稳态通量锥的极射线相吻合[46]。在这项工作中采用的可逆约定中，我们引入了一个分割表示，并在分割空间中计算极射线，然后将其投影回原始反应坐标。这些算法及其随网络大小缩放的详细描述在附录B中提供。值得注意的是，虽然对于小型到中型CRN，这些算法可以提供一个实用工具来识别EFMs，但更大的网络将需要更有针对性的计算任务。在下文中，我们将讨论所呈现的结果如何指导这种有针对性的搜索的开发。

C. 守恒定律与反应空间
完全描述一个CRN的另一个关键组成部分是识别守恒定律。它们封装在\(\mathbf{M}\)的左零空间中，称为\(\mathbf{M}\)的余核，或等价地，\(\mathbf{M}\)的核[4,27]。我们用\(\mathcal{L}\)表示所有守恒定律的集合，这在代谢网络的背景下也被称为“代谢池”[47]。因此，对于CRN的任何守恒定律\(L\)，都有\(\mathbf{L}(\mathbf{M})\)成立。在一个开放的CRN中，引入化学计量器会打破一些守恒定律，特别是那些涉及受化学计量器控制的物种的定律[29]。因此，\(\mathbf{M}\)的一个子空间代表了开放系统中的守恒量。这个子空间，用\(\mathcal{K}\)表示，正好是\(\mathbf{M}\)的余核。对于任何\(\mathcal{L}\)，我们可以识别一个相应的守恒量\(L(\mathbf{L}(i)\)：
\[ \mathcal{L}(i) = \mathbf{M}(\mathbf{M}_i) \]
图2(c)报告了图2(a)中示例反应网络相关的内部守恒定律。尽管在存在受化学计量器控制的物种的情况下一些守恒定律被违反，但整个集合\(\mathcal{L}\)仍然是识别给定CRN中允许的反应的宝贵工具。这些是原则上可以发生的反应，因为它们与\(\mathbf{L}\)确定的守恒定律兼容，即使它们在调查的CRN中不可能。类似于EFMs的概念，我们引入了基本守恒定律（ECLs）的集合。每个ECL都有一个最小且不可简约的支持集。因此，任何在给定的开放CRN中不存在的可行反应，都由化学计量矩阵中的列来表示，这些列满足对所有i的约束条件。我们将所有与完整化学计量矩阵相关的守恒定律兼容的CRN中可能发生的反应的集合称为“反应空间”。反应空间将是确定一般CRNs热力学空间的基础点。实际上，引入一个虚构的测试反应可以让我们通过探索在CRN中可能发生但实际不存在的化学转变来限定浓度范围，即那些属于某些条件的反应。CRNs的全局热力学特性主要由LDB条件（方程式（8）来表征。该条件将反应通量的不可逆性与吉布斯自由能的变化联系起来，吉布斯自由能的变化等于（负）反应的亲和力，并且代表了一种局部关系。此外，CRNs还具有与其拓扑结构相关的全局热力学特性，这些特性由EFMs（有效流模式）捕获。对于任何给定的EFM，我们可以通过以下方式评估一个循环亲和力：（14）这是一个考虑了参与该循环的所有反应的净亲和力。注意，方程式（14）对于稳态和非稳态都成立[48,49]。根据定义，沿着EFM进行操作不会改变内部物种的浓度；因此，循环亲和力仅由化学调控物种的浓度变化产生，因为这些变化被外部通量所补偿[见方程式（6）]。因此，循环亲和力等于执行相应反应循环的化学库中吉布斯自由能变化的负值。这是CRN的一个纯粹的结构特性，与内部物种的浓度无关[38]。图2(d)显示了图2(a)中涉及反应4的网络的EFMs。我们还报告了与每个EFM相关的循环亲和力。理解这一特性如何反映CRN的平衡和非平衡状态至关重要。在平衡状态下，所有净通量都消失，即所有亲和力都等于零。这意味着，使用方程式（9）中亲和力向量的定义，平衡条件等价于要求所有循环亲和力为零。封闭的CRNs总是能达到平衡状态。相反，开放的CRNs往往无法达到平衡，因为它们的性质通常与方程式（14）中的循环亲和力条件不相容。事实上，化学调控物种的存在通常使开放的CRNs保持非平衡状态，根据动力学情况，可能会导致出现各种吸引子（例如，固定点、极限环、奇异吸引子）。然而，如果所有化学调控物种都保持在它们的平衡浓度，即对于每个EFM都满足条件，则开放的CRNs会放松到平衡状态[27,50]。

IV. 化学反应网络的通用热力学界限
A. 亲和力的热力学可达区域
沿着EFMs的循环亲和力提供了对CRNs热力学特性的全局视角，并量化了它们与平衡条件的偏离。然而，确定整个网络中所有反应亲和力的稳态值仍然具有挑战性。实际上，亲和力分布源于外部化学调控物种的存在，其确定需要求解具有详细动力学的系统。大多数非平衡CRNs的固有非线性通常排除了解析解的可能性，因此必须使用数值方法来表征非平衡稳态的热力学特性。尽管亲和力分布依赖于动力学细节，但我们在这里展示了如何仅使用CRN的全局热力学和拓扑特性来建立反应亲和力的热力学界限，这些界限是局部量。这种方法为理解非平衡系统的行为提供了宝贵的见解，而不需要完全了解其动力学。此外，它为理解热力学驱动力如何在特定范围内有效增强（或降低）反应亲和力提供了一个有前景的视角，这一效应在各种化学系统中被证明是相关的[14,51,52]。这部分的主要结果总结在第二节A中，下面我们将介绍其推导过程。

1. CRN中反应的推导
为了推导这些通用热力学界限，我们首先假设CRN存在一个或多个具有非零物种浓度的稳态解，这些解由方程式（6）的固定点表示。我们回顾第三节B的内容，由于CRN中的每条边都允许共形分解，因此沿边的每个非零通量——等价于每个反应通量——必须至少包含在一个沿着该边方向相同的EFM中，即一个共形EFM。此外，对于一个给定的反应，通量的分量与亲和力具有相同的符号。因此，共形EFM也与稳态亲和力向量对齐。然后考虑一个被维持在非平衡状态的反应（即，和具有相同的符号）。在任何共形分解中，至少会出现一个包含该反应的EFM。设这样的EFM为，因此和。我们可以重写方程式（14）来单独得到该反应的亲和力：（15）由于与共形，求和中的所有项都是正的，所以（16）这个界限依赖于在稳态通量的共形分解中识别一个EFM，而这又需要知道和，因此，需要知道整个动力学。为了进一步推导出纯热力学的界限，我们最大化所有沿着具有与相同符号的分量的EFM，即所有这样的EFM。因此，（17）由于的符号取决于未知的稳态通量，我们区分了两种可能的情况，和。通过定义集合为所有具有分量的EFM，我们有（18）（19），这里我们利用了当时，方程式（17）最右边的值变为，利用了的负性。由于和代表具有相反全局方向的相同EFM，并且比率在任何符号反转下都是不变的，因此在取EFM的极值时可以使用任一方向。将方程式（18）和（19）结合起来，对于任意的，我们最终得到（20），其中极值是在包含的所有EFM上取得的。区分和只是为了跟踪的分量的两种可能符号，而比率在EFM的全局反转下是不变的。因此，最终界限不需要知道稳态通量的符号或识别哪些EFM与其共形，从而消除了指定整个动力学的必要性。在这里，我们在不等式两边都包括了等号，因为注意到在平衡状态下应该成立等式。因此，尽管推导是针对非平衡系统进行的，但它也可以很容易地扩展到平衡系统。实际上，在没有由化学调控物种引起的非平衡驱动的系统中，所有循环亲和力都消失（对于所有循环）。这个条件确保了在稳态下，所有反应亲和力也都消失；因此，稳态是一个平衡状态，等式成立。然而，在存在非平衡驱动的情况下，方程式（20）揭示了每个反应亲和量偏离平衡的热力学极限，这完全基于全局的非平衡特性。我们在这里指出，当非平衡驱动导致净通量为零时，方程式（20）中的不等式也会变成等式，从而亲和力也为零。系统实际上表现为一个平衡系统，等式如上所述得以恢复。此外，当不参与任何循环的反应时，方程式（20）中的等式也成立，因为它满足详细平衡，因此表现出零亲和力。一个重要的观察是，这种公式化将反应亲和力这种局部量与不需要知道与它共形的EFMs的全局热力学特性联系起来。这个结果适用于CRN中的任何反应。

此外，一个反应的亲和力可以分解为来自内部和外部物种的贡献，即[见方程式（10）]。将移到方程式（20）的右边，我们得到了内部物种亲和力的一个上限：（21）这里，我们将定义为一组反应路径，这些路径来自包含的EFMs，它们将内部物种从反应的右边转换到左边，如负号所示。更准确地说，中的每个元素由表示。这个集合包含了反应本身及其亲和力（取负号），以及所有由消除了反应本身后的EFMs构建的路径。实际上，方程式（21）中的最后一个等式可以通过注意到当时，方程式（17）的最右边变得等于来容易理解，这里利用了的负性。由于和代表具有相反全局方向的相同EFM，并且比率在任何符号反转下都是不变的，因此在取EFM的极值时可以使用任一方向。将方程式（18）和（19）结合起来，对于任意的，我们最终得到（20），其中极值是在包含的所有EFM上取得的。区分和只是为了跟踪的分量的两种可能符号，而比率在EFM的全局反转下是不变的。因此，最终的界限不需要知道稳态通量的符号或识别哪些EFM与之共形，从而无需指定整个动力学。这里，我们在两个不等式上都包括了等号，因为注意到在平衡状态下应该成立等式。因此，尽管推导是为非平衡系统进行的，但它也可以很容易地扩展到平衡系统。实际上，在没有由化学调控物种引起的非平衡驱动的系统中，所有循环亲和力都消失（对于所有循环）。这个条件确保了在稳态下，所有反应亲和量也消失；因此，稳态是一个平衡状态，等式成立。然而，在存在非平衡驱动的情况下，方程式（20）揭示了每个反应亲和量偏离平衡的热力学极限，这完全基于全局的非平衡特性。我们在这里指出，当非平衡驱动导致净通量为零时，方程式（20）中的不等式也会变成等式，因此亲和力也为零。系统有效地表现为一个平衡系统，等式如上所述得以恢复。此外，当是一个不参与任何循环的反应时，方程式（20）中的等式也成立，因为它满足详细平衡，因此表现出零亲和力。一个关键的观察是，这种公式化将具有全局热力学特性的局部量作为反应亲和量进行了界定，而不需要知道与它共形的EFMs。这个结果适用于CRN中的任何反应。

此外，一个反应的亲和力可以分解为来自内部和外部物种的贡献，即[见方程式（10）]。将移到方程式（20）的右边，我们得到了内部物种亲和力的一个上限：（21）这里，我们将定义为一组反应路径，这些路径来自包含的EFMs，它们将内部物种从反应的右边转换到左边，如负号所示。更准确地说，中的每个元素由表示。这个集合包含了反应本身及其亲和力（取负号），以及所有从消除了反应本身的EFMs构建的路径。实际上，方程式（21）中的最后一个等式可以通过注意到所有使得的EFMs都将从的右边转换到左边的来容易理解。从中构建的路径的亲和力是，这里我们除以，因为在一些EFMs中可能会有更多的通过的转换。在消除之前，我们称之为的路径集合在反向中包含了本身。从这些观察中，方程式（21）直接得出。通过采用类似的方法，通过对中的所有路径的最大亲和力的负值进行取值，我们可以得到的一个下限。

2. 对于有效反应的推导
尽管方程式（20）和（21）具有普遍性，但它们仅适用于CRN中存在的反应亲和力。然而，可以克服这一限制，为原本不在CRN中的反应的有效亲和力推导出上下热力学界限，而这些反应原则上与CRN兼容。我们引入了“化学探针”的概念，即一个属于反应空间但不属于原始网络的测试反应。由此产生的扩展化学计量矩阵是，其中新引入的反应由列表示。我们假设进行得如此缓慢，以至于其净通量可以忽略不计，即。然而，它仍然保持一个有限的亲和力，这取决于参与反应的物种的化学势。为了确定与这个测试反应相关的亲和力，我们考虑一个包含的共形EFM。根据之前的关系，我们有（22）由于由化学探针引起的反应属于原始CRN的反应空间，添加它至少闭合了一个涉及探针和原始网络中路径的循环，其中是沿着探针反应的单位向量。此外，因为探针是以无限小的通量引入的，所以原始反应的稳态亲和力保持不变。除了包含并且与和对齐的EFM之外，因此，可以通过向包含探针的循环中添加原始CRN中的一个通量模式来构建第二个通量模式，该模式通过相同的路径但方向相反，并由不同的路径闭合。这个操作保持探针的方向不变，取消了的贡献，并产生了一个其包含的原始反应上的分量与对应的稳态亲和力符号相反的通量模式。然后，这个通量模式的共形分解包含了一个具有相同符号结构的EFM，我们将其表示为。因此，我们有和。在图3中，我们用一个简单的CRN来说明这个想法，以建立直观的理解。我们还指出，在原始CRN中参与循环的所有边都存在反对齐。在所有其他情况下，不属于任何循环的边可以忽略，因为它们不会对循环亲和力做出贡献。通过将这些考虑纳入方程式（15），我们得到（23），通过使用所有反应中的反对齐。我们可以任意将的符号设为正，根据构造，的符号也将是正的。因此，通过将方程式（22）和（23）结合起来，我们有（24），这里我们使用了之前描述的平衡条件，即和，因为亲和量和都改变了符号。很容易检查，如果，这个三明治不等式仍然成立。至关重要的是，方程（24）相对于方程（20）提供了一个更严格的界限，因为我们使用了描述化学探针的额外信息。这个夹心不等式适用于任何属于反应空间的反应，并且在给定决定非平衡条件的能量预算的情况下，确定了任何有效亲和力的热力学极限。因此，方程（20）和（24）构成了这项研究的第一个主要结果。

图3. 化学探针的示例。(a) 显示了一个线性CRN，所有反应都用黑色（和灰色）线条表示，探针用虚线红色线条表示。箭头表示稳态反应亲和力。与探针相关的稳态通量消失，但其亲和力仍然有限。用绿色显示了与修改后的网络（即原始CRN加上探针）的亲和力相符的EFM。我们突出显示了参与此EFM的节点和边。(b) 修改后的网络的EFM与探针的亲和力对齐，但与原始CRN的稳态亲和力反向对齐。通过将亲和力分解为内部和外部物种的贡献，我们可以将上界重写为：（25）对于下界也是如此。这个表达式与原始CRN中不存在的反应的亲和力相比有一个主要区别。在第一种情况下，由直接反应给出的项不会出现，因此必须考虑将反应从右向左转换的路径。这一性质与这里引入的符号完全一致，也符合这种化学探针仅仅是被用来研究系统的数学工具，不包含任何物理信息的事实。最后需要指出的是，这样的化学探针可以类比为电路中的电压表。两者都以最小的干扰来测量电位差：探针逐渐引入，导致网络的边缘电流无限缓慢，而理想的电压表具有无限大的电阻，导致流量为零。因此，探针利用其消失的化学通量来测量网络节点之间的化学势差，类似于电压表测量电压的方式。正如我们在下一节中看到的，这种方法对于确定CRN的整个热力学空间至关重要。

B. CRNs的热力学空间：物种浓度的可访问区域
从反应亲和力的界限中，我们可以推导出物种浓度的上界和下界。这些限制决定了在给定CRN的能量预算下每种化学物种的可访问区域。我们将这个区域称为“热力学空间”。这个概念已经在参考文献[22]中提出，适用于线性和催化CRNs。此外，在各种情况下，如代谢网络[47]、信号放大[36]和自组装[40]中，给定外部非平衡约束的热力学可访问空间的概念也具有根本重要性。这部分的结果在第II B节中进行了总结，而其推导过程在下面给出。

1. 热力学空间的推导
尽管稳态浓度不仅仅由热力学决定，因为它们通常取决于CRN的动力学性质，但它们通过化学势进入热力学描述。我们回顾一下，亲和力也可以用化学势来表示。通过指定方程（9）对于，我们有（26），其中表示沿的反应的标准吉布斯自由能变化，表示沿反应路径的吉布斯自由能变化。通过结合方程（21）和（26），并利用亲和力的对偶性，即对于任何，我们有（27），我们得到了参与反应的内部浓度的界限：（28）。这里，是将化学物种从反应的左侧转换到右侧的反应路径的有效平衡常数。实际上，如果CRN只有一个路径，例如，系统将达到平衡，涉及物种的浓度乘积与化学计量指数将恰好等于。如方程（28）所示，浓度被所有这些路径上最小化的平衡常数所下限。按照同样的程序，我们还可以推导出一个上界，得到所有参与的反应的物种的可访问浓度范围：（29）。方程（29）中的界限是为通过反应连接的物种的浓度推导出来的。然而，从方程（25）开始，我们可以看到，对于任何通过想象的与守恒定律兼容的反应连接的物种，也有类似的界限；即，它属于CRN的反应空间。这些界限与方程（29）相同，只需将替换为，并且路径必须从原始CRN构建（因为化学探针从未存在），这与系统的物理内容一致。因此，方程（29）构成了这项研究的第二个主要结果。

从计算角度来看，确定这些热力学界限并不一定需要完整的EFM枚举，如上所述并在附录B中展示的，这对于小型到中等规模的CRN是可行的。实际上，只需要识别涉及驱动反应的路径，因为需要确定方程（29）中热力学贡献极值的路径来确定热力学空间。因此，目前的公式自然指向一个包含预定反应的定向EFM搜索问题。最近已经提出了这方面的现有算法[53]，而开发 refined 的定向方法以减轻与完整EFM确定相关的组合爆炸仍然是研究的活跃领域[54,55]，超出了本工作的范围。

我们注意到，方程（29）中的每个界限都可能在由单一路径动力学主导的CRN中达到非平衡条件下的饱和，这些条件下的通量消失。这样的路径决定了自由能差异的极端贡献，即导致方程（28）中的下界。饱和条件与参考文献[22,56]中讨论的线性系统中的条件类似。

让我们研究平衡状态和非平衡状态之间的差异，以及热力学空间可以提供的关于系统的信息。在热力学平衡状态下，所有驱动力都消失，因此平衡常数对所有路径都相同，等于内部物种的标准吉布斯自由能差异，正如预期的那样。因此，方程（29）两边变为等式；因此，所有浓度的稳态值趋同于由它们的自由能和CRN的守恒定律决定的唯一值。相比之下，非平衡状态的CRNs表现出一系列可行的稳态浓度。这构成了物种浓度整个相空间的一个子空间，并且它取决于网络的全局非平衡热力学性质。一般来说，非零的热力学空间使得复杂动态成为可能，并且是复杂行为出现的先决条件，如自由能转换、双稳态、模式形成和运动响应状态。热力学空间，编码在方程（29）的界限中，量化了非平衡条件与任何CRN中出现的复杂性之间的密切关系，因此提供了复杂生化过程潜在功能多样性和效率的度量。

C. 测试候选反应拓扑的热力学空间
上述推导出的界限需要指定CRN的拓扑结构。然而，它们可以被转化为一个用来询问反应网络本身的工具。因为由方程（29）定义的热力学空间包含了任何与给定拓扑兼容的动力学可达到的每一个稳态，所以超出其限制的测量浓度无法通过调整速率常数来解释。这样的结果将 falsify 与所提出的拓扑一致的所有动力学类别。因此，热力学空间是一种无需动力学的候选反应集的一致性测试，正如我们现在对超图展示的，违规甚至可以指出缺失的具体反应。

考虑图2(b)中的网络[也在图4(a)中报告，为了完整性并表示为]，一个二聚体可以通过两个活性单体（）和一个非活性单体（）或两个非活性单体（）组装而成。二聚体放大，其中，其上限由包含燃料驱动反应的基本通量模式决定。特别是模式，在组装之前激活两个单体，因此经历了两次驱动反应，从而将上限设置为（30）。在图4(b)中，上限及其对应的模式用红色显示。

图4. 使用热力学空间测试反应拓扑。(a) 图2中网络的两个候选拓扑：包含反应1、2、3和4，而移除了反应2。（b）相对于其平衡值，二聚体放大的可访问范围，作为驱动的函数。平衡值（绿线）是共享的下限；模式只经历一次驱动反应，并将上限设置为斜率1。这定义了热力学空间的上限（绿色阴影区域）；相反，只有在反应2存在时才可访问的模式，使用了两次驱动反应，并将上限提高到斜率2。这设置了热力学空间的上限（绿色阴影区域）；两个上限之间的区域（红色阴影区域）是反应2的推断窗口。点是通过随机动力学速率和总浓度的模拟得出的稳态状态：（灰色的）保持在下限斜率1以下，而（蓝色的）占据了窗口，因此，在推断窗口中的测量值不可能来自，并揭示了反应的存在。

现在想象移除反应2。相应的网络显示在图4(a)中，并表示为。这种消除将破坏反应模式，因此，剩余的模式将设置二聚体放大的上限。然而，在中，驱动反应只出现一次，从而将上限降低到（31）。在图4(b)中，这个上限及其对应的模式用黄色显示。

因此，在这个简单的例子中，添加一个额外的超边使得可实现的放大倍数加倍，导致了一个推断窗口，如果没有活性单体组装反应，这个窗口是无法捕捉到的，并揭示了这种隐藏反应的存在，这与速率常数不符。这种论证主要作为候选拓扑的 falsification 测试，并可以扩展到所有可以通过添加特定超边多次使用驱动反应的情况。显然，在线性网络中，驱动反应的重复使用是没有类似情况的，这是超图结构的一个真实特征。这种推断方法是一面的且是健壮的：超过上限会 falsify 拓扑，而未重复使用驱动反应的缺失候选反应则保持界限不变。因此，我们的方法原则上能够识别可以发现的缺失路径类别，从而补充了第VI节中讨论的数据驱动的驱动力推断。

D. 与先前结果的联系
在这项工作中，我们得出了两个主要结果：稳态下反应亲和力的热力学约束，即方程（20）和（24）；以及物种浓度的界限，即（伪）平衡常数，即方程（29），这些界限定义了CRN的热力学空间。类似的结果之前已经为线性和催化CRNs等特殊情况进行过研究。因此，我们的结果将这些发现推广并扩展到了通用开放化学系统的具有挑战性的情况，引入了热力学空间的概念，并提供了一种方法来为任何CRN找到它。在线性和催化CRNs中，化学计量矩阵的每一列的条目被限制为只有内部物种的一个和，有效地使CRN成为一个图而不是超图。对于这样的图网络，适用一个简单的守恒定律。此外，矩阵树定理允许通过生成树表示稳态解，从而便于推导出稳态可观测量上的热力学界限。先前的结果主要在随机动力学的背景下讨论，用概率比率表示界限。这些概率可以解释为相对于守恒定律的内部物种的标准化浓度。这种表述方式允许系统在确定性和随机性描述之间建立直接联系：（32）这可以直接从方程（29）得出，通过选择一种化学探针，使得和。上述结果等价的形式最早由Maes和Neto?ny在研究线性网络时发现[33]，他们利用最小和最大耗散路径的热耗散来表达上下界，即将（伪）平衡常数重新表述为自由能的变化。Lin通过将电路与线性CRN进行类比，也得出了类似的结果[57]。随后，我们之前的工作证明这些界限也适用于非线性催化网络[22]，在这些网络中，催化效应不会改变化学计量矩阵的元素，并保持CRN的有效图结构而非超图结构。此外，Liang和Pigolotti采用基于轨迹的方法证明了这些界限也适用于线性网络，并将其分析扩展到限制稳态通量不对称性[58]：（33）这个结果可以看作是方程（20）在线性CRN中的一个特例，它利用LDB条件（方程（8）表达了正向和反向半通量的比率。最终，我们指出这里展示的结果包含了所有之前的工作，将他们的发现（关于亲和力和物种浓度）推广到了一般的CRN。这种扩展构成了对化学系统操作范围的首次数学推导，具有至关重要的意义。实际上，许多中心生物学过程，从代谢和信号传导到自我组装[3,4,40]，以及像多稳定性和模式形成这样的复杂现象[8,10,22]，都从根本上依赖于这类通用网络中的多分子相互作用，因此需要一个超越简单模型的框架。因此，我们的结果可以揭示各种不同系统的普遍特性。它们可以揭示由于非平衡条件而产生的内在潜力和限制，并在从校对到生命起源等多个领域提供普遍的答案，而迄今为止只获得了特定系统的解决方案。

V. 应用
A. Schl?gl模型与双稳态
非平衡状态下的开放CRN可能会表现出复杂的涌现现象，例如多稳态，其中反应动力学允许多个稳定固定点。这种非平衡多稳态对于各种生物学功能至关重要，包括基因调控、信号转导和细胞命运决策[59–61]。然而，理解热动力学驱动施加的定量约束仍然是一个具有挑战性的问题。Schl?gl模型是一个用于多稳态的最简模型，它可以在特定的参数范围内维持两个稳定的静止状态[8,62]。分析这样一个典型模型中的这一特征对于阐明非平衡热力学驱动力在限制和塑造化学反应系统涌现出的双稳态行为中的作用至关重要。
Schl?gl模型由以下反应和化学计量矩阵定义：（34），其中是一个内部物种，和是两个被化学计量控制的外部物种。这些反应可以用图5(a)中显示的反应超图来表示。该网络包含一个EFM（反应流图），如图5(b)所示。沿着这个EFM，一个外部物种被转化为另一个外部物种，伴随着相应的自由能变化。因此，其循环亲和力为。应用方程（20），我们推导出内部物种两种反应的亲和力界限：（35）为了限制内部物种的浓度，我们将方程（29）应用于第一个反应中涉及的所有物种。在这种特殊情况下，通过考虑第二个反应也会得到相同的结果，因为我们可以将内部物种视为。实际上，反应亲和力的界限也是一致的[见方程（35）]。我们首先注意到，通过两条路径可以实现内部物种的净生产：（第一个反应）和（第二个反应的逆反应）。然后我们确定了考虑这些路径时的最小和最大平衡常数，并得到内部物种稳态浓度的以下界限：（36）稳态浓度的范围限制了反应动力学的固定点。在图5(c)中，我们展示了这个区域包含了稳定和不稳定的固定点。此外，在图5(d)中，我们针对不同的循环亲和力值限制了整个分岔图。方程（36）的上下界在图5中分别用蓝线和红线表示，在时间尺度分离极限下可以渐近地接近。实际上，当第一个反应进行的速度远快于第二个反应时，系统达到了上界，从而达到了准平衡状态。相反，当第二个反应显著更快，主导系统动力学时，就接近了下界。这个时间尺度分离原则揭示了动力学参数如何将系统推向其热力学允许范围内的不同极端，突出了热力学约束和动力学因素在确定稳态浓度方面的相互作用。
图5. Schl?gl模型的热力学空间。(a) Schl?gl模型的反应超图，包含一个内部物种和两个被化学计量控制的外部物种。(b) Schl?gl模型的EFM，其中内部物种从环境中产生（用表示），然后通过反应1和2降解回环境中。(c) 模型的零线。固定点（导数的零点）被限制在热力学空间内（绿色阴影区域）。(d) 以循环亲和力为控制参数的分岔图。稳定（实线）和不稳定（虚线）解被热力学空间限制（绿色阴影区域）。(c)和(d)中的蓝线和红线分别代表从方程（36）导出的热力学空间的下界和上界。图(c)和(d)中的数值结果是用，和获得的，而驱动力通过调整被化学计量控制的外部浓度来变化。(c)中的。

这种反应网络的复杂双稳态行为源于第一个反应的高阶自催化性质。这种自催化作用基于物种的浓度调节反应速率，导致两个稳定固定点接近由方程（36）定义的热力学空间的上下界。这些界限有效地限制了稳定固定点的范围。在平衡系统中，当时，上下界重合，热力学空间萎缩，排除了双稳态解的可能性。这种萎缩阐明了非平衡驱动力对于CRN中复杂行为出现的重要性，即使对于简单模型也是如此。尽管我们提出了一个用于多稳态的最简模型，我们的热力学约束适用于任意CRN，并且可以潜在地扩展到细胞系统中的复杂生化过程，揭示了基因调控或其他生物功能中双稳态模式的起源，其鲁棒性在于多稳态的存在。此外，由于我们的结果不需要任何详细的动力学信息，它提供了一个工具，用于理解在物理限制内化学和生物系统中复杂行为出现的热力学约束。

B. 手性对称性破缺
手性对称性破缺（CSB）对于理解生命的起源至关重要，其中生物分子的纯手性被认为是由自发的对称性破缺机制产生的。虽然动态模型依赖于自催化和循环来解释CSB的出现[63–65]，但它们在区分能量相同的对映体时严重依赖于非平衡条件[66]。尽管已经认识到耗散的必要性，但对于非平衡驱动强度与CSB程度之间的定量关系仍不清楚。在这里，我们应用我们的热力学框架来阐明这种联系。
我们分析了一个热力学上一致的、可逆的Frank模型[65]版本，这是一个已知会表现出CSB的最简系统。该模型描述了两种手性异构体，和，它们从一个非手性前体自催化产生，随后通过二聚化形成一个复合体[见图6(a)]。所有反应都被认为是可逆的，并遵守局部详细平衡[见方程（8）]。物种和被化学计量控制，它们的化学势建立了非平衡驱动。反应和化学计量矩阵由（37）给出。系统的EFM是，对应于净转化。相关的循环亲和力量化了热力学驱动力。
图6. 手性对称性破缺的热力学分析。(a) 由手性物种和，非手性前体以及二聚体组成的热力学一致的Frank模型。和（绿色方块）被化学计量控制。虚线红色边缘表示化学探针。(b) 关键的EFM，涉及到探针的向量表示。(c) 从(b)中的EFM派生出的分子转化路径，其净效应是将物种从探针的左侧转化为右侧。(d)，(e) 在浓度平面上的相图，显示了零线（橙色；，蓝色），固定点（稳定，•；不稳定，），以及系统轨迹（流线）。绿色阴影区域是由方程（29）定义的热力学可达空间。(d)中的蓝线和红线分别代表热力学空间的下界和上界。在(c)和(d)中。图(c)和(d)中的数值结果是用，，和获得的，而驱动力是通过调整被化学计量控制的外部浓度来变化的。在(c)中，。
这种反应网络的复杂双稳态行为源于第一个反应的高阶自催化性质。这种自催化作用基于物种的浓度调节反应速率，导致两个稳定固定点接近由方程（36）定义的热力学空间的上下界。这些界限有效地限制了稳定固定点的范围。在平衡系统中，当时，上下界重合，热力学空间萎缩，排除了双稳态解的可能性。这种萎缩阐明了非平衡驱动力对于CRN中复杂行为出现的必要性，即使对于简单模型也是如此。
虽然我们提出了一个用于多稳态的最简模型，但我们的热力学约束适用于任意CRN，并且可以潜在地扩展到细胞系统中的复杂生化过程，揭示了基因调控或其他生物功能中双稳态模式的起源，其鲁棒性在于多稳态的存在。此外，由于我们的结果不需要任何详细的动力学信息，它提供了一个工具，用于理解在物理限制内化学和生物系统中复杂行为出现的热力学约束。

B. 手性对称性破缺
手性对称性破缺（CSB）对于理解生命的起源至关重要，其中生物分子的纯手性被认为是由自发的对称性破缺机制产生的。虽然动态模型依赖于自催化和循环来解释CSB的出现[63–65]，但它们在区分能量相同的对映体时严重依赖于非平衡条件[66]。尽管已经认识到耗散的必要性，但对于非平衡驱动的强度与CSB程度之间的定量理解仍然不清楚。在这里，我们应用我们的热力学框架来阐明这种联系。
我们分析了一个热力学上一致、可逆的Frank模型[65]版本，这是一个已知会表现出CSB的最简系统。该模型描述了两种手性异构体，和，它们从一个非手性前体自催化产生，随后通过二聚化形成一个复合体[见图6(a)]。所有反应都被认为是可逆的，并遵循局部详细平衡[见方程（8）]。物种和被化学计量控制，它们的化学势建立了非平衡驱动。反应和化学计量矩阵由（37）给出。系统的EFM是，对应于净转化。相关的循环亲和力量化了热力学驱动力。在平衡状态下，高能量或动力学上稳定的结构虽然通常是人们所期望的，但由于系统倾向于最小化自由能，这些结构的出现概率可能会很低。将组装过程推向远离平衡的状态可以显著提高这些状态的产率和特异性[67,68]。然而，预测这种非平衡条件所施加的性能限制仍然是一个重大挑战。需要一个严谨的框架来定量地将能量输入与给定目标组装的可实现丰度和复杂性联系起来。在这个应用中，我们将在最小自组装模型上运用我们的热力学空间形式主义。我们展示了它定义复杂结构放大倍数和可达浓度严格界限的能力。系统如图7(a)所示，包含三种内部物种：单体、二聚体和三聚体，以及三种反应：一种将燃料转化为废物的反应驱动二聚体的形成，二聚体和单体之间的自发结合形成三聚体，以及三聚体自发分解为三个单体。这里，燃料和废物被视为化学计量的外部物种，使系统保持非平衡状态。详细的反应方案在方程(41)中给出，同时给出了化学计量矩阵，而相应的反应超图在图7(a)中用黑色超边表示：(41)。对于这种自组装CRN，在反应超图中有一个EFM（能量流动机制），其循环亲和力源于燃料到废物的转化驱动的二聚体形成。我们进一步假设在模型中这种化学势差是严格正的。应用我们的框架，根据方程(20)，每个反应在稳态下的亲和力被限制如下：(42)。使用方程(29)，我们还可以限制参与这三种反应的内部物种的浓度比：(43)，其中我们使用、、和来分别表示内部物种的标准化学势。浓度限制定义了图7(d)中所示的热力学空间。

图7. 自组装的热力学空间。(a) 自组装模型的反应超图，显示了三种内部物种（单体、二聚体和三聚体），两种外部物种（和），以及一个化学探针（虚线红色边）。(b) 包括化学探针的三个EFM，用向量表示。(c) 通过原始网络对应于(b)中三个EFM的反应序列，将探针反应两侧的物种进行转换。(d) 热力学空间的图示（绿色）。它与守恒律平面（黄色）相交，定义了内部物种可能的稳态浓度。(e) 从方程(43)和(45)计算出的热力学空间（绿色阴影区域），其中绿色边界代表平衡比率，而虚线红色边界确立了非平衡放大的限制。点表示使用随机动力学速率和固定热力学驱动以及标准化学势下的稳态浓度。

此外，为了确定内部物种浓度的可访问区域，我们还必须考虑现有的守恒定律[见图7(d)]。该模型有一个由右零空间派生的ECL（能量转换限制），由给出。这意味着单体单元总数的守恒。这个守恒定律允许一个额外的有效反应：(44)，该反应在原始CRN中并没有明确给出，但它属于反应空间。因此，我们可以通过化学探针引入这个反应[图7(a)和7(b)]。这相当于添加了一个新的化学计量列，根据定义，它与守恒定律正交，即。扩展的化学计量矩阵展示了三个EFM，它们包含了由化学探针给出的额外反应：如图7(b)所示。由此，我们可以使用方程(29)找到化学探针两侧的物种浓度界限。首先，我们找到将二聚体转化为三聚体的路径，然后评估它们相应的最大和最小（伪）平衡常数：(45)。类似地，我们可以使用方程(24)来限制这个额外反应的有效亲和力：(46)。

图7(e) 明确展示了各种物种浓度组合下的热力学空间，如方程(43)和(45)所示。在这个例子中，边界的一侧代表通过标准化学势差确定的平衡值，而另一侧则规定了给定驱动下的非平衡放大的限制。散点表示在固定驱动下，使用随机动力学参数和总浓度进行数值模拟得到的稳态浓度。所有点都落在预测的界限内，验证了我们的理论预测。尽管这个例子使用了一个最小模型，但它清楚地展示了一个关键原理：非平衡驱动（通过）与网络拓扑（反映在EFMs和守恒定律中）之间的相互作用，为自组装系统中的稳态浓度定义了一个明确的可访问热力学空间。这一框架进一步提供了关于复杂结构（例如这个例子中的三聚体）的丰度可以在其平衡值之外被放大的程度的严格界限。

我们的方法可以直接应用于涉及多个步骤和中间状态的更复杂的自组装模型，从纯粹的拓扑和热力学的角度提供了一个关于耗散性自组装能力的通用界限。这为理解自然自组装系统的结构原理提供了有力工具，并指导合成自组装系统的设计，特别关注于预测如何利用能量输入来实现特定的结构或功能结果。

反应-扩散模式是耗散性自组织的典型例子，量化它们的热力学成本是非平衡热力学的核心挑战[20,21,69]。在这个最终应用中，我们将我们的热力学空间框架从混合良好的系统扩展到这些空间扩展的现象。虽然我们之前的工作为形成模式的催化网络提供了界限[22]，但现在我们将这种分析推广到任意反应拓扑。这使我们能够建立关于模式对比度的通用热力学约束，揭示了网络化学计量结构如何从根本上限制宏观空间组织。我们研究了一个通过最小两物种反应网络展示非平衡模式形成的反应-扩散系统。该系统具有两种相互竞争的转化途径：(47)第一种反应代表内部物种之间的自发四阶转化。第二种描述了一个燃料驱动的催化过程，其中的分子作为催化剂，促进额外的分子转化为，并且反之亦然。这种催化途径通过化学计量的燃料物种和废物物种之间的转换使系统保持非平衡状态。

我们关注一个没有边界流的质量守恒反应-扩散系统。通过对方程(48)中的稳态反应-扩散方程进行适当的化学计量线性组合，反应项相互抵消，得到(50)。由于是一个满足连通域上均匀诺伊曼边界条件的谐波函数，它必须是常数。因此，系统遵守以下守恒定律，该定律限制了空间分布：(51)，这在浓度空间中定义了一个流量平衡子空间。对于一维系统，方程(51)可以通过直接积分得到。如图8(a)所示，并遵循相空间几何方法[23]，这个子空间与反应零线（reactive nullcline）的交集决定了稳态模式探索的极端浓度。我们的框架揭示了反应零线，即与CRN在方程(47)中相关的反应项消失的浓度集，从根本上受到以下热力学界限的约束：(52)。通过识别流量平衡子空间[方程(51)]与热力学空间[方程(52)]的边界之间的交点，我们可以找到两种内部物种和的稳态模式的可达浓度界限。在图8(a)中，它们分别用和表示物种。然而，这些界限仍然依赖于扩散系数，因为流量平衡子空间的斜率由决定。为了获得与扩散无关的界限，我们进一步在流量平衡子空间的所有可能斜率上最大化可达范围。由于这个斜率是非正的并且在极限时趋近于零，因此在这个极限下获得了最宽的间隔。流量平衡子空间变成了一条水平线，稳态浓度的界限取最大值。特别是，从方程(52)我们得到(53)。

因此，在这个极限下，模式对比度满足(54)。

图8. 反应-扩散模式的热力学约束。(a) 浓度的相空间显示反应零线嵌入在热力学空间中，其边界由和决定[见方程(52)]。热力学空间与流量平衡子空间的交集给出了稳态模式的浓度范围。这个说明性面板使用了和，通过将设置在方程(49)中获得反应零线。(b) 模式对比度与循环亲和度。理论上的最大对比度是黑色曲线，而数值模拟（圆圈）得到的对比度是针对不同的扩散系数比率。理论线接近于。(c) 在固定扩散系数比率下，对于三个循环亲和度得到的稳态浓度剖面，对应于(b)中的方点，展示了模式幅度随着驱动力的增加而增加，始终受到热力学约束的限制。(d) 在固定循环亲和度下，对于三个扩散系数比率[图(b)中的五边形]得到的稳态浓度剖面。在(b)–(d)的数值模拟中，我们设置了，，和，同时变化和。束缚的数值验证显示在图8(b)–8(d)中，我们取和。与理论界限的偏差与两种物种的扩散系数比率有关——当的扩散远快于时，几乎所有的非平衡驱动都分配在物种上，方程(54)中的界限在极限下饱和。

在本节中，我们将热力学空间框架应用于两个实验数据：在第一个应用中，我们测试了我们的界限在描述蛋白质构象变化的具体CRN中的有效性；在第二个例子中，我们展示了如何直接从实验数据中提取有用的热力学信息。

A. 通过伴侣蛋白稳定天然蛋白质
蛋白质在其天然构象下才能发挥功能[70]。然而，由于冲击或环境压力的影响，它们可能会发生变性，落入自由能景观的深层局部最小值，从中很难自发逃逸，并最终形成稳定的聚集体，如β-淀粉样蛋白和α-突触核蛋白纤维[70–72]。由于错误折叠的蛋白质和聚集构象被确定为几种神经退行性疾病的主要原因，它们最近引起了大量关注[73,74]。在这种情况下，伴侣蛋白起着主导作用，能够利用ATP的能量重新折叠错误折叠的蛋白质，从而增加天然状态的非平衡稳定性[71]。

在图9(a)中，我们展示了GroEL（伴侣蛋白）对MDH蛋白（苹果酸脱氢酶）的作用的CRN（计算网络），忽略了伴侣蛋白存在的情况下的聚集[71]。通过从参考文献[71]中提取实验数据和参数值，我们计算了不同浓度对的热力学界限，[ND]（天然二聚体）与（天然单体），[ND]与（中间体），以及[ND]与（错误折叠的蛋白质）。在图9(b)和9(c)中，我们展示了在不同总蛋白质浓度下但这些ATP量相同的情况下的这些比率；也就是说，我们改变了守恒定律施加的约束，同时保持非平衡驱动不变。总蛋白质浓度来自参考文献[71]，并在粗粒度上考虑了在伴侣蛋白引入之前由于聚集导致的可溶性蛋白质的损失。然后我们添加了固定量的伴侣蛋白。我们展示了和无法从它们的平衡值被置换[图9(b)]，而天然二聚体与错误折叠蛋白质的比率远高于平衡状态[图9(c)]，所有数据都落在热力学空间预测的范围内。这个结果表明，ATP水解的唯一作用是通过伴侣蛋白机制有效增加错误折叠蛋白质的自由能，其增加的量受到能量预算的限制。这一观察结果与参考文献[71]中概述的情景完全一致。

图9. 热力学空间的数据驱动应用。(a) 描述伴侣蛋白（GroEL）对MDH蛋白（苹果酸脱氢酶）作用的CRN。天然MDH二聚体可以通过中间状态错误折叠，而伴侣蛋白使用ATP来稳定天然折叠构象。灰色盒子中包含了一个简化版本的反应图。（b）天然二聚体（[ND]）的浓度与天然单体（[NM]）和中间体（[I]）的浓度之间的关系。即使在ATP存在的情况下，实验数据（点）仍然位于平衡线上。（c）错误折叠蛋白质（[M]）的浓度偏离了与[ND]的平衡状态，位于热力学空间（绿色区域）内。蓝线和红线分别表示下限和上限，以及相应的路径（彩色盒子）。参数来源于参考文献[71]。（d）CRN描述了一种光驱动的光酸。在光照下（黄色箭头），MCH向其顺式异构体的转化速率加快，导致溶液中[某物质]的净变化。（e）化合物1的MCH浓度（[MCH]）与[某物质]的关系，即含有长度为[某值]的化学桥的MCH。参考文献[75]中的数据与平衡线不吻合，但它们位于为增加的非平衡驱动[某值]而绘制的热力学空间（绿色区域）内。彩色盒子代表通往上下限的路径。（f）[某值]的估计精度表示为[某值]，其中[某值]是实验值与推断出的[某值]之间的相对误差，该误差与观察到的浓度相符。不同的[某值]表示不同的化合物。（g）从[某值]我们可以看到每种化合物的功率依赖性量子产率的下限（黑线）。参数来源于参考文献[75]。B. 光驱动的光酸

第二个例子涉及水中的质子化孟氰蓝的动态。这是一种光酸，它会经历一系列结构和酸碱转化，如图9(d)所示。从质子化形式（MCH）开始，系统可以通过反应速率常数[某值]去质子化生成孟氰蓝（MC），并释放[某物质]到溶液中。通过环闭合反应，MC可以转化为螺吡喃（SP），其反应速率常数为[某值]，SP再次可以被质子化形成MCH的顺式异构体（[某值]），其反应速率常数为[某值]。这种物种最终会通过反应速率常数[某值]重新转化为MCH，从而完成反应循环。与逆循环相关的反应速率常数为[某值]和[某值]。在光照下，MCH转化为[某物质]的反应加速，而在无光条件下则极其缓慢。根据参考文献[75]，为了简化问题，我们假设光照使MCH向其顺式异构体的反应速率加快了[某值]的光化学速率，使系统进入一个光稳态，导致溶液中[某物质]的可测量变化。由于这一特性，这些化学物种可以作为可调[某物质]的缓冲剂[76]。这一系列循环操作首先在参考文献[75]中通过实验和理论进行了表征，研究了不同长度化学桥的MCH家族，以及不同功率的可见光源。这种表征提供了一个可靠的基础，我们可以用它来验证我们从界限中得到的预测。我们提到，为了进行这样的比较，如上所述，我们遵循了文献中引入的模型，因此没有将无光子和光辅助的路径视为独立的实体。因此，得到的结果应被解释为在这个有效的粗粒度模型内的热力学一致性界限。进一步的工作可以基于实验数据，对光驱动的光酸进行扩展分析和明确的热力学描述。

在参考文献[75]的实验中，我们在图9(e)中报告了在不同光稳态条件下[某物质]与[某物质]的关系，每种条件都与不同的激光功率相关联。我们展示了这些浓度受到将[某物质]转化为[某物质]的两个反应途径相关的平衡常数的限制，根据方程（29）。平衡线与实验数据不吻合，但随着非平衡驱动[某值]的增加，可进入的空间扩大了。

1. 非平衡驱动的推断

想象进行一个不同的、更简单的实验。我们观察到溶液中的光酸，无论是[某物质]、2还是3，暴露在未知功率的可见光源下。通过选择性波长的UV-Vis吸光度，可以容易地推断出溶液中所有物种的稳态浓度比率[某值]、[某值]、[某值]和[某值]。从方程（29）出发，为了简单起见，我们关注一个单一比率，得到（55），其中上标[某值]表示观测到的量。在这种情况下，我们无法通过[某值]直接测量[某物质]来验证界限。然而，[某值]是[某值]和[某值]的函数。通过先前独立进行的校准实验（例如，参见参考文献[75]），我们通过拟合现有数据来推断函数[某值]——从而推断出[某值]。将[某值]代入方程（55）并结合[某物质]反应的平衡估计[某值]，我们得到了功率的最小值[某值]，从而使上下限成立。特别是，我们发现上限[某值]提供了更有信息的约束，允许我们在未知功率下确定与观测结果相符的最小非平衡驱动[某值]。在受控功率条件下的实验还提供了每种化合物的平均量子产率的估计。由于[某值]（[某值]量化了每个分子的光子能量[75]），我们将这个最小功率映射到一个与我们的观测结果相符的[某值]的最小值。

为了测试[某值]的推断，我们将其应用于参考文献[75]中的[某值]和[某值]。在这种情况下，我们知道所有参数的正确值，特别是[某值]和[某值]的基准值，因此我们能够量化我们方法的稳健性。在图9(f)中，我们绘制了三种化合物在不同激光功率下的[某值]与[某值]之间的相对误差。我们注意到推断方法的性能非常好，高达75%-85%的非平衡驱动被方程（55）中的热力学界限所捕捉。实验上，在不同功率下发现了略有不同的量子产率。显然，[某值]的估计为每种化合物的功率依赖性[某值]提供了下限[图9(g)中的黑色实线]。这些结果表明，我们的热力学界限原则上可以用来推断CRNs中的非平衡驱动信息。

VII. 讨论

在这项工作中，我们推导了通用CRNs在稳态下的基本热力学约束，引入了热力学空间的新概念来描述与系统能量预算相容的物种浓度范围。类似的界限也可以用于反应亲和力，并用来估计CRN中任意两个物种之间有效亲和力的限制。我们的方法基于CRN超图几何，并利用新颖的化学探针技术来研究那些不是通过现有反应直接连接的物种。因此，它显著扩展并推广了之前关于线性和催化网络的发现，适用于任意网络拓扑的挑战性案例。这种推广对于复杂的多分子CRNs的应用至关重要，例如多稳态生化网络和自组装过程，在生物系统中很常见，但这些研究的热力学分析还不够全面。我们相信，我们引入的结果将对于仅从全局热力学性质出发来表征CRNs至关重要，而不需要依赖对整个动力学的详细了解。

我们通过四个示例展示了我们框架的适用性：Schl?gl模型用于双稳态，一个热力学上一致的手性对称性破坏模型，一个简单的自组装放大，以及反应-扩散模式的开始。这些案例研究说明了我们的界限如何阐明生物学和生化系统中各种复杂行为出现的热力学前提条件。我们的方法桥接了非平衡热力学、网络理论和系统生物学的概念，为分析生物系统提供了一个强大的跨学科工具。我们还开发了一个开源的python包TACOS，用于确定任意小型到中型CRNs的热力学空间（见附录C），以便于广泛应用所呈现的结果。

我们框架的通用性为超出此处探索的应用开辟了途径。例如，它可以扩展到限制多分子CRNs中对动力学参数扰动的响应，建立在之前对于线性系统的研究基础上[31]。我们的方法还有潜力推广到不施加伪单分子约束的分子模板网络中的信息传播研究[56,77]。此外，构建热力学空间的想法之前已在代谢网络的背景下作为受限优化问题进行过讨论[78]。我们的框架提供了这个空间的确切形式，从而将优化动力学速率常数的非线性、高维问题转化为最优反应路径的结构热力学识别。此外，推导出的界限不仅仅是简单地划定稳态解：它们揭示了实现最优性的化学路径。因此，这些结果可能有助于识别基于涉及自由能和外部化学调控物种的代谢活动操作范围的基础构建块。未来的工作还可以探索生化系统中的其他非线性现象的应用，例如生化振荡。一个令人兴奋的方向将是研究各种吸引子类型（例如，极限环、奇异吸引子）的热力学界限，可能允许人们将热力学空间的概念扩展到稳态之外。这可能会提供关于生物系统中复杂动态行为（如昼夜节律和钙振荡[10,79]）的热力学要求的洞察。

总之，这项工作为分析生物和合成系统中的非平衡现象提供了一个通用框架，有助于更深入地理解控制复杂化学和生物过程的原则。通过建立CRNs运行的基本热力学约束，我们的方法为设计具有生命特性的合成分子系统提供了新的视角，例如通过复杂的自组装过程出现的信息序列。最终，这项工作为理解复杂、远离平衡的生物和生化系统的热力学理论做出了重要贡献。