![]()
導語
長期以來,統計物理中的變分原理為平衡態體系提供了堅實的理論基礎,但非平衡態體系始終缺乏一個同等普適且自洽的理論框架。隨著研究對象從宏觀熱學系統拓展至少顆粒、生物與信息系統,傳統假設逐漸失效。本文系統綜述最大口徑原理這一基于路徑熵最大化的非平衡統計方法,闡明其理論基礎、數學形式與操作流程,并通過擴散、電流分配及基因電路等實例,展示其在復雜少顆粒動力學建模中的獨特優勢與廣泛應用前景。
關鍵詞:非平衡統計物理、最大口徑原理(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.