Robustness, Stability and Efficiency of Phage lambda Gene Regulatory Network: Dynamical Structure Analysis

噬菌体 λ 基因调控网络的稳健性、稳定性和效率：动态结构分析

https://arxiv.org/ftp/q-bio/papers/0403/0403016.pdf

https://arxiv.org/pdf/q-bio/0403016

摘要

基于最近发展的复杂网络动态结构理论和Shea和Ackers在1980年代的开创性工作，我们为控制噬菌体λ生命周期的基因调控网络构建了一个透明且简洁的数学框架，该框架自然地包含了随机效应。动态结构理论指出，一个复杂网络的动态由其四个基本组成部分决定：耗散（类似于降解）和随机力、由势能决定的驱动力以及横向力。势能可以被解释为噬菌体发育的景观，包括吸引盆地、鞍点、峰值和谷值。耗散使得噬菌体在由势能定义的景观中具有适应性：噬菌体总是倾向于接近附近吸引盆地的底部。随机波动赋予噬菌体通过鞍点在势能景观中搜索的能力。在我们的模型中，分子参数主要由野生型噬菌体的实验数据确定，并由一种突变体的数据补充，我们对突变体的计算结果在蛋白质数量、溶菌频率和溶菌培养中的裂解频率方面与对其他突变体的可用实验观察定量一致。计算结果再现了噬菌体λ遗传开关的观察到的鲁棒性。这是第一个成功代表如此广泛的主要实验现象的数学描述。具体来说，我们发现：1) 噬菌体λ开关的稳定性和效率的解释是鞍点穿越率对势垒高度的指数依赖性，这是在景观中的随机运动的结果；2) cI抑制基因转录的正反馈，通过CI二聚体协同结合增强，是噬菌体λ遗传开关对突变和参数波动的鲁棒性的关键。

连载标题:噬菌体λ调控网络的鲁棒性建模

关键词：噬菌体λ、调控网络、鲁棒性、随机效应、发育景观

I.导言

噬菌体λ是一种生长在细菌上的病毒。它是最简单的生物体之一。噬菌体λ的染色体由一个被蛋白质外壳包裹的DNA分子组成。在感染宿主大肠杆菌细胞时，噬菌体λ将其染色体注入细菌内部，并将蛋白质外壳留在外部。在细菌内部，它选择两种生长模式中的一种。它可以指导细胞产生新的λ噬菌体颗粒，导致细胞裂解。或者，它可以在溶菌状态下建立休眠驻留，将其基因组整合到宿主的DNA中，并作为宿主染色体的一部分进行复制。在这两种不同的生命周期中，由于分子相互作用，表达不同的噬菌体基因集。对这样一个过程的鲁棒性和稳定性进行现实建模一直是生物计算中的主要挑战。

噬菌体λ的发育过程和遗传控制被视为决定发育途径的分子机制的典范。人们相信，类似的分子相互作用可能潜在地支配着许多发育过程。通过研究噬菌体λ的基因调控网络，人们希望在分子水平上对主要生物学功能的调控获得深入的理解。这些功能之一是表观遗传状态的编程：噬菌体决定它将遵循溶菌还是裂解生长状态的方式。在过去的50年中，广泛的生物学调查在这方面提供了相当好的定性图景。有一个合理的场景可以指导对实验观察的理解。

基因调控网络执行的另一个功能是遗传开关的维持和操作。在溶菌状态下生长的噬菌体保持休眠，除非被激发。当发送信号激活RecA蛋白时，它们会切割CI单体，使噬菌体进入裂解生长。观察到噬菌体λ开关既高度稳定又高效。当噬菌体在溶菌状态下生长时，它会在许多代中保持休眠。自发诱导发生在不到百万分之一的细胞分裂中。一旦噬菌体暴露于适当的信号，它几乎以100%的速率转变为裂解状态。噬菌体λ中遗传开关的稳定性和效率仍然是一个谜。我们从生物学或物理化学方面都没有很好的理解。

对噬菌体λ进行建模的数学和数值活动一直在持续进行中。其基本原理相当直接：生物学功能应该作为基于噬菌体调控元件的分子机制和它们独立测量参数的系统属性从模型中显现出来。Shea和Ackers提出的噬菌体λ基因调控的优雅物理化学模型已成为后来研究的基础。然而，不久之后，Reinitz和Vaisnys指出理论结果与实验数据之间的不一致可能表明存在额外的协同作用。Arkin等人对噬菌体λ的发育进行了随机模拟，以决定早期的溶菌状态。最近，Aurell和Sneppen使用基于Onsager-Machlup函数的方法分析了噬菌体λ遗传开关的鲁棒性，并发现他们的理论分析无法再现噬菌体λ开关的鲁棒性。

遗传开关的稳定性和切换效率的共存是一个明显的困境，原因如下。溶菌状态异常稳定。生长环境中的波动和基因调控网络的内在波动不会偶然翻转开关。那么当噬菌体受到威胁时，切换过程如何能在如此少的外部干预下变得如此完整？这些模型内部的一致性问题自然产生：Shea和Ackers的工作中容易操作的诱导或高效率的切换是否是以牺牲基因调控网络的鲁棒性为代价的。换言之，如果一个模型被构建得能够忠实地再现观察到的基因调控网络的鲁棒性，它是否会失去开关的效率。

毫无疑问，一个可信的噬菌体λ模型应该同时再现遗传开关的鲁棒性、稳定性和效率属性。从这样的模型中，我们还应该能够计算噬菌体发育的观察量，如蛋白质数量和溶菌频率。在噬菌体λ的实验数据基础上建立这样一个数学框架是本工作的主要目标。

我们的程序首先是构建一个噬菌体λ基因调控网络的最小模型，以定量再现尽可能多的实验结果，并避免与实验数据的任何定性不一致。特别是，我们要求这个模型的计算结果与测量遗传开关鲁棒性的实验一致。在成功获得这样一个模型后，我们接着研究遗传开关的鲁棒性、稳定性和效率之间的关系。我们在模型中包括了蛋白质数量的随机波动。之所以有必要包括这种波动，不仅是因为它们无处不在，还因为它们是影响遗传开关内在稳定性和切换效率的关键量。

通过结合新开发的强大的非线性动力学分析方法，该方法考虑了随机力，并以前建立的物理化学模型为基础，我们构建了一个数学框架，以计算表观遗传状态和发展路径的数量：以野生型为参考，计算突变体的蛋白质数量、蛋白质分布、每种状态的寿命和溶菌频率。我们发现我们的计算结果与可用的实验数据定量一致。

我们解决数学模型的方法具有独特的优势。它有助于将基因调控网络结构的隐藏层带到表面。在数学分析中，我们计算基因调控网络的势能作为蛋白质数量的函数。这种势能可以被解释为噬菌体发育的景观图，包括吸引盆地、鞍点、峰值和谷值。这种发育和分化的景观概念在生物学中早已被构想出来。在这个势能景观中进化的噬菌体类似于在电磁场中移动的带电粒子。类似于电动力学中的案例，势能景观提供了基因调控网络动态结构的可视化。势能的最小值给出了可能的表观遗传状态的位置。鞍点的高度给出了这些表观遗传状态之间的时间尺度。

有了这个可用的模型，我们可以通过组织基因元素并探索它们之间的关系，以数学上容易但实验上不能或尚未完成的方式来从数学建模中对基因调控网络的生物学运作投以新的光线。通过这种操作，我们可能能够洞察基因调控网络的生物学结构。鼓励来自于这样一个事实，即我们的建模确实首次展示了遗传开关所需的鲁棒性。通过数学上消除反馈机制，然后逐一将它们放回，调整它们的相对强度，我们发现CI生产的正反馈，通过CI二聚体在操纵子上的协同结合增强，似乎是开关结构鲁棒性的关键。

我们的表述为稳定性-效率困境提供了一个直观吸引人的解决方案。遗传开关的数学等价物是一个在双井势中受随机力影响的带电粒子。在势能最小值中的寿命，代表表观遗传状态的寿命，对势垒高度有指数依赖性。通过降低势垒高度，噬菌体在溶菌状态中的寿命急剧下降。噬菌体迅速进化到更有利的生长状态，即裂解状态。正是这种指数依赖性负责开关的效率和表观遗传状态的稳定性。

我们还将讨论我们表述的实际相关性。对于任何新引入的理论量，一个关键的问题是如何实验上探测其结构。可以从物理学和化学中的压力概念中得出类比。除了其热力学意义和微观起源外，压力也是一个可以直接宏观测量的量。我们证明，本工作中引入的基因调控网络的势能和摩擦也与实验有这样透明的联系。从实验中确定动态结构的一个主要问题，即所谓的“逆向工程”，是无法唯一确定微观动态参数。我们的方法可能提供了一个解决这个问题所需的工具，通过设计更相关的实验来探索给定描述水平上可测量量之间的关系。

本文的其余部分安排如下。第二节总结了关键的生物学实验研究。第三节总结了关键的先前生化建模要素。第二节和第三节的目的是为我们的建模提供一个简洁、连贯的生物学、化学和物理背景。第四节将在噬菌体λ的最小模型内讨论动态结构分析方法。

计算结果、与生物数据的比较以及讨论将在第五节中详细呈现。在第六节中，我们将对当前工作进行评述。

II. λ Switch 

控制和维持噬菌体λ功能的遗传开关由两个调控基因cI和cro，以及λ DNA上的调控区域OR和OL组成。通过蛋白CI维持的已建立的溶菌状态通过阻断操作子OR和OL，阻止所有裂解基因的转录，包括cro。在溶菌状态下，CI的数量作为细菌状态的指示器：如果DNA受到损伤，比如紫外线照射，RecA的蛋白酶活性被激活，导致CI降解。CI数量的减少允许裂解基因的转录，首先是cro，其产物是蛋白Cro。

决策或切换集中在操作子OR上，由三个结合位点OR1、OR2和OR3组成，每个位点可以被Cro二聚体或CI二聚体占据。如图1所示，这三个结合位点控制两个启动子PRM和PR的活动，分别对应cI和cro的转录。cro的转录从PR开始，部分与OR1和OR2重叠。cI的转录从PRM开始，与OR3重叠。RNA聚合酶对两个启动子的亲和力，以及随后两种蛋白的生产，取决于Cro和CI在三个操作子位点上的结合模式，从而建立溶菌状态，每个细菌约有500个CI分子。然而，如果CI的数量变得足够小，Cro的增加生产将开关翻转到裂解状态。

对噬菌体λ开关稳定性的定量实验研究有着悠久的历史。最近，三个独立小组报告了缺失recA基因菌株中自发诱导的频率，这一结果由Aurell等人进行了综述。他们都证实了两个早期的基本观察结果：存在开关行为，且开关是稳定的。此外，尽管使用了不同的菌株背景，在不同的大陆和不同的时间进行，他们都获得了一致的开关频率数值。计算尝试去定量理解这种行为并未成功，即使允许野生型可能更稳定的假设。

更近期的数据表明，野生型可能比之前在文献[15]中观察到的值稳定两个数量级：转换到裂解状态的开关率可能小于每分钟。除了需要更多的实验研究之外，这也使得理论建模面临更大的挑战。在接下来的部分，我们将使用这个新的野生型数据作为我们模型的主要输入，以完全确定我们的模型。我们还将讨论以前的数据，以说明我们的模型中明显的指数敏感性。

III. Physical Chemical Model 物理化学模型

细胞内的CI和Cro蛋白分子被假设处于稳态平衡。在任何特定时间，并不总是有相同数量的CI和Cro二聚体结合在操作子上。这些数量是波动的，而平衡假设应该给出这些波动的大小。关键输入是CI和Cro二聚体化常数以及它们结合到三个操作子位点OR1、OR2和OR3的吉布斯自由能。

遵循Ackers等人的研究25和Aurell等人的研究16，我们用三个数字（i, j, k）来编码CI和/或Cro结合到OR的状态s，分别对应OR3、OR2和OR1。编码规则是：如果相应的位点是空的，则为0；如果位点被CI二聚体占据，则为1；如果位点被Cro二聚体占据，则为2。在Shea和Ackers的大正则方法6中，有is个CI二聚体和js个Cro二聚体结合到OR的状态s的概率是

RNA聚合酶（RNAp）占据OR1和OR2，或者OR2和OR3，而不是其他配置。我们进一步简化pR(s)的表达式，注意到OR的控制是由CI和Cro蛋白操作的，而不是RNAp2,3。如果OR1和OR2没有被CI或Cro占据，RNAp将以由RNAp结合能决定的概率与它们结合。RNAp首先结合到OR1和OR2，然后阻止CI和Cro结合的情况被排除，这是基于只有CI和Cro控制调节行为的假设。除了实验观察之外，如果CI和Cro结合的时间尺度比RNAp结合的时间尺度短，这个假设是合理的。除了一个总体常数，我们将其纳入转录速率中，RNAp结合不再相关。因此，我们将其从pR(s)的表达式中取出。

总状态数减少到27。这种简化首次由Aurell和Sneppen使用9。我们将不再使用下标R来表示结合概率pR。我们应该指出，先前的实验和理论结果已经被Aurell等人简洁地回顾了16，我们将遵循他们的惯例。

二聚体和单体浓度由二聚体的形成和解离决定，这给出了二聚体浓度与蛋白质总浓度之间的关系：

我们使用了Aurell和Sneppen9中提出的这四个参数，它们是根据溶原状态和裂解状态下得到的蛋白质数量推导出来的。

自由能∆G(s)是通过体外研究确定的。然而，体内条件可能会有所不同。所测得的蛋白质-DNA亲和力可能会敏感地依赖于缓冲溶液中存在的离子以及其他因素。这一观察结果在我们对比理论结果与实验数据时将十分重要。另一方面，体内因这些变化产生的影响应该得到补偿，例如，通过腐胺26等物质来平衡氯化钾浓度的变化，以及通过其他离子和拥挤效应27来进行调节。

我们注意到Record等人27已经观察到体内和体外分子参数之间可能存在显著差异。Darling等人引用的数据是在200mM KCl浓度下获得的，这类似于体内条件。因此，尽管我们预计体内和体外数据之间会有差异，但差异可能不会很大。

图1中描述遗传调控的数学模型是一组耦合方程，用于描述细胞内CI和Cro数量的时间变化率：

IV.随机效应与动力结构分析

IV.1随机动力学

如果CI和Cro的数量在宏观上很大，那么方程（1）将是动态的完全准确描述，因为数量的波动是N的1/2阶，而修正是1/N的1/2阶。然而，这些数量只在数百的范围内。因此，波动不可忽视。实际的蛋白质生产过程受到许多偶然事件的影响，例如CI或Cro在溶液中找到自由操作子位点所需的时间，或者RNA聚合酶分子找到并附着到可用启动子所需的时间，这表明有更多的随机来源。因此，作为一个具有有限N噪声的网络的最小模型，我们考虑以下两个耦合的随机微分方程系统，有两个独立的标准高斯和白噪声源：

方程（9）定义了一个2×2的扩散矩阵D。噪声强度可能包含来自产生和衰减率的贡献，假设每个都由一个单一独立反应主导，正如Aurell和Sneppen9所采用的。这样的噪声可以被称为“内在”噪声。还存在其他噪声源，“外在”噪声32。我们处理噪声以包含内在和外在两种来源：所有这些都假定为高斯和白噪声。这一假设的一致性应该通过实验来测试，如下文所述。然而，某些概率事件在建模的背景下可能不会表现为高斯和白噪声，这可以通过单独的生物实验来确定，如下文将讨论的pRM240突变17。

已经证明11,12，存在一种独特的分解，使得随机微分方程，方程（8），可以转化为以下形式：

从方程（8）和方程（9）到方程（10）和方程（12）的分解由以下一组方程确定：

可以从方程（13）和方程（14）中用F和D求解Λ和Ω。实际上，这是可以形式上完成的。一旦知道了Λ和Ω，要求方程（10）可以简化为方程（8）就给出了（Λ + Ω）F = -∇U(N)，这被用来获得U。一般来说，这种分解是一个复杂的数学和数值任务。摩擦矩阵的简化进一步简化了问题。通常，扩散矩阵D在生物学上是未知的：没有足够的测量来明确地确定噪声。因此，我们可以将半正定对称矩阵Λ视为需要通过实验确定的参数。在我们的计算中，我们假设D是对角矩阵。根据方程（14），Λ在二维情况下是对角矩阵。实验测量的recA-1溶菌体转换为裂解状态的比例被用来确定Λ的元素。

在这里，我们想对数学过程给出一个直观的解释。方程（8）对应于一个虚构的无质量粒子在由两种蛋白质数量NCI和NCro形成的二维空间中的动态，既有确定性也有随机力。很容易检查，在一般情况下∇×F(r) ≠ 0和∇•F(r) ≠ 0。因此，由于运动方向的横向力和摩擦力，F(r)不能简单地由标量势的梯度表示。当存在横向力和摩擦时，二维运动的最简单情况是一个带电粒子在磁场和电场中运动，这正是方程（10）的形式。

从方程（10）出发，我们可以将半正定对称矩阵Λ解释为摩擦矩阵，将反对称矩阵Ω解释为“磁场”的结果。摩擦矩阵代表物理学中的耗散。它类似于生物学中的降解。标量函数U扮演势函数的角色，将决定噬菌体的最终稳定分布。当最终分布函数由以下给出时，将达到全局平衡：

势能U，即系统的景观，在图2中描绘（参见图4）。

噬菌体在势能景观中看到两个最小值和一个鞍点。这两个最小值对应于裂解状态和溶菌状态。一旦噬菌体处于其中一个最小值，它移动到另一个最小值的概率率由Kramers速率公式给出，形式如下33,34：

在势能垒高度的情境下，该高度是鞍点与初始最小值之间的势能差，同时考虑到时间尺度以及由摩擦力、鞍点周围的曲率决定的尝试频率ω0。我们在此指出，尝试频率通常是方程(10)中动力学量的复杂函数。在本文中，其具体形式将通过经验来确定。关于一般性的数学讨论，我们建议读者参考文献[34]。

IV.2 基因调控网络的动态结构

方程（10）以四个组成部分的形式给出了基因调控网络的动态结构：摩擦、驱动力的势梯度、横向力和随机力。这种动态结构分类服务于两个主要目的。它本身提供了对基因调控网络主要特征的简洁描述，并且它提供了比较不同基因调控网络的定量措施，例如野生型噬菌体与其突变体之间。

势能可以被解释为噬菌体发育的景观图。每一个表观遗传状态由一个势能最小值表示，其周围区域形成一个吸引盆地。摩擦所代表的耗散使得噬菌体在由势能定义的景观中具有适应性：噬菌体总是倾向于接近附近吸引盆地的底部。势能在最小值附近的变化，连同摩擦，给出了松弛的时间尺度：在表观遗传状态受到扰动后达到平衡所需的时间。一旦我们知道了摩擦和势能在最小值附近的情况，我们就能很好地掌握松弛时间，τ = η/U″，这里η是摩擦的强度，U″是势能的二阶导数，都是在沿着相关轴的一维近似中。当U″是一个常数时，松弛时间与势能最小值附近扰动的幅度无关。

这里需要做两点说明。1. 摩擦矩阵的含义与力学中的含义相同：如果没有外部驱动力，系统倾向于在其附近的最小位置停止。在生物学中最接近的概念是“降解”：在给定条件下，总是存在一个自然的蛋白质状态。2. 结果表明，横向力在当前的类似开关的行为中并不是一个主导因素。然而，它的存在是振荡生物行为的必要条件，这在本文中将不再进一步讨论。

势能提供的另一个时间尺度是表观遗传状态的寿命，这是通过Kramers速率公式（方程（16））通过势垒高度给出的。这样的尺度衡量了在波动环境中表观遗传状态的稳定性。在噬菌体λ的情况下，溶菌状态的寿命非常长，除非噬菌体在其操作子位点发生突变。当噬菌体受到刺激时，分隔溶菌和裂解状态的势垒高度降低。由于其对势垒高度的指数依赖性，溶菌状态的寿命急剧缩短，发生切换。从不同的角度来看，随机力赋予噬菌体通过穿越鞍点在势能景观中搜索的能力，并推动切换事件。Kramers速率公式是这种优化能力的定量度量。

V. Results 结果

V. 1 Determining in vivo parameters. 确定体内参数。

首先，我们需要确定在理论模型中使用的自由能。迄今为止，所有测量的噬菌体λ的结合能量都是通过体外研究确定的。体内条件与体外条件之间的差异可能包括缓冲溶液中的离子浓度和基因组DNA的空间配置，例如环状结构35,36,37。表I中从体外到体内合作能量的相对较大的改变可能部分是由于环状效应，尽管当前模型中没有直接考虑环状效应。我们注意到，在体内条件下，所有操作子都处于同一种环境中，包括离子条件和DNA配置。后者的原因是操作子在基因组中彼此紧密相邻。如果基因组DNA的弯曲增加了或减少了DNA-蛋白质结合，这些紧密且短的操作子位点很可能经历相同数量的变化。因此，我们假设除了体外DNA-蛋白质结合能量之外，总体结合能量差异分别加到所有的CI和Cro蛋白上：。

为了确定∆GCI（∆GCro），我们需要的实验输入不仅仅是体外测量。为了避免不必要的不确定性，在构建的模型中，我们尽量包含最少数量的参数。包括了两个CI二聚体之间的协同结合。两个Cro二聚体之间的协同结合、CI和Cro二聚体之间的协同结合、CI和Cro的非特异性结合并未包括在内。我们后续的计算验证了CI协同结合对遗传开关特性是必不可少的，而我们忽略的结合对计算结果没有显著影响。我们需要调整的有三个参数，即CI（∆GCI）的体内与体外结合能量之差，Cro（∆GCro）的，以及CI二聚体协同作用（∆G(cooperative)）的。

我们首先使用野生型和突变体λOR121的CI数量来确定∆GCI，然后我们通过要求野生型的裂解和溶菌状态同样稳定来确定∆G(cooperative)和∆GCro，这是通过Kramers速率公式计算得出的。调整后的体内结合能量和其他我们用于建模的参数在表I中给出。使用这些调整后的参数，我们再现了噬菌体遗传开关的鲁棒性（如图3所示）。

由Little等人研究的突变体λOR3'23'[15]，Hochschild等人[38]对其进行了特征化，特别是针对其与OR3'的结合能力。为了产生所需的蛋白质水平，我们发现OR3'与Cro蛋白之间的结合能比OR3与Cro蛋白之间的结合能小1.8 kcal/mol，这与Hochschild等人的结果一致。从OR3到OR3'的CI结合能略有增加，为1 kcal/mol，这也与测量结果相符。

我们假设摩擦矩阵Λ是一个对角常数矩阵。与Aurell和Sneppen9类似，我们假设方程（2）中的随机波动与蛋白质数量的平方根除以松弛时间成比例：DCI = Const×τCI/NCI（溶菌），DCro = Const×τCro/NCro（裂解），其中NCI（溶菌）是溶菌状态下的CI数量，NCro（裂解）是裂解状态下的Cro数量。Const需要通过实验来确定。在方程（7）中，我们注意到如果反对称矩阵Ω很小，那么Λ就是D的逆矩阵。我们假设Λ = D^-1来计算Ω，并发现确实在关注的区域，即通过鞍点连接两个势能最小值的势能谷中，A是可以忽略不计的。我们最终使用的具体参数是：

V.2噬菌体λ基因调控网络的动态结构

原始问题，由方程（8）描述，可以被解释为一组二维微分方程，用以描述一个粒子的运动，如果我们将蛋白质数量NCI和NCro视为坐标，并将粒子位置视为）在时间t。有一个确定性力和一个随机力作用于这样的粒子。这个确定性力同时具有摩擦力、势能力和横向力的特征。我们之前讨论的分解，方程（10）允许我们将这些组成部分分开。我们在这里分别讨论它们。

野生型λ噬菌体及其一些突变体在势能景观中看到两个最小值和一个鞍点（图4）。这两个最小值对应于裂解和溶菌状态（参见图2）。势能最小值的位置给出了裂解和溶菌状态的平均蛋白质数量。有一个相对狭窄的山谷连接这两个最小值。沿这条山谷的最高点是鞍点。由于具有大势能的区域不易接近，而低势能区域形成了一个山谷，我们可以沿着山谷可视化势能，并在一维图中展示它，如图2和图12所示。

反对称矩阵Ω可以由沿z方向的单个标量B表示为ΩF = B z× F，假设x=NCI，y=NCro。野生型噬菌体的横向场B通过数值求解方程（13）获得，结果绘制在图5中。这个场在除了沿两轴的区域外都是小的。沿两轴，横向B场没有效果，因为运动是由陡峭的势能引导到山谷中的。一旦噬菌体从原点演化开来，当CI和Cro的数量都很小的时候，在后期发展中，横向力可以从方程（10）中取出，而不会改变噬菌体的动态。在计算松弛时间和溶菌状态的寿命时，由于上述原因，我们可以忽略横向力。

V. 3 溶菌和裂解状态的蛋白质数量

势能最小值的位置给出了裂解和溶菌状态的平均蛋白质数量。在溶菌（裂解）状态中，实际上只存在CI（Cro）蛋白。通过改变模型中的参数，我们还可以了解当实验条件变化时蛋白质数量如何变化（如表II所示）。当温度升高时，裂解和溶菌状态中的蛋白质数量都会增加。如果CI蛋白的降解时间减少，溶菌状态中CI蛋白的数量就会减少。我们将在后面展示，这种减少会降低溶菌状态的稳定性，使噬菌体倾向于裂解生长。

V.4 溶菌状态的蛋白质分布

势能最小值附近给出了平均值周围的蛋白质分布情况。这样的分布可以通过实验测量，尽管对于噬菌体λ尚未进行此类测量。最有可能的测量是溶菌培养中CI蛋白的分布。我们在图6中给出了分布的形状，这是使用方程（15）计算得出的。突变体的势能限制较浅，导致蛋白质数量的分布更广泛。

V.5 溶菌状态的松弛时间

当蛋白质数量在受到扰动后偏离其在势能最小值处的值时，它们倾向于松弛到势能最小值处的值（如图7所示）。松弛时间由势能最小值附近的势能和摩擦强度决定：τ = η/U″，其中η是摩擦强度，U″是势能的二阶导数，η和U″都是沿着扰动方向的。我们在表III中给出了根据Little等人的数据测量的野生型和突变体的松弛时间。

V.6 溶菌和裂解状态的稳定性

首先我们讨论在recA被禁用时溶菌状态的内在稳定性。当噬菌体在溶菌培养中生长时，偶尔的波动可能会使噬菌体λ转向裂解生长。在我们解出U（势能）之后，计算Kramers逃逸率所需的势能差是已知的。我们仍然需要尝试频率。尝试频率的计算通常较为复杂，并且它对逃逸率的影响不如能量差那么显著。因此我们简单地将其视为一个拟合参数，得出 /分钟。Little等人的实验中突变体的势能差如下，表IV所示。为了完整性，我们还包括了裂解状态的寿命。这是一个噬菌体在裂解状态下生长可能转变为溶菌状态的时间尺度。由于野生型噬菌体在一个小时左右就会裂解细胞，这个时间尺度在野生型噬菌体中是不可观察的。只有在裂解被抑制的情况下，才能观察到裂解状态的内在寿命[43], [44]。

对于recA+噬菌体，CI单体可能被RecA蛋白切割。当噬菌体没有被激发且SOS系统没有被激活时，可以合理假设这种CI单体的切割只是偶尔发生。这种切割对噬菌体的主要影响是增加了CI蛋白的噪声水平。我们发现，如果噪声水平翻倍，势垒高度就会减半。对于这样的噪声水平，已经转换到裂解状态的溶菌体的比例与实验结果吻合得很好（如表II所示）。

V.7 遗传开关的鲁棒性

遗传开关的鲁棒性是指其对某些参数变化不敏感，并且能够容忍某些遗传多态性。λ开关对突变的稳定性由Little等人[15]测量，并与我们的计算结果一起总结在表II中。Aurell等人[16]的计算未能复制出这种鲁棒性。他们发现，同时复制野生型recA-溶菌体的测量稳定性以及突变体323中稳定溶菌体的存在是不可能的。

我们基于新的参数值集对开关进行了直接测试。由于存在两个势能最小值和一个鞍点表明了开关的稳定性，我们寻找配置的范围。我们发现，对于平均CI单体数量从相对于野生型值，通过将方程（8）中的NCI替换为αNCI，以及Cro单体数量从相对于野生型值，开关是稳定的。我们也让CI降解时间τCI从0.1变化到15相对于其野生型值，以及τCro从相对于其野生型值，再次，开关是稳定的。请注意，其中一个参数已经变化了9个数量级。通过以上所有对稳定性的测试，我们得出结论，这个噬菌体λ开关非常稳定，对参数的小变化不敏感（图8）。

通过系统地改变结合能和其他参数，我们试图探究负责开关鲁棒性的关键要素。我们发现CI二聚体的协同结合是最重要的因素。如果我们将CI二聚体的协同结合减少一半，野生型不会受到显著影响。然而，被观察到稳定的突变体λOR323，在计算中变得不稳定（如图9所示）。我们可以将协同结合的效果理解为增强正反馈的一种方法。有了这种协同性，当CI二聚体结合到OR1或OR2时，正反馈就会被激活。

V.8 开关效率

对噬菌体λ遗传开关的鲁棒性分析表明，其表观遗传状态对参数变化，包括蛋白质数量变化，具有鲁棒性。那么开关是如何实现转换的呢？从理论角度来看，噬菌体从溶菌生长转换到裂解生长有两个途径。实际上，噬菌体似乎同时使用了这两种策略。为了清晰起见，我们首先分别讨论这两个途径。

第一个诱导途径是在保持所有其他条件不变的情况下增加CI蛋白数量的噪声水平。从数学上讲，这意味着增加方程（8）中的ζCI和方程（9）中的DCI，同时保持方程（8）和方程（9）中的所有其他项不变。通过分解过程改变摩擦矩阵Λ。结果，势能U也会发生变化。因此，对于不同的噪声水平，噬菌体在不同的势能景观中移动。这种噪声水平的变化有剧烈的效应。它通过使溶菌势能井变浅来改变其最小值。作为一个好的近似，溶菌势能井的势垒高度与噪声强度成反比。噪声水平翻倍，势垒高度减半。结果，增加的噪声水平显著降低了溶菌状态的寿命，如图10所示。另一方面，裂解状态的寿命保持不变。势能景观中的这两种变化的结合有效地使噬菌体转向裂解生长。

第二个诱导通道是通过方程（8）中的确定性项。对于确定性项来说，例如，引入CI单体的裂解相当于将方程（8）中的NCI替换为αNCI，这里的α是一个代表裂解强度的因子（图11）。如果α小于0.02，我们发现溶菌状态不再稳定，即不再是势能最小值。这样的小α值的解释是几乎每个CI单体都被裂解了。如果α较小，比如说0.1，意味着90%的CI单体被裂解，溶菌状态仍然稳定，其寿命几乎不变。显然，仅通过均匀裂解CI单体而不引入额外的CI水平噪声，并不是一个有效的诱导方式。

噬菌体可能同时使用了这两个通道。第二个通道显然被使用了，因为RecA蛋白会裂解CI单体。第一个通道也被使用的迹象强烈，这来自于观察结果：在没有激发的情况下，recA+噬菌体显示出比recA-噬菌体短得多的溶菌状态寿命。在没有激活RecA蛋白的可观察规模下，溶菌寿命的这种显著减少可以通过CI噪声水平翻倍来解释。图12给出了切换过程的示意图。

在Shea和Ackers6的早期工作中，没有包括随机效应。在他们的模型中，即使溶菌势能最小值很浅，也会限制噬菌体继续在溶菌状态下生长。只有在溶菌势能最小值完全消失时，才会发生切换。对于他们使用的参数，他们发现当20%的CI单体被裂解时，这样的切换就会发生。正如Aurell等人[16]所指出的，在这些早期工作中，遗传开关建模结果没有显示出观察到的鲁棒性。在我们要求遗传开关应该展现出观察到的鲁棒性之后，溶菌势能最小值的消失被推到了2%。然而，实际的切换发生在溶菌势能最小值消失之前，当溶菌势能最小值太浅而无法限制波动时。如果10%的CI单体被裂解，我们至少期望由于CI单体数量的减少，DCI增加10倍。溶菌状态的势垒降低到小于1，因此变得太浅，无法允许持续的溶菌生长。

V.9 与Little等人实验的定量比较

我们总结了与Little等人[15]的测量相关的计算结果在表V中，因为他们的数据是最更新和系统的。在他们的实验中，他们测量了recA+和recA-噬菌体每个溶菌细胞的自由噬菌体数量，但没有将recA+转换为转换到裂解状态的溶菌体的比例。如果我们假设recA+和recA-噬菌体的爆发大小相似，我们对RecA+蛋白的计算与他们的测量结果在定量上是一致的。

在表V中，我们假设了噬菌体λ的基因开关具有双稳态性，并计算了裂解状态下的蛋白质水平。这当然不适用于野生型，因此提出了一个测试计算出的Cro水平的实验问题。实现双稳态的一种方法可能是通过抑制裂解作用，达到所谓的抗免疫表型[43][44]。

正如在当前模型的公式化中所讨论的，我们做出了简化的假设，将所有偶然或概率事件视为高斯白噪声。这个假设影响了两个可测试的生物学量：溶菌状态的寿命[方程（16）]，以及溶菌状态下CI数量分布的形状[方程（15）和图6]。同时测量这两者可以用作对高斯白噪声假设的一致性检验。我们将recA+对切换动态的影响视为高斯白噪声，以简化我们的计算，这与本文中最小建模类型方法的相同推理。实际上，我们已经假设有了recA+，总噪声强度翻倍，与图10和11中讨论的相符。使用这个假设，我们计算了裂解状态的切换率，由表V的最后一列表示。带有recA+的CI分布应该比带有recA-的情况宽两倍。这两个结果都受到实验测试的制约。

可能有一些偶然事件不能在当前公式中被视为高斯白噪声。一个例子已经在生物实验中被提出[17]，pRM240突变大大削弱了启动子，因此削弱了产生CI的能力。这种突变使得溶菌体几乎不稳定，并估计至少对野生型中观察到的裂解切换负责99%。我们已经将这个输入用于表I和V。然而，我们可以重新计算所有菌株裂解状态的切换率，假设相同的最小模型，具有相同的切换率函数形式，但使用先前的实验数据[15]。以这种方式获得的切换率是：野生型(λ+)，；。实际上，野生型的稳定性降低了超过2个数量级。对于野生型，总体噪声强度增加了60%，导致溶菌状态下CI分布更宽。其他量，如蛋白质水平，没有明显变化。分子参数中唯一值得注意的总体变化是体内协同能量，从-6.4 kcal/mol变为-6.7 kcal/mol。模型与实验之间存在良好的总体定量一致性。

确实，在自然科学中的任何数学建模都需要经验输入来完全确定其数学结构。对于噬菌体λ的建模，已经有大量的分子数据使我们能够几乎确定我们的模型。我们的参数中额外的自由度是通过野生型的数据来固定的，比如切换频率。上述对我们数学结构的不太敏感的频率，即分子参数的几个百分点的变化可以导致频率两个数量级的变化，这是我们建模内部一致性的显著证明。它证明了切换对某些分子参数具有指数级的敏感性。除了需要更多的理论努力来超越我们目前的最小建模之外，显然需要更多的实验来测试当前模型：我们模型中的精确体内分子参数以及蛋白质数量的分布和时间相关性应该被视为预测。

V.11 从直接实验确定动态结构

我们为基因调控网络引入了四个动态量：摩擦、势能、横向力和随机力。摩擦和随机力的强度是相关的。对于遗传开关来说，横向力与动态属性无关。因此，遗传开关的两个关键量是摩擦和势能。这些量可以从更微观的模型中，用分子参数来计算。目前的定量成功在于允许体内和体外的差异，以及各种噪声贡献的存在。

然而，这四个量可以直接通过生物学方法确定。

有三个实验可以确定遗传开关的摩擦和势能（图13）。第一个实验是每个表观遗传状态周围的蛋白质分布，给出的公式为：

第三个实验是表观遗传状态的寿命。噬菌体从一个生长的表观遗传状态演变到另一个状态的概率由Kramers速率公式给出，即我们的方程（16）。

其中 ∆Ub 是势垒高度，ω0 是尝试频率。ω0 由摩擦和势能势垒的曲率决定。势能的曲率与势能势垒的高度以及势能最小值附近的形状有关。因此，可以从其表观遗传状态的寿命中确定∆Ub：∆Ub = ln(ω0) – ln(P)。

V.11 动态结构分析作为一种现象学描述

噬菌体λ的基因调控网络是一个复杂的动态系统。收集这个数学建模所需的参数需要数十年的巧妙实验研究和辛勤工作。对于一个更复杂的系统，资源和时间可能限制了我们以如此细致的方式研究每个分子元素的能力。一种对细节要求较少但仍然能够捕捉主要特征的方法是极具兴趣的。动态结构理论分析提供了构建这样一个现象学模型的指导。

动态结构理论中引入的数量，摩擦、势能梯度、横向力以及与摩擦相关的随机力，都是在给定描述水平上可测量的数量。这些数量类似于热力学中的温度、压力、自由能，它们由微观细节决定，但可以独立于细节进行测量。一旦建立了这些数量之间的关系，如方程（10）所示，我们就准备好了写下网络的有效运动方程，而无需依赖于细节。

VI. Perspective远景

十年前，在发现细菌噬菌体λ戏剧性切换行为的四十年后，人们承认[2]我们现在应该“感谢那些早期工作者的洞察力，他们认识到这种病毒的生长——特别是它以两种不同模式生长的能力——是基因表达的一个启示性例子”。这些研究的重要性在于其背景，即控制基因表达的分子机制也适用于许多其他生物调节过程。

可以将噬菌体λ在分子生物学中的研究与氢原子在现代物理学中的研究进行比较。在这两种情况下，研究对象都是同类中最简单的：噬菌体λ是最简单的生物体之一，氢原子是最简单的原子。在理解它们的过程中，出现了新的理论框架和实验方法。与氢原子研究相关的量子力学的发展，深刻地改变了科学和技术的世界。尽管我们无法预见噬菌体λ的研究最终将带给我们什么，但其重要性是毫无疑问的。

我们认为，当前的研究是一个展示定量方法力量的例子。它表明，数学框架既可以作为定量描述工具，也可以作为生物发现工具。通过数学建模讨论了被研究的生物学属性，噬菌体遗传开关的鲁棒性、开关的稳定性和效率，以及允许它们共存的机制。通过改变模型中的参数并研究结果，我们发现协同作用在遗传开关的鲁棒性中的重要性，以及波动在开关效率中的作用。生物体内的基因调控系统是一个复杂的网络。其动态方程涉及至少数十个分子参数。这些参数可以通过实验以合理的精度确定。正如我们在这项工作中所展示的，即使是像噬菌体λ这样简单的生物体，其数十个分子参数确实已经被测量，我们仍然可以允许体外测量的参数和体内参数之间的差异。因此，对于比噬菌体λ更复杂的生物体来说，进行全面的建模，包含所有测量的分子参数，似乎是一项艰巨的任务。对于这些复杂的系统，我们在前一节讨论的最后部分，从直接实验中建立现象学模型，可能提供实际帮助。然后，显然重要的是，尽可能多地对像噬菌体λ这样的简单系统进行全面建模，以获得洞察力，以便为更复杂的系统建立现象学模型。

尽管当前取得了定量成功，但基因调控网络的研究才刚刚开始。我们还没有触及深刻的生物学问题，例如为什么噬菌体选择这样的结构，或者是什么进化原则指导了这种选择[45]。最重要的是，即使在我们已经获得了噬菌体λ基因调控网络的彻底定量描述之后，我们仍然需要找出如何改进我们已经发展的知识和技术，以便将其应用于更复杂的高等生物。一旦这些问题得到回答，噬菌体λ将引导我们理解比我们现在能想象的更深刻、更宏大的知识，这并非不可想象。

https://arxiv.org/pdf/q-bio/0403016