反向工程基因网络：将遗传扰动与动力学建模相结合

摘要

基因组项目正在迅速生成控制细胞行为的基因和蛋白质的广泛列表，对这些列表的分析正在提供大量临床相关信息。同时，在描述许多细胞系统的调节机制方面取得了令人瞩目的进展(1). 转录调控，细胞用于控制基因表达(2,三)当调节蛋白通过增强或阻断启动子区聚合酶结合的生化反应增加或降低转录速率时，就会发生。由于许多基因编码调控蛋白，可以激活或抑制其他基因，因此出现了一个复杂的相互作用的基因和蛋白质网络或电路。阐明遗传水平的亚细胞过程如何在表型水平的宏观现象中表现出来，将是后基因组研究的主要目标。

许多细胞过程在遗传水平上通过类似复杂电路的图表进行描述(4)最近，人们对与这种基因组电路有关的两大研究途径产生了兴趣。光谱的一端是量化基因调控的基本规律。在电路类比的背景下，这个问题涉及到一组介观方程的推导，这些方程忠实地量化了遗传电路中包含的信息。自然的攻击计划是使用前沿工程方法，根据基础生物化学生成的一组方程设计和测试相对简单的电路。这一领域的最新研究成功地将动力系统分析与相对简单的遗传回路的构建相结合，例如自动调节单基因网络（参考文献。5; F.Isaacs、J.H.、C.R.Cantor和J.J.C.，未发表的工作），遗传开关(6)和遗传振荡器(7)。

另一个极端是推断自然发生的大规模网络中基因的连接性的项目。这项工作是由最近的技术进步推动的，这些技术进步允许同时测量数千个基因的表达水平。这种微阵列技术可以快速生成大量的基因活动模式目录，这突出了对系统工具的需求，以识别潜在基因网络的结构和动力学。这里，系统识别问题(8)自然属于逆向工程; 一个复杂的遗传网络是大量表达数据的基础，其任务是推断遗传电路的连通性。

逆向工程方法需要大量数据集和大量计算资源。通常有大量的网络体系结构与给定的表达式数据集兼容，而这样的映射问题最初使推导特定网络的任务看起来很艰巨。因此，有几项研究使用遗传算法、非线性模型、时间序列分析和贝叶斯模型以小型网络为目标(9–15)但尚不清楚这些技术是否适用于大型网络（>100个基因）。从全基因组网络分析大数据集的技术包括聚类分析和系统搜索与某些病理状态相关的基因表达特征模式(16–20)通常只提供有关网络结构的间接信息。

如上所述，在数学建模的背景下，已经构建并研究了几种新颖的小规模设计师基因网络（参考文献。5–7和21; F.Isaacs、J.H.、C.R.Cantor和J.J.C.，未出版作品）。在目前的工作中，我们探索了这种设计者基因网络在大规模网络逆向工程中的应用。这些小型设计器网络可以插入细胞，并用于为基因表达实验提供受控的扰动机制。我们的理论基础是，mRNA水平的变化提供了有关网络拓扑的间接信息。我们的逆向工程方案旨在通过分析一系列合理构建、设计扰动的微阵列实验生成的数据，为实验人员提供一个稳健的方法来推断网络拓扑。

方法

在这里，我们讨论了如何构建一个捕获基因网络结构的动态模型，以及如何设计一个对噪声鲁棒的逆向工程方案，只使用基因表达的稳态变化，同时使用现实的先验的基因网络性质的统计约束。由于蛋白质和代谢物浓度的可靠大规模测量尚不可行，我们将重点放在mRNA动力学上。

虽然我们的动机来自于微阵列实验中单个mRNA浓度的测量，但在这里我们将通过利用以数字表示生成模拟微阵列数据的模型。我们将这些模型指定为数据生成模型，我们将通过扰动这些模型，然后使用生成的数据推断模型的潜在连通性来证明我们的过程的有效性。

我们可以将网络的动力学投影到一个通用的线性映射模型上，因为我们可以围绕稳态传递扰动。我们认为N个基因，具有典型的时间尺度τ₁, τ₂, . . . , τ_N个我们将基因的mRNA表达水平表示为x个₁,x个₂, . . . ,x个_N个。在没有任何交互的情况下，我们让我th mRNA物种以一定速率降解γ_我然而，mRNAjth可能间接影响另一mRNA的动态我th通过蛋白质和代谢物等中间产物，从而改变其转录速率。我们用有效的基因-基因耦合系数来表示这一点w个_ij公司.我们通过使用斜坡函数来扰动基因对_我(囊性纤维变性。图。​图11 一和b). 线性化的映射模型x个₁=一₁, . . . ,x个_N个=一_N个，是

1

对于我= 1, . . . ,N个，带2N个+N个²未知参数(N个γ，N个τ’s，以及N个² w个的）。

在单独的窗口中打开

图1

(一)三基因网络示意图。箭头和实心圆圈分别表示激活和抑制的大小w个_ij公司来自基因x个_j到x个_我.基因切换扰动x个₁通过使用斜坡功能。(b)在模拟基因表达实验中，拨动（虚线）对基因动力学（实线）的影响。之前基因活动的平均值(一_我在时间0至10 a.u.（粗线）和之后（时间60至80 a.u.，粗线）估计扰动。

图中的三基因网络说明了我们在本研究中进行的实验和分析的一般性质。​图1。1使用遗传拨动开关（图。​（图11一; 裁判。6)，我们可以选择性地干扰给定基因的活性（此处选择x个₁). 监测持续扰动之前、期间和之后基因表达的时间进程（图。​（图11b). 一个瞬态刺激将开关从较低状态切换到较高状态，然后开关保持在较高状态（图中的虚线）。​图11b). 其他两个基因的活性(x个₂,x个_三)因为他们与x个₁因此，测量扰动前后的基因表达水平可以为我们提供有关网络结构的信息。

因为逆向工程方法必须对微阵列中固有噪声或实验变异性引起的瞬态波动具有鲁棒性，所以我们只测量基因表达的平均稳态值（图中的粗线。​图11b). 因此，等式。1变为0=w个′_我1(x个₁−一₁)++w个′_我,我−1(x个_我−1−一_我−1) − (x个_我−一_我) +w个′_我,我+1(x个_我+1−一_我+1) + . . . +w个′_国际号码(x个_N个−一_N个)，我们在那里吸收γ_我进入之内w个_ii（ii）并重新缩放耦合参数，即，w个′_ij公司=w个_ij公司/(w个_ii（ii）− γ_我)，只留下N个²参数。逆向工程问题是推断所有未知的参数w个′_ij公司，构成矩阵W公司，来自基因表达的诱导变化Δ_我=x个_我−一_我如果我们必须搜索所有可能性，那么大型密集网络的推理问题在计算上是很难解决的。这是因为与数据一致的可能解决方案数量之大令人望而却步。来自细胞网络的数据，包括蛋白质-蛋白质相互作用(22)和代谢网络(23)，建议稀疏拓扑，因为输入的最大数量(k个_最大值)到一个单位是k个_最大值≪N个此约束减少了反向工程算法中的搜索空间和计算次数。

以数字表示实验和结果

在本节中，我们用两个例子来说明我们的逆向工程算法的有效性以数字表示实验。我们证明，通过有选择地迭代扰动稳态周围的基因网络，我们可以确定数据生成模型的底层架构。在第一个以数字表示实验中，我们考虑了一个与映射模型等价的线性模型。因为在这种情况下映射是精确的，所以我们能够得出与映射引起的误差无关的结论。在第二次以数字表示实验中，我们考虑了一个先前报道的非线性数据生成模型，该模型描述了果蝇属分段网络(24). 在这种情况下，我们能够证明，使用线性映射可以正确推导出基础非线性模型的连通性。

例1：线性基因网络。

为了说明我们的方法，我们构造了一个随机网络W公司具有n个=40和k个_最大值= 3. 我们任意选择10个基因具有三个输入连接，20个基因具有两个输入连接和10个基因有一个输入连接。然后我们随机分配这些连接。在这个假设的40个基因网络中k个_最大值=3，有N个(N个!/(k个_最大值!(N个−k个_最大值)!))2^{k个_最大值}≈ 10⁵给定基因的可能输入集。可能的矩阵数W公司因此数量级为10²⁰⁰我们的逆向工程算法的目标是将连通矩阵的数量减少到一个，下面我们将展示如何使用合理选择的扰动序列中的表达式数据来推断唯一的W公司.

作为一个例子，我们在这里着重于识别特定基因的输入(17)我们任意选择干扰基因1。然后我们检查了基因表达的变化，Δ_我（图。​（图22 顶部)，以及一组可能的解决方案（Var[(M（M）₁₇(:,j))]; 图。​图22 中间). 诱导的基因表达变化从小到大不等（图。​（图22 顶部). 我们跟踪了以下可能的输入一致的用表达式数据生成第一个扰动。直觉上，我们预计特定基因的表达水平会发生微小变化(j)对基因是否j影响另一个基因(我). 给定基因的拟议输入变量之间的预期关系j到基因我数值实验证实了诱导表达式变化的幅度（图。​（图22 底部). 该观察结果是选择在一致输入中具有最小表达变化和最大变化的基因进行后续扰动的基础。确定W公司，我们对每个步骤中的所有基因执行此计算。

在单独的窗口中打开

图2

(顶部)基因表达的变化（Δ_我=x个_我−一_我)对于所有基因(x个轴，n个=40）当x个₁感到不安。(中间)对于每一个基因，都有来自其他基因的几个不同的可能输入（解决方案），这些输入与一个微扰实验生成的表达数据一致。这里我们绘制了标准偏差σ_{我,j,e（电子）}属于M（M）₁₇(:,j)（钻石；年轴）不同基因(x个轴，j)和平均重量〈w个〉 [〈(M（M）₁₇(:,j))〉; 开环]跨越所有可能的基因输入解决方案x个₁₇在第一次微扰实验之后(e（电子）= 1). 两个实际的输入连接（来自基因14和18）用+表示。(底部)表达式水平的变化与建议解决方案中的变化之间的相关性。

我们对几个随机产生的40个基因网络进行了逆向工程(k个_最大值= 3). 在图中。​图3，三，所有基因的一致可能解的平均数k个_最大值＝3相对于扰动实验的数量绘制。为了进行比较，我们研究了W公司当所有单基因扰动被随机选择而没有重复时（图中的圆圈）。​图3）。三). 这里，需要10个扰动实验来确定网络中所有基因的连通性。按照我们的方案，更明智地选择基因更有效。它只需要七个扰动（图中的三角形）就可以得到正确的网络结构。​图3）。三). 因为W公司从定义上讲是稀疏的，单个基因的扰动通常导致整个网络中的表达活性几乎没有变化。因此，通过在每个实验中为不同的基因引入多个扰动可以获得效率（图中的平方。​图3）。三)。

在单独的窗口中打开

图3

一致可能解的数量（对数标度）作为扰动实验数量的函数，采用三种不同的选择程序：圆圈表示每个实验扰动一个随机选择的（RP）基因；三角形对应于每个实验中扰动的一个算法选择（AP）基因；和方块代表了每个实验中扰动的四个算法选择的基因。三角形和正方形都表明了N个_CPS公司当我们通过算法选择基因时，（一致可能解的数量）与实验数量（虚线）成对数下降(n个= 40). 注意，因为第一个扰动是随机选择的，所以这两个选择方案之间的差异被揭示出来了之后第一个扰动。

推断网络结构所需的扰动数量不仅取决于由N个和k个_最大值，以及可能的错误来源，包括(我)实验分辨率(ii（ii）)映射错误，例如非线性引起的错误，以及(三)有限模型分辨率（即w个_ij公司). 在这里，我们检查了实验的临界数量E类_C类取决于误差ɛ和Δ之间的比率，其中Δ是基因表达的绝对变化|x个_我−一_我|所有基因和所有实验的平均值（图。​（图44一). 因为两者都是N个和k个_最大值确定网络的组合复杂性，我们考虑比率之间E类_C类和日志(N个_CPS公司)，其中N个_CPS公司是一致的可能解决方案的数量。

在单独的窗口中打开

图4

缩放行为所需数量的扰动实验，以确定基因调控网络(一)增加误差会增加关键实验的数量(E类_C类). 通过使用相同的模拟表达式数据集，实验分辨率降低（增加ɛ）；因此，信号（Δ）不变。估计的斜率，所需的实验次数E类_C类按日志缩放(N个_CPS公司)，来自图。​图3三绘制（装箱）以供参考（RP，随机选择的扰动；AP，算法选择的扰动）。为了便于参考，绘制了一条具有单位斜率的虚线(n个= 40). (b)随着网络规模的增加，扩展行为(N个). 使用多个开关（四到八个），ɛ/Δ=0.5。填充框表示n个=40例，每个实验使用四个算法选择的扰动。(c（c）)实验临界数的理论估计，E类_C类，作为的函数N个，最大输入数k个_最大值，以及ɛ/Δ项的两个不同值。(d日)小幅增加E类_C类随着分辨率的增加w个（即。，秒)对于ɛ/Δ项的两个不同值(n个= 2,000).

为了参考，我们绘制了图。​图3三在图中。​图44一（装箱）。显然，与每个实验只扰动一个基因（圆圈）相比，每个实验扰动多个基因（正方形）可以更有效地减少可能的解决方案数量。此外，增加同一组表达式数据的误差ɛ（即不同情况下Δ不变）需要更多的实验来解决所有情况下的网络问题（图。​（图44一). 选择算法产生½的斜率，而E类_C类/日志(N个_CPS公司)使用随机选择方案时（三角形）。因此，我们期望我们的基因选择算法在ɛ和Δ之间的比率较大时特别有用。

我们还研究了E类_C类和日志(N个_CPS公司)取决于基因的数量N个（图。​（图44b). 使用多重扰动，我们发现实验的临界数，E类_C类，很好地近似于对数(N个_CPS公司)对于不同的N个这一发现使我们能够表达E类_C类就网络参数而言

2

哪里秒是的不同可能值的数量w个_ij公司和参数米和第页是通过拟合找到的常数。在这里米≈0.75和第页≈0.5，对于每个实验四个算法选择的扰动（图。​（图44一)，而第页≈1和米对于随机选择的扰动，≈0.75。在N个»k个_最大值限制，E类_C类按日志缩放N个（图。​（图44c（c）). 图。​图44c（c）还显示了如何E类_C类取决于两者k个_最大值和ɛ/Δ。注意，通过展开日志Nk公司_最大值对于主要订单，我们可以使用½k个_最大值[1+ɛ/Δ]对数[N个]估计E类_C类.

对E类_C类在秒，不同值的数量w个_ij公司，相对较弱（图。​（图44d日)，因为秒在对数范围内（等式。三). 因此，在w个_ij公司，网格尺寸Δw个= 2/(秒+1），不会显著增加实验的临界数量。反向工程类似于酵母中蛋白质网络的网络(22)，图。​图44d日说明我们的方法通过25–100个实验（ɛ/Δ=1−10，k个_最大值= 5,秒= 10,n个=2000），而50-300次实验(k个_最大值=15）需要恢复所有连接。通常，我们有E类_C类≪N个尽管解决密度更大的网络需要更多的实验。

例2：非线性基因-蛋白质网络。

尽管前面例子中的线性数据生成模型为我们的方案提供了一个系统基准，但自然发生的基因调控网络可能包含显著的非线性(21). 在这个以数字表示实验中，我们探索了我们的方案在先前报道的描述分段极性网络的数据生成模型中的应用果蝇属(24). 尽管该模型包含强烈的非线性，并且有一些我们认为无法直接测量的蛋白质-蛋白质相互作用，但我们发现可以使用线性模型恢复有效的基因-基因相互作用。

这个果蝇属模型由描述基因和蛋白质时间演化的10个非线性方程控制（参见支持文本作为支持信息发布在PNAS网站上，网址：www.pnas.org). 方程的形式由基因的信使核糖核酸的进化给出x个,dx公司/日期=年₁对^v（v）₁/κ₁^v（v）₁+年₁对^v（v）₁，其中对= 1 −年₂^v（v）₂/（κ₂^v（v）₂+年₂^v（v）₂). 这个年₁蛋白质激活基因x个和年₂蛋白质起阻遏作用。半最大活化系数由κ决定₁和κ₂分别为和v（v）₁和v（v）₂是希尔系数。完整模型（图。​（图55一)有几种蛋白质-蛋白质相互作用，例如PTC、CI和CN之间的相互作用。然而，当我们对分割网络进行反向工程时，我们只监测和考虑mRNA表达的变化，试图恢复有效的基因-基因相互作用（图。​（图55b)，在缺乏关于蛋白质动力学的信息的情况下。

在单独的窗口中打开

图5

(一)的接线图果蝇属分割细胞模型。基因显示为圆形，蛋白质显示为方形。(b)有效的基因到基因接线图。此处的链接表明至少存在一条连接两个基因的蛋白质通路。箭头和实心圆圈分别表示激活和抑制。

如线性部分所述，扰动网络足以识别主要连接。我们的逆向工程算法发现ptc公司基因与自身之间存在显著差异，表明存在强有效的负耦合。该算法还表明ci公司到ptc公司（图中的虚线。​图6），6)，而所有其他连接ptc公司建议为零。这与原始的分割网络很好地比较。从图中可以看出。​图5，5，有两条途径ptc公司基因抑制自己的表达。第一条途径涉及ptc公司基因激活蛋白PTC，PTC刺激蛋白CN的产生，CN随后抑制ptc公司在其他净负性途径中，PTC蛋白也抑制CID蛋白的产生，CID蛋白增强ptc公司表达式。负面消息ptc公司因此，我们的算法检测到的自我作用对应于这两条通路的联合抑制效应。对果蝇属模型。通过我们的反向工程方案（图。​（图6），6)是基因-蛋白质途径网络中有效基因-基因相互作用的准确重建（图。​（图55b). 所有具有统计意义的连接（图。​图6）6)我们的算法检测到的是正确的，因为它们与图中显示的有效基因到基因连接具有相同的符号。​图55b.

在单独的窗口中打开

图6

恢复的接线图果蝇属分割细胞模型。具有统计意义的连接显示为实线；虚线表示不具有统计意义的建议连接。星形图表明，我们的算法给出的连接符号与从相应固定点的分割细胞模型方程中数值计算雅可比矩阵时的连接符号相同。箭头和实心圆圈分别表示激活和抑制。

作为对我们方案的独立测试，我们对果蝇属分割模型。由此发现的基因与基因之间的联系如图中的星号所示。​图6。6所有这些连接都通过反向工程程序进行识别（图。​（图6）。6). 有趣的是，原始路径图中的符号（图。​（图55 一和b)不是中间反应级联净强度的可靠预测值。事实上，我们的方案表明重量（wg）自我互动和ptc公司-至-小时投影功能较弱。到的三个输入中小时，来自ptc公司基因具有最小的雅可比项，这与我们从逆向工程方案中发现的一致。我们的算法还发现了完整基因-蛋白质网络中并行路径的净效应ci公司-至-重量（wg）相互作用被确定为净正弱连接。将两条蛋白质途径中的雅可比项相加，可以确认总和为正但很小（图中的星号）。​图6）。6). 这个英语自耦合是没有相应雅可比项的单个连接。最后，我们注意到，如果我们可以扰动和监测任何基因或蛋白质的活动，那么我们可以重建与图对应的完整网络。​图55一（未显示数据）。这个问题就相当于重建一个没有蛋白质的基因网络。

讨论

我们开发了一种适用于重建基因调控网络的迭代逆向工程方法。通过使用最小线性模型，并通过选择性地迭代扰动基因网络中的基因，我们可以通过少量实验恢复网络拓扑。我们已经通过一系列在数字中实验，包括非线性测试果蝇属基因-蛋白质网络。

我们的算法以行为基础对网络进行逆向工程，因此其效率取决于k个_最大值，输入边数的界限，但它与输出边数无关。我们的算法需要O（运行）(k个_最大值（ɛ/Δ）对数[N个])确定每个基因最多包含一个网络的扰动实验k个_最大值输入。这个k个_最大值日志[N个]因子与之前基于布尔模型缩放行为的推测一致(16). 我们还发现，当误差较小且基因表达的变化较大时，我们直觉上需要较少的微扰实验。我们观察到，使用多个开关进行微扰比使用单基因微扰更有效。信号（Δ）随着多次扰动而变大，因为我们不太可能遇到只有一小部分基因活动发生改变的情况。这与个体扰动较大的情况不同，个体扰动在增加信号的同时，具有引起显著非线性效应的不希望的、可能的后果。因此，使用多重遗传扰动是增加基因表达平均变化而不产生非线性效应的有效方法。我们的方法不同于其他每次只干扰单个基因的方案(13,25). 在这些早期的研究中，基因敲除被用于诱导基因表达的变化；然而，基因敲除可能导致二次代偿性改变，从而改变连接矩阵W公司相反，用合成基因网络（如拨动开关）扰乱基因动力学(6)，提供了一种快速且可控的方法来诱导基因表达的变化而不改变W公司.

就给定网络所需的实验数量而言，我们的反向工程算法非常有效。然而，这并不一定意味着计算效率。实际上，在所有可能的解决方案中搜索大值k个_最大值和N个因为这样的计划需要O（运行）(N个^{k个_最大值+1})初始化步骤中的计算。为了缓解这种计算效率低下的问题，我们的选择算法可以很容易地与其他计算方法结合，以搜索可能的解决方案。因为我们使用基因表达的微小变化作为初始选择标准（我们算法的第2步），所以我们可以执行几个初始扰动实验，而无需检查全部的一致的可能解决方案，从而减少N个_CPS公司在初始化步骤中。许多方法，如奇异值分解(26)和动态编程技术(27)，然后可用于从此类数据集构造一组可能的初始解（我们算法的步骤3.1）。通过使用精心选择的扰动，可以进一步减少剩余解的数量（我们算法的步骤4）。因此，我们可以根据网络的复杂性和可用的实验资源灵活调整计算效率。

我们的逆向工程算法的实用性取决于控制良好的基因表达技术的应用程度和基因表达测量的质量。扰动的大小受两个约束约束。一方面，如果扰动过大，基因网络将被推到其线性响应范围之外，这将削弱算法的线性假设。因此，预测网络中的误差将增加。另一方面，如果扰动过小，基因网络的响应将被仪器和生物噪声掩盖。可接受扰动的确切大小必须通过实验来确定。尽管扰动基因的表达是增加还是减少并不重要，但调控过表达技术更可取，因为其实现更简单，性能更可靠。除了遗传开关外，还有多种原核和真核表达系统可供使用，它们提供了可忽略的基线和所需的可控表达，包括Pbad系统(28,29)，ARGENT系统(30,31)和Pip系统(32). 定量实时PCR提供了可靠测量小扰动引起的基因表达细微变化所需的精度。实时PCR还提供了足够的吞吐量来有效测量100个或更多基因的反应。

有了基因调控网络的地图，我们可以从基因回路的角度提出有关疾病的新问题，并尝试操纵基因网络的功能输出。此外，我们可以有效地识别和验证特定的药物靶点，并消除化合物的影响。这可能导致基于细胞动力学网络方法的新型药物的开发。

补充材料

致谢

脚注

工具书类