蛋白质的磷酸化和去磷酸化
这篇笔记是要把前面的理论用到分子生物学与细胞生物学中。
生物系统中常见的信号开关如磷酸化-去磷酸化环（PdPC）与G蛋白偶联信号系统。它们能够产生超灵敏度。下面考察PdPC系统。
在一定的简化假设下（激酶和磷酸酯酶的浓度保持不变），PdPC的动力学可以用以下ODE简化表示：
其中x表示激活（磷酸化）的信号蛋白的浓度。这一动力系统会呈现出鞍结点分叉，其分岔图为（下面都把k_1作为分岔参数）
也就是说改变参数时，这一系统不仅是“超灵敏”，甚至出现了双稳态，直接从低态跳向高态。这是生物系统中非常常见而又不平凡的现象。
不过，如果仅有这样一个ODE的话，我们并不知道跳跃的这个相变点在什么地方，也不知道具体的分布是怎样的。这就必须要请出随机的化学主方程模型。不过我们更加喜欢的是Ito扩散过程，所以不如直接把化学主方程作Kramers-Moyal展开变成随机微分方程比较好（很多Remarks：Kramers-Moyal展开到二阶，当然有Gillespie的一个基于化学计量学矩阵的直接的公式，不过Kramers-Moyal总是可以用的；必须要展开到二阶才有Fokker-Planck方程的形式；也只能展开到二阶或者无穷多阶， Pawula theorem说明其他阶都无法保证合理的概率密度函数；展开的时候把R跟p放在一起展开，出来的结果自然是向前Kolmogorov方程，如果出来的是向后Kolmogorov方程，就说明展开搞错了）。近似的SDE结果为：
其中X_t表示浓度。V正好就是大偏差参数。这个式子同时吻合涨落的量级O(N^{-1/2})。不过需要注意的是它不是常系数的噪声，分析起来会有点麻烦，比如没法直接套用exit problem的结论。
下面考察其平稳分布的大偏差（它跟escape time的大偏差是热力学和动力学的两个方面）。求解Fokker-Planck方程得到平稳分布为
这已经是一个LDP的形式，速率函数称为非平衡态景观函数（但是这里其实是平衡态？不过下面的结果会说明它对非平衡态也成立），它有两个稳定点，在V\rightarrow\infty的极限下，只有其中比较小的那个留下来。所以由这两个点的高低可以确定一阶相变的相变点，这是ODE模型做不到的。
然而问题在于，这个方式得到的rate function跟后面直接基于离散状态Markov链（生灭过程）得到的rate function形式上是不一样的。虽然很容易证明它们有同样的极值点，但是机值点处的值并不相同，最严重的情况是双稳态的两个最低点的高低关系都有可能不同（存在这样的例子）。这里的问题在于Kramers-Moyal的扩散近似中丢失了一些东西，参考文献中说“Van Kampen repeatedly emphasized that the Fokker–Planck approximation can only be obtained for master equations with small individual jumps”，如果想要更精确地刻画这个体系，还是要从最原始的CME或者叫生灭过程出发。这是一个一维的生灭过程，解析地做起来还是非常方便的。而且下面的分析对于一般的一维生灭过程也都适用。首先求出不变分布：
它可以直接近似为LDP形式
其中的rate function为
它的确与前面得到的速率函数不一样。我们来看看这个速率函数（非平衡态景观函数）有什么性质。最重要的一点，它就是这个化学反应体系的“势”，即一个Lyapunov函数。一个确定性的反应必然沿着这个Lyapunov函数向下走。证明如下：
只有在landscape的极值点，即R^+=R^-的时候，它才会停止，系统处于定态。对于现在这个磷酸化-去磷酸化环，可以画出某一参数下的landscape的图像为
从这个landscape的角度理解这个化学反应体系的动力学性质。在短时间内，系统在某一个稳定点附近形成一个亚稳定的Gauss分布（弛豫）。在很长的时间尺度内，可以看到两个稳定点之间的跳跃（Markov，就像经典的离子通道信号一样），在更长的时间尺度内，系统处于不变分布。
对于这个信号开关系统来说，上面的性质说的也就是：
信号蛋白（的浓度）有激活态（大量激活）和未激活态（少量激活）两个状态，二者是分隔开来的，不能通过连续的变化把二者连接起来；
短时间内，系统处于某一个稳定态附近，并且有一个Gaussian的涨落，这是CLT的一种体现，也和常识相符；同时也是rate function展开到二阶的结果；
改变某个参数（上面是k_1），一开始蛋白只能在未激活态，经过一次鞍结点分叉，激活态和未激活态都是稳定的，此时系统在长时间尺度内会在激活态和未激活态之间不断跳跃；但是在V趋于无穷时，只有一个态能够保留下来；
随着参数的改变，最终保留下来的态会有变化，这是一个相变；
两个态之间的跳跃是一个Markov过程，可以用两状态jump process来刻画，这是对这个体系在长时间尺度内的约化；
实际上我们很自然的会选择用Markov链给后面这个约简过程建模。而前面的论证说明了这一建模从底层来说就是合理的。
下面考察在两个稳定点之间跃迁的速率。我们不套用Freidlin-Wentzell理论的结果，而是从Laplace方法推导。这也是一个非常有大偏差意味的工具。
所谓Laplace方法就是说，对于一个很大的M（可以想象成大偏差中的参数），可以对大偏差形式的积分做一个近似（x_0为唯一全局极大值点，不在边界上）：
右边的指数部分其实相当于大偏差中收缩原理的思想；而指前因子相当于Eyring-Kramers理论中的指前因子（形式也长得差不多）。
Laplace方法的一个经典应用就是证明Stirling公式：
右边直接套用Laplace公式得到
Laplace方法有一个简单的推广，在下面会用到
回到原来的问题。对于一个jump process，其首达时的期望值是可以通过列线性方程直接解析求解出来的（其实就是Dynkin方程的离散情况，可以写成Au=-1，也可以按照离散时间马氏链的直观“跳步法”理解）。考虑前面的生灭过程，从n_1^*跑到n_2^*。列出线性方程：
强行递归得到（用平稳分布简化表达）
根据前面的连续近似的p(x)表达式，可将其化简为
它可以用Laplace的思想化简。不严格地说，外面这个积分大概在势能最高点x_3^*处贡献最大，所以可以限定积分在x_1^*-x_3^*，而里面这个积分在这个范围内，有唯一一个极小值点x_1^*，所以里面可以先用一次Laplace，外面再用一次Laplace，最终结果为
这其实是类似于重新推倒了一遍上一篇文章的公式，不过是对于jump process来推的，所以是之前没有的新的结果（这个形式一看就是在x_1和x_3用了两次Laplace）。而且这种Laplace方法是一般的，可以用于任意这种形式的近似。
上面这个速率公式有着非常本质的意义。下面作分析：
这是以前从来无法想象的结果。如果遵循认知发展的道路：从传统的化学家角度看，一个化学反应体系，如果是闭的，就不可能有双稳态；然而生物是开化学反应体系，于是出现了双稳态、极限环。对于熟知ODE理论的人来说，双稳态是确定的，落入哪一个稳态只取决于初始位置，所以这个图像中不可能有两个稳态之间的跳跃；人们也可能觉得这种跳跃的rare events是在太过稀少，没有实际的物理意义。然而从随机的角度看，我们现在能非常定量地刻画这种跳跃：知道了它具体满足怎样的指数分布。
抛开这个具体的例子，我们一般地去想象一个化学反应，非平衡态景观函数到底做了些什么？系统经过一定时间弛豫到一个稳定状态，这个状态的浓度有一个分布，从这个分布可以定义出一个landscape，也叫做非平衡态景观函数。乍看起来它只是一个平衡态的东西，类似于某种势能，浓度按照这个势能有一个Boltzmann分布；它作为一个大偏差速率函数，在最低点的展开能够给出稳态分布的中心极限定理，这也跟普通的“误差累加造成Gauss分布”的常识吻合。但是这个非平衡态景观函数出乎意料地展现了很多非平衡态的性质。首先，它是确定性动力学的Lyapunov函数，也就是说，它的确是一种“势”，而不是强行定义出来的。我在上一篇文章中讨论的问题是dX_t=-\nabla V(X_t)+\sqrt{2\epsilon}dB_t，这是一个明确的“势场”中的随机游动；可是对于化学反应的CME来说，并没有先天存在的这样一个V(x)（即使在Kramers-Moyal展开下可以积分强行构造出一个V(x)，但是这个V(x)不满足Boltzmann分布也不能用于Eyring方程，说到底是因为那个dB_t部分不是常数的影响，如果强行构造的话就是丢失了化学反应intrinsic的噪声特性），这是个很大的不同。可是我们发现平稳分布的非平衡态景观函数就担当了这个V(x)的角色。接下来跟上一篇文章中一样，这个非平衡态景观函数还恰好可以用来讨论双稳态之间的跳跃，并且有一个精确的Arrhenius-type（Eyring-type）的公式。
进一步，虽然没有证明，但是根据前面的Freidlin-Wentzell理论，我们至少可以猜测这个分布是指数分布。所以我们可以写出长时间尺度约化的马氏链的Q矩阵来。这就是一个从化学诞生之日就开始用，却一直到现在才被阐明的东西。
数值模拟
下面用这个PdPC做数值模拟，验证一下上面的公式。首先理一理需要做的事情：
画出样本轨道，跟确定性动力学对比；
在弛豫时间尺度上，看其分布是否是Gauss；
在长时间尺度上，看其分布是否遵循landscape；
统计跳跃时间，看其均值与分布是否与预测相同；
用标准Monte-Carlo方法与Doob-Gillespie方法做数值模拟。
还是上面的模型，k_1取40。则可以具体写出两个速率为
其两个稳定不动点为0.0609与0.5923，中间的不稳定不动点为0.3466。
其非平衡态景观函数的图像为：
左边这个不动点能量更低。在V\rightarrow\infty的时候，只有左边这个点留下来。
在x_1处景观函数值为-0.0525，二阶导数值为9.9709。在x_2处景观函数值为-0.0223，二阶导数值为0.8818。在x_3处景观函数值为-0.0139，二阶导数值为-0.8102。在x_3处的正/逆向流量为3.4664。
对于不同的V，其分布的样子为：
所以说，即使是双稳态的系统，在热力学极限下也只会留下一个稳定点。然而，这个结论是误导性的：因为，此时能垒也变成无穷大，所以两个稳定点之间根本不能相互跃迁。结论就是：热力学极限下，完全按照确定性动力学跑，看初始态决定落在什么稳定点上。这个地方会出现极限交换不成立的情况：
若先t\rightarrow\infty，则达到平稳分布，即双峰，再把V\rightarrow\infty，就只留下一个稳定点；
若先V\rightarrow\infty，则在固定的时间内，系统会弛豫到某一个稳定点；再把t\rightarrow\infty，因为势垒已经是无穷大，没法跳过去了，所以还是保留在这个稳定点上面，哪个稳定点完全取决于初始条件。
这就是物理中实际存在的热力学极限中极限交换的问题。所以有人说物理里根本不用考虑极限交换的条件的时候，这就是一个反例。（用于打脸物理系看不起数学分析？）
确定性动力学的模拟轨道为：
这是个典型的双稳态。
固定V=100，其Gillespie模拟得到的随机轨道为：
从低态跳向高态的平均时间的理论计算结果为30.2702。从高态跳向低态的平均时间的理论计算结果为4.9674。这个统计懒得做了，从图上来看至少是差不多的。
蛋白质折叠的ZBS化学主方程模型
跳开来随便写点别的。不过还是化学主方程以及Q过程的首达时相关的东西。
首先是Levinthal悖论。假设一条多肽有100个氨基酸，每个氨基酸有3种状态等概率地选择，则所有的可能性有3^100种。假设每遍历一种的时间为10^{-9}s，则最长的总时间将达到10^33年（平均到达时间算起来更复杂，要用马氏链，见下面的计算）。而宇宙年龄只有10^10年，所以这种折叠方式从根本上是行不通的。
Zwanzig,Szabo与Bagchi用化学主方程提供了一种解决Levinthal悖论的方法。关键在于外界提供能量导致错误与正确之间转换的速率并不平衡，这一稍许的不平衡将导致期望时间极大缩短。设正确到错误的速率为k_0，错误到正确的速率为k_1。定义平衡常数K=k_0/k_1\leqslant 1。考虑从完全错误出发到完全正确的期望时间，也是用强马氏性列线性方程，求解得到的精确结果为
它可近似为
假设N=1000，可以画出折叠时间关于提供能量的关系
可见时间随能量的下降是极其惊人的，如果不提供能量（Levinthal悖论的情况），所需时间根本不是宇宙年龄比得上的；但是只要提供5k_BT的能量，就可以把时间缩短到1s以内。这只是一个极度简略的模型，只为了提供Levinthal悖论的一个解决思路。
关于还原论与构建论以及复杂系统的哲学
昨天在关于神经电位传播的cable equation的话题下面看到很多学生物的人坚定地认为生物系统是不可能服从Maxwell方程的，生物电就是生物电，跟物理电在本质上不同，它只能从生物角度解释，用方程理解生物都属于邪教或者民科。且不说这种人不知是受到了高中老师的什么关于“学科隔离”的神奇熏陶（毕竟高考也闹过“化学的酶不包括核酶”的笑话），或者本身神经生物学学歪了，这种想法在50年前可能还是非常普遍的（虽然已经有Hodgkin-Huxley了？）。现在不管是为了混经费，发文章还是各种各样的目的，所谓系统生物学、生物物理、计算生物学、计算系统生物学之类的各种名字都成为一门学科，关于复杂系统的哲学应该更加深入人心了（虽然我完全不认为给几个蛋白质列几个ODE数值模拟个极限环出来的水文章能对理解生命有什么贡献；所以我总是自我标榜“理论生物学”这个说法）。不管怎样，Anderson的More is different一文还是值得重新读一下的。
所有有机物和无机物的运行机制, 就我们所知而言, 都被认为受同一组基本定律所支配。对于这一组基本定律, 我们相信, 除了某些极端情形之外, 我们已经有了很好的理解。然而这并不意味着，只有那些研究真正是基础的东西的科学家才是探索这些定律的人（天体物理学家, 基本粒子物理学家, 逻辑学家和数学家等）。纵观20世纪科学的发展, 人们可以看到两种潮流：”内涵性(intensive)研究”和”外延性(extensive)研究”。简言之: 内涵性研究探求基本定律, 而外延性研究致力于按照已知的基本定律来解释现象。外延性科学并不意味着它等同于工程学，这是因为还原论并不蕴涵着构建论。大型和复杂的基本粒子集合体的行为, 并不能按照少数基本粒子性质的简单外推来理解.。事实上, 在复杂性的每一个层次, 都会有崭新的性质出现。为理解这些新行为所进行的研究, 本质上是同样基础性的（实际上统计力学做的就是这样一件事情，而且本质上是数学的。这就是为什么我觉得四大力学之中统计力学是涵义最深的一门，理论力学的Lagrangian formulation与Hamiltonian formulation，量子力学的薛定谔方程的表述的路径积分表述，电动力学由Maxwell方程或者作用量出发构建，它们都是“内涵性”的部分；而统计力学是连接“内涵性”和“外延性”的部分，这部分人们的认知实际上还非常粗糙。但是人们对统计力学更多的是一种习惯性的操作，内部深层次的东西[并不是说他们看不起的遍历理论，而是概率与大偏差]，却并没有得到重视）。心理学不是应用生物学, 生物学也不是应用化学，往上走一层都需要新的规律。
多嘴一句，上世纪八十年代曾经很火的“灾变论”，作为“系统科学”或者“复杂性科学”的一门，现在已经销声匿迹，连google一下intro都不好找了。这东西（似乎？）就是ODE分岔理论的扩展，但是被过度引申了，在社会学、政治学、历史学里面都有杂七杂八的胡说八道式的“阐述”。这属于步子太大扯着蛋，人类的数学还不足以把这样的复杂系统考虑在内，非要强行这么做，就堕入民科了。
参考文献
Vellela M, Qian H. Stochastic dynamics and non-equilibrium thermodynamics of a bistable chemical system: the Schlögl model revisited[J]. Journal of The Royal Society Interface, 2008, 6(39): 925-940.
Qian H. Nonlinear stochastic dynamics of mesoscopic homogeneous biochemical reaction systems—an analytical theory[J]. Nonlinearity, 2011, 24(6): R19.
Zhang Y, Ge H, Qian H. van't Hoff-Arrhenius Analysis of Mesoscopic and Macroscopic Dynamics of Simple Biochemical Systems: Stochastic vs. Nonlinear Bistabilities[J]. arXiv preprint arXiv:1011.2554, 2010.
