![]()
导语
长期以来,统计物理中的变分原理为平衡态体系提供了坚实的理论基础,但非平衡态体系始终缺乏一个同等普适且自洽的理论框架。随着研究对象从宏观热学系统拓展至少颗粒、生物与信息系统,传统假设逐渐失效。本文系统综述最大口径原理这一基于路径熵最大化的非平衡统计方法,阐明其理论基础、数学形式与操作流程,并通过扩散、电流分配及基因电路等实例,展示其在复杂少颗粒动力学建模中的独特优势与广泛应用前景。
关键词:非平衡统计物理、最大口径原理(Maximum Caliber)、路径熵(Path Entropy)、少颗粒体系、动力学变分原理
Kingshuk Ghosh等丨作者
张博文丨译者
赵思怡丨审校
![]()
论文题目:The Maximum Caliber Variational Principle for Nonequilibria 论文链接:https://pmc.ncbi.nlm.nih.gov/articles/PMC9827727/ 发表时间:2020年2月19日 论文来源:Annual Review of Physical Chemistry
摘要
自1865 年克劳修斯(Clausius)与 1877 年玻尔兹曼(Boltzmann)的开创性工作开始,熵及其最大化的原理便成为了从微观性质推导宏观物质平衡态的理论基石。然而,尽管相关研究已开展颇丰,科学界至今仍未建立一个同样令人满意的、适用于非平衡态体系的普适变分原理。直到 1980 年,E.T. Jaynes Shore与Johnson的研究为该领域开辟了新的方向。本文在此综述最大口径原理(Maximum Caliber)——一种类似最大熵原理,能够在给定动力学约束的条件下推断路径上的流分布的方法。该方法正为复杂体系的研究提供全新视角,尤其适用于少粒子复杂系统,例如基因电路、蛋白质构象的反应坐标、网络流量、鸟群集群行为、细胞迁移以及神经元放电等。
非平衡统计物理学最近的研究方向
非平衡态物理(Nonequilibrium Physics, NEP)的研究核心是流—— 通常为物质流、热流或电流。传统模型多聚焦于宏观尺度:在此尺度下,分子或粒子的数量足够庞大,其分布可被表示为空间 x 与时间 t 的连续可微函数(如密度或浓度c(x,t) )且其涨落效应可忽略不计。典型的宏观流模型包括纳维 - 斯托克斯流体力学方程、欧姆电流定律、基于菲克定律(Fick’s law)的粒子梯度流,以及基于傅里叶定律的热流。
有趣的是,当前正有三大因素推动该领域发展产生新的研究成果:
a. 非平衡态物理的研究范围已突破传统热学材料的范畴拓展至多个领域,包括以下:互联网中的信息流 (1)、城市间的人口流动 (2)、股票市场中的资产流 (3)、鸟群的集群行为 (4)、科学论文的引用流 (5)、细胞网络内生化分子与蛋白质的运输及信号传导、细胞内基因与蛋白质的进化动力学 (6−9)、大脑中的神经信号 (10)、生物的进化与发育过程等其他诸多应用场景。这些体系均不涉及过去的传统非平衡态物理的核心研究对象 —— 热力学浴、物理功或粒子碰撞。
b. 得益于单粒子与少粒子实验技术的出现,科研人员对微观尺度的研究兴趣日益浓厚,此类研究正好属于非平衡态统计物理的研究范畴,系统的涨落效应、速率分布与路径分布是该类研究的核心关键因素:
涨落效应:少数粒子体系中,微观状态的随机波动不可忽略;
速率分布:同一反应中,不同粒子的反应速率并非均一,而是呈现一定的分布;
路径分布:微观粒子的动态过程并非只有一条固定路径,而是存在多种可能的路径,且不同路径具有不同的概率。
c. 科学界仍在持续探索一种普适变分原理作为非平衡态统计物理理论基石,就像热力学第二定律与玻尔兹曼分布作为理论基石支撑着平衡态物质统计物理。
![]()
图1 非平衡态物理领域的发展历程简要时间表,最上行为发明者,中间行为非平衡过程与对应模型,最下行为平衡态。蓝色描述的是研究和承认的原理。
图 1 描述了非平衡态物理领域的简明发展历程。19 世纪中叶,一系列唯象模型相继建立,包括牛顿粘性流体模型、菲克粒子流定律、傅里叶热流定律与欧姆电流定律。工业革命推动了人们对蒸汽机中功与热的理解。克劳修斯于 1865 年左右发现,平衡态的预测可基于一个物理量的最大化趋势——他将该物理量命名为熵,其定义为 :
![]()
此后,热力学第二定律这一变分原理应运而生。不久之后,气体动理论、统计热力学与玻尔兹曼 - 吉布斯分布定律相继被提出。这些理论极大地增强了热力学第二定律的解释力与丰富度,具体体现在两个方面:(a) 建立了通过分子与材料的微观性质解释宏观平衡态的模型;(b) 将我们对热力学第二定律中熵最大化的理解,建立在微观态的概率分布之上。然而,这些原理存在一个关键局限——它们仅适用于平衡态或近平衡态体系。
科学界始终在探索适用于非平衡态的变分原理。研究要回答的核心问题是能量耗散率 (dU/dt) 或熵产生率 (dS/dt) 是否存在极大值或极小值的趋势。相关研究实例包括 Onsager的最小能量耗散原理 (11)、Prigogine的最小熵产生原理 (12) 与最大熵产生原理 (13)。然而,这些候选原理始终未能完全令人满意。其局限性在于:它们仅适用于近平衡态体系,需要引入局部平衡或与热浴弱耦合等假设;且主要应用于宏观尺度,通常是预设而非预测唯象关系。建立普适非平衡态统计物理变分原理的一个主要障碍是,目前尚未发现一个能够与平衡态中的 (具有同等核心地位的实验关系)。
如今,关于流的模型已十分丰富,包括扩散方程、朗之万方程、主方程、随机游走、玻尔兹曼输运理论等。那么,我们为何还需要一个底层的变分原理?原因其实很简单,当前的模型均依赖一定的简化假设而非普适的底层原理。例如,粒子或时间步的独立性、近平衡态、线性力-流关系、高斯噪声、大粒子数,或聚焦于热浴与粒子碰撞。对于更具挑战性的非线性少粒子复杂动力学(传统简化假设不再成立),以及热物理之外的应用场景,我们需要一个普适变分原理作为研究的指导框架。本文综述的证据表明,最大口径原理或许正是这样一种原理。
最大熵与玻尔兹曼分布:
热力学第二定律的微观基础
在讨论非平衡态前,我们先回顾最大熵原理与玻尔兹曼分布。玻尔兹曼提出的著名公式 S = k\ln W 揭示了一个核心结论:克劳修斯提出的宏观热力学原理,其本质源于系统可能的微观排布方式的数量 W 。这一公式同时确立了玻尔兹曼指数分布定律是热力学第二定律变分原理的微观表现形式。更具体地描述,对于任意离散选项 i=1,2,3,… 上的概率分布 {pi}=p1,p2,p3,…,我们可以定义该分布的数学熵(mathematical entropy)为:
![]()
该物理量可对任意概率分布进行计算,但需要注意的是,这种数学熵 Smath 与克劳修斯提出的物理熵 (physical entropy)SClausius 并非同一概念。同时,Smath也不是我们构建物理平衡态理论模型所需的熵;对于平衡态模型,我们需要的是态熵 Sstate ,其定义如下:
首先,将概率分布 pi 限定为系统微观态上的分布;其次,我们认为:只有能使熵达到最大值的唯一特定分布 所对应的熵,即满足:
![]()
这个关系的熵才与热力学第二定律对平衡态物理行为的预测相关。式(2)中,kB为玻尔兹曼常量。需要强调的是,kB是针对某一特定分布Sstate定义的,而非任意数学分布。对于正则系综(即系统与一个热浴相接触,热浴可固定系统的平均能量),其预测过程为:在平均能量约束下,对概率分布 pi 最大化数学熵Smath,即
![]()
其中,Ei为微观态 i 的能量,〈E〉为系统的平均能量。同时,概率分布需满足归一化条件:,
此时上述约束优化问题的解即为吉布斯-玻尔兹曼分布(Gibbs–Boltzmann distribution):
![]()
其中,为满足上述所有约束条件的微观态概率。式(4)中,T 为系统的温度。该分布是平衡态统计物理的核心,可用于预测模型中所有微观态 i=1,2,3,… 的平衡态布居数。
对熵的不同含义缺乏严谨区分,是造成相关概念混淆的主要原因:
SClausius 仅能预测宏观平衡态热力学行为,例如:热量倾向于从高温物体流向低温物体、粒子倾向于向低浓度区域扩散、材料内部的密度趋于均匀等。它无法提供任何关于微观尺度或分布函数的信息。
Sstate 是我们基于平衡态的微观物理模型计算得到的物理量。为了将平衡态模型与对应的宏观实验结果相关联,我们通常采用等效关系:Sstate = SClausius。
区别于与上述介绍的热力学相关的熵,在下文中,我们将介绍另一种与动力学过程相关的熵,即路径熵( Path Entropy )Spath。
最大口径原理:适用于动力学过程的变分原理
E.T. Jaynes四十年前在Annual Review of Physical Chemistry中首次提出了最大口径原理(14)。本文将综述该原理作为非平衡态普适原理的当前研究现状,及其部分应用。与其他非平衡态变分原理相比,最大口径原理的独特之处体现在四个方面:
其理论基础为粒子的轨迹,而非宏观的浓度;
其优化目标为路径熵的最大化,而非态熵的最大化;
基于实验数据约束推导微观模型时,其存在的逻辑矛盾与混淆更少;
其具有概率论的公理化基础 (15)。
![]()
图2 路径熵用于量化不同通路间的流量分布均匀性。线条粗细代表流量密度,即通路的出现概率。
以下是该原理的简要概述:最大口径原理之于动力学,正如最大熵原理之于平衡态。最大熵原理关注的是微观态的概率分布,而最大口径原理关注的是系统演化的路径或轨迹的概率分布。设 {Γ} 为系统在时间演化过程中所有可能的轨迹集合。 Γ 可描述多种类型的动力学过程:例如,考虑系统从初始时刻 Ti 的初态演化至终末时刻 Tf 的终态的轨迹(见图2),或处于稳态的系统的演化轨迹。对于前者,单条轨迹可表示 ,即系统在时间点 Ti 至 Tf 之间可能经历的所有状态序列。其他类型的轨迹将在下文讨论。设 PΓ 为定义在轨迹系综 {Γ} 上的概率分布。
设 J(Γ) 为定义在轨迹空间上的泛函。J 的具体例子包括:流所携带的质量/热通量、沿轨迹的平均耗散量,或沿轨迹的平均能量。若我们希望在约束平均泛函值固定的前提下,推断轨迹空间上的概率分布 PΓ,则该约束条件可表示为:
![]()
在这里以及本文的后续部分中,我们将统一使用大写字母 P 表示轨迹上的概率分布,而使用小写字母 p 表示通用概率;当使用小写 p 时,我们会明确说明其具体含义。
需要注意的是,满足上述约束条件的概率分布 PΓ 可能有无限多种。与平衡态的情况类似,我们在此也采用熵最大化的策略,但此时的优化对象是所有可能的轨迹,而非状态。对应的路径熵定义为:
![]()
我们的优化目标为:在式(5)的平均泛函约束与概率归一化约束下,最大化上述路径熵(见图2)。其中,qΓ 为轨迹空间上的参考分布或先验分布。
该约束最大化问题可通过引入拉格朗日乘子来求解。我们定义一个无约束优化函数,将其称为口径(Caliber),记为 C:
![]()
在式(7)中, γ 是用于将系综平均 〈J〉 固定为给定值的拉格朗日乘子,而 α 是用于保证概率分布归一化的拉格朗日乘子。对口径 C 关于 PΓ 求最大值后,我们得到最优轨迹分布:
![]()
其中,
![]()
是对所有轨迹的权重求和,称为动力学配分函数(dynamic partition function)。其在动力学中的作用与平衡态统计物理中的配分函数完全对应。
最大口径原理的一个核心结论是:可测量的平均速率类物理量与模型的动力学配分函数之间存在明确的解析关系。该关系可表示为:
![]()
在实际应用中,最大口径原理数学框架的操作流程如下:
明确与当前研究问题相关的轨迹类型;
基于相关的实验约束条件,构建轨迹空间上的概率分布(式9)。每条轨迹的概率由其包含的所有演化步骤的统计权重决定;
利用式9,可进行与平衡态统计物理类似的理论预测。例如,将所有轨迹的权重求和得到动力学配分函数 Qd;
利用式11,可计算出与给定的平均泛函值〈J〉 及其他可能的约束条件相一致的所有统计权重与轨迹概率。
为了进一步阐释其数学框架与应用方法,我们将在下文给出最大口径原理的若干具体应用实例。
最大口径原理提供新的研究视角
两态动力学:基于平均速率预测通路分布
![]()
图3 采用最大口径方法结合马尔可夫模型表征 A ↔ B 两态动力学。
(a) 不同的轨迹由不同的势能面表示。最大口径方法对所有轨迹进行枚举,其权重为初始未知的路径权重。
(b) 实验测得的平均量(例如粒子滞留于态 A 的平均次数 ⟨Naa⟩)可用于确定这些路径权重,进而获得所有通路的相对概率。
(c) 基于平均值预测的方差与实验结果高度吻合。
考虑一个可在两个势阱 A 与 B 之间跳跃的单胶体粒子(见图3)。Phillips 及其团队通过双激光光镊构建这一两态体系,并开展了相关实验研究(16)。图3展示了作者在不同光镊条件下观测到的粒子轨迹实例。本研究采用最大口径原理结合两态马尔可夫模型对该体系的动力学行为进行建模。图中枚举了所有可能的粒子通路,并标注了各通路对应的统计权重。单个时间步长内的统计权重定义如下:上跳(A→B)的权重为 u,滞留于态A的权重为 a,滞留于态 B 的权重为 b,下跳(B→A)的权重为 d。上述四个权重中仅有两个是独立的,原因在于:从任意初始态出发的所有跃迁路径的权重之和必须为1,且到达任意终态的所有跃迁路径的权重之和也必须为1。
最大口径方法的核心流程为:首先从实验轨迹中提取两个独立的权重参数(例如 u 与 d),随后将其代入图3所示的最大口径公式中。在这两个独立参数确定后,即可获得完整的通路概率分布。例如,最大口径方法能够基于平均速率预测速率分布的方差,且该理论预测的方差值与实验测定结果高度吻合。
这一马尔可夫动力学与最大口径方法相结合的研究框架,已被进一步应用于单分子三态循环体系(A ↔ B ↔ C ↔ A)的研究中(17)。研究结果表明:平自选速率的提升会使噪声以更快的速度衰减。此外,该研究还证实:基于轨迹的建模方法能够捕捉所有可能的粒子轨迹,而非仅局限于平衡态附近的轨迹,因此其适用范围可拓展至远离平衡态的体系。
基于超越菲克定律平均形式的变分原理推导唯象定律
![]()
图4 扩散方程建模将浓度 c(x,t) 视为连续且可微的物理量 —— 例如,菲克定律的表达式即为 ⟨J⟩=−D∂x∂c。
本小节将通过一个简单的微观模型,从最大口径这一变分原理出发推导菲克扩散定律(),以阐释最大口径方法的核心思想。非平衡态物理(NEP)中的传统扩散方程模型通常假设:体系中的粒子数量足够多,因此其密度或浓度可被表示为连续且可微的函数(见图4)。
![]()
图5 显微镜载玻片上胶体颗粒的微观尺度扩散行为(18)。
(a) 用于测量少量胶体微球自由扩散过程、以量化其涨落的微流控装置。
(b) 用于追踪微球随时间运动轨迹的视频快照。
(c) 利用该装置测得的三条典型浓度分布曲线。
这些分布曲线表明,本实验所测量的此类少颗粒体系的扩散流存在显著涨落。而菲克定律这类唯象理论仅能描述大颗粒数体系的平均行为,无法解释本文所观测到的涨落现象。
然而,近年来的胶体实验已开始探索小数量粒子体系的扩散行为(见图5)。小数量粒子扩散具有诸多值得关注的特性(例如通量分布的宽度〈J2〉),这些特性对于模型构建具有重要意义。利用最大口径方法的轨迹分析框架,结合本研究中提出的狗-跳蚤模型(将不同的空间位点视为“狗”,将位于这些位点上的粒子视为可在“狗”之间跳跃的“跳蚤”),可便捷地计算出这些特性。
![]()
图6 (a) 浓度梯度 c(x),标注了相邻的两个统计区间 i 与 i+1,直观显示两区间内的粒子数 N 存在差异。
(b) 一条可能的粒子运动轨迹,其中各步均标注了对应的统计权重。
图6a展示了沿空间坐标 x 分布的浓度梯度 c(x) 。在位置 x i 处存在 N i 个粒子。在狗-跳蚤模型的表述中(18), N i 代表位于空间位点 x i 处的“狗”身上的“跳蚤”数量。在一个时间间隔 Δ t 内,有 r 只“跳蚤”从该“狗”身上向右跳跃,另有 ℓ 只“跳蚤”向左跳跃。本研究设定两个约束条件:平均向右跳跃的跳蚤数为 ⟨r⟩ ,平均向左跳跃的跳蚤数为 ⟨ℓ⟩ 。那么,在单次观测中恰好出现 r 次右跳与 ℓ 次左跳的概率 P i 是多少?
我们通过在已知平均速率的约束下最大化路径熵来计算概率 P i ,并采用拉格朗日乘子 λ r 与 λ ℓ 来引入这些约束条件,最终得到:
![]()
其中,我们通过定义 σ= exp (λ) 对公式进行了简化;而 Q i 为动力学配分函数,其表达式为:
![]()
式中,简并度 ,该表达式的推导基于粒子跳跃行为相互独立的假设。向右的平均通量 ⟨r⟩ 可表示为:
![]()
类似地,相邻的下一个空间列(包含 N i+1 个粒子)向左的平均通量为:
![]()
接下来,我们假设左右跳跃具有对称性(即体系无漂移运动),此时可得 λ r = λ ℓ 。我们定义 ,则净向右通量可表示为:
![]()
我们将该粒子数差值转换为通量(即单位时间、单位面积内跳跃的粒子数),并利用浓度替代粒子数: Δ N=−A Δ c Δ x 。其中, A Δ x 为包含这些粒子的体积;负号的含义为:我们定义的正方向是从位点 i 到位点 i+1 ,而浓度梯度 Δ c= c i+1 − c i 的定义方向恰好相反。代入后可得:
![]()
公式(16)即为适用于两个相邻粒子平面的菲克扩散定律。根据该模型,我们可得到扩散系数的表达式 。这一模型揭示了粒子的平均流动与浓度差成正比的根本原因:高浓度区域具有更多的粒子流动通路。
![]()
图7 基于狗 - 蚤模型,利用最大口径原理推导得到的简单扩散理论,成功预测了少颗粒体系的实验结果(19)。该理论的核心结论包含以下五点:
1. 它从变分原理出发,推导出了菲克第一定律,即 ⟨J⟩=−D ∂ x ∂ c ;
2. 它证明了菲克定律在少颗粒极限下依然成立;
3. 它能正确预测完整的速率分布;
4. 它计算了一个类似麦克斯韦妖的物理量 ——反向流(wrong-way flows)的占比(由 ϕbadactor 定量描述),结果表明,随着净通量的增大,该反向流占比会变得可以忽略;
5. 它能准确给出通量涨落关系。
尽管存在多种推导菲克定律的方法,但本研究采用的最大口径方法具有以下三个核心优势:
(a) 该方法基于一个普适的底层变分原理推导菲克定律,具有更强的理论普适性;
(b) 该方法不仅能得到平均通量 ⟨J⟩ (公式16),还能给出完整的速率分布 P i (见图7);
(c) 最大口径方法的教学直观性强,易于理解和推广。
基尔霍夫电流定律的推导:含两个电阻的节点处电流如何分配?
![]()
图8 最大口径原理可推导出基尔霍夫电流定律 —— 即节点处的电流分配比例,与各支路的流动电阻成反比。
基尔霍夫电流关系指出:粒子流或流体在节点处的分配比例,与节点下游各支路的电阻成反比。本文总结了Jaynes的相关论证(14):最大口径原理可正确推导出这一规律,而一个假想的最小熵产生率原理则无法实现。尽管基尔霍夫原理适用于任意类型的流动,但为简化分析,我们在图8中以两个并联电阻 R l (左电阻)与 R r 右电阻)为例进行说明。
根据欧姆定律,通过电阻的电流可表示为 i= V 0/ R ,其中 V 0 为外加电压, R 为电阻值。我们假设两个电阻分别与温度为 T l 和 T r 的热浴相连。由此引出核心问题:当两个电阻并联时,总电流 I= I l + I r 。如何在左支路(流经 R l 的电流 I l )与右支路(流经 R r 的电流 I r ) 之间分配?
在最小熵产生原理的框架下,电流的分配比例由热熵产生率的极值条件决定:即对电流在两条支路中的分配方式求熵产生率 的极值。我们将熵产生率写作:
![]()
其中, J l 与 J r 分别为左、右两个电阻上的热流速率。
式(17)中,电阻的耗散项采用了经典形式:力×通量(即 V 0 ×I )。同时,我们也引入了欧姆定律的关系: V=IR 。在总电流约束 I= I l + I r 下,对 I l 和 I r 求 的最小值,可得到:
![]()
这一结果显然是错误的。原因在于:对于外加电压为 V 0 的电阻,其电流大小仅由电阻值 R 决定;温度 T 本不应出现在该表达式中。当然,电阻值本身可能与温度相关(即 R=R(T) ),但对于基尔霍夫定律而言,唯一的决定因素是电阻的实际取值,而非温度本身。
接下来,我们采用与前文相同的方法,将狗-跳蚤模型与最大口径原理相结合,推导正确的电流分配规律。假设存在 N 个电子(对应模型中的“跳蚤”),在一个给定的时间间隔 Δ t 内,有 ℓ 个电子可以跃迁至左支路, r 个电子可以跃迁至右支路。我们将 ℓ 与 r 作为两个约束条件引入模型——这两个约束本质上对应于两条支路的固有属性。通过引入两个拉格朗日乘子 λ l 和 λ r 来体现这些约束,口径可表示为:
其中, P j 代表某一微观轨迹的概率,该微观轨迹由 ℓ 与 r 的一种特定分配方式定义(归一化条件 ∑ P j =1 ) 已被隐含假设)。对口径进行最大化,可得到微观轨迹的概率分布: P j ∝ exp (− λ l ℓ− λ r r) 。
随后,我们将微观轨迹的概率分布转化为宏观轨迹的概率分布 P M (ℓ,r) ——宏观轨迹仅由 ℓ 与 r 的取值定义。这一转化过程中会引入一个组合因子: 。利用该组合因子,我们可以得到 P M (ℓ,r) 服从多项分布,并进一步推导出: ⟨ℓ⟩=N P l , ⟨r⟩=N P r 。其中,左、右支路的跃迁概率分别为:
![]()
由此,我们可以得到一个关键关系:
![]()
下面我们将这一模型结果与欧姆定律建立联系,从而推导出基尔霍夫电流定律:
1. 欧姆定律中的电流是宏观统计平均值,因此在本模型中,支路电流的比值满足 i l / i r =⟨ℓ⟩ / ⟨r⟩ ;
2. 模型中的固有属性(跃迁概率的比值)与欧姆定律中的电阻值相对应,即 P l / P r = R r / R l 。
结合以上两点,我们即可推导出基尔霍夫电流定律——平均流量的分配比例与支路电阻的比值成反比。
这一推导过程具有两层重要意义:
第一,最大口径原理仅通过约束条件即可完成推导,不会像最小熵产生原理那样,错误地引入温度等无关因素;
第二,与仅适用于平均通量的基尔霍夫定律不同,最大口径原理与狗-跳蚤模型的结合,还能给出节点处流量的完整速率分布。
最大口径原理对少颗粒复杂系统与基因电路的建模
化学反应、生化网络及基因电路通常并非简单的线性系统;它们可能包含非线性元件、负反馈或正反馈环路、振荡器、开关、门控及类逻辑元件。这些系统的底层作用细节往往是未知的。此外,此类系统的建模挑战还会因数据的高噪声特性而进一步加剧——这是基因表达过程中涉及的粒子数量极少所导致的固有问题。更重要的是,实验通常仅能测量少数带荧光标记的蛋白质的动力学行为,而其背后的基因表达过程可能由多个调控因子共同驱动。因此,仅通过实验数据往往无法推导出所有底层相互作用及其速率常数。那么,我们能否换一种思路:构建一个与实验数据自洽的有效动力学模型,并利用该模型进行可靠预测?
最大口径原理是构建此类模型的理想工具——它能够基于有限的信息(如少数物种的实验数据),建立少颗粒系统的有效模型。最大口径原理所构建的模型具有最少的参数数量,同时这些参数能够直接从实验数据中最大化地确定。下文将通过三个实例,展示最大口径原理如何捕捉少颗粒基因电路的复杂动力学行为。在基因电路中,蛋白质由DNA转录翻译生成,同时这些蛋白质又可以结合到DNA上,调控自身或其他蛋白质的合成速率。
单基因电路中自激活行为的建模
![]()
图9 最大口径原理可预测自激活基因电路的动力学行为。(a) 基因α负责合成蛋白质 A。当 A 蛋白的二聚体A2结合至启动子区域时,A 蛋白的合成速率会显著加快(20)。需注意,负反馈电路可采用完全相同的方法进行建模;区别仅在于:本图中被称为启动子的浅蓝色区域将被替换为阻遏物结合区域,且阻遏作用的效果是减慢而非加快 A 蛋白的合成。(b) 实验测得的具有随机性的开关型时间轨迹。(c) 以该轨迹为输入,最大口径原理可预测出 A 蛋白在正常合成状态(速率g)与加速合成状态(速率g∗)下的合成速率,以及其降解速率d。
图9展示了单基因调控电路中自激活行为的实现机制(20-22)。基因 α 负责合成A型蛋白质(A蛋白),该蛋白质以降解速率 d 发生降解。基因 α 的DNA序列两侧带有启动子区域。当两个A蛋白分子形成二聚体 A 2 并结合到该启动子区域时,基因 α 合成A蛋白的速率将显著高于基础水平。
实验数据的表现形式为:单位时间内A蛋白分子数量 N A (t) 随时间变化的、具有噪声的开关型时间轨迹(见图9b)。仅通过这些数据,我们无法推导出所有的微观速率参数。因此,在对该电路进行建模时,我们的目标是:从完整的随机时间轨迹中(而非仅对轨迹进行平均)提取最大量的信息,以构建一个有效模型。我们利用该随机时间轨迹,来推断三个核心参数:蛋白质的基础合成速率 g 、降解速率 d ,以及启动子被激活时的有效加速合成速率 g ∗ 。我们通常无法先验地知道这些速率常数、 A 2 与启动子的结合亲和力,或其他类似的机制性变量。其他建模方法可能会明确引入这些变量,但这通常需要增加额外的参数,而这些参数的取值往往缺乏实验依据。
最大口径原理能够直接从随机实验数据中提取出这三个核心物理量(20, 22)。此外,最大口径原理还能预测一个有效反馈参数 K ,该参数用于量化蛋白质与其自身启动子之间的耦合强度。该模型还能自然地产生双峰分布——这是模型能够成功描述开关行为的必要条件。更进一步,最大口径模型还能为我们提供深刻的洞见:如何通过调节反馈参数来改变这些分布特征(20, 22)。
下面我们具体说明最大口径原理如何应用于这个简单的基因电路。我们将口径定义为:
![]()
其中, ℓ α 为合成状态变量,其取值为0到某个预设最大值 M 之间的整数; ℓ A 为降解状态变量,用于描述在一个时间间隔结束时,仍未降解的、先前已存在的蛋白质分子的数量。这两个约束条件对应的拉格朗日乘子分别为 b α 和 b A ; 定义为观测到一组特定的 ℓ α 与 ℓ A 组合的概率。
因此,式(21)中的第一项为路径熵;第二项和第三项分别为对平均合成速率与平均降解速率施加的约束;最后一项通过拉格朗日乘子 K A 对 ℓ α ℓ A 的平均值施加约束,其物理意义是强制实现蛋白质分子数量 N A 与A蛋白合成速率之间的正相关关系。这是捕捉反馈本质所必须施加的、两个变量之间耦合的最低阶项。
基于该口径的定义,我们可以得到对应的口径最大化的路径概率分布:
![]()
其中,归一化常数 Q d 的表达式为:
![]()
实验观测到的轨迹的似然函数 L 可以用路径概率 来表示。通过最大化该似然函数,我们可以确定四个关键参数: M 、 b α 、 b A 与 K A 。在确定这些拉格朗日乘子之后,我们可以进一步利用它们来推断电路的不同速率常数、反馈参数及其他特征(20)。这些参数均无法通过实验直接测量。定义反馈参数为蛋白质合成状态变量 ℓ α 与当前存在的蛋白质状态变量 ℓ A 之间的皮尔逊相关系数。该反馈参数的取值是一个有效度量,用于描述A蛋白的存在对其自身合成速率的影响程度。在分子水平上,该度量可能会受到多种变量的影响,其中包括二聚体蛋白分子与启动子位点的结合常数。
基因双稳态开关电路
![]()
图10 最大口径原理(Max Cal)可给出双稳态开关基因电路中的速率分布。(a) 基因α合成蛋白质 A,基因β合成蛋白质 B。蛋白质 A 的结合会阻遏 B 的合成,蛋白质 B 的结合则会阻遏 A 的合成。(b) 该调控模式的最终效应为双稳态(赢家通吃):当 A 或 B 任意一种蛋白质的含量出现少量过剩时,其丰度会进一步提升,最终完全主导整个体系。图中为实验测得的随机时间轨迹。图 b 下半部分经授权改编自参考文献 22。(c) 无需其他任何先验信息,Max Cal 仅以两种蛋白质的随机时间轨迹为输入,即可推断出不同的速率参数(基础合成速率g、阻遏态合成速率g∗、降解速率d等)
图10展示了一个基因双稳态开关电路,该电路最初由Gardner等人设计(23)。基因 α 以速率 g ∗ 合成A蛋白,基因 β 以速率 g ∗ 合成B蛋白。每个基因的DNA序列两侧均带有阻遏物结合区域。当一个B蛋白分子结合到基因 α 的阻遏物区域时,会减缓A蛋白的合成速率;同理,当A蛋白结合到基因 β 的阻遏物区域时,也会减缓B蛋白的合成速率。
对于该电路,实验数据为两种蛋白质的分子数量 N A (t) 与 N B (t) 随时间变化的时间轨迹(见图10b)。最大口径模型能够产生交替输出的行为——这是双稳态的标志性特征:当A蛋白的分子数量开始超过B蛋白时,系统会进入“赢家通吃”的状态,A蛋白成为主导;反之,当B蛋白的数量占优时,系统也会切换到B蛋白主导的状态。
值得注意的是,最大口径原理仅需利用原始的输入数据,就能准确预测出蛋白质在基础状态和阻遏状态下的合成速率、降解速率,以及A蛋白与B蛋白之间的相关关系(22, 24)。
阻滞振荡器——一种振荡型基因电路
![]()
图11 最大口径原理可描述阻遏振荡器基因电路的少颗粒动力学行为。(a) 基因α、β、γ分别合成蛋白质 A、B、C。蛋白质 A 的结合会阻遏 B 的合成,蛋白质 B 的结合会阻遏 C 的合成,蛋白质 C 的结合则会阻遏 A 的合成。(b) 该调控模式产生的效应为振荡型时间轨迹。实验中常规测得的 A、B、C 三种蛋白质的丰度分布,其包含的信息要少于随机时间轨迹。(c) 无需其他任何先验信息,Max Cal 仅以三种蛋白质的随机时间轨迹为输入,即可推断出不同的速率参数(基础合成速率g、阻遏态合成速率g∗、降解速率d)及反馈强度K。图 c 经授权改编自参考文献 26。
图11展示了由Elowitz与Leibler设计的阻遏振荡器基因电路(25)。如图c所示,这是一个由三种蛋白质构成的环形电路:A、B和C。基因 α 合成A蛋白,基因 β 合成B蛋白,基因 γ 合成C蛋白。三种蛋白质的基础合成速率均为 g 。A蛋白能够结合到B基因的启动子区域,从而减缓B蛋白的合成速率;同理,B蛋白能够减缓C蛋白的合成速率,而C蛋白又能够减缓A蛋白的合成速率。蛋白质在阻滞状态下的合成速率为 g ∗ ,该速率远低于基础合成速率 g 。三种蛋白质的降解速率相同,均为 d 。
实验的原始数据为每种蛋白质的分子数量随时间变化的轨迹(见图11b)。与前文的建模方法相同,我们通过最大化观测到的带噪声振荡轨迹的似然函数,来推断底层的有效速率常数: g 、 g ∗ 与 d 。这一过程充分挖掘了隐藏在动力学数据中的全部信息(见图11)。通过该方法推断出的底层参数,与用于生成模拟数据的原始模型的参数具有很好的一致性(26)。最大口径原理还能预测有效反馈强度 K 。
当详细模型未知时,最大口径原理是建模的优选方法
![]()
图12 基因网络的传统动力学模型。(a) 质量作用模型(MA):通过任意非线性函数 f 描述平均行为(⟨A⟩),以实现对反馈的建模(kd 为降解速率)。(b) 质量作用 - 随机噪声耦合模型(MA + 随机噪声):在质量作用方程的基础上引入随机噪声,得到朗之万型方程。(c) 化学主方程模型(CME):基于跃迁概率 W 描述概率分布 P 的时间演化;而跃迁概率的确定,需要引入一组辅助物种 {Y},这类物种在实验中通常无法观测。(d) 化学主方程 - 质量作用耦合模型(CME + MA):一种粗粒化模型,其核心是将质量作用模型中使用的唯象函数 f 替代跃迁概率 W,以此描述概率的时间演化。
传统建模方法在处理少颗粒动力学与复杂系统时存在固有限制(图 12)。质量作用模型(Mass-Action Model, MA)仅能描述质量作用效应(即宏观尺度的体相平均动力学),无法刻画动力学涨落或速率分布。质量作用 - 噪声耦合模型(MA + Noise Models)试图弥补这一缺陷:其通过预设的分布形式引入时间涨落,朗之万方程便是典型代表(27)。然而,这类模型所预设的涨落分布并非在所有场景下都成立。例如,此类噪声模型仅适用于动力学行为服从简单单峰分布的系统,无法处理双稳态开关电路中存在的双峰分布问题。此外,该类模型通常通过引入非线性函数 f 来描述复杂动力学过程 —— 如希尔形式的函数 。这类函数往往是人为特设的,其形式无法通过独立实验验证,且需要引入多个可调参数。
化学主方程模型(Chemical Master Equation, CME)(28)能够明确且恰当地描述动力学涨落。但构建化学主方程模型的前提是,研究人员必须掌握系统完整的反应网络细节 —— 即明确多物种体系下的所有反应状态,以及状态之间的跃迁关系(对应图 12 中的 Y)。然而,这类多物种反应的细节往往无法通过实验进行验证。这就导致化学主方程模型存在两个关键问题:其一,模型包含的参数数量过多;其二,模型的相空间规模过于庞大。即便是双稳态开关电路或阻遏振荡器这类相对简单的系统,其相空间的计算量也可能达到难以处理的程度。尽管有限状态投影法(Finite State Projection, FSP)(29)能够在一定程度上压缩相空间的规模,但多物种体系带来的组合爆炸问题仍然是一个巨大的挑战。化学主方程 - 质量作用耦合模型(CME + MA Models)(30)虽然能够描述系统的涨落行为,但同样需要引入类似希尔函数 fHill 的非线性函数,且模型中的参数往往无法反映系统的真实底层作用机制。
当关于底层电路的可用信息极其有限时,最大口径原理便是建模的优选方法。它为推导观测轨迹的概率分布提供了一套严谨的理论框架 —— 即便是针对复杂动力学系统,也无需对非线性函数的形式作出常规的预设。得益于其自上而下的建模本质,该方法能够将相空间的规模控制在最小限度,从而解决了化学主方程面临的计算难题;同时,它又避免了化学主方程 - 质量作用耦合模型中人为特设的假设。最大口径原理会充分利用轨迹中包含的全部信息,而非仅对轨迹进行平均化处理 —— 即便是在实验无法精确测量目标物理量的情况下,这一优势依然存在。因此,将最大口径原理与用于参数求解的最大似然法相结合,能够妥善处理以荧光信号形式呈现的实验数据,而无需直接获得分子数量的精确值 —— 这正是实验中最常见的情况。通过引入荧光强度 - 分子数量转换的分布函数 ,最大口径原理可以直接构建原始荧光轨迹的观测似然函数。借助这一方法,该原理能够将单颗粒荧光信号的固有不确定性与少颗粒体系的本征涨落(由电路自身的结构特征决定)分离开来,并基于此构建电路的细节模型。已有研究证实,将最大口径原理、最大似然参数估计法与荧光 - 分子数转换模型相结合的建模框架,在多种不同的基因电路中均取得了成功的应用(20, 22, 26)。
在此,我们针对分子动力学建模的相关问题,提出一个具有普适性的观点。动力学教材中所描述的反应机制,通常是指在反应物转化为产物的主导反应路径上,可能存在的中间步骤。这类机制的研究对象是平均化的主导反应路径,因为相关实验通常是在包含大量分子的烧杯中进行的。而本文旨在解决一个截然不同的挑战:即如何利用路径分布所提供的额外信息(以随机时间轨迹的形式呈现),而非仅依赖平均化的结果,来研究少颗粒动力学系统中的反应机制。对于少颗粒流动系统而言,反应路径的分布本身也能提供丰富的机制性信息(即前文所述的拉格朗日乘子相关物理量)。但需要强调的是,最大口径原理与中间态建模并非相互排斥。我们可以轻松地对上述最大口径模型进行扩展,引入任何描述中间态所需的额外变量;这些变量的加入,将为中间态附近的反应路径分布提供更为深入的洞见。换言之,只要有相应的实验数据支持,最大口径原理始终可以纳入更多的信息。而在缺乏数据的情况下,该方法会在避免引入不必要假设的同时,构建出一个最小且自洽的有效模型。
最大口径原理可以进行单细胞数据驱动的网络参数的分布推断
细胞内部存在由速率系数 k 构成的生化网络,这些速率系数在时间尺度上近似保持恒定(31)。生化物种的丰度 x(k,t) 随时间的涨落服从泊松少颗粒统计规律 ,这类涨落被称为内源噪声(intrinsic noise)。该内源变异性的强度与物种平均丰度 N 成反比(标度关系为 ∝1 / N ),其中 N 为某物种的典型平均分子数。
另一种变异性被称为外源噪声(extrinsic noise)(32),其产生的原因是:速率参数 k 本身会在不同细胞之间存在差异。这种变异性会导致我们在不同细胞中观测到彼此不同的随机时间轨迹。那么,我们能否推断出包含外源噪声影响的物种丰度轨迹的分布 P(Γ) —— 即描述丰度 x(k,t) 整体分布的模型?
这一问题面临的核心挑战在于:流式细胞术、免疫荧光染色(33)及单细胞 RNA 测序(34)等实验技术,无法直接获取单个细胞内部的随机时间轨迹。这些技术仅能在某一特定的时间快照下,提供由细胞间差异所导致的生化物种丰度分布。最大口径原理能够基于此类时间快照数据,实现对外源变异性的定量分析(35)。Dixit 及其团队近期开发了一种新方法(35, 36),可从单细胞快照数据中直接推断出参数的分布 P(k) 。此外,最大口径原理还能进一步推断出细胞群体中物种丰度轨迹的分布 P(Γ) 。
最大口径原理可实现网络上输运动力学的估计
![]()
图14 最大口径原理可快速实现网络输运流量的估计。若已知在网络中流动的可移动单元的定态布居数,最大口径方程(式 25)便能给出完整的跃迁速率矩阵 —— 该矩阵即为使路径熵最大化的唯一解。
我们考虑如下具体问题:某一生物分子存在多种不同的亚稳态构象。我们希望计算得到完整的跃迁速率矩阵,其中包含任意两个构象态 a 与 b 之间的所有跃迁速率 k ab 。当采用分子动力学模拟方法解决这一问题时,计算过程往往既缓慢又具有极大挑战性——其根源在于,构象跃迁属于稀有事件,涉及高自由能的中间态,而这类中间态在模拟中难以被充分采样。然而,最大口径原理提供了一种快速的解决方案:只需输入易于快速获取的有限信息,即可对该速率矩阵进行近似估计。若计算机模拟能够在构象空间中进行充分搜索,以识别并采样得到各亚稳态的布居数 p a ;且我们已知一个或两个全局速率量(例如,整体跃迁过程的发生速率),那么最大口径原理便可预测出使路径熵最大化的跃迁速率矩阵(图13)。
具体而言,对于一个包含 N 个状态 {a,b,…} 的马尔可夫系统,其路径熵可由定态分布 {pa} 与跃迁概率 {kab} 明确表示。此时路径熵的可以表示为(37):
![]()
在有限的速率信息约束下,最大口径原理通过最大化式(24)中的路径熵函数,实现对跃迁概率矩阵的估计。
此处的有限信息可包含以下几类:
a .完整的定态分布 p a (例如,参见文献38, 39);
b. 定态的平均值 E=∑ p a E a ;
c. 动力学物理量的路径系综平均值。
根据所施加约束条件的不同,该熵最大化问题可通过解析或半解析的方式求解。具体而言,Dixit 等人(40)证明,最大口径原理所得到的跃迁速率 r ab 可表示为:
![]()
其中, γ 为拉格朗日乘子,其作用是强制动力学平均值 J 满足预设的约束值。
最大口径原理所构建的马尔可夫过程,已在多个研究领域得到成功应用,包括:生物分子动力学的解析(37–46)、生化反应网络的建模(47)、决策理论(48)以及机器学习(49)。在此,我们将选取其中两个应用实例,进行详细的阐述。
基于分子模拟的构象变化动力学估计
蛋白质分子的分子模拟常被用于获取构象变化的速率与路径,因为这些过程往往是决定生物机制的关键。如前所述,这类模拟面临的核心挑战在于:分子动力学模拟对稀有动力学跃迁(高自由能态)的识别与采样能力,远弱于其对稳定态与亚稳态的识别和探索能力。而最大口径方程(式25)为解决这一问题提供了一种实用、简便且高效的方案——若已知稳定态与亚稳态的布居数,且掌握一个或两个全局平均速率量,即可通过该方程估计出这些构象态之间的所有跃迁速率(见图13)。已有研究以一个七残基丙氨酸肽为模型系统,验证了该公式的估计精度:该系统的完整速率矩阵已通过大量精确计算获得,而最大口径理论的预测结果与之一致(40)。
式25也被成功应用于更大、更复杂的蛋白质构象变化过程——即G蛋白偶联受体(GPCR)在配体激活动力学中发生的别构跃迁(50)。应用式25的一个主要问题在于:如何从无偏动力学系综中确定动力学平均值 ⟨J⟩ ,并进而求解拉格朗日乘子 γ 。大多数用于确定平衡态能量景观的计算采样技术,均采用有偏系综,因此无法用于估计动力学量。近期,Meral等人(50)提出了一种巧妙的解决方案。他们利用了以下关键事实:在元动力学模拟中,可通过集合变量坐标(collective variable coordinate,CV),从有偏系综中估计出无偏的动力学平均值(51)。基于此,他们在马尔可夫态模型(Markov State Mode, MSM)中,同时实现了两项关键任务:估计平衡态分布,以及获得多个物理量的无偏动力学平均值。随后,他们将这些结果代入式25,成功估计出了构象跃迁速率。
分子模拟中反应坐标的确定
![]()
图14 最大口径原理可实现优质反应坐标(RC)的识别。复杂的势能景观(图 a)可被投影至任意反应坐标上;ΔG代表自由能(图 b、c)。最大口径原理使我们能够快速估计沿任意反应坐标的近似动力学,并识别出能实现时间尺度最大程度分离的反应坐标(图 d、e)。在本简单实例中,反应坐标 1(RC1)是优于反应坐标 2(RC2)的优质反应坐标。
生物分子变化的计算机模拟中,一个重要的挑战是找到主导反应坐标。构象空间具有高维性,模拟对其的采样往往十分稀疏;即便已知一个优质的反应坐标,若没有足够的采样以获得收敛的布居数,也无法得到沿该坐标的跃迁速率。近期,Tiwary及其团队(41, 42, 45)开发了一种巧妙的方法,利用最大口径方程(式25)来确定反应坐标。他们采用元动力学对目标过程进行模拟——该方法首先需要选择一些与当前研究问题相关的集体变量。其具体步骤如下:
第一步:估计沿任意一组集体变量线性组合的自由能剖面。
第二步:对于任意给定的反应坐标线性组合,利用最大口径方程(式25),在网格上估计出近似的速率矩阵。
第三步:将使速率矩阵的能隙(谱隙)达到最大值的线性组合,选为最优反应坐标(42)(见图14)。
值得注意的是,对线性组合的优化无需额外的模拟计算。这是因为,元动力学等增强采样模拟技术具有一项关键优势:若已知沿某一组线性组合的自由能剖面,即可直接估计出沿任意其他线性组合的自由能剖面(42)。
该方法所基于的核心原理如下:反应坐标通常描述的是大尺度运动,其速率慢于小尺度运动(如侧链旋转、溶剂分子的微小位移等)。因此,具有清晰分离的慢运动过程的路径,是优质反应坐标的理想候选。该方法已被扩展至多维反应坐标的确定(45),并在多个实例中得到了成功应用(43, 44)。
基于新数据的最大口径原理:马尔可夫模型的修正与更新
![]()
图15 生长因子激活通路的模型。受体的四个状态(①~④)定义如下:① 细胞表面的无配体结合受体(绿色)可被配体(黄色)结合,进入状态②;随后发生磷酸化,进入状态③。④ 所有状态的受体均可被内吞并降解,不同状态下的内吞降解速率存在差异。图中箭头代表跃迁速率;最大口径原理预测的变化最显著的跃迁速率 ,在图 (a) 中以灰色标注,在图 (b) 中以粉色标注。
我们考虑如下一个常见的实际问题:某模型网络的所有微观速率均已完成估计。但在后续研究中,可能会出现模型预测与实验数据不一致的情况——其原因可能是实验系统受到了扰动,或是初始的速率估计存在误差。例如,在蛋白质折叠模拟中,由于全原子力场的不精确性,计算得到的折叠速率可能与实验结果不符。采样不足也可能导致此类误差,这是生物分子模拟中一个众所周知的问题(52)。那么,如何对原有的速率矩阵进行修正呢?在大多数情况下,这个问题并没有唯一的解。而最大口径原理为我们提供了一种最优方案:它可以对完整的微观速率矩阵进行修正,使模型与新的实验数据达成一致(47, 53)。
若原始的计算速率为 { q ab } ,则修正后的、与新数据自洽的速率可表示为 { k ab } 。这些修正速率的求解方式为:最大化如下的相对熵
![]()
再举一个具体的实例(47):一种生长因子膜受体蛋白会经历一个四态生化循环,其状态包括:未结合配体态、配体结合态、磷酸化激活态以及降解态(见图15)。对于正常的野生型蛋白,这些状态之间的跃迁速率是已知的。
而当该受体发生某些突变时(例如,在某些癌症中出现的突变(54)),会导致激活态的布居数 p act 出现可观测的增加。我们希望利用这一单一的观测结果,来更新四态模型的预测布居数与跃迁速率。为实现这一目标,我们在激活态布居数的新值作为约束条件下,最大化相对路径熵。值得注意的是,尽管存在无穷多种更新速率矩阵以拟合该观测数据的方式,但该方法的预测结果显示:激活态布居数的增加,最有可能是通过降低受体的内吞速率实现的(47)。这一结论与生长因子信号通路中一个已被充分证实的异常现象完全一致(54, 55)。
最大口径原理的其他应用
最大口径原理已成功应用于鸟类群集行为(56)、细胞迁移运动(57)及神经元放电活动(10, 58)等研究领域。针对鸟类群集行为,Cavagna 等人(56)构建了一套最大口径理论框架,该框架可实现对连续两个时间步长下观测变量的关联分析。在利用模拟数据进行基准测试时,该模型的性能显著优于仅基于静态信息构建的模型。Tweedy 等人(57)则利用最大口径原理,建立了基于细胞形态的细胞迁移运动模型。他们通过分析细胞形态轨迹的时间演化,推断出对应的拉格朗日乘子。研究表明,这些推断得到的拉格朗日乘子能够有效区分经药物处理的细胞与未处理的对照细胞,以及基因修饰细胞与未修饰的亲本细胞。此外,最大口径原理还被用于捕捉神经元群体中复杂的临界动力学行为(10, 58)。当视网膜神经元暴露于自然图像刺激时,会表现出全或无的放电响应。已有研究利用最大熵原理(59),证实了神经元群体放电的统计特性中存在临界行为。这些研究发现,最大熵模型中被精细调节至临界态的拉格朗日乘子,其对应的物理模型为自旋玻璃模型。Mora 等人(10)将这些观测结果扩展至动力学范畴:他们通过对跨时间的神经元间关联施加约束,构建了神经元群体行为的动力学模型。研究结果表明,利用最大口径原理纳入神经元放电的动力学信息后,模型可预测出视网膜神经元同样处于动力学临界态。
最大口径原理自然导出了非平衡物理学的经典结果
为何马尔可夫模型在非平衡物理中如此普遍?
马尔可夫建模对广泛的动力学过程均具有良好的适用性。在一阶马尔可夫模型中,其核心假设为:仅需获知某一给定状态及其相邻动力学状态(即紧邻的前序与后续状态)的性质,即可对各状态的布居数及状态间的跃迁速率进行充分近似。该模型忽略了任何更长期的记忆效应。那么,为何马尔可夫模型在自然系统的建模中如此普遍且实用?最大口径原理为此提供了一种合理的解释。在各类建模方法中,最大口径原理是最大程度的数据驱动方法——其仅使用直接测量得到的可观测量。此处的核心观点在于:某一特定实验所能提供的数据性质,决定了能够捕捉这些数据的最优模型。例如,对于一个两态过程 A↔B ,若我们仅获知四个物理量——从 A 到 B 、 A 到 A 、 B 到 B 及 B 到 A 的跃迁频率——那么,使口径最大化的模型类别即为一阶马尔可夫模型(60-63)。除非我们拥有更多的信息,否则引入更多参数的模型(例如高阶马尔可夫模型)均缺乏合理的依据。通过在上述给定约束集下最大化口径,可直接推导出如下结论:连续两个时间步长内,两状态间的跃迁速率仅依赖于前一个时间步的状态。当数据涉及两个以上的时间步(即系统存在记忆效应)时,Lee与Pressé(62)给出了严格的数学描述。
近平衡统计物理的已知结论可由最大口径原理推导
最大口径原理可推导出近平衡统计物理中一系列著名的结论,这些结论建立了通量、流与熵产生之间的关联。其中包括格林-久保关系(Green–Kubo relationship)、昂萨格倒易关系(Onsager’s reciprocal relationship) 以及普里戈金最小熵产生原理(Prigogine’s minimum entropy production principle)(64)。例如,继Thomson关于电流与热流耦合的实验之后,昂萨格考虑了两种流的耦合问题——即两种流(记为 a 与 b )在两个热库之间的流动。若通量 J a 与 J b 与驱动力( λ a , λ b )呈线性关系,则其耦合形式可表示为:
![]()
昂萨格基于微观理论提出了倒易关系的论证,即 L ab = L ba 。而利用最大口径原理,该关系的推导过程可被极大简化。首先,若对于任意给定轨迹 Γ ,任意时刻的粒子通量为 j a ( Γ ,t) 与 j b ( Γ ,t) ,则最大口径原理预测的轨迹分布为:
![]()
其中, Q d 为动力学配分函数, λ a (t) 与 λ b (t) 为拉格朗日乘子, p( Γ ) 为平衡态分布。由于最大口径原理是一种基于配分函数的方法,我们可直接得到其一级与二级偏导数:
![]()
与平衡态热力学中的麦克斯韦关系(Maxwell relations)类似,此处得到的二级偏导数与求导顺序无关。由此,昂萨格倒易关系可直接推导得出(推导细节参见参考文献64)。
这些耦合关系的应用范围十分广泛,不仅适用于热电学或功与热的转换,还可用于描述:生化能量源如何驱动化学反应、分子马达与离子泵的运行;能量源如何提高生物分子钟与校对机制的精度;以及光伏材料如何实现光与电流或热的耦合。
本综述的问题、批评、展望与局限性
本节将提供更广阔的研究背景。但首先需要说明的是,受篇幅限制,我们无法对非平衡统计物理中其他重要、活跃且相关的研究方向进行综述,例如随机热力学(stochastic thermodynamics)(65)与大偏差理论(large-deviation theory)(66)。
最大口径原理能否正确处理耗散过程?
最大口径原理能否正确处理耗散过程?(67, 68)若 Γ 表示一条轨迹,且仅施加单一约束 ⟨J⟩ ,则最大口径原理预测的路径布居数为:
![]()
但我们可以设想一个具体的实例:一个粒子在粘性流体中运动。其平均速度可通过两种不同的方式实现:在高粘度介质中施加一个较大的力,或在低粘度介质中施加一个较小的力。我们有理由认为,这两种情况下的路径分布应当是不同的。然而,式(32)却暗示这两种路径分布不存在差异。这是否意味着最大口径原理的失效?答案是否定的。事实上,这两种情况下的路径分布确实是不同的。但对于此类耗散系统,仅施加单一的 ⟨J⟩ 约束是不充分的——在这类系统中,能量输入与输出的速率 同样是至关重要的。对于这种情况,我们需要引入如下的分布形式(69):
正如在一个精确模型中所展示的那样,这一额外的约束将导致正确的速率分布(69)。然而,需要注意的是,最大口径原理本身并不会指定某一特定问题需要哪些约束。这是建模者需要做出的决策,其取决于该特定问题中哪些物理量是相关的。这一挑战与平衡态热力学中的情况类似:在体积与浓度同时发生变化的情况下,仅指定温度是不充分的;此时,还必须同时指定压力与化学势。
从数据中学习拉格朗日乘子的数值问题
将最大口径原理应用于含噪声的生物数据时,往往会面临数值计算的挑战。利用数据同时拟合多个拉格朗日乘子的计算成本可能极高(70, 71)。数据拟合过程可能需要同时求解 N 个非线性方程,而这些方程通常并非相互独立。此外, ⟨J⟩ 如何依赖于统计权重的解析表达式往往是不可得的。同时,实验数据通常也会包含实验误差。因此,从数据中获取拉格朗日乘子的过程,可能需要对耦合的非线性方程组进行随机采样。这导致我们有时无法将约束条件确定为精确的固定值,而只能将其表示为拉格朗日乘子的可能分布(35)。最大口径原理应用中的另一项挑战在于状态空间的确定。正如在基因网络或生物分子的马尔可夫模型中所展示的那样,我们通常需要先验地指定一个状态空间。而更复杂的模型则无需预先指定该空间(72, 73)。
最大口径原理是一种推断原理,还是一种物理原理?
综上所述,最大口径原理是一种通用的方法,用于在动力学过程的模型中推断速率与路径的分布。给定一个模型,以及有限的数据集——例如几个平均速率或其他任意矩——最大口径原理可预测出与该模型、数据及概率论规则均自洽的分布。我们认为,这一思想与平衡态统计力学推断分布的方式在本质上是完全相同的。我们将整个统计力学视为对物理模型进行推理的过程。
一个相关的问题是:统计力学能否从力学中推导出来以及其是否比单纯的推理更为深刻?我们的观点是,在实际应用中,这是无法实现的。虽然热力学第一定律关乎能量与力学,其建立在物理量的基础之上;但热力学第二定律则关乎布居数——因此,其本质是关于推理或概率论的。我们无法从第一定律推导出第二定律。并非所有的物理规律都能仅从纯粹的力学中推导出来。借用E.T. Jaynes的表述:玻尔兹曼的过人之处在于,他认识到尽管气体的行为在原则上可以通过追踪所有弹子球式的碰撞过程来确定,但计算气体性质的唯一实用方法,是用统计描述来替代详细的力学过程——这正是“统计力学”一词的由来(74)。一旦我们接受了玻尔兹曼的这一智识飞跃,并采用熵的表达式 S=k ln W 或 S=−∑ p i ln p i ,我们就必然将熵的变化视为从数据中对模型进行推理的任务。
尽管如此,物理推理与非物理推理之间仍存在差异。若我们仅获知平均通量,最大口径原理可以推断出其分布,但无法提供更多的信息。然而,若我们在最大口径原理中使用的模型更具物理性,我们则可以获得更多的信息——例如,速率如何依赖于粒子的性质或流网络的结构。此外,在哈密顿动力学适用的情况下,其可以提供额外的机制性洞见,将力与流与底层分子的性质关联起来。而且,在各类不同的约束中,温度在统计热力学中占据着特殊的地位。热平衡态具有一个动力学所不具备的实验真值,即克劳修斯关系(Clausius relationship) S Clausius = q/ T = ⟨U⟩/ T ,该关系可通过平均能量确定熵。但克劳修斯关系仅在平衡态下成立,且还受到其他多种限制性条件的约束。目前,对于远离平衡的动力学过程,我们尚不知晓其对应的实验真值基础。
最大口径原理如何解决非平衡物理的关键难题?
在此,我们简要总结最大口径原理如何解决非平衡物理中的典型挑战与问题。
最大口径原理作用于轨迹(路径),而非状态
最大口径原理模型可以包含所有相关的路径,包括那些远离平衡态的路径。Shore与Johnson(15)提出了一个重要但未被充分重视的观点:函数 −∑ p i ln p i 具有两个关键性质:
a. 它是唯一能保证概率论规则自洽性的函数;
b. 它仅对使自身最大化的单一分布函数 具有有效的预测能力。
最大口径原理中所引入的唯一路径熵为: ,
其中 为使该熵最大化的路径概率。我们无需考虑该分布的任何偏差;因此,它满足Shore-Johnson关于自洽概率推理的判据(15, 63, 75)。
其适用范围不限于近平衡态
最大口径原理不需要以连续函数为起点——例如,状态熵 S state =S(U,V,N) ,该熵是空间与时间的广延变量的连续可微函数,且仅通过克劳修斯关系在平衡态下被严格定义。为了保证这种光滑性,局域平衡假设是必然的推论。这一假设具有很强的限制性,它将研究范围局限于那些仅发生微小步骤、并在过程中不断达到平衡的过程(76)。
其适用于非热学系统
最大口径原理不仅限于热学过程或温度热库,因此,它可轻松应用于广泛的流动问题。它对模型或数据的来源不做任何假设。其适用范围不仅限于材料、分子、分子碰撞、哈密顿系统、刘维尔定理、热库或温度。最大口径原理具有更广泛的普适性,可应用于随机动力学系统、基因电路、网络及其他各类系统。
最大口径原理为非平衡物理的不同类别提供了合理化的依据
非平衡过程通常被划分为若干广泛的类别,例如:平衡态、近平衡弛豫过程(例如,由施加于反应网络的条件变化或从高浓度到低浓度的扩散所导致的弛豫)、近平衡定态(通过电阻的恒定欧姆电流,或球体在粘性液体中的缓慢拖曳),或远离平衡态(驱动球体在流体中产生湍流)。根据在最大口径原理的构建中哪些约束是相关的,这些分类可以被自然地定义。当然,平衡态需要满足细致平衡条件,且不施加任何额外的速率约束。近平衡过程是一类耗散过程,但其耗散与固定外力下的某一流率 ⟨J⟩ 呈线性比例关系。例如,欧姆热耗散与电压×电流成正比;另一个例子是斯托克斯定律,即球体在粘性液体中的耗散与力×速度成正比。在这些情况下,近平衡态可仅通过指定相关的速率 ⟨J⟩ 来定义——因为耗散与流的线性比例关系,无需施加额外的耗散约束。近平衡态具有以下一项或多项相应的特征:(a)线性的力-流关系;(b)熵产生与耗散的对应关系;(c)局域平衡假设的适用性;以及(d)格林-久保关系、昂萨格倒易关系与普里戈金最小熵产生原理的成立(64)。相比之下,远离平衡的过程可能需要额外的信息来同时描述其宏观与微观性质——例如,用于描述超 ⟨J⟩ 所指定耗散的额外耗散模型,或无法被解释为热力学强度量简单梯度的力约束。
总结:最大口径原理
是一种适用于动力学与路径的通用变分原理
最大口径原理是一种通用原理,用于在给定有限数据的情况下,推断动力学模型中速率与路径的分布。它可以推导出著名的近平衡态结论,并能区分近平衡态与远离平衡态,同时其在远离平衡态的情况下依然有效。它克服了传统非平衡物理所面临的诸多问题。结合模型,它可以推导出唯象定律,例如菲克定律。它适用于少粒子复杂系统,例如基因电路。其核心逻辑对于构建数据所允许的最简模型具有显著的优势。尽管其在所有情境下的普适性尚未被完全证明,但所有的迹象均表明,它是一种适用于非平衡态的通用原理。
参考文献
Broder A 2000. Graph structure in the web. Comput. Netw. 33:309–320
Zipf GK 1949. Human Behavior and the Principle of Least Effort. Cambridge, MA: Addison-Wesley
Marsili M, Maslov S, Zhang YC 1998. Dynamical optimization theory of a diversified portfolio. Phys. A Stat. Mech. Appl. 253:403–418
Cavagna A, Queirós SD, Giardina I, Stefanini F, Viale M 2013. Diffusion of individual birds in starling flocks. Proc. R. Soc. B Biol. Sci. 280:20122484
Peterson J, Dixit P, Dill KA 2013. A maximum entropy framework for nonexponential distributions. PNAS 110:20380–20385
Zeldovich K, Chen P, Shakhnovich E 2007. Protein stability imposes limits on organism complexity and speed of molecular evolution. PNAS 104:16152–16157
Zou T, Williams N, Ozkan S, Ghosh K 2014. Proteome folding kinetics is limited by protein halflife. PLOS ONE 9:e112701
Dixit PD, Maslov S 2013. Evolutionary capacitance and control of protein stability in protein-protein interaction networks. PLOS Comput. Biol. 9:e1003023
Dixit PD, Pang TY, Maslov S 2017. Recombination-driven genome evolution and stability of bacterial species. Genetics 207:281–295
Mora T, Deny S, Marre O 2015. Dynamical criticality in the collective activity of a population of retinal neurons. Phys. Rev. Lett. 114:078105
Onsager L 1931. Reciprocal relations in irreversible processes. Phys. Rev. 37:405
Prigogine I, Kondepudi D 1998. Modern Thermodynamics: From Heat Engines to Dissipative Structures. New York: John Wiley
Martyushev LM, Seleznev VD 2006. Maximum entropy production principle in physics, chemistry and biology. Phys. Rep. 426:1–45
Jaynes ET 1980. The minimum entropy production principle. Annu. Rev. Phys. Chem. 31:579–601
Shore J, Johnson R 1980. Axiomatic derivation of the principle of maximum entropy and the principle of minimum cross-entropy. IEEE Trans. Inf. Theory 26:26–37
Wu D, Ghosh K, Inamdar M, Lee H, Fraser S et al. 2009. Trajectory approach to two-state kinetics of single particles on sculpted energy landscapes. Phys. Rev. Lett. 103:050603
Pressé S, Ghosh K, Phillips R, Dill KA 2010. Dynamical fluctuations in biochemical reactions and cycles. Phys. Rev. E 82:031905
Ghosh K, Dill K, Inamdar M, Seitaridou E, Phillips R 2006. Teaching the principles of statistical dynamics. Am. J. Phys. 74:123–133
Seitaridou E, Inamdar M, Phillips R, Ghosh K, Dill K 2007. Measuring flux distributions for diffusion in the small-numbers limit. J. Phys. Chem. B 111:2288–2292
Firman T, Balazsi G, Ghosh K 2017. Building predictive models of genetic circuits using the principle of maximum caliber. Biophys. J. 113:2121–2130
Nevozhay D, Adams R, Itallie EV, Bennett M, Balázsi G 2012. Mapping the environmental fitness landscape of a synthetic gene circuit. PLOS Comput. Biol. 8:e1002480
Firman T, Wedekind S, McMorrow TJ, Ghosh K 2018. Maximum caliber can characterize genetic switches with multiple hidden species. J. Phys. Chem. B 122:5666–5677
Gardner T, Cantor C, Collins J 2000. Construction of a genetic toggle switch in Escherichia coli. Nature 403:339–342
Pressé S, Ghosh K, Dill K 2011. Modeling stochastic dynamics in biochemical systems with feedback using maximum caliber. J. Phys. Chem. B 115:6202–6212
Elowitz M, Leibler S 2000. A synthetic oscillatory network of transcriptional regulators. Nature 403:335–338
Firman T, Amgalan A, Ghosh K 2019. Maximum caliber can build and infer models of oscillation in a three-gene feedback network. J. Phys. Chem. B 123:343–355
Van Kampen NG 1992. Stochastic Processes in Physics and Chemistry (Vol. 1). Amsterdam: Elsevier
Gillespie D 1992. A rigorous derivation of the chemical master equation. Phys. A Stat. Mech. Appl. 188:404–425
Munsky B, Khammash M 2006. The finite state projection algorithm for the solution of the chemical master equation. J. Chem. Phys. 124:044104
Frigola D, Casanellas L, Sancho J, Ibanes M 2012. Asymmetric stochastic switching driven by intrinsic molecular noise. PLOS ONE 7:e31407
Llamosi A, Gonzalez-Vargas AM, Versari C, Cinquemani E, Ferrari-Trecate G et al. 2016. What population reveals about individual cell identity: single-cell parameter estimation of models of gene expression in yeast. PLOS Comput. Biol. 12:e1004706
Raj A, van Oudenaarden A 2008. Nature, nurture, or chance: stochastic gene expression and its consequences. Cell 135:216–226
Wu M, Singh AK 2012. Single-cell protein analysis. Curr. Opin. Biotechnol. 23:83–88
Saliba AE, Westermann AJ, Gorski SA, Vogel J 2014. Single-cell RNA-seq: advances and future challenges. Nucleic Acids Res. 42:8845–8860
Dixit PD, Lyashenko E, Niepel M, Vitkup D 2019. Maximum entropy framework for predictive inference of cell population heterogeneity and responses in signaling networks. Cell Syst. https://doi.org/10.1016/j.cels.2019.11.010 [Crossref] [Medline] [Web of Science]
Dixit PD 2013. Quantifying extrinsic noise in gene expression using the maximum entropy framework. Biophys. J. 104:2743–2750
Dixit PD, Dill KA 2014. Inferring microscopic kinetic rates from stationary state distributions. J. Chem. Theory Comput. 10:3002–3005
Wan H, Zhou G, Voelz VA 2016. A maximum-caliber approach to predicting perturbed folding kinetics due to mutations. J. Chem. Theory Comput. 12:5768–5776
Zhou G, Pantelopulos GA, Mukherjee S, Voelz VA 2017. Bridging microscopic and macroscopic mechanisms of p53-MDM2 binding with kinetic network models. Biophys. J. 113:785–793
Dixit PD, Jain A, Stock G, Dill KA 2015. Inferring transition rates of networks from populations in continuous-time Markov processes. J. Chem. Theory Comput. 11:5464–5472
Tiwary P, Berne B 2016. How wet should be the reaction coordinate for ligand unbinding?. J. Chem. Phys. 145:054113
Tiwary P, Berne B 2016. Spectral gap optimization of order parameters for sampling complex molecular systems. PNAS 113:2839–2844
Tiwary P 2017. Molecular determinants and bottlenecks in the dissociation dynamics of biotin–streptavidin. J. Phys. Chem. B 121:10841–10849
Hovan L, Comitani F, Gervasio FL 2018. Defining an optimal metric for the path collective variables. J. Chem. Theory Comput. 15:25–32
Smith Z, Pramanik D, Tsai ST, Tiwary P 2018. Multi-dimensional spectral gap optimization of order parameters (SGOOP) through conditional probability factorization. J. Chem. Phys. 149:234105
Dixit PD, Dill KA 2019. Building Markov state models using optimal transport theory. J. Chem. Phys. 150:054105
Dixit PD 2018. Communication: Introducing prescribed biases in out-of-equilibrium Markov models. J. Chem. Phys. 148:091101
Dixit PD 2018. Entropy production rate as a criterion for inconsistency in decision theory. J. Stat. Mech. Theory Exp. 2018:053408
Dixit PD 2019. Introducing user-prescribed constraints in Markov chains for nonlinear dimensionality reduction. Neural Comput. 31:980–997
Meral D, Provasi D, Filizola M 2018. An efficient strategy to estimate thermodynamics and kinetics of G protein-coupled receptor activation using metadynamics and maximum caliber. J. Chem. Phys. 149:224101
Tiwary P, Parrinello M 2014. A time-independent free energy estimator for metadynamics. J. Phys. Chem. B 119:736–742
Sawle L, Ghosh K 2016. Convergence of molecular dynamics simulation of protein native states: feasibility versus self-consistency dilemma. J. Chem. Theory Comput. 12:861–869
Dixit PD, Dill KA 2018. Caliber corrected Markov modeling (C2M2): correcting equilibrium Markov models. J. Chem. Theory Comput. 14:1111–1119
Herbst RS 2004. Review of epidermal growth factor receptor biology. Int. J. Radiat. Oncol. Biol. Phys. 59:S21–S26
Huang HJS, Nagane M, Klingbeil CK, Lin H, Nishikawa R et al. 1997. The enhanced tumorigenic activity of a mutant epidermal growth factor receptor common in human cancers is mediated by threshold levels of constitutive tyrosine phosphorylation and unattenuated signaling. J. Biol. Chem. 272:2927–2935
Cavagna A, Giardina I, Ginelli F, Mora T, Piovanni D et al. 2014. Dynamical maximum entropy approach to flocking. Phys. Rev. E 89:042707
Tweedy L, Witzel P, Heinrich D, Insall RH, Endres RG 2019. Screening by changes in stereotypical behavior during cell motility. Sci. Rep. 9:8784
Vasquez JC, Marre O, Palacios AG, Berry MJ II, Cessac B 2012. Gibbs distribution analysis of temporal correlations structure in retina ganglion cells. J. Physiol. Paris 106:120–127
Tkačik G, Mora T, Marre O, Amodei D, Palmer SE et al. 2015. Thermodynamics and signatures of criticality in a network of neurons. PNAS 112:11508–11513
Stock G, Ghosh K, Dill K 2008. Maximum caliber: a variational approach applied to two-state dynamics. J. Chem. Phys. 128:194102
Ge H, Pressé S, Ghosh K, Dill K 2012. Markov processes follow from the principle of maximum caliber. J. Chem. Phys. 134:064108
Lee J, Pressé S 2012. A derivation of the master equation from path entropy maximization. J. Chem. Phys. 137:074103
Pressé S, Ghosh K, Lee J, Dill K 2013. Principle of maximum entropy and maximum caliber in statistical physics. Rev. Mod. Phys. 85:1115–1141
Hazoglou MJ, Walther V, Dixit PD, Dill KA 2015. Communication: Maximum caliber is a general variational principle for nonequilibrium sta...
特别声明:以上内容(如有图片或视频亦包括在内)为自媒体平台“网易号”用户上传并发布,本平台仅提供信息存储服务。
Notice: The content above (including the pictures and videos if any) is uploaded and posted by a user of NetEase Hao, which is a social media platform and only provides information storage services.