厦门论文发期刊杂志发表投稿-同位素辅助代谢通量分析作为等式约束的非线性程序，可提高可扩展性和鲁棒性

丹尼尔?卢格，加内什·斯里拉姆

出版日期： 2022年03月24日

抽象

稳定同位素辅助代谢通量分析（MFA）是估计代谢网络中碳流量和分配的强大方法。MFA的核心是一个参数估计问题，其中通量和代谢物池大小是模型参数，通过优化进行估计，以考虑稳态或同位素非平稳同位素标记模式的测量。随着MFA问题的规模化发展，它们需要高效的计算方法来实现快速而稳健的收敛。MFA 问题的结构使其能够被转换为相等约束的非线性程序 （NLP），其中相等约束是从 MFA 模型方程构造的，目标函数定义为模型预测和一组标记测量值之间的平方残差之和 （SSR）。这种NLP可以通过使用代数建模语言（AML）来解决，该语言提供了最先进的优化求解器，以实现鲁棒的参数估计和对大型网络的卓越可扩展性。以这种方式实现时，执行优化时不区分状态变量和模型参数。在此类优化的每次迭代期间，系统状态都会更新，而不是从头开始显式计算，并且这与模型参数估计值的改进同时发生。这种优化方法与传统的"射击"方法形成鲜明对比，在传统的"射击"方法中，状态变量和模型参数保持不同，并且在逐步优化的每次迭代期间重新计算系统状态。我们的NLP公式使用Wiechert等人的MFA建模框架[1]，该框架适合将模型方程合并到NLP中。NLP约束由基本代谢物单元（EMU）或cumomers上的平衡组成。在这个公式中，稳态和同位素非稳态MFA（inst-MFA）问题都可以作为NLP来解决。对于 inst-MFA 情况，使用并置将描述标记动态的常微分方程 （ODE） 系统转录为 NLP 的代数约束系统。可以使用在AML上实现的NLP求解器有效地解决这种大规模的NLP。在我们的实现中，我们使用了在通用代数建模系统（GAMS）中实现的简化梯度求解器CONOPT。NLP框架对于inst-MFA特别有利，可以很好地扩展到具有许多自由参数的大型网络，并且与在每次迭代时计算系统状态和灵敏度的拍摄方法相比，具有更强大的收敛特性。此外，当使用cumomer框架时，这种NLP方法支持将串联MS数据用于稳态和inst-MFA。我们组装了一个软件eiFlux，它是用Python和GAMS编写的，它使用NLP方法，同时支持稳态和inst-MFA。我们在几个例子中证明了NLP配方的有效性，包括基因组规模的inst-MFA模型，以突出这种方法的可扩展性和稳健性。除了典型的inst-MFA应用程序外，我们预计该框架和我们的相关软件eiFlux对于将inst-MFA应用于复杂的MFA模型特别有用，例如为真核生物（例如藻类）开发的模型以及与多种细胞类型的共培养物。

作者简介

同位素辅助代谢通量分析（MFA）是一个计算密集型参数估计问题。同位素非稳态 MFA （inst-MFA） 代表了计算最繁重的 MFA 应用程序。我们将稳态和inst-MFA问题的公式呈现为等约束非线性程序（NLP），由代数建模语言中实现的最先进的求解器解决。我们表明，与传统方法相比，这种公式导致鲁棒的收敛特性，特别是对于inst-MFA。我们开发了一个软件eiFlux，它使用NLP公式来执行稳态和inst-MFA。我们在几个例子中展示了eiFlux的应用，包括基因组规模的inst-MFA模型，并表明即使从对参数的非常糟糕的初始猜测开始，它也具有鲁棒的最佳收敛性。eiFlux是使用Python编程语言和通用代数建模系统（GAMS）实现的，使用CONOPT求解器。eiFlux可根据要求提供，等待机构批准，并免费供学术使用。

引文： Lugar DJ，Sriram G（2022）同位素辅助代谢通量分析作为平等约束的非线性程序，用于提高可扩展性和鲁棒性。PLoS Comput Biol 18（3）：e1009831。https://doi.org/10.1371/journal.pcbi.1009831

编辑 器： Costas D. Maranas，宾夕法尼亚州立大学，美国

收到： 2021年6月3日;接受： 一月 12， 2022;发表： 三月 24， 2022

版权所有： ? 2022 卢格，斯里拉姆。这是一篇根据知识共享署名许可条款分发的开放获取文章，该许可允许在任何媒体上不受限制地使用，分发和复制，前提是注明原始作者和来源。

数据可用性： 数据可用性 重现我们的结果和 eiFlux 演示版本所需的所有数据和模型都可以在手稿或随附的补充文件中找到。有关更多详细信息，请参阅 S1 文件的标题。代码可用性 我们在 S1 文件中提供了 eiFlux 的演示版本，该版本专门运行本文中讨论的所有模型以及所有必要的输入文件。这可以作为编译的Python代码，eiFlux_Limited.pyc。在运行时，此编译的代码打印 NLP 中使用的模型参数（例如，描述问题的数组），以便完全可验证结果。此外，eiFlux 使用的 GAMS 脚本位于 S1 文件中GAMS_Model文件夹中。这些脚本包含用于稳态和内部 MFA 问题的常规 NLP。这些脚本仅供参考。编译后的 Python eiFlux_Limited.pyc 在后台使用这些脚本。最后，作者正在与他们的机构合作，授权完整版的eiFlux，使用户能够定义自己的模型。在许可过程完成后，此完整版本将免费提供给学术和非营利研究人员。以下 GitHub 存储库将指示许可过程的状态，并列出使用 eiFlux 的说明：https://github.com/SriramLabUMD

资金： 这项工作由美国国家科学基金会部分资助（授予GS的编号MCB-1517671）。这项工作还得到了国家需求领域研究生援助（GAANN）奖学金（奖项编号P200A180093），马里兰大学大脑与行为研究所的大脑和行为倡议种子赠款以及美国农业部国家食品和农业研究所的农业和食品研究计划竞争性赠款No. 2019-67015-29412的部分支持。资助者在研究设计，数据收集和分析，出版决定或手稿准备方面没有任何作用。

相互竞争的利益： 作者宣布不存在相互竞争的利益。

这是一篇PLOS计算生物学方法论文。

介绍

稳定同位素辅助代谢通量分析（MFA）是一种用于实验量化体内代谢通量的强大方法[2-4]。这种方法涉及两个部分 - 实验数据收集和计算通量估计。实验上，将同位素标记的底物送入细胞培养物。然后测量取决于通量的细胞内代谢物的标记测量，通常使用质谱（MS）进行测量。然后将碳原子重排的代谢网络模型拟合到测量的标记模式（称为质量同位素分布（MID），以估计最能解释这些测量的通量分布。标记测量值可以在稳态标记条件或非稳态（同位素非稳态）标记条件下收集，其中时间序列测量值是在稳态方法上收集的。在稳态情况下，通量为未知参数的代数模型适合于标记数据。在同位素非平稳的情况下，该模型由常微分方程组（ODE）组成，其中通量和代谢物池大小是未知参数。这是一个计算密集型参数估计问题，需要使用专用软件。

几种建模公式，特别是那些涉及累积体[1]和基本代谢物单元（EMU）[5]的公式，已经应用于通过以分段线性形式求解模型方程来减轻计算负担。其他试图通过以替代形式表达状态变量来减轻计算负担的公式包括键聚体[6，7]（对于涉及单个，均匀标记的碳源的特定类型的MFA很有用）和通聚体[8]。诸如此类的公式导致稳态MFA问题变得可以有效解决。EMU方法向同位素非稳态MFA（inst-MFA）的扩展使得非稳态问题在计算上易于处理[9]。

MFA 应用程序通过求解参数估计问题来估计通量，该参数估计问题可最大限度地减少实验测量的标记状态（或 inst-MFA 情况下的动态标记状态 [DLS]）之间的平方残差 （SSR） 之和。这种最小化需要一种优化算法来估计最能解释实验DLS的最佳通量（以及inst-MFA的代谢物池大小）。据我们所知，所有可用的inst-MFA软件都采用"拍摄方法"来确定最佳效果。在该方法中，通过顺序算法执行参数估计，其中选择一组模型参数，然后使用该组参数值的数值积分方案评估DLS。在每次迭代时，使用高斯-牛顿方法（例如Levenberg-Marquardt）或无导数方法（例如模拟退火[10，11]）对参数估计进行细化，并重新计算DLS并将其与测量值进行比较。重复这些步骤，直到达到最佳状态。

参数估计过程（如 MFA）的理想质量是鲁棒性。稳健的估计将是准确（收敛于全局最优的邻域）和可靠的（从任意起点收敛到相同的最优）[12]。在每次迭代时重新计算DLS的拍摄方法并不可靠，因为它受到一个小的收敛区域[13]的影响，特别是对于具有大量自由参数的问题。这通常是由于存在许多局部最优，通常远离全局最优解决方案。在这种情况下，对于一般网络，向全局或接近全局最优的鲁棒收敛通常需要对模型参数进行非常好的初始猜测。但是，在 inst-MFA 中，池大小和通量通常难以先验估计，并且对于此问题中的大量可用参数，通常没有良好的初始猜测。因此，优化通常必须从可行参数空间中的各个点重新启动多次，以增加找到全局最优解的机会。

为了解决这些局限性，我们将 MFA 问题的公式作为等约束非线性程序 （NLP） 提出。在NLP公式中，EMU或cumomer平衡构成了一个相等约束系统，必须保持其可行性，同时最大限度地减少测量和模型预测之间的误差。对于 inst-MFA，我们使用并置 [14] 将常微分方程组 （ODE） 离散化为适合于最先进的 NLP 求解器的形式。这种离散化通常被称为转录[14]。在NLP方法中，参数估计和状态预测在优化期间同时发生，这与使用每次迭代的参数值从头开始计算状态的传统方法形成鲜明对比。使用在通用代数建模系统（GAMS）[17]（一种代数建模语言（AML））中实现的NLP求解器CONOPT [15，16]可以有效地解决NLP问题。我们采用Wiechert等人提出的cumomer建模公式[1]，非常适合将模型方程合并到NLP框架中。此外，我们还使用此建模公式来开发EMU平衡方程。在这个公式中，EMU分数丰度是状态变量，而不是传统上的EMU向量[5，9]。使用在 AML 上实现的求解器（如 GAMS）求解 NLP 可以高效、可靠地解决 MFA 和内部 MFA 问题。AML具有自动微分和有效使用稀疏线性代数等功能。因此，它们有效地集成了最先进的优化算法，并使其适用于大规模问题。

我们表明，NLP公式对于inst-MFA特别有益，其中并置用于转录ODE系统，从而产生鲁棒收敛和相对较少的局部最优，即使对于大型网络模型也是如此。这与使用拍摄方法进行内部MFA参数估计的传统方法形成鲜明对比，并且往往会受到小区域收敛的影响。NLP公式还可以很好地扩展到具有大量独立参数的问题[10]，为具有许多独立通量和池大小的大型网络模型的通量估计提供了一个计算有效的框架，例如真核系统模型（例如藻类），共培养条件或基因组规模模型中遇到的通量。此外，由于同位素和聚单体很容易映射到串联MS测量[18]，因此该公式支持将串联MS数据用于稳态和内联MFA。

我们开发了一个软件eiFlux（用于isotope辅助代谢通量分析的相等约束非线性编程），它使用NLP公式并支持稳态和inst-MFA（图1）。使用GAMS Python API，eiFlux使用Python编程语言为给定代谢网络的模型方程组装矩阵和向量，然后使用NLP求解器CONOPT解决GAMS中的MFA问题，将模型方程强制执行为相等约束。代谢网络和碳原子重排以类似于其他软件的用户友好方式指定，使模型易于在eiFlux和现有软件之间转移。

thumbnail 下载：

个人电脑幻灯片

巴新放大图片

断续器原始图像

图 1. 有几个因素使eiFlux在当前可用的MFA软件中独一无二。

eiFlux 使用相等约束的非线性规划来解决稳态和内部 MFA 问题。在NLP配方中使用搭配处理ODE系统可产生鲁棒的收敛特性。支持cumomer建模方法的eiFlux可以直接拟合串联MS数据。

https://doi.org/10.1371/journal.pcbi.1009831.g001

方法

作为等约束非线性程序的稳态 MFA

为了在任意代谢网络上执行稳态MFA，应将积液或EMU天平，测量映射矩阵和SSR目标函数组装成一般模型结构。我们首先介绍如何为cumomer天平开发此模型结构，然后讨论EMU框架如何适应此结构。我们采用了Wiechert等人开发的系统程序[1]，该程序使代谢网络及其相关的碳原子重排能够有效地转化为cumomer平衡系统。M?llney等人[19]描述了如何使用映射矩阵将同位素和cumomer分布映射到实验数据集。基于此公式的稳态 MFA 问题的 NLP 如 （Eq 1） 所示。

(1)

NLP的目标函数（Eq 1）包括（i）测量的同位素标记数据和模拟数据之间的残差项，其中n是测量指数;（ii）任何测量通量和计算通量之间的残差的术语，其中r是通量指数。在这里，σn和σv，r分别是同位素标记测量和通量测量的标准偏差。相等约束是 （i） cumomer 或 EMU 平衡，其中 i、j 和 k 是 cumomer 或 EMU 指数，以及 Qrkij、Prkj，并且是定义网络碳原子重排平衡的数组，由Wiechert等人开发[1];（ii） 化学计量通量天平，S二维码是化学计量矩阵;（iii） 对（可选）测量比例因子[19]的约束，以等于其相应的片段比例因子;（iv） 将积木或动车组映射到测量值;（v）如果使用积木框架，则定义零阶积木。不等式约束包括（i）通量上的任何线性不等式约束，包括可逆反应的净通量的边界;（ii）通量的上限和下限。所有通量值必须是非负的，因此对于可逆反应，正向反应和反向反应必须各自视为独立的。

如果用户拟合完整且规范化的 MID，则这些比例因子应固定为等于 1。矩阵 U 的元素数控在此约束中定义在 （Eq 2） 中。

(2)

这个NLP也可以使用EMU框架来定义。在这种情况下，EMU 分数丰度是状态变量，而不是积木分数。在 Antoniewicz 等人提出的 EMU 分解之后 [5]，模型数组，Prkj、和 M新泽西州可以有类似的定义。将动车组存储在复合向量 x（通过连接动车组向量形成）中，定义数组 Q：

(3)

这里是代谢物的化学计量系数，对应于指数k与指数r反应的EMU分数。数组 Prkj同样可以定义为：

(4)

阵列，其元素对于代谢物进入网络的反应不为零，定义如下：

(5)

测量映射矩阵 M新泽西州定义为：

(6)

最后， x因普定义为输入复合动车组矢量。对于馈送的同位素标记化合物，根据同位素富集的模式定义相应的EMU元素。其余的EMU元素是根据天然同位素丰度来定义的。

在我们的实现中，这个NLP是使用CONOPT降低梯度算法求解的，该算法在GAMS中实现[15，16，20]。CONOPT 算法适用于相等约束总数与变量总数相似的优化问题。在简化梯度算法的每次迭代中，变量被划分为两组：基本和非基本。此分区由求解器自动执行，无需区分参数和状态变量。在算法的每次迭代中，非基本变量被传播到一个更好的点（一个改进目标函数的点），并且基本变量被传播到新的点，使得它们满足相等约束[16，20]。重要的是，基本变量的数量与相等约束的数量相同。因此，基本变量的值由这些约束确定，并且任何约束冲突都将在给定的 CONOPT 求解器迭代中使用牛顿方法的迭代进行更正。非基本变量的数量与自由参数（通量、池大小和比例因子（如果使用））大致相同，因此在以这种方式制定的典型 MFA 问题中，通常只有 O（100） 个非基本变量。在 CONOPT 中，仅计算用于优化的非基本变量的约简黑森值形式的二阶导数信息。

因此，对于相等约束数与变量数相似的非线性程序，非基本变量较少，并且算法的行为效率更高。降低梯度算法非常适合 MFA 问题，其中典型网络具有数百到数千个累积体或 EMU 分数丰度，每个累积体都贡献一个相等约束。有关将通用简化梯度算法用作简单示例的简单分步示例，请参阅 S1 文本。重要的是，一般的NLP求解器（如CONOPT）不会将模型参数（即通量）与模型变量（即累积体和模拟测量值）区分开来，并将两者都视为优化变量。

因此，Cumomer或EMU余额不会在每次迭代中从头开始求解，也不会以前面提到的分段线性方式求解。该算法以随机或用户输入的可行通量分布进行初始化。从这里开始，算法在一个方向上从一个点移动到另一个点，使得约束（其中积木或EMU余额是一个大子集）在统计上仍然可行（到一个小公差）。随着迭代的进行，变量值将更新以保持约束可行性。我们用一个示例NLP来说明这种方法，其中迭代轨迹可以很容易地可视化到最优（图2）。S1 文本提供了有关此示例中解决的 NLP 的详细信息。在这里，g1(x） 和 g2(x） 是非线性约束。该算法从初始可行点 x 开始0.从这里开始，它以这样的方式通过连续的候选点进行，使得约束g1(x） 和 g2(x）保持可行，并改进目标函数，直到达到最优。最优值是满足这些约束的目标函数的最小值，通常不同于无约束的最小值。

thumbnail 下载：

个人电脑幻灯片

巴新放大图片

断续器原始图像

图 2. 对于一个简单的说明性示例，显示了通用简化梯度算法的前六次迭代（有关详细信息，请参阅 S1 文本）。

在减少梯度算法的每次迭代中，变量都会传播到约束仍然可行的改进点。在此示例中，可行区域（以蓝色阴影显示）位于约束 g 之上1和低于约束 g2.在这种情况下，迭代将沿着 g 进行1非线性约束边界。

https://doi.org/10.1371/journal.pcbi.1009831.g002

CONOPT 能够在优化的每次迭代期间有效地处理非线性约束，同时采用相当大的步长。这是因为作为减少的梯度求解器，CONOPT 使用相等约束来显著减小有效问题大小。回想一下，相等约束由 EMU 或 cumomer 天平产生，参数数为 O（102） 对于典型模型。因此，NLP中的变量总数仅比相等约束的总数多几百个。CONOPT 求解器选择一组 O（102）变量（非基本），其余变量（基本）根据相等约束精确确定。这有效地将问题减少到仅在非基本变量中为一个，其中只有O（102).基本变量的值（其数量与约束相同）在每次迭代时使用牛顿迭代进行更新，这些迭代涉及约束雅可比的反函数，直到在公差内满足相等约束。CONOPT 中最大和最小可行性公差的默认值为 1x10-7和 4x10-10分别。在 eiFlux 中求解的模型使用这些默认值。通过这种方式，CONOPT求解器能够在有效大小为O（10）的优化问题的每次迭代期间保持大量相等约束的可行性。2） 变量。在实践中，使用 CONOPT 求解的 MFA 模型在优化迭代期间很少会失去可行性。在失去可行性的迭代中，不可行性通常非常小，CONOPT能够通过移动到附近的可行点来快速恢复可行性。有关 CONOPT 保持约束可行性的方法的数学详细信息，请参见 S2 文本。

作为相等约束非线性程序的 Inst-MFA

在inst-MFA中，代谢通量使用细胞内代谢物中同位素标记的时间序列测量来估计。在这种情况下，系统不处于同位素标记稳态，并且从瞬态标记动力学中提取通量信息[3，4]。这与稳态 MFA 不同，稳态 MFA 从稳态标记数据中提取通量估计值。虽然目标是相同的 - 使用同位素标记数据准确估计细胞内通量 - 但inst-MFA需要将ODE系统拟合到数据中，而不是像稳态情况下那样的代数方程。这带来了额外的计算负担。此外，通量和代谢物池大小都是inst-MFA中的模型参数。这与稳态 MFA 形成鲜明对比，后者对池大小不敏感。

在各种情况下，Inst-MFA 可能比稳态 MFA 更可取。对于稳态标记实验，细胞倍增时间可以被认为是在细胞内代谢物中达到稳态标记所需时间的下限[4]。当倍增时间很长时，特别是对于哺乳动物细胞来说很常见，在这个时间间隔内保持代谢稳定状态可能很困难。由于inst-MFA仅要求系统保持代谢稳定状态，因此实验的时间长度可能要短得多，从而克服了一些实验困难[3]。

光自养代谢是一种重要情况，其中inst-MFA不仅优于稳态MFA，而且是必需的。在这种情况下，稳态标记谱不提供代谢通量的信息[3，4]。所有通量信息都包含在首次提供细胞培养物之间的标记瞬变中13C标记的碳源及其达到同位素标记稳定状态的时间。在这些系统上成功实施inst-MFA可以为光合生物（如藻类）提供重要的工程见解。

对inst-MFA的第一次数学处理[21，22]提出了一种使用cumomer框架模拟DLS实验和参数灵敏度的方法。随后的工作[23]表明，通过使用伴随方法计算测量残差的梯度，可以加速这些计算。目前，有两个软件能够执行inst-MFA——INCA [24]和OpenMebius [25]。所有这些方法都在优化例程的每次迭代中模拟DLS。例如，INCA使用同位素非平稳动车组方法[9，24]在Levenberg-Marquardt算法中进行DLS仿真以进行参数估计。将 inst-MFA 问题视为受相等约束的 NLP 将提供许多好处，包括健壮性和可扩展性。但是，这要求将 inst-MFA 中的 ODE 系统转换为代数方程。

为了开发这样的公式，我们构建了一个框架，该框架使用并置[14]将任何inst-MFA问题转录为一组代数约束，以便合并到一般的NLP框架中。正如Shin等人所描述的[10]，这种方法对于动态系统的大规模参数估计问题有几个重要的好处。具体而言，这种方法允许大型模型由AML有效处理，例如GAMS，它采用自动微分方案来计算导数，以及稀疏线性代数算法，可有效处理通过OTE转录获得的代数方程的结构化系统[10]。此外，这种方法可以很好地扩展到具有大量自由参数的问题[10]。由于池大小和自由通量都是inst-MFA中的模型参数，因此在考虑分区的真核系统（如藻类或基因组规模模型）时，参数总数很容易超过100。基于传统拍摄方法的方法涉及灵敏度方程的仿真，不能有效地扩展到具有大量自由参数的问题，从而阻碍了它们对此类问题的可用性[10，26]。在这种情况下，NLP框架提供了一个明显的优势。

除了这些考虑因素之外，使用高斯-牛顿型（例如Levenberg-Marquardt）优化算法的拍摄方法在为参数值提供较差的初始猜测时没有鲁棒收敛性[13]。在这种情况下，这种方法往往会收敛到远离全局最优的较差局部最优。在 inst-MFA 中，对于模型中的大多数或全部通量和池大小，通常没有良好的初始猜测。在测量残差较大的情况下，这种方法也可能具有缓慢的收敛[27]。虽然在非凸优化问题（例如 inst-MFA）上使用任何凸优化方案时必须考虑局部最优，但使用并置的 NLP 框架具有相对稳健的收敛属性 [14]，即使对参数的初始猜测很差。要搜索全局最优值，通常从多个（随机）初始起点重新启动优化。这种强大的收敛意味着与拍摄方法相比，用户将需要执行更少的随机重启来找到全局（或接近全局）的最佳状态。诸如无导数搜索方法之类的替代方案，其中包括模拟退火和遗传算法，可以避免局部最优，但不能很好地扩展到具有大量自由参数的参数估计问题[10]，从而阻碍了它们对大型MFA模型的适用性。

描述积木或动车组动力学的ODE系统如（等式7所示）。这类似于（Eq 1）中的第一个约束，它由稳态下的余额组成。每个积木或动车组分数（指数 k）的动力学由 ODE 描述。在此等式中，pk是对应于积液或EMU级分x的代谢物的池大小k.所有余额的集合共同构成了耦合的非线性 ODE 系统。

(7)

在NLP公式中，搭配用于转录ODE[14]。使用并置方法，对时域进行划分，多项式在每个时间间隔内近似于 ODE 解。使用并置将ODE转录为代数方程组允许使用通用NLP求解器（例如CONOPT）解决参数估计问题。在eiFlux中，我们选择使用Radau IIA正交搭配方法，因为它具有刚性稳定性[28，29]。搭配方法还具有相应的完全隐式 Runge-Kutta 方法 [28，29]。有关搭配方法的推导及其与Runge-Kutta方法联系的详细说明，请参见Huynh [29]。使用 5 的简单示例千-订单 Radau IIA 搭配方法在 S3 文本中提供。对于形式为（Eq 8）的ODE系统，相应的完全隐式Runge-Kutta方法需要在每个时间间隔内求解（Eq 9）和（Eq 10）。

(8)

(9)

(10)

在 （公式 9） 和 （公式 10） 中，a伊杰指的是 Runge-Kutta 矩阵的元素，对于完全隐式方法（如 Radau IIA 方法），每个元素都可能不为零。c 的值我是区间 [0，1] 和 b 上的并置节点j是权重。s 的值是方法中的阶段数。表 1 总结了用于转录 inst-MFA 问题的设置索引。

thumbnail 下载：

个人电脑幻灯片

巴新放大图片

断续器原始图像

表 1. 用于开发 inst-MFA 问题的 NLP 的索引集。

希腊字母是用于离散化ODE系统的索引。

https://doi.org/10.1371/journal.pcbi.1009831.t001

选择一组时间节点，一般的Runge-Kutta方案，将整个时域离散化如（Eq 11）和（Eq 12）。在这里，hγ = tγ?tγ?1.

(11)

(12)

方程 （Eq 11） 可以通过引入变量 K 分为两个方程αγi.

(13)

(14)

将 （Eq 12）、（Eq 13） 和 （Eq 14） 应用于 （Eq 7） 中的稳态平衡，可为转录的 ODE 系统提供以下代数约束系统：

(15)

要近似计算积木或 EMU 分数值，请 xμi，在某个时间点 tμ如果在γ的时间间隔内进行测量，多项式基函数 Bα(τ）在（Eq 16）中定义的搭配方法，用于开发（Eq 17）。这些基函数是通过积分拉格朗日多项式找到的。对于 Radau 并置方法，这些拉格朗日多项式是从相应的 Radau 多项式 c 的零构造的1...cs [28，29]。有关使用 Radau 并置来近似于 ODE 系统的解决方案的示例，请参阅 S3 文本。S3 文本中的公式 S2.7 列出了基函数 Bα(τ） 显式用于 3 阶段方法。

(16)

(17)

（Eq 15）和（Eq 17）一起用于生成具有转录ODE系统（Eq 18）的NLP。作为参考，表 1 中定义了开发 NLP 中使用的所有索引集。为清楚起见，用于转录 ODE 系统的集合使用希腊字母表示。

(18)

在这里，基质D用于约束相同代谢物的所有累积体的池大小相等（Eq 19）。因此，代谢物池大小的向量映射到积木或EMU馏分池大小的向量，pk (等式19）。

(19)

初始条件 x0i被指定为输入积木或动车组分数。在 eiFlux 中，进纸标签以 x 为单位指定因普.未指定标记模式的代谢物被假定为每个原子具有天然同位素丰度。那些指定了标记模式的原子被假定具有未指定为标记的原子的天然同位素丰度。初始条件的这种定义允许代谢物指定其初始标记。这允许灵活地定义同位素标签的脉冲输入（作为通常步骤输入的替代方案），这对于将标签推注添加到培养物中的建模场景可能很有用。测量时间点处的积液或EMU值（具有指数μ）使用测量映射矩阵M映射到代谢物标记测量值mnμ.池大小受用户定义的上限和下限的约束。目标函数还包含一个术语，用于表示任何测量池大小的测量残差。

支持串联 MS 数据。

串联MS（MS/MS）是传统"单"MS的扩展，它比单个MS提供更多关于标记状态的信息[30，30-32]，因此有可能提高通量估计的准确性。在MS / MS（S1图）中，分析物串联通过两个质量分析仪，并在两个分析仪之间发生碎片。第一个质量分析仪选择母离子的特定质量，然后将其碎片化以产生具有不同质量的子离子。这些子离子在第二个质量分析仪中进行分析。与单个MS相比，母离子和子离子MID的组合解释提供了有关标记状态的更丰富的信息。eiFlux 的一个新特性是，它支持稳态 MFA 和 inst-MFA 的串联 MS 数据。此功能具有优势，因为MS/MS数据对MFA特别有用[32]。

尽管之前大多数MFA报告都是基于气相色谱（GC）-MS分析，但液相色谱（LC）-MS / MS分析的使用正变得越来越普遍。LC-MS / MS的两个明显好处是（i）它不需要分析物的衍生化，以及（ii）它通常可以比GC-MS测量更多的中枢碳代谢中间体，特别是那些衍生效率低下的人。使用LC-MS / MS进行inst-MFA的研究人员通常仅使用母体离子测量（MS数据），而不是母体和子离子测量（MS / MS数据），因为目前没有inst-MFA软件来适应MS / MS数据。

结果

玩具网示例

我们给出了一个简单的示例来说明 eiFlux 如何使用并置进行 inst-MFA 通量和池大小估计。对于此示例，我们模拟了一组给定参数（通量和池大小）值的测量值，因此最佳估计问题的确切解是先验已知的。在这个玩具网络中（图3），代谢物A可以通过三种不同的途径转化为代谢物B，其中一种途径是可逆的。每条途径都有一个独特的碳原子重排。假设在t = 0时，将细胞培养物喂入含有100%U-13断续器xt.细胞内代谢物A，B，C，D，E和F将逐渐富集，随着t→∞而完全富集。在这种情况下，有关磁通量和池大小的所有信息都包含在标记瞬态中。因此，稳态 MFA 不可用，并且需要 inst-MFA 来估计这些参数。

thumbnail 下载：

个人电脑幻灯片

巴新放大图片

断续器原始图像

图 3. 玩具网络显示了具有化学计量通量约束的反应和碳原子重排。

灰色虚线表示单元格边界。由于 v0假设是确切已知的，基于化学计量约束有三个独立的通量。还有六个细胞内池大小参数。

https://doi.org/10.1371/journal.pcbi.1009831.g003

这个玩具网络非常简单，使用MATLAB（MathWorks，Natick，MA）可以很容易地配制和数值积分同位素非平稳的积木天平。因此，为了演示INST-MFA的NLP方法，我们通过对积量天平与设定的参数值进行数值积分来生成准确的合成数据。对于此演示，所选参数值显示在 S4 文本中。合成数据是为这个玩具网络生成的，假设没有自然13C丰度和100%13C A 的富集xt在 t = 0 时。[0 20] 所述综合测量值取自合成数据，时间间隔为2。这些合成测量值用于解决逆问题 - 从片段B[12]，B[23]和B[123]的合成时间序列测量中估计通量和池大小。有关用于模拟数据的通量和池大小、eiFlux 收敛的解决方案以及模型拟合数据的详细信息，请参阅 S4 文本。S1 数据包含拟合到此模型的模拟数据。

由此产生的代谢物B的时间依赖性标记曲线如图4所示。该 5千-顺序 Radau IIA 并置点在每个时间间隔内显示为标记。在此示例中，时域被拆分为五个区间，每个区间中有三个并置点。t = 0 处的点表示初始条件。使用inst-MFA的NLP公式，合成测量值精确拟合到模型中。此外，由 5 生成的连续近似千-订单 Radau IIA 方法准确匹配通过求解 MATLAB 中的积木天平生成的合成数据。

thumbnail 下载：

个人电脑幻灯片

巴新放大图片

断续器原始图像

图 4.

（顶部）通过拟合综合测量值生成的DLS的连续近似值。标记指示 5 的并置点千-排序 Radau IIA 方法，较大的标记指示每个时间间隔的终点（也是 Radau IIA 方法中的搭配点）。垂直虚线表示每个时间间隔的边界。（底部）图中显示了与合成数据的拟合，每个合成测量值的标准差任意分配 0.01。合成数据是通过使用4以数值方式求解具有已知通量和池大小分布的积木平衡而生成的。千-订单显式 Runge-Kutta 方法和时间步长 0.02。由 5 生成的连续近似千-订单 Radau IIA 搭配方法拟合合成测量值叠加以进行比较。显然，使用搭配方法，合成测量值被模型精确拟合，连续近似值与合成数据准确匹配。

https://doi.org/10.1371/journal.pcbi.1009831.g004

内部MFA问题由一个ODE系统组成，其中许多是僵硬的。对于典型的标记实验，问题的性质知道刚性区域位于时间跨度的开始。在实验时间跨度的早期分配较短的时间间隔非常重要，这样可以准确地捕获刚性动力学。如果指定包含 ODE 刚度区域的时间间隔太宽，则用户可能会注意到此时间间隔内多项式近似中的数值伪影。然而，由于Radau IIA方法的刚性稳定性，数值解在随后的间隔内保持稳定。

为了证明Radau IIA方法的重要稳定性，我们接下来在动力学的刚性区域中使用大于允许的步长来求解玩具网络（图5）。在本例中，我们将时域划分为两个连续的区间：[0， 10] 和 [10， 20]，并使用 5 执行参数估计千-订购方法。正如预期的那样，我们发现该解决方案在[0，10]区间内近似较差。该区域的高测量残差表明了这一点。然而，后续区间[10，20]中的解没有受到前一个区间中近似值差的显着影响，如该区间中的小测量残差所示。此外，尽管在第一个时间间隔内近似值较差，但解的数值近似逐渐接近系统的真实稳态解。在这种情况下，稳态解决方案是所有代谢物都完全富集，这与数值近似一致。这种数值稳定性是适用于刚性ODE系统的方法的一个重要特征，例如Radau IIA。

在选择时间点将时域划分为间隔时，用户必须注意典型同位素标记实验的动态。在典型的实验中，上游代谢物的标记状态在将同位素标记的化合物引入细胞培养物后立即迅速变化。在这种初始的快速瞬变之后，随着系统向稳态标记分布松弛，动态变慢。其动力学在两个（或多个）不同时间尺度下运行的系统通常被称为僵硬。同位素标记状态经常表现出僵硬的动力学，动力学的刚性，快速变化的区域发生在开始时。因此，用户应分配时间点，以便在早期时间间隔较短，以便在快速变化的区域中准确捕获ODE解决方案。以后可能会分配更宽的时间间隔，因为解决方案预计会缓慢变化。

我们发现，估计 inst-MFA 模型的通量的最有效方法是使用 3 个初始可行通量和池大小分布（有关 eiFlux 使用的方法，请参阅 S5 文本）最初执行一组优化重启。第三-订购拉道 IIA 搭配方法。随后应将这些重新启动的最佳解决方案用作使用 9 进行优化的初始点千-订购 Radau IIA 方法以提高解决方案的准确性。eiFlux允许用户指定此细化，以便它自动发生。

在选择时间点和实现搭配方法时，我们建议遵循以下准则：

应分配时间点，以便在早期时间产生较短的时间间隔，以准确捕获该区域中快速变化的动态，并且在解决方案缓慢变化并且分离为太多时间间隔时，在以后的时间间隔更长的时间间隔将在计算上浪费。

为了最大限度地提高计算效率，应使用三阶 Radau IIA 并置方法从随机初始参数值多次重新启动优化，以搜索全局最优值。这些重启的最佳解决方案是作为使用高阶方法进行优化的起点，例如 5 阶或 9 阶 Radau IIA 方法（首选 9 阶）。

如果用户观察到在绘制时解的近似值很差，如时间间隔之间导数的明显不连续性所表明的那样，则应在近似值较差的时间间隔内添加一个额外的时间点，以将此间隔一分为二并获得更好的近似值。

thumbnail 下载：

个人电脑幻灯片

巴新放大图片

断续器原始图像

图 5. 在刚度区域中为玩具网络分配大于允许步长的结果。

当时域分为两个区间时，eiFlux 收敛的总测量残差：绘制 [0， 10] 和 [10， 20] 图，表明第一个时间区间的近似值较差不会显著影响后一个时间区间的解。为 B 的所有测量片段绘制解。显然，该解在第一个时间间隔内近似较差，但第二个区间中的解保持精确近似，并且渐近接近真正的稳态解。

https://doi.org/10.1371/journal.pcbi.1009831.g005

大肠杆菌串联在MFA示例

我们在一个现实的例子上演示了eiFlux，即E。大肠杆菌中枢碳代谢。该网络基于Young等人提出的网络[9]。它总共含有58种具有5152种同位素（5094种含糖体）的代谢物和具有26种独立通量的84种反应。对于给定的通量分布和随机分配的代谢物池大小分布，我们模拟了来自假设实验的瞬态串联MS测量，其中E.大肠杆菌均匀喂养99%13C标记的葡萄糖。模拟的串联MS数据包括Rühl等人描述的代谢物和子片段[32]。在这样的实验中，稳态同位素标记分布没有提供有关代谢途径通量的信息。因此，必须从瞬态（同位素非平稳）标记数据中提取此信息。总共有2338个测量值适合该模型。基于 95% 置信区间的可接受 SSR 为 2368。为了适应这种数据类型，我们使用cumomer建模形式主义来构建NLP。

我们使用随机正态分布误差对模拟数据进行扰动以模拟实验误差，然后使用eiFlux将此数据拟合到模型中以估计通量和池大小。由于模拟数据集由完整的 MID 组成，因此未使用测量比例因子。为了寻找最佳拟合解决方案，使用 3 从随机初始可行通量和池大小分布重新启动优化 15 次第三-订购 Radau IIA 搭配方法来处理 ODE。这些起点与最优值相去甚远，每个起点都有一个大约为 10 的初始目标函数值6.所有重新启动都收敛到一个狭窄的 SSR 范围 — 2294（3 次）、2293（2 次）和 2288（10 次）。从非常糟糕的起点紧密收敛的优化证明了并置NLP方法对inst-MFA的鲁棒收敛特性，并展示了与使用拍摄方法的传统方法相比的明显优势。这 15 次重新启动的最佳解决方案被用作使用 9 次优化的初始点千-订购Radau IIA方法，提高溶液的精度，实现了2212的SSR。在单个处理器内核上运行的英特尔酷睿 i7-7700 CPU @ 3.60 GHz 上，这 15 次重新启动总共花费了 144.3 分钟，即每次重启平均需要 9.62 分钟（中位数：8.78 分钟）。使用 9 的细化步骤千-订购方法花了5.80分钟。

理想情况下，我们将恢复生成模拟数据的通量分布。然而，由于数据受到人为不确定性的干扰，原始的通量分布可能无法准确恢复。有利地，通过发酵醋酸酯合成的助焊剂收敛到接近其实际值，尽管不是测量的助焊剂。这表明，使用同位素非平稳串联MS数据可以潜在地估计一些代谢物流入和流出，而无需直接测量它们。通量图和模型拟合到一些代表性数据如图6所示。为了估计通量置信区间，我们使用自举蒙特卡罗[33]。这种方法以前曾用于MFA研究;有关示例和简要说明，请参见 [34]。实际磁通量值的子集与估计通量值之间的比较如表2所示。显然，eiFlux通过这些重要反应准确地估计了通量，尽管测量受到模拟测量误差的扰动。

thumbnail 下载：

个人电脑幻灯片

巴新放大图片

断续器原始图像

图 6.

（左）使用埃舍尔生成的大肠杆菌模型的中心碳代谢反应。未显示产生生物质代谢物（氨基酸）的外周反应。有关网络反应及其相关碳原子重排的完整列表，请参见S4文本。（右）一组代表性模型拟合到几种代谢物的模拟串联MS测量数据。有关适合所有模拟数据的完整模型集，请参阅 S4 文本。

https://doi.org/10.1371/journal.pcbi.1009831.g006

thumbnail 下载：

个人电脑幻灯片

巴新放大图片

断续器原始图像

表 2. 将 E 反应子集的 True 和估计通量值进行比较。大肠杆菌中心碳代谢网络。

计算真实通量和估计通量之间的差值。显然，eiFlux准确地估计了受模拟测量误差扰动的一组模拟测量的通量。缩写：磷酸戊糖途径（PPP），Entner-Doudoroff途径（EDP），三羧酸（TCA），葡萄糖6-磷酸（G6P），果糖6-磷酸（F6P），甘油3-磷酸（GAP），1，3-双磷酸甘油酯（BPG），6-磷酸葡萄糖酸盐（6PG），2-酮-3-脱氧-6-磷酸葡萄糖酸盐（KDPG），草酰乙酸酯（OAA），乙酰辅酶A（AcCoA），柠檬酸盐（Cit），异柠檬酸盐（iCit），乙醛酸盐（Gox），琥珀酸盐（Suc），苹果酸盐（Mal），丙酮酸盐（Pyr），磷酸烯醇丙酮酸盐（PEP），乙酸酯（AcCoA）。

https://doi.org/10.1371/journal.pcbi.1009831.t002

运行eiFlux后，会自动生成一个包含网络反应的JSON文件。此JSON文件可由Escher读取，Escher是用于可视化代谢网络和通量分布的应用程序[35]。这使用户能够可视化其网络模型并快速轻松地评估通量。这种快速可视化可以帮助模型开发和分析。E. 的通量图。图6中的大肠杆菌网络是使用埃舍尔生成的。eiFlux 还会自动绘制模型拟合到数据以进行视觉评估，并有一些代表性的模型拟合到 E。大肠杆菌串联MS示例数据如图6所示。完整的模型拟合集显示在 S4 文本中。S1 数据包含拟合到此模型的模拟数据。

与现有软件的比较

为了直接比较NLP方法与拍摄方法，我们使用eiFlux和MFA软件INCA（版本2.0）[24]（使用拍摄方法）求解了相同的模型。我们解决了E。之前使用合成的单MS时间序列数据集（S1数据）描述的大肠杆菌模型。为了开发这个合成数据集，我们使用与生成模拟串联MS数据时相同的通量和池大小分布，为[36]中列出的GC-MS片段集生成了模拟测量值。该数据是为喂养80%1，2-的假设实验而生成的13C和20%均匀13C标记的葡萄糖。然后，我们用模拟的测量误差扰动这些测量。为了能够与INCA直接比较，本例中在eiFlux中使用了测量比例因子。我们使用每个软件从随机初始可行通量和池大小分布中重新启动优化 20 次。这些初始通量和池大小分布是在INCA中生成的，方法是将给定重新启动的"扰动初始值的日志数"选项设置为4或5，从而确保起点远离最佳值。这些值被转移到eiFlux中，作为优化的起点。这保证了每次重启时，INCA和eiFlux中初始可行通量和池大小分布相同，从而确保了公平的比较。使用eiFlux，这20次重启都是使用3次执行的第三-订购 Radau IIA 方法，然后使用 9 进行改进千-订购拉道IIA方法，以提高解决方案的准确性。为了比较这两种方法，我们报告了每次重启的平均收敛时间，以及每个平台在最佳计算的最佳值的一个小邻域内收敛的重启比例。在此示例中，eiFlux 和 INCA 均在个人计算机上运行，使用在单个处理器内核上运行的 Intel Core i5-8250U CPU @1.60 GHz 处理器。

使用 eiFlux，使用 3 次重新启动的 20 次重新启动中计算出的最佳解决方案第三-订单Radau IIA方法为660.39，使用9进行细化后降至644.14千-订购方法。在这 20 次重启中，有 16 次收敛在最佳计算最佳值的 2% 以内，定义为在 660.39 – 673.60 范围内与 SSR 收敛的重新启动，使用 3第三-订购拉道IIA方法。其余 4 次重新启动收敛到比此最佳值高出 10% 的较差解决方案。使用 3 进行优化第三-订购 Radau IIA 方法需要每次重新启动的平均时间为 2.55 分钟。9号酒店千-订单细化平均需要2.65分钟才能收敛。使用INCA，从20次重新启动中计算出的最佳解为646.34。在这 20 次重新启动中，6 次收敛在最佳计算最佳值的 2% 以内（SSR 在 646.34 – 659.27 范围内）。总共有 13 个收敛到一个比此最佳值高出 10% 以上的较差解。这些结果总结在图 7 中。使用INCA的优化每次重启的平均时间为7.57分钟。这种比较表明，对于核心网络，NLP方法的运行时间与拍摄方法相当，但它为最佳解决方案提供了更强大的收敛。

thumbnail 下载：

个人电脑幻灯片

巴新放大图片

断续器原始图像

图 7.

（左）20 次重启中的 20 次重启，收敛于 NLP 方法和拍摄方法的最佳计算最佳值附近 （百分比）。对于NLP方法，这20次重新启动是使用3阶Radau IIA方法执行的，然后使用9阶方法进行细化。最佳最优成绩为3阶法660.39，9阶细化后为644.14。对于拍摄方法，我们使用软件INCA，最佳达到最佳值为646.34。（右）基因组规模的Synechocystis示例的100次重新启动的结果。使用NLP方法，大多数优化从随机初始点收敛到最佳计算最优解的10%以内，很大一部分在2%以内。这凸显了 NLP 方法对内部 MFA 应用程序的鲁棒性。

https://doi.org/10.1371/journal.pcbi.1009831.g007

基因组规模的胞囊肿 inst-MFA 示例

为了证明这种方法的可扩展性，我们解决了Gopalakrishnan等人提出的Synechocystis PCC 6803基因组规模的inst-MFA模型[37]。该模型适合Young等人提出的数据集[38]。在Gopalakrishnan等人的分析中[37]，作者使用拍摄方法进行了inst-MFA，以估计Synechocystis PCC 6803（以下简称Synechocystis）基因组规模模型中的通量和池大小。在这种方法中，作者首先执行FVA以消除模型中在研究的光自养生长条件下不能携带通量的反应，这导致模型尺寸减小。然而，该模型仍然相当大，需要8.4×10的仿真。5Levenberg-Marquardt 优化算法每次迭代时的 ODE。这些方程的数量包括模拟动车组动力学所需的状态方程及其对拟合参数的灵敏度。

使用此方法时，必须在优化算法的每次迭代中从头开始求解状态和灵敏度方程。在Levenberg-Marquardt算法中，灵敏度方程用于在每次迭代期间确定良好的步进方向以进行参数估计。在 inst-MFA 情况下，状态和灵敏度方程都是 ODE。灵敏度方程的数量缩放为n×p，其中n是EMU质量分数的数量，p是自由参数的数量[23]。这种扩展意味着随着网络规模的增加，灵敏度方程的数量会变得非常大。作者报告了851个动车组，有2311个动车组质量分数，每个都贡献了一个平衡。有 367 个拟合参数。灵敏度方程的数量是通过将这两个值相乘来计算的，得到报告的 8.4×105在Levenberg-Marquardt算法的每次迭代中必须从头开始求解的ODE。这就是拍摄方法的可扩展性差的原因。

虽然作者没有报告求解该模型所需的计算时间和内存，但考虑到ODE系统的大小，它可能相当大。作者发现，核心模型和基因组尺度模型之间的许多预测通量差异很大。这表明，在开发核心模型时所做的建模假设（例如通路集总）可能会明显偏倚通量估计值。这些结果表明，使用inst-MFA时，更大，更详细的模型更适合精确阐明通量。这突出表明，如果大规模模型的此类分析要成为常规，则需要对MFA进行可扩展的方法。

我们使用 eiFlux 实现 EMU 建模方法，使用基因组规模的 Synechocystis 模型重复了通量和池大小估计。由于拟合的 MID 不完整且未校正，因此比例因子在优化中作为自由参数处于活动状态。通过将额外原子的效应分解到测量映射矩阵M中来解释附加原子的天然同位素丰度新泽西州，在 NLP 中。这使我们能够适应本研究中提供的数据。我们使用 3 从 100 个随机初始可行通量和池大小分布重新开始优化第三-订购拉道IIA方法。每次重启的平均收敛时间为15.89分钟。这种基因组规模模型的短收敛时间表明，NLP方法可以很好地扩展到大型，详细的网络，这些网络代表了inst-MFA应用的当前前沿。在这 100 次重启中，94 次收敛在 10% 以内，19 次收敛在最佳计算最佳 （568.23） 的 2% 以内，当运行 3第三-订购拉道IIA方法。这表明优化并不经常收敛到远离全局最优解的较差局部最优，并且即使对于基因组规模模型，这种方法仍然稳健有效。这些结果总结在图 7 中。为了提高解决方案的准确性，这 100 次重新启动中的最佳解决方案被用作使用 9 进行优化的初始点千-订购拉道IIA方法。这 9千-阶细化收敛到最佳 SSR 为 513.30 的解决方案。

对于使用 3 的优化第三-Order Radau IIA方法，NLP在CONOPT的模型预处理步骤之后包括约80，000个代数等式约束和变量。对于 9千-序法，NLP包括~182，000个代数相等约束和变量。与8.4×10相比，这是对问题的显着简化。5在拍摄方法方法中，在每次迭代Levenberg-Marquardt算法时必须从头开始求解的ODE。这种简化的结果是，灵敏度方程不需要在NLP方法中构建和求解。

该优化还具有内存效率，即使对于9千-订单细化步骤。CONOPT 所需的大部分内存用于存储约束的雅可比矩阵及其 LU 分解。虽然像这样的大型模型可能具有10的数量级5相等约束和变量，雅可比矩阵非常稀疏。当使用 9 运行细化时千-订购拉道IIA方法，雅可比人只有9.7×105CONOPT 的模型预处理步骤之后的非零元素。GAMS 和 CONOPT 使用稀疏线性代数例程来加速计算并保持内存效率。值得注意的是，二阶导数信息仅以简化的黑森变量的形式计算非基本变量。这种减少的黑森州很小，与雅可比及其LU分解相比，对CONOPT算法的内存使用的贡献可以忽略不计。此示例表明，可以使用单个处理器内核在标准个人计算机上有效地求解大规模模型。

讨论

MFA，特别是inst-MFA，是一个计算密集型参数估计问题，涉及包含O（1000）ODE和O（100-1000）测量值的模型中的O（10-100）自由参数（通量和代谢物池大小）。目前为inst-MFA开发的可用算法使用了顺序优化算法，其中DLS在每次优化迭代时都显式地从头开始重新计算。这种拍摄方法可能效率低下，并且可能无法很好地扩展到大型网络。众所周知，它遭受小区域收敛[13]，如果没有良好的初始猜测，通常无法稳健地收敛到全局最优。在这份手稿中，我们遵循了一种不同的方法，并将MFA问题表述为一个相等约束的非线性程序（NLP）。在这种方法中，DLS 不是在每次迭代时从头开始求解，而是更新以满足相等约束。换句话说，DLS 是隐式处理的，而不区分状态变量和模型参数。这种方法显著提高了优化的效率，特别是当DLS由O（100-1000）标记瞬变组成时。当 inst-MFA 模型包含大量状态变量和参数时，它也可以很好地缩放。NLP方法使cumomer或EMU框架能够得到有效和高效的使用，从而实现强大的收敛和快速的运行时间。

我们通过使用Wiechert等人提出的框架将MFA问题表述为平等约束的NLP[1]。这种NLP可以使用在最先进的AML（如GAMS）中实现的优化算法来解决。为了解决 inst-MFA 参数估计问题，我们使用并置将 ODE 系统转录为 NLP 的一组代数约束。我们发现，在GAMS中实现的大规模降低梯度算法CONOPT在解决稳态和内部MFA问题方面特别快速和健壮。但是，原则上，可以使用任何通用的NLP求解器。对于inst-MFA问题，与使用拍摄方法的传统方法相比，这种方法可以很好地扩展到具有许多自由参数的大型网络，并且具有更强大的收敛特性。此外，由于 MFA 问题可以使用 cumomer 框架来表述，因此对于稳态和 inst-MFA 都可以轻松支持串联 MS 数据。我们组装了一个软件包eiFlux，它使用Python和GAMS构建，它采用NLP公式，并使用cumomer或EMU框架支持稳态和inst-MFA。

在用玩具网络演示概念验证后，我们在逼真的E上使用了我们的配方和eiFlux。大肠杆菌串联MS示例。在这里，我们证明了我们的NLP公式在从较差和随机的初始点执行的每次重启时都会稳健地收敛到狭窄的SSR范围。计算速度很快，使用cumomer框架进行inst-MFA平均收敛到最佳值平均不到10分钟。在此示例中，在 CONOPT 的模型预处理步骤之后，基础 inst-MFA NLP（包括转录的 ODE 系统）中的约束和变量数量分别约为 180，000。

接下来，我们提供了使用NLP公式的eiFlux和使用核心E上拍摄方法的INCA之间的直接比较。大肠杆菌网络拟合模拟单MS数据。我们发现两个平台的计算速度相当，但eiFlux收敛到最优时更稳健，这可以从随机初始猜测的重启分数中得到证明，这些猜测收敛在最佳计算最优的一个小社区内。最后，我们使用eiFlux求解了Gopalakrishnan等人提出的Synechocystis代谢的大型基因组规模inst-MFA模型[37]由Young等人提出的拟合时间序列测量[38]。这个例子表明，NLP公式可以很好地扩展到非常大的网络，并且使用这种公式，可以在个人计算机上有效地求解这样的网络。这种方法的可扩展性使大规模 inst-MFA 模型的常规解决方案成为可能，并显著扩展了将来可能应用 inst-MFA 的范围。我们预计这种方法将允许复杂的真核网络（例如藻类）和共培养代谢的大规模模型通过inst-MFA有效地处理。

过去，提高稳态和inst-MFA中参数估计的速度和效率，强调在优化算法的每次迭代中快速求解标记状态。在这里，我们专注于在改进的优化框架中制定问题，该框架允许将MFA问题有效地扩展到更大，更详细的网络。我们发现这种方法对于inst-MFA特别有用，可以有效地解决基因组规模模型。最后，NLP 方法的灵活性允许在磁通量和池大小上加入额外的线性和非线性约束，这些约束不容易纳入当前的 MFA 或 inst-MFA 应用程序中。这些类型的约束可以纳入 MFA 的未来扩展和使用，并代表这项工作的未来可能方向。

支持信息

Example of a reduced gradient nonlinear programming algorithm.

Showing 1/9: pcbi.1009831.s001.pdf

Skip to figshare navigation

1简化梯度非线性规划示例1算法2普遍平等-约束非线性规划问题如所示（S1.1).这里()fx是3目标函数最小化，这可能是非线性的。向量函数()gx是4相等约束，通常可能是非线性的。变量，x，可以添加在限定下5下限和上限，磅和ub分别。6最小化s.t.()()fg磅ubxx0x（S1.1)7拉格朗日王朝（S1.2)可由模型组件引入（S1.1).8T( , )( )( )Lfg三十（S1.2)9在等式中（S1.2),( , )Lx是拉格朗日量，并且是对应于10平等约束，()gx.每()我gx具有相关的拉格朗日乘数我.优化平等-11受约束程序英沃尔韦斯查找拉格朗日量的静止点（S1.3).静止点是12被标识为其中( , )Lx0.13T( , )( )( )Lfgxxx0（S1.3)14在（S1.3),( , )Lx是拉格朗日量相对于模型变量的梯度，x.此渐变15在最大值、最小值或鞍点处可以为零的拉格朗日量。术语()gx是一个矩阵16其行由约束的渐变组成（S1.4)，也被称为Ja科比安的17约束()Jx.18TTT12()()( )( )()mgg吉杰gxx二十x（S1.4)19用符号替换（S1.4),（S1.3)成为（S1.5).20T( , )( ) ( )LfJxxx0（S1.5)2122中所示的表单（S1.5)用于构造简化梯度算法。23广义降低梯度 （GRG） 算法适用于非线性选择中的 imization 问题24其中目标函数和约束是平滑和可微的[1,2].特别好-25适用于等式约束总数与总约束数相似的优化问题26模型变量的数量。这是因为 GRG 算法涉及将变量划分为27集。一套组成s 的基本变量，bx，以及另一个非基本变量，铌x.在 每次迭代时28GRG算法，非基本变量被传播到一个更好的点（一个改善29目标函数），a基本变量被传播到新点，使得它们满足30平等约束[1,3].以下通用 GRG 算法中的步骤列表改编自 Drud31[3]:321.找到一个可行的起点，0x.此时，计算目标函数0()fx，的渐变33目标函数0()fx和雅可比安约束0()Jx.342.选择一组n基本变量，bx，使得bJ，基本列的子矩阵J是35非奇异。36

23.求解乘数T断续器杰夫.374.计算降低的梯度T断续器r.385.如果停止标准已满足。396.查找搜索方向，d，对于基于r并计算切线方向。407.沿方向执行搜索d.对于每个步骤，调整bx满足( , )b铌gxx0使用伪-41牛顿过程。428.重复这些ps 与当前点。43接下来，将在一个简单的示例中说明此过程。此示例 GRG 算法仅用于44说明目的。大-scale GRG 算法，如 CONOPT[2,3],要复杂得多，45使它们更加有效和高效。因此，以下示例仅供说明之用，并且是46不代表CONOPT的实施。47平等的例子-约束非线性程序显示在（S1.6).48最小化s.t.22121 212221131221 2143521( ) 200360240168( )0.5 0( )2 秒( )00fxxx x xxxgxxxxgx x xxxgxxxxxxxx（S1.6)49这个 exa之所以选择 mple，是因为它由具有非线性约束的非线性目标函数组成。50具体说来3()gx是线性的，而1()gx和2()gx都是非线性的。另请注意，（S1.6)在一般情况下51格式显示于（S1.1).在这里，目标函数只是两个模型变量的函数，1x和522x.这种结构的问题至少很常见-方块模型拟合。这个问题有一个总数53的 5 个变量和 3 个相等约束。54操作 GRG 算法需要重复使用目标函数的梯度和55约束的雅可比。此处列出的示例问题供参考：56122113214002003604002002402 1 1 2 0 00( ), ( )10 1 001 1 0 0 10二十二十二十fJxx二十（S1.7)5758步骤 1：使用可行的起点进行初始化59从满足的点开始s 约束()gx0.计算目标函数，梯度60目标函数，以及此时的雅可比函数。我们将从任意选择的可行点开始61T00 0.5 0 1 0.5x.选择这个起点仅仅是因为它满足了约束。这62目标函数0()fx、目标函数梯度0()fx和雅可比语00()婉婷x在这一点上是63显示于（S1.8).下标，0，是当前迭代的索引。64000260401 1 0 0 00( ) 98.0,( ),1.5 0 0 1 001 1 0 0 10ffJ二十（S1.8)6566

1 / 9

Download

figshare

S1 文本。 简化梯度非线性规划算法的示例。

https://doi.org/10.1371/journal.pcbi.1009831.s001

（英文）

S2 文本。 CONOPT 处理相等约束的数学细节。

https://doi.org/10.1371/journal.pcbi.1009831.s002

（英文）

S3 文本。 简单的同位素网络来说明搭配。

https://doi.org/10.1371/journal.pcbi.1009831.s003

（英文）

S4 文本。 玩具网和E.大肠杆菌串联质谱法以MFA为例。

https://doi.org/10.1371/journal.pcbi.1009831.s004

（英文）

S5 文本。 在 eiFlux 中生成初始可行通量和池大小分布的内部方法。

https://doi.org/10.1371/journal.pcbi.1009831.s005

（英文）

S6 文本。 符号列表。

https://doi.org/10.1371/journal.pcbi.1009831.s006

（英文）

S1 图 通过串联MS测量的同位素信息。

https://doi.org/10.1371/journal.pcbi.1009831.s007

（英文）

S1 数据。 用于玩具网络和E的模拟数据。大肠杆菌的例子。

https://doi.org/10.1371/journal.pcbi.1009831.s008

（XLSX）

S1 文件。 包含以下内容的 Zip 文件：本文中提供的所有示例的 eiFlux 模型。

这些文件夹中的文本文件用作 eiFlux 的输入。文本文件 modelfile.txt 用于指定要运行的模型。eiFlux 使用的 GAMS 脚本在文件夹中GAMS_Model。eiFlux的演示版本专门运行上面列出的两个模型，从而完全复制了本手稿中的结果。这是一个编译的Python文件，eiFlux_Limited.pyc。此演示版本仅运行上面列出的两个模型。用户可以更改模型参数，例如时间节点，但不能更改模型本身。按照Installation_Instructions.pdf文件中给出的说明进行操作后，用户可以运行eiFlux_Limited.pyc文件来运行该软件。

https://doi.org/10.1371/journal.pcbi.1009831.s009

（邮编）

确认

作者感谢丹麦ARKI Consulting & Development A/S的Arne Stolbjerg Drud博士分享了对CONOPT算法内部工作的解释，并感谢马里兰大学的Stephen （Will） Carroll对eiFlux初始开发的贡献。

引用

1.Wiechert W， M?llney M， Isermann N， Wurzel M， de Graaf AA.代谢网络中的双向反应步骤：III。同位素标记系统的显式解决方案和分析。生物技术生物。1999;66: 69–85.pmid：10567066

查看文章PubMed/NCBI谷歌学术搜索

2.安东尼耶维奇先生代谢工程中代谢通量分析指南：方法、工具和应用。Metab Eng. 2021;63： 2–12.下午：33157225

查看文章PubMed/NCBI谷歌学术搜索

3.谢烨，年轻的京东。同位素非平稳代谢通量分析（INST-MFA）：将理论付诸实践。库尔·奥宾生物技术。2018;54: 80–87.下午：29522915

查看文章PubMed/NCBI谷歌学术搜索

4.Wiechert W， N?h K. 同位素非平稳代谢通量分析：复杂但信息量大。库尔·奥宾生物技术。2013;24: 979–986.pmid：23623747

查看文章PubMed/NCBI谷歌学术搜索

5.Antoniewicz MR，Kelleher JK，Stephanopoulos G.基本代谢物单位（EMU）：用于模拟同位素分布的新框架。Metab Eng. 2007;9： 68–86.pmid：17088092

查看文章PubMed/NCBI谷歌学术搜索

6.Sriram G， Shanks JV.使用碳键标记实验改进代谢通量分析：键聚体平衡和布尔函数映射。Metab Eng. 2004;6： 116–132.下午：15113565

查看文章PubMed/NCBI谷歌学术搜索

7.van Winden WA， Heijnen JJ， Verheijen PJT.累积键聚体：二维通量分析的新概念 [13C,1H] COSY NMR 数据。生物技术生物。2002;80: 731–745.pmid：12402319

查看文章PubMed/NCBI谷歌学术搜索

8.Srour O， Young JD， Eldar YC.通量单体：一种新方法13C代谢通量分析。BMC 系统生物学 2011;5： 1–14.下午：21194489

查看文章PubMed/NCBI谷歌学术搜索

9.年轻的JD，Walther JL，Antoniewicz MR，Yoo H，Stephanopoulos G.一种基于基本代谢物单元（EMU）的同位素非平稳通量分析方法。生物技术生物。2008;99: 686–699.pmid：17787013

查看文章PubMed/NCBI谷歌学术搜索

10.Shin S， Venturelli OS， Zavala VM.用于动态生物系统模型中参数估计的可扩展非线性规划框架。PLOS Comput Biol. 2019;15： e1006828.pmid：30908479

查看文章PubMed/NCBI谷歌学术搜索

11.柯克帕特里克S，Gelatt CD，Vecchi MP。通过模拟退火进行优化。科学。1983;220: 671–680.pmid：17813860

查看文章PubMed/NCBI谷歌学术搜索

12.罗德里格斯-费尔南德斯 M， 门德斯 P， 小班加一种混合方法，用于在生化途径中进行高效而稳健的参数估计。生物系统。2006;83: 248–265.pmid：16236429

查看文章PubMed/NCBI谷歌学术搜索

13.Peifer M， Timmer J. 使用多次拍摄方法对生化过程的常微分方程中的参数估计。IET 系统生物学 2007;1： 78–88.pmid：17441551

查看文章PubMed/NCBI谷歌学术搜索

14.Betts JT. 使用非线性规划的最优控制和估计的实用方法，第2版。工业与应用数学学会;2010.

15.Drud A. CONOPT：用于大型稀疏动态非线性优化问题的GRG代码。数学课程。1985;31: 153–191.

查看文章谷歌学术搜索

16.Drud A. CONOPT - 一个大规模的GRG代码。ORSA J Comput.1994;6: 207–216.

查看文章谷歌学术搜索

17.通用代数建模系统 （GAMS）。2751 Prosperity Ave， Suite 210， Fairfax VA 22031： GAMS Development Corporation.弗吉尼亚州费尔法克斯， 美国;2021.

18.Choi J，Antoniewicz MR.串联质谱：代谢通量分析的新方法。Metab Eng. 2011;13： 225–233.pmid：21134484

查看文章PubMed/NCBI谷歌学术搜索

19.M?llney M， Wiechert W， Kownatzki D， de Graaf AA.代谢网络中的双向反应步骤：IV。同位素标记实验的优化设计。生物技术生物。1999;66: 86–103.下午：10567067

查看文章PubMed/NCBI谷歌学术搜索

1000万Abadie J， Carpentier J. Optimization.学术出版社;1969.

21.N?h K， Wahl A， Wiechert W. 同位素稳态计算工具13代谢稳态条件下的C标记实验。Metab Eng. 2006;8： 554–577.pmid：16890470

查看文章PubMed/NCBI谷歌学术搜索

22.维彻特 W， 诺赫 K.从静止到静止代谢通量分析。Adv Biochem Eng Biotechnol.2005;92: 145–172.pmid：15791936

查看文章PubMed/NCBI谷歌学术搜索

23.Mottelet S，Gaullier G，Sadaka G.使用伴随方法在同位素标记实验中的代谢通量分析。IEEE/ACM Trans Comput Biol Bioinform.2017;14: 491–497.下午：28113867

查看文章PubMed/NCBI谷歌学术搜索

24.年轻的京东。INCA：用于同位素非平稳代谢通量分析的计算平台。生物信息学。2014;30: 1333–1335.pmid：24413674

查看文章PubMed/NCBI谷歌学术搜索

25.Kajihata S， Furusawa C， Matsuda F， Shimizu H. OpenMebius： An Open Source Software for Isotopically Non stationary13基于C的代谢通量分析。BioMed Res Int. 2014;2014： e627014.pmid：25006579

查看文章PubMed/NCBI谷歌学术搜索

26.Fr?hlich F， Kaltenbacher B， Theis FJ， Hasenauer J. 基因组规模生化反应网络的可扩展参数估计。PLOS Comput Biol. 2017;13： e1005331.pmid：28114351

查看文章PubMed/NCBI谷歌学术搜索

27.Tjoa IB， Biegler LT. 微分代数方程组参数估计的同步解和优化策略.工业化学研究. 1991;30： 376–385.

查看文章谷歌学术搜索

28.Hairer E， Wanner G. 由 Radau 方法求解的刚性微分方程。J Comput Appl Math. 1999;111： 93–111.

查看文章谷歌学术搜索

29.Huynh HT. 搭配和 Galerkin 时间步进方法。第19届AIAA计算流体动力学。德克萨斯州圣安东尼奥：美国航空航天研究所;2009. https://doi.org/10/gg4rxz

30.Antoniewicz MR. 串联质谱法，用于测量稳定同位素标记。库尔·奥宾生物技术。2013;24: 48–53.pmid：23142542

查看文章PubMed/NCBI谷歌学术搜索

31.Choi J， Grossbach MT， Antoniewicz MR. 使用气相色谱/串联质谱法测量天冬氨酸的完整同位素分布。肛门化学. 2012;84： 4628–4632.pmid：22510303

查看文章PubMed/NCBI谷歌学术搜索

32.Rühl M， Rupp B， N?h K， Wiechert W， Sauer U， Zamboni N. LC-MS/MS 中中心碳代谢物的碰撞片段化提高了13C代谢通量分析。生物技术生物。2012;109: 763–771.pmid：22012626

查看文章PubMed/NCBI谷歌学术搜索

33.Press WH， Teukolsky SA， Vetterling WT， Flannery BP. Numeric Recipes 3rd Edition： The Art of Scientific Computing.第3版. 剑桥大学出版社;2007.

34.Nargund S，Misra A，Zhang X，D. Coleman G，Sriram G.通量和反流：植物悬浮细胞中的代谢物反流及其对同位素辅助代谢通量分析的意义。摩尔生物系统。2014;10: 1496–1508.下午：24675729

查看文章PubMed/NCBI谷歌学术搜索

35.国王ZA，Dr?ger A，Ebrahim A，Sonnenschein N，Lewis NE，Palsson BO。埃舍尔：一个Web应用程序，用于构建，共享和嵌入生物途径的数据丰富的可视化。PLOS Comput Biol. 2015;11： e1004321.pmid：26313928

查看文章PubMed/NCBI谷歌学术搜索

36.年轻的JD，Allen DK，Morgan JA。代谢通量分析中的同位素物测量技术II：质谱。植物代谢：方法和协议。纽约：施普林格科学+商业媒体有限责任公司;2014.

37.Gopalakrishnan S， Pakrasi HB， Maranas CD.使用基因组尺度碳图谱模型阐明了Synechocystis PCC 6803中的光自养碳通量拓扑。Metab Eng. 2018;47： 190–199.下午：29526818

查看文章PubMed/NCBI谷歌学术搜索

38.年轻的JD，Shastri AA，Stephanopoulos G，Morgan JA。使用同位素非平稳13C通量分析绘制光自养代谢图。Metab Eng. 2011;13： 656–665.pmid：21907300

查看文章PubMed/NCBI谷歌学术搜索