基于时间过程数据和稳态通量曲线的代谢途径幂律建模的简化方法

摘要

背景

为了提高对代谢系统的理解，人们尝试从时间进程构建S系统模型。由于时间序列参数估计的非线性特性，通常使用非线性曲线拟合算法进行建模。然而，所需的大量迭代计算阻碍了大规模代谢途径模型的发展。为了解决这一问题，我们提出了一种新的方法，包括根据目标系统的雅可比矩阵对代谢途径进行幂律建模，并通过S系统的线性化对稳态通量剖面进行建模。

结果

两个分别模拟直线和分支路径的案例研究的结果表明，我们的方法减少了需要估计的未知参数的数量。传统动力学模型模拟的时间过程和我们方法描述的时间过程在代谢物浓度的大范围扰动下表现类似。

结论

该方法降低了计算复杂度，有利于构建代谢途径的大规模S系统模型，实现了目标系统雅可比矩阵和稳态通量分布的动态模拟模型逆向工程的实际应用。

背景

系统建模已成为理解代谢系统数学特性的有力工具。代谢测量技术的快速发展推动了建模的进步，特别是使用代谢产物浓度扰动影响的数据，这些数据包含有关代谢途径结构和调节的宝贵信息[1]. 用幂律近似表示酶促反应，称为生化系统理论，是理解代谢系统的有效方法[2,三]. 广义质量作用（GMA）和S系统[4,5]通常用作幂律建模方法，具有广泛的表示空间，允许充分表达酶动力学[6]尽管形式简单固定。此外，由于S系统形式的参数数量少于GMA形式，因此S系统是一个合适的建模框架。从给定的实验数据推导S系统模型不仅是理解非线性特性的有力工具，而且也是确定系统调控结构的有力工具[7,8].

由于S系统的非线性特性，从时间过程数据进行S系统建模通常很困难。非线性拟合算法，如遗传算法或人工神经网络，已被用于解决此问题[9–14]. 尽管这些方法可以应用于代谢途径，但在涉及许多紧密相连、未指定参数的目标模型的情况下，需要巨大的计算能力[13]. 此外，目标代谢途径的范围越广，由于参数搜索空间的扩展，出现局部极小值的可能性越大。网络结构分割方法可以减少遗传网络建模中的总参数搜索范围[15]但由于可逆反应之间的密切关系和变构调节，很难将其应用于代谢途径。Diaz-Serra和Fairén提出了一种基于稳态假设的方法，允许从系统的雅可比矩阵构造S系统模型[16]. 由于雅可比约束优化阶段未指定参数的搜索范围，这些作者的方法允许有效的参数估计。然而，需要估计的参数数量过多的问题仍未解决。

我们提出了一种基于目标系统雅可比矩阵和稳态通量曲线的代谢途径幂律建模方法，并对S系统进行了线性化。这减少了未指定参数的数量。两个数值实验表明，该方法生成的S系统模型描述了与传统动力学模型相似的动力学行为。

方法

从时间历程数据中检索雅可比矩阵

作为第一步，必须从代谢时间过程数据中获得雅可比矩阵。在本节中，我们总结了Sorribas的方法等[17]，我们在这项工作中使用。

在生物化学系统中，雅可比矩阵可以定义为：

哪里J型是雅可比矩阵，和δX（X）表示一个小扰动，包含浓度X_我作为其元素。雅可比矩阵的元素可以通过线性最小二乘拟合从扰动的时间历程中获得[17,18]. 该方法基于这样一个事实，即瞬态在稳态条件下对小扰动产生线性响应[5]. 其数学基础是，线性表示构成泰勒级数展开式的一阶项，在这种情况下，它足够精确。

S系统动力学级数的测定

S系统是由两个术语构成的幂律表示：生产率和降解率：

其中α_我和β_我是速率常数，g_ij公司和h_ij公司是动力学级数，和x_我表示化合物的浓度。

在稳态条件下，它可以简单地表示为：

通量=V（V）⁺=V（V） ^-    (3)

其中，通量是稳态通量与x的总和_我、和V⁺和V^-分别是生产和消费术语。

在稳态条件下，等式中的生产项和消耗项(2)因此可以表示为：

其中x_{j、 0个}表示x的稳态浓度_j个.

因此α_我和β_我产量：

定义X_我作为：

=年 _我 (x个₁, ...,x个 _n个 )     (我= 1, ...,n个)     (7)

等式的雅可比(2)可以表示为：

在稳态条件下，方程式(8)可以通过替换公式(5)带有α_我和等式(6)带有β_我，提供：

因此：

曾经的雅可比J型_ij公司，x个_{j、 0个}，和通量，公式(10)是确定动力学阶数g的约束条件_ij公司和h_ij公司Savageau将S系统的线性化描述为稳定性分析的“F因子”[19]. 我们使用这种表示来估计参数值。

在大多数情况下，g_ij公司和/或h_ij公司可从代谢途径的结构中获得；然而，在缺乏已知动力学阶数的情况下，需要进行参数估计来确定它们。在这种情况下，方程式(10)被用作限制参数搜索范围。

结果

案例研究1：线性生化途径

在第一个案例研究中，我们将此方法应用于已发布的酵母半乳糖代谢生化模型，如图所示1[20]. 该模型由五种代谢物和四种酶反应组成，主要用Michaelis-Menten方程描述。附录1中列出了动力学方程、系统参数和初始条件。

图1

半乳糖代谢途径[20]物质（X_1–5)代表参考文件中的以下代谢物：X₁，外半乳糖（GA_E类); X（X）₂，内部半乳糖（GA_我); X（X）_三，半乳糖1-磷酸（G1P）；X（X）₄，UDP-葡萄糖（U_{德国劳埃德船级社}); 和X₅，UDP-半乳糖（U_{通用航空公司}). 反应为V_信托收据半乳糖转运蛋白；V（V）_GK公司半乳糖激酶催化；V（V）_燃气轮机半乳糖-1-磷酸尿苷转移酶催化；V（V）_欧洲药典，由UDP-半乳糖4-差向异构酶催化。X（X）_三是速率为V的抑制剂_GK公司该模型的动力学方程、系统参数和初始条件如附录1所示。

在这种情况下，无需估计即可确定S系统模型的所有参数。生化模型可以转换为以下S系统形式。

自X起₁是一个自变量，在列出的速率方程中省略了它。

路径结构提供了以下约束：h₂₂=克₃₂，小时₂₃=克₃₃，小时₃₃=小时₄₃=克₅₃，小时₃₄=小时₄₄=克₅₄，小时₅₄=克₄₄，小时₅₅=克₄₅, β₂= α_三, β_三= β₄= α₅、和α₄= β₅稳态浓度为：X_1,0=0.50 mM（定值），X_2,0=0.146毫米，X_3,0=0.007毫米，X_4,0=0.817 mM，和X_5,0=0.243 mM。稳态通量为：V₁=V₂=V_三=V₄=0.081 mM/min。雅可比系数为：

我们的方法能够根据稳态代谢物浓度、稳态通量及其雅可比矩阵生成线性路径的S系统模型。例如，g₂₁和g₃₂由以下因素决定：

因为h₂₂等于g₃₂,

在此建模过程中，所有参数均通过此类简单计算确定。由此产生的S系统模型为：

案例研究2：分支生化途径

分支生化途径如图所示2在第二个示例中进行了测试。该途径具有生化途径的典型特征：分支、反馈调节和产物抑制（见附录2）。X（X）_三是两个V的抑制剂_三和V₄.

大多数参数都是通过我们的方法获得的。然而，一个参数无法通过计算确定，因此使用了参数估计方法。模型可以转换为以下S系统形式：

图2

分支路径.X（.X）_三是两种V的抑制剂_三和V₄.

路径结构提供了以下约束：h₁₁=克₂₁.稳态浓度为：X_1,0=0.067，X_2,0=0.049，X_3,0=0.081，X_4,0=0.041，稳态通量为：V₁=伏₂=0.1，伏_三=V₅=0.043和V₄=V₆= 0.057. 雅各布斯是：

一半以上的参数是通过使用稳态浓度、稳态通量和雅可比矩阵进行简单计算得到的。以下四个参数无法确定：α_三，克₃₃, β_三、和h₃₃例如，h₁₁是一个参数，可以通过以下方式确定：

因为g₃₃是一个无法通过计算确定的参数，需要通过某种约束提供估计。Asα_三，小时₃₃、和β_三可以被视为依赖参数，一旦g₃₃得到所有四个参数：

此外，g的搜索范围₃₃可能受到h定义的限制₃₃自X起为正值_三是反应V的底物_三.g的搜索范围₃₃设置为0到-0.95。当未指定参数的个数是搜索范围有限的参数时，可以使用线性拟合算法来拟合原始动力学模型中的时程数据。在这个模型中，g₃₃通过线性优化方法确定。因此，生成了以下S系统模型：

上述计算中使用的参数列表如表所示1.

表1案例研究2中S系统模型的参数

扰动后的瞬态动力响应比较

为了评估该方法的通用性，通过计算平均相对误差（MRE），将S系统模型响应各种扰动的瞬态动力学与原始Michaelis-Menten模型的瞬态动力学进行了比较：

其中n是模型中代谢物的数量，m是时间过程中采样点的数量，x'_我（t） 表示根据创建的S系统模型计算的时间进程，以及x_我（t） ，根据原始Michaelis-Menten模型计算的时间进程。本实验共有10个采样点，采样点间距为0.5，X₂案例研究1和X₁在案例研究2中，扰动的目标范围从0%到200%。代谢物对100%扰动的响应时间过程如图所示三在这两个示例中，S系统模型和参考模型中观察到类似的动态行为，作为对扰动的响应。

图3

一种代谢物扰动100%时代谢物的变化.（a）和（b）表示案例研究1中X时代谢物的变化₂是分别使用参考模型和S系统模型导出的扰动目标。（c） 和（d）是X时案例研究2中的变化₁是分别由参考模型和S系统模型给出的扰动目标。

MRE随扰动范围的变化如图所示4无扰动的初始MRE在1%以内。尽管在这两个案例研究中观察到MRE略有增加，但在200%的扰动范围内，它们仍保持在4%以内。

图4

平均相对误差（MRE）的变化对目标代谢物进行0%至200%的扰动。（a） 案例研究1中MRE的演变，其中X₂是扰动目标。（b） 案例研究2中MRE的演变，其中X₁是扰动目标。

讨论

根据代谢物浓度的时间进程数据进行幂律建模通常需要根据目标代谢网络的大小进行参数估计。特别是在开发大规模代谢模型时，这需要针对具体问题进行简化。我们的方法可以通过使用稳态通量剖面和目标系统的雅可比矩阵来减少未指定参数的数量，从而适用于大规模幂律建模。为了验证我们的方法，我们使用了两个现有的生化途径模型，即由动力学方程描述的直通路和分支通路模型。在S系统中，用12个参数对线性代谢途径进行建模（案例研究1），我们的方法准确地确定了所有参数，而Diaz-Serra和Fairen开发的方法[16]至少有四个参数未指定，其中包括错误校正参数。在具有16个参数的分支代谢途径的情况下（案例研究2），我们的方法确定了除一个参数外的所有参数，而Diaz-Siera和Fairén的方法有12个未指定的参数，其中包括误差校正参数。案例研究表明，我们从时间序列中建立S系统模型的方法可以比以前报告的方法更有效地减少未指定参数的数量[16]. 此外，扰动响应实验表明，我们的方法创建的模型可以再现类似于参考模型的动力学，因为当扰动范围为200%时，MRE约为3%（图。4). 用我们的方法生成的S系统模型可以在很宽的稳态范围内提供精确的仿真。与许多幂律代谢模型一样，这种限制并不妨碍代谢途径的建模和分析[4,5].

对实验噪声的鲁棒性是建模方法实际应用的重要要求。由于雅可比矩阵对时间序列中的噪声相当敏感[18]，对雅可比矩阵中包含的数值误差作出稳健的校正响应对于开发涉及其使用的模型非常重要。采用Diaz-Serra和Fairén的方法[16]在S系统模型中加入误差修正参数，以减小实验噪声的影响。然而，这种方法可能会导致未指定参数的数量增加。作为降低时间历程数据中实验噪声的另一种方法，使用非线性神经网络模型估计适当的斜率是有效的[10]，并可以提供一个误差受控的时间过程，从而能够高精度地获得雅可比矩阵。已经讨论了获得一般类型扰动具体影响的精确估计的方法[21]并且可以对在扰动下获得的时间尺度代谢物浓度等数据进行更精确的分析。

为了评估我们的方法对实验噪声的鲁棒性，我们测量了手动将数值误差插入雅可比通量或稳态通量时的计算误差。表2总结了S系统和原始Michaelis-Menten模型的模拟时间轨迹的MRE。显然，雅可比矩阵对数值误差非常敏感，因为它直接影响雅可比函数的动力学阶数。此外，稳态通量分布数据可能包含实验噪声。然而，用我们的方法建立的模型对雅可比通量和稳态通量的误差相对稳健，表明这些方法在实践中是有用的。

表2作为实验误差函数的平均相对误差（MRE）。案例研究2中测量了S系统模型和原始模型的MRE。从模型中获得的时间历程数据包括状态变量的扰动，以检查S系统模型和原始Michaelis-Menten模型之间的动态响应差异。雅可比通量和稳态通量的再现具有100%的数值误差。雅可比矩阵的数值误差平均地插入到雅可比函数的所有元素中。X（X）₁是扰动的目标。时间进程数据来自S系统模型，其中X₁对50%的增长感到不安。为计算抽取了10个时间点，它们之间的间隔为0.5 s。MRE是根据S系统模型和原始Michaelis-Menten模型中的时间进程数据计算得出的。在分支生物化学途径的情况下（案例研究2），雅可比增加了100%，并且三个稳态通量J_1–2、J_3–5、和J_4–6表示通过V的通量₁和V₂，通过V的通量_三和V₅和通过V的通量₄和V₆分别为。

索里巴斯（Sorribas）和卡斯坎特（Cascante）提出了一种使用一组给定的对数增益测量值来识别调节模式的方法[7]. 在他们的论文中，他们建议选择可能模式的一种策略是进行扰动实验并测量相应的动态响应。对于这种方法的实际应用，开发适当的方法来执行所需的实验至关重要；然而，由于缺乏详尽的测量方法，测量各种稳态通量产生的对数增益是不现实的。我们的方法可以通过使用一组合适的稳态通量曲线来开发相当准确的代谢途径模型。尽管稳态通量剖面数据仍然难以直接和全面测量，但已经开发了几种使用同位素测量稳态通量剖面的方法[22,23]. 我们的方法假设获得的雅可比矩阵是准确的。因此，在雅可比矩阵以线性方式表现的适当小扰动后，获得反映瞬态动力学的时程数据是很重要的[17].

随着代谢测量技术的进步，综合代谢组数据无疑会积累起来[24–26]. 在我们的实验室中，我们开发了一种使用毛细管电泳质谱的高通量技术，该技术可提供涉及数百种离子代谢物的有效时间过程数据[27,28]. 我们的方法有望通过利用累积代谢组和稳态通量剖面数据以及代谢组测量技术的预期发展，为大规模代谢途径提供高通量建模。

结论

我们的方法提供了代谢途径的稳定和高通量S系统建模，因为它通过使用稳态通量剖面和从时间过程数据中检索的雅可比矩阵，大大减少了未指定的参数。该方法生成的S系统模型可以在稳态点附近的大范围内提供精确的仿真。结合代谢组测量技术，它应该允许对大规模代谢途径进行高通量建模。

附录

附录1：酵母半乳糖代谢生化模型

Atauri构建了酵母半乳糖利用途径模型等. [20]，其反应图如图所示1.

模型的速率方程为：

其中速率表达式为：

使用的参数可参考[20]. 计算的稳态条件为：X₁; 0.5 mM（固定），X₂; 0.146毫米，X_三; 0.00703毫米，X₄; 0.817毫米，X₅; 0.243 mM，以及通过路径的流量；0.0081毫摩尔·秒^-1.

附录2：分支生化途径

分支模型的速率方程，其反应方案如图所示2是：

=V（V）₁-V（V）₂

=V（V）₂-V（V）_三-V（V）₄

=V（V）_三-V（V）₅

=V（V）₄-V（V）₆

其中速率表达式为：

V（V）= 0.1