基于马尔可夫状态模型（MSMs）实现化学主方程（CME）谱特性分析和可视化

创作时间:

作者:

@小白创作中心

基于马尔可夫状态模型（MSMs）实现化学主方程（CME）谱特性分析和可视化

引用

CSDN

https://blog.csdn.net/j_jinger/article/details/145330248

化学动力学领域的核心挑战之一在于理解和预测复杂反应网络中发生的随机行为。传统的确定性模型，如常微分方程（ODE），虽然在描述平均行为方面取得了显著成功，但却无法捕捉由分子随机性引起的涨落和多稳态现象。为了更准确地刻画此类复杂体系，化学主方程（CME）应运而生。然而，CME的求解，尤其对于具有大量物种和状态的大规模系统而言，常常面临着计算上的巨大挑战。近年来，基于马尔可夫状态模型（MSMs）的建模方法为CME的谱特性分析和可视化提供了新的视角，使得理解随机多稳态和状态转换成为可能。本文旨在深入探讨基于MSMs的CME分析方法，并阐述其在揭示复杂化学系统动力学特征方面的优势。

化学主方程（CME）与随机动力学

化学主方程（CME）是一种描述离散且随机的化学反应网络的偏微分方程。它关注的是系统处于特定状态的概率随时间变化的动态，而非物种浓度的平均值。CME的表达式通常涉及到状态概率的变化率，其由一系列由反应速率决定的概率转移项构成。具体而言，CME的推导基于以下假设：反应是离散的，物种数量只能取非负整数值，且系统的演化是马尔可夫过程，即当前状态的概率只依赖于上一时刻的状态，而与更早的历史无关。

CME的强大之处在于它能够自然地捕捉由分子涨落带来的随机性效应，尤其是在低分子数的情况下。然而，CME的求解本身就具有很高的计算复杂度。对于具有多个物种和多种反应的系统，状态空间会呈指数级增长，使得CME的精确求解变得不可行。因此，寻找有效的近似方法来分析和理解CME的动力学行为至关重要。

马尔可夫状态模型（MSMs）及其构建

马尔可夫状态模型（MSMs）提供了一种将复杂的连续状态空间离散化的方法。MSMs的核心思想是将状态空间划分为若干个离散的宏观状态，并通过估计这些状态之间的转移概率来构建一个马尔可夫链。与CME直接描述微观状态的演化不同，MSMs关注的是宏观状态之间的转移，从而大大降低了模型的复杂度。

构建MSMs通常需要以下几个步骤：

状态空间采样：首先，需要从系统中采样大量的状态轨迹。这可以通过分子动力学模拟、随机模拟（如Gillespie算法）或其他方法来实现。
状态聚类：将采样得到的轨迹聚类到若干个离散的宏观状态。聚类算法的选择取决于具体的系统和需求，常用的方法包括k-means聚类、谱聚类等。
转移概率估计：基于采样轨迹，估计不同宏观状态之间的转移概率。通常采用最大似然估计或贝叶斯方法来完成这项任务。
模型验证：构建的MSMs需要进行验证，以确保其能够准确地描述原始系统的动力学行为。这可以通过比较MSMs预测与原始模拟结果来完成。

基于MSMs的CME谱特性分析

通过构建MSMs，我们可以将CME的分析转化为对离散马尔可夫链的分析。对于马尔可夫链，我们可以研究其转移矩阵的谱特性，如特征值和特征向量。这些谱信息揭示了系统动力学的关键特征。

特征值：MSMs的转移矩阵的特征值与系统的弛豫时间尺度密切相关。模值接近于1的特征值对应于缓慢的动力学过程，如多稳态之间的状态转换，而模值远离1的特征值则对应于快速的动力学过程，如同一宏观状态内的涨落。特别是，特征值为1对应的特征向量代表稳态分布，揭示了系统在长时间尺度上的平均行为。
特征向量：与特征值相对应的特征向量携带了关于状态之间转移路径和概率分布的信息。例如，与最大特征值（1）对应的特征向量描述了系统的稳态概率分布，而与次大特征值对应的特征向量则可能揭示了不同稳态之间的转换路径。

通过分析MSMs的谱特性，我们可以提取出CME中隐藏的动力学信息，从而更好地理解系统的随机行为。

基于MSMs的可视化随机多稳态和状态转换

除了谱特性分析之外，MSMs还可以用于可视化随机多稳态和状态转换。

稳态分布可视化：通过MSMs的稳态分布（与特征值为1对应的特征向量），我们可以可视化系统在不同宏观状态下的概率分布。在多稳态系统中，稳态分布会呈现出多个峰值，每个峰值对应于一个稳态，并且峰值的高度反映了该稳态的稳定性。
状态转换路径可视化：基于转移矩阵和特征向量，我们可以构建状态转换网络，并可视化不同状态之间的转移路径。例如，我们可以根据转移概率的大小绘制状态之间的连线，并根据特征向量的信息标记出状态的稳定性和过渡状态。
自由能景观重构：利用MSMs的稳态分布，可以通过Boltzmann反演方法重构系统的自由能景观。自由能景观能够直观地展现多稳态系统的势能面，从而理解状态转换的动力学过程。

部分代码

%SSA sampling for model A (toggle switch)
%generate N_rows SSA trajectories of lengths #cells * N_vert * tau
%and compute some statistics like the weights of the clusters
N_ssa=0;
tau=5000;
N_rows=10;
N_vert=500;
num_cells=60;
%num_trajectories = num_cells * N_vert;
%tf=num_trajectories*tau;
UL=zeros(1,N_rows);
LR=zeros(1,N_rows);
UL_rel=zeros(1,N_rows);
LR_rel=zeros(1,N_rows);
diagonal=zeros(1,N_rows);
crossings=zeros(1,N_rows);
out_count=zeros(1,N_rows);
num_hor=zeros(1,N_rows);
Pii=zeros(1,N_rows);
cell=0;
voronoi_table=[];
N0=[];
for n=1:N_rows
n
x0=[1;1];
X_total=[];
%break sampling down in smaller parts to save some intermediate results
for k=1:15
%save_name = sprintf('SSA_results/toggle_reference_SSA_results');
%save(save_name,'UL','num_hor','out_count','crossings')
figure(1)
boxplot([UL_rel',LR_rel'])
figure(2)
boxplot([Pii',Pii'])