非线性动力学模型广泛用于描述控制复杂生物途径系统的过程。在过去十年中，由于使用分子生物学方法通过高通量实验收集的数据，这些模型的验证和进一步开发成为可能。虽然这些数据非常有益，但它们通常是不完整和有噪声的，这使得复杂动态模型的参数值推断具有挑战性。幸运的是，许多生物系统都嵌入了可以利用的线性数学特征，从而改进拟合并导致优化算法更好地收敛。在本文中，我们使用一种新的方法来探索动态模型的推理选项可分离非线性最小二乘优化并将其性能与传统的非线性最小二乘法进行了比较。大量模拟的数值结果表明，所提出的方法至少与传统的非线性最小二乘法一样准确，但通常更优越，同时也大大减少了计算时间。
1.简介
非线性动力学模型广泛用于描述控制复杂生物途径系统的过程。在这种情况下，人们特别感兴趣的是所谓的规范格式，它们在可能的响应方面非常灵活，但涉及功能形式的非常有限的领域。在线性系统之外，最著名的规范格式是Lotka–Volterra（LV）模型[1–4]在生物化学系统理论（BST）的框架内使用二项式项和幂律系统，后者专门使用幂函数的乘积。BST最初设计用于分析生化和基因调节系统，但后来在各种生物医学和其他领域得到了更广泛的应用[5，6]. 虽然很容易以符号形式建立复杂生物系统的LV或BST模型，但确定最佳参数值仍然是一个重大挑战。因此，估计常微分方程组（ODE）的参数仍然是一个活跃的研究领域，吸引了各种科学领域的贡献（例如[7–12]. 事实上，近年来已经提出了许多用于ODE模型的优化方法，但没有一种在广泛的应用中特别有效，原因是生物数据的内在问题（稀疏性、不确定性、噪声等）跨越了整个频谱技术和计算问题（大量局部极小值、不可识别性、草率等）。基于斜率的估计等方法（例如[13])和动态通量估计[14–16]缓解这些问题，但不是万能药。
在这里，我们重新审视并实现了早期的想法[17]将估计任务分为线性和非线性两个方面。然而，我们的主要关注点是不真的是一种新的估计方法就其本身而言相反，我们感兴趣的是关于参数估计的更一般和更高层次的观点，而不是技术文章中通常提出的观点。具体来说，本文讨论了动态模型的参数估计，其数学格式包含线性特征，允许参数和系统状态的自然分离。一个简单的例子是线性常微分方程，其中向量场在参数中是线性的θ，使用表示…的导数关于t吨作为一个更有趣的例子，ODE向量场在参数中可能是部分线性的，就像BST中所谓的S系统模型一样[5].
示例1。S系统[18]定义为在这里为非负速率常数，而是反映变量对给定流入或流出影响的强度和方向性的实值动力学顺序。非正式地，人们可以将这个系统视为一个回归方程，其中“协变量”是变量在右手边，而“响应”变量是导数在左侧。注意，矢量场在速率常数中是线性的，但在动力学级是非线性的.
利用动态模型中参数和系统状态的可分性的估计方法历史悠久；参见[19]对于特殊情况。然而，直到最近才对这种方法进行了严格的统计分析[20]. 在一篇关于动力学模型推理的经典论文中，Varah[17]顺便提到“可以使用可分性或可变投影的概念（参见[21]或[22])在隐式求解线性参数的情况下，对非线性参数求解所产生的（完全）非线性最小二乘问题，然后利用它们对非线性参数的表示获得线性参数。由于这减少了要解决的非线性最小二乘问题的规模，因此值得一试。”令人惊讶的是，鉴于常微分方程的参数估计通常被认为是建模动态过程的计算瓶颈，瓦拉的建议在实践中没有得到广泛的遵循。事实上，在致力于动态模型参数拟合技术的大量文献中，我们只知道两个相关的参考文献：使用直接积分方法，Dattner等人[23]将可分离非线性最小二乘技术应用于非均匀环境中捕食者-食饵系统参数的推断，而Wu等人[24]利用可分性估计高维线性常微分方程系统的参数。此外，瓦拉利用可分性估计ODE参数的想法最近才在一个公开可用的软件包中实现[25]. 该软件的相关细节将在后面的章节中讨论。
本文的分析旨在说服读者，瓦拉的想法确实值得追求。为了支持这一观点，我们探索并比较了两种动态模型的通用数据拟合方法：传统的非线性最小二乘法（NLS）和提出的可分离非线性最小二乘法。通过对具有代表性的复杂模型进行广泛的Monte-Carlo模拟，我们识别并量化了通过这种分离方法获得的显著统计和计算收益。我们最终会得出这样的结论：模型可分性是非常有益的，对于任何具有显著线性特征的复杂动态系统，都应该考虑SLS方法。
论文组织如下。在节中2，我们将在动态模型的背景下详细介绍SLS方法。章节3描述了模拟设置，量化了我们用于比较SLS和NLS性能的统计度量，并给出了数值结果。在节中4，我们指出了未来的研究方向，而结论在第节中提供5.
2.可分离非线性最小二乘法（SLS）和瓦拉的思想
2.1. 概述
根据Varah在动态模型推理上下文中的原始想法，利用可分性进行参数估计的主要优点如下[26]:（i）优化所需的初始猜测更少（ii）优化问题的条件更好（iii）融合速度更快
这些优势在一些出版物中得到了令人信服的证明。例如，请参见马伦[27]物理和化学的实施和应用；Chung&Nagy公司[28]对于高维情况，其中参数数量远远大于观测数量；Gan等人[29]他比较了SLS问题的几种算法的性能；和Erichson等人[30]who通过变量投影研究了稀疏主成分分析。可分离模型具有广泛的实际适用性，因此形成了积极的理论和应用研究课题。例如，当分析可分离结构的“简化”非线性优化问题时，与原始NLS问题（例如[20，31]).
在下文中，我们将重点介绍作为常微分方程系统（例如[32]). 具体来说，考虑以下公式给出的方程组哪里在中获取值， ，和.为了我们的目的，我们明确地将函数中的线性分量与非线性分量分开F类通过设置哪里，和符号表示矩阵转置（参见[20]). 在这里，大小向量，表示与状态变量不可分离的“非线性”参数x个，同时，大小的矢量，包含“线性”参数；请注意.作为中的向量场(3)在线性参数向量中是可分离的，我们将相应的ODE系统称为参数中的线性（参考线性回归模型的情况），尽管解决方案θ对系统可能是高度非线性的，甚至是隐式的。
示例2。让
然后，人们看到了这个等式(1)是的特例(2)–(3).
2.2. 解决方案战略
让， ，是初值问题的解(2). 我们假设噪声测量在系统上的时间点收集.这种情况的常见统计公式是
这里是随机变量是不可观测的独立测量误差（不一定是高斯误差），平均值为零，方差有限。
Varah的ODE模型参数估计方法的工作原理如下。让表示使用样条或局部多项式获得的更平滑的数据（参见例如[33，34]和[35]用于处理各种平滑方法和大量参考书目）。此平滑器近似于解决方案至ODE(2). 瓦拉建议在方程式中插入平滑器(2)，现在只满足近似值，并将由此产生的参数差异最小化ξ和.最小化则为.这个想法在布鲁内尔建立了坚实的统计基础[36]古古斯维利和克拉森[37]. Varah的原始方法需要使用导数作为，这是一个缺点，因为众所周知，从噪声和稀疏数据估计导数可能相当不准确；参见例如，Vilela等人[38]还有周和沃伊特[7]或者更常见的是Fan和Gijbels[33]. 最近的研究[20，23，39–46]已经表明，将瓦拉的思想移植到方程的解水平上更有成效(2). 为了实现这一转变，我们定义了一个积分准则函数因为它在基于积分的估计方法中是典型的（参见上面的参考资料）。在这里，是欧几里德规范。最小化(5)超过产生一个参数估计器，其特征通常与基于微分方程本身的估计器略有不同。在实践中，积分被离散化并用和代替，因此可以使用典型的非线性最小二乘法来执行最小化，例如fmin搜索在Matlab中。离散格式为
NLS解决方案没有考虑ODE的特定线性形式(3)，但在(2).
正是在这个阶段，瓦拉建议利用可分性，而没有实际研究这种方法。在这里，我们提供了必要的细节（参见[20]). 表示
然后，与保持固定，将(5)由给定
哪里表示单位矩阵。符号和强调解对非线性参数的依赖性.这个解决方案重新插入(5)，得出简化的积分准则函数（参见[23]):
一次最小化超过和解决方案
获得，估计值ξ和θ立即跟进，并通过以下方式给出（轻度滥用矩阵转置符号）分别是。方程式(7)和(8)由瓦拉驾驶[17]上述建议。事实上，注意非线性优化仅用于估计非线性参数，与NLS方法相比，它可以大大降低非线性优化问题的维数。
从以上推导可以清楚地看出，SLS问题是一类特殊的NLS问题，对于不同的变量集具有线性和非线性目标函数。虽然Lawton和Sylvestre已经提出了利用可分性改进参数估计的想法[47]，随后的许多文献似乎都是基于Golub和Pereyra提出的变量投影方法[21]. Golub和Pereyra[26]回顾了30年来对此问题的研究。
3.仿真框架和结果
为了研究SLS和NLS的相对性能，我们设计并执行了一个大型Monte-Carlo仿真，其结果将在本节中介绍。
所有计算均在R（右）在Amazon EC2 m5a.4xlarge实例上使用西蒙德Yaari和Dattner软件包[25]（常微分方程的可分离积分匹配）。该包中应用的统计方法使用积分匹配标准函数的平滑和最小化，利用了微分方程的数学结构，如参数与方程的可分性。与基于导数的方法相比，光滑化和基于积分的方法在常微分方程参数估计中的应用可以获得更准确和稳定的结果[20]. 这里，我们在中使用了默认平滑和优化设置西蒙德在这方面，SLS和NLS都得到了平等的待遇。明确地，西蒙德使用交叉验证（参见，例如[35])以确定最佳的平滑量。有关使用该软件包的详细指南，请参阅Yaari和Dattner[25]. 再现数值结果的代码可以在GitHub上访问（请参阅https://github.com/haroldship/complexity-2019-code/tree/master/Final代码首次提交）。我们依靠ggplot2中的包R（右）见威克姆[48].
3.1. 蒙特卡洛研究设计
我们选择了几个在各种科学学科中产生的具有代表性和挑战性的ODE模型。那些是（i）用于模拟传染病传播的SIR模型（ii）具有正弦季节调整的Lotka–Volterra人口模型（iii）BST内的广义质量作用（GMA）系统，例如代谢途径系统（iv）沿神经元轴突的FitzHugh–Nagumo动作电位系统
下面给出了关于这些系统和我们使用的具体实验装置的更多数学细节。
在每种情况下，我们通过对系统进行数值积分，然后向时间进程中添加独立的高斯噪声来生成观测值，如(4). 我们考虑了各种参数设置、样本大小和噪声水平，如下所述。ODE参数通过NLS和SLS进行估算，如方程式所示(6)和(8)分别是。
作为性能标准，使用了执行优化所需的时间和所得参数估计的准确性。虽然比较计算时间很简单，但有许多选项可用于比较准确性。我们重点讨论了两种优化方案之间的主要区别，即它们处理线性参数估计的方式。SLS不需要对这些参数进行初始猜测。相比之下，NLS确实需要对每个线性参数进行良好的初始猜测；否则，它可能会偏离或陷入局部极小值。因此，寻找非线性优化问题的“好”解通常需要在参数空间中进行“好”的初始猜测。显然，关于这些参数的一些“先验信息”对于优化目的至关重要。关键的见解是，这些先验信息被封装在ODE本身的数学形式中，例如(3). 重要的是，虽然NLS没有考虑ODE的特殊数学特性并以统一的方式处理所有参数，但SLS并非如此。因此，人们可能会先验的当关于线性参数的先验信息质量较低时，预计SLS比NLS更有效，可能更准确。另一方面，当一个人拥有关于线性参数的高质量先验信息时，我们预计SLS和NLS的表现将类似。人们可能会注意到，几乎所有GMA和S-系统中的非线性参数都是非常有界的，通常在-1和+2之间，它们的符号通常是已知的，而线性参数在GMA系统中是无界的，在S-系统中是非负的，并且对它们的大小一无所知（见[18]). 因此，不需要SLS中线性参数的先验信息可能是一个巨大的优势。
对于Monte-Carlo研究，我们通过对参数值使用高、中、低质量的初始猜测来改变先验信息。这里，更高的质量意味着最初的猜测更接近事实。更具体地说，NLS使用的线性参数的初始猜测是以真实参数值为中心的高斯随机变量，其标准偏差等于真实参数乘以先验信息值（换句话说，先验信息值也可以理解为“事先分发”）。诚然，“高”、“中”和“低”的具体量化有点主观，并且在不同的ODE模型中有所不同，如下所述。为了使优化算法（尤其是NLS）更快更好地收敛，将非线性参数约束在给定的范围内，无论我们如何改变线性参数的先验信息，该范围都是相同的。此外，在每个Monte-Carlo迭代中，我们对NLS和SLS的非线性参数使用了完全相同的（伪随机）初始猜测。因此，就非线性参数的信息而言，无论线性参数的先验值如何，每个基准模型都保持不变。因此，这两种算法都获得了关于非线性参数的相同先验信息，并且都没有得到优先处理。
我们使用的噪声级（信噪比，SNR）定义如下。对于给定的解决方案对于ODE方程，我们计算标准偏差.然后SNR，比如，和由给定，和，分别，其中σ是高斯测量误差的标准偏差ϵ，定义见方程式(4). 我们将这些SNR分别称为“低噪声”和“高噪声”（参见[49]; 尽管背景不同）。然后，我们比较了所得参数估计值的均方误差（MSE），从而在统计上可识别的ODE模型中进行了有效的比较（参见[20]相关定义和结果）。作为另一种精度测量，我们使用了标准(5)和(7)以最佳参数值进行评估。我们提出的两个标准虽然合理，但却是不同的。因此，在每一个实验装置中，预计它们不会一致。然而，他们在第节中得出的全球结论5连贯且有利于SLS。
我们现在提供有关模型和实验装置的数学细节。
3.1.1. SIR年龄组
感兴趣的系统是SIR型（已恢复的易感感染）的流行病学模型，包括年龄组和季节性成分（例如[50]). 各年龄组的感染过程和每个季节使用两个易感比例方程进行描述(S公司)并被感染(我)种群中的个体（恢复个体的比例由):
模型的参数为传输矩阵β、回收率γ，以及，表示季节性流行的流感病毒株等的相对传染性与第一季相比(用作参考并固定在1）。如Yaari等人[46]，考虑到该模型的可分性特征有利于数据拟合。具体而言(9)初始值为非线性，它们通常是未知的，必须进行估计。就我们的目的而言，考虑一个具有一个年龄组和一个季节的模型就足够了。使用了以下参数设置：.我们考虑了两个样本大小，18和36，以及两个噪声级，和.之前使用的信息是，分别对应于高、中、低质量。蒙特卡洛研究的规模为500个模拟。
3.1.2. Lotka–Volterra，季节性强制
作为另一个基准，我们考虑了经典捕食者-食饵模型的扩展，即Lotka–Volterra模型，该模型包括捕食率的季节强迫，使用两个控制幅度的额外参数()和相位(ω)强制：
非线性参数为ϵ和ω.我们考虑了时间间隔内的动态.参数设置如下所示
初始值为.研究了四种实验场景，对应于100和200的样本量，以及和.先前的信息值为，分别对应于高、中、低质量。蒙特卡洛研究的规模为500个模拟。
3.1.3. GMA系统
我们分析的GMA系统由三个变量中的三个微分方程组成([18]; 第84–85页）。他们是
这里的线性参数是速率常数γ而非线性的是指数动力学阶（f）注意，参数（f）被允许变成负值，他们的符号可能是已知的，也可能是未知的。我们考虑了时间间隔内系统的动态.参数设置如Voit中所示[18]; 即，
初始值为.研究了四种实验场景：样本大小为100和200，SNR为和.先前的信息值为，分别对应于高、中、低质量。蒙特卡洛研究的规模为500个模拟。GMA系统的参数估计被认为是一项具有挑战性的数值任务[18].
3.1.4. FitzHugh–Nagumo系统
FitzHugh–Nagumo（FHN）系统[51–53]模型在神经元传递中激发动作电位。它是由
这个带有两个状态变量的系统是为了简化霍奇金和赫胥黎的研究中提出的更复杂的模型[54]用于研究和模拟巨型鱿鱼轴突的神经功能。FHN模型在神经生理学中用作观察到的动作电位的近似值。
系统(10)参数为线性，和b，但非线性c（c）。我们考虑了两种样本大小，和，和两个SNR和.参数设置为.初始值为.在时间间隔内获得了真实的解.此处使用的先前信息是，分别对应于高质量、中等质量和低质量（参数的初始猜测被保证为正）。蒙特卡洛研究的规模为500个模拟。许多研究人员研究了FHN模型的参数估计问题。特别是Ramsay等人[55]坎贝尔和斯蒂尔[56]还有拉姆齐和胡克[11]指出了估算ODE系统参数的几个困难。
3.2. 蒙特卡罗分析结果
我们的发现通过图表和表格呈现。主要总结为表格1和2，其中我们报告了线性参数估计的均方误差（蒙特卡罗模拟的平均平方误差）比率（对于非线性参数，请参阅本节末尾的讨论）。从表格中可以得出几个结论。（i）给定高质量的先验信息，NLS和SLS的精度是可比的，并且在各种实验设置中都不优越（从表中的原始数字中看到的至少一些差异可能归因于蒙特卡洛模拟误差，因此似乎微不足道）。（ii）当先验信息质量下降到中等或较低时，SLS的性能在大多数情况下都明显优于NLS（在一定程度上取决于具体的实验设置）。（iii）对于固定噪声级，随着样本量的增加，SLS的优势变得更加明显。（iv）对于固定样本量，随着噪声水平的增加，SLS仍优于NLS，但程度较小。
这些结果也可以通过简单统计图的组合进行可视化。因此，图1显示了在几个实验设置下比较两种方法的MSE的折线图。而表中的数字1和2是MSE的比率，此处的数字表示MSE的绝对值。从图中可以看出，SLS相对于NLS的优势在于先验信息不太理想。请注意，在这个特定的设置中，SLS在获得高质量的先验信息方面的表现比NLS差。一个似是而非的解释可能如下：然而，在我们的实验设置下，SLS通过(3)在整个仿真过程中是固定的，NLS原则上可以接收对线性参数的任意精确的初始猜测。因此，人们可以设想一个阈值，即使用后一种信息比使用结构关系的好处更大(3). 然而，除了观察到它似乎在对可能的参数值有很好了解的场景中表现出来外，几乎不可能对这种现象进行精确量化。实际上，这种理想的先验信息是罕见的。
图的面板（a）1进一步表明，在我们报告的特定场景中，SLS随着噪声水平的降低而提高，这与同一图中的NLS不同。
图2是NLS和SLS损耗的散点图(5)和(7)（对数标度）在最佳参数估计值下进行评估。该图以另一种方式强调了NLS先验信息的重要性：很明显，后者的性能受到初始参数猜测质量的强烈影响。同样，当先验信息质量较高时，NLS和SLS的表现类似。然而，当先验信息的质量不如理想时，就像在大多数应用中一样，NLS会比SLS差得多。散点图还快速反映了估算结果的可变性。
我们从图中得出的结论2由图的面板（a）确认3，显示根据标准测量的NLS和SLS损失箱线图（对数标度）(5)和(7). 模式很清楚：SLS比NLS表现更好，而且NLS的劣势随着先验信息的降低而变得更加明显。
图的面板（b）3总结了计算时间。SLS通常要快得多。NLS的执行时间受先验信息质量的影响，有趣的是，随着这种质量的提高，NLS的运行时间也会增加。所有其他模型（和噪声级）的结果相似，因此省略。
最后，表格3和4提供有关非线性参数的信息。在NLS的情况下，可以观察到关于线性参数的先验知识如何传播到非线性参数的估计精度中。特别是，对于线性参数的先验信息不太理想的情况，即使在非线性参数的情况下，SLS也比NLS具有明显的优势。
4.展望
复杂动力系统中的数据拟合仍然是一个具有挑战性的问题，即使利用了可分性，也无法以轻松的方式进行处理。例如，为了揭示第节中的模式3在这项工作中，我们必须仔细设计实验研究，因为否则模拟可能不会收敛或收敛到较差的解。NLS和SLS都是如此，但无论何时观察到它们，NLS的收敛问题都要严重得多；他们对FitzHugh–Nagumo系统的情况特别敏感。这一结果突出了先前信息关于用于优化的参数或不同表达的初始参数猜测的质量。这里我们主要关注关于线性参数的先验信息的影响。然而，也很明显，非线性参数的先验信息对于优化目的具有同样重要的作用，这对于NLS和SLS都是如此（数据未显示）。
作为我们探索性工作的结果，我们展望了未来的以下有希望的研究方向。
4.1. 动态系统SLS的数值实现
我们分析中的所有计算都是在R（右）使用西蒙德Yaari和Dattner软件包[25]. 然而，使用ODE的可分性属性的想法与特定的编程语言无关，也可以在其他软件包中实现。事实上，自从Golub和Pereyra首次引入变量投影法以来，在变量投影法的背景下已经做了很多工作[21]. 在非线性回归背景下，Golub和Pereyra的变量投影法[21]在中实现R（右）在中国家统计局命令；请参阅Venables和Ripley[58]以及第218-220页的应用示例。此外，我们意识到TIMP公司马伦和范斯托库姆套餐[59]; 它实现了变量投影方法。因此，下一步可以是结合不同包的优势，例如。，西蒙德和TIMP公司为了开发动态系统中变量投影的高级软件。
4.2。特定类复杂动力系统的定制算法
众所周知，优化方案的性能关键取决于用于描述数据的底层数学模型。因此，不同类别的动态模型似乎需要根据其特性定制特定算法。例如，GMA系统的参数估计与使用SIR时遇到的挑战不同（参见第节3). 我们预计，将未来的研究重点放在特定类型的模型上，并为其参数估计开发稳定的算法，会有很多收获。
4.3. 动态系统背景下SLS的理论性质
古古斯维利和克拉森[37]在平滑的一般背景下研究了NLS的统计特性，而Dattner和Klaassen[20]专门针对参数线性（函数）的ODE系统。人们可能会认为古古斯维利和克拉森使用的一些假设[37]当问题接近达特纳和克拉森所考虑的问题时，可以放松[20].
4.4. 部分观测、高维和未指定动态系统的扩展
最近处理高维ODE模型推断的工作表明，利用参数中的线性对开发成功的估算方法至关重要（参见[24，39]). 更一般地说，使用变量投影方法来研究部分观测、高维和可能指定错误的动态系统的情况是很有趣的。这项工作可能还需要使用高维正则化技术（例如[39])用于平衡数据和模型，并特别考虑到潜在的模型错误指定（请参见[55]).
5.结论
在这项工作中，我们设计了一项广泛的模拟研究，以探索动态系统中两种推理优化方案的相对统计和计算性能：典型的非线性最小二乘法（NLS）和一种新的可分离非线性最小二乘法（SLS）。作为基准，我们考虑了多种生物领域中广泛使用的ODE模型。我们通过估计值的均方误差（MSE）测量了这两种方法的统计性能。作为另一个性能标准，我们在最佳参数估计值处使用了损失函数值。通过完成每个优化所需的执行时间，还比较了这些方法的计算性能。
我们的总体建议如下：每当复杂动态系统包含大量线性参数时，其参数的估计应采用可分离非线性最小二乘法，而不是更常用的通用非线性最小二乘法。从我们的研究中得出的一般模式是，SLS的性能至少与NLS一样好，并且通常比NLS更好，尤其是在有关系统的先验信息不理想的情况下，这是实践中的典型情况。这一说法在所有测试模型中都是一致的。
数据可用性
再现数值结果的代码可以访问https://github.com/haroldship/complexity-2019-code/tree/master/Final_code_First_Submission网站.
披露
文章的预印本可在ArXiv上找到，网址为https://arxiv.org/abs/1908.03717.
利益冲突
提交人声明他们没有利益冲突。
致谢
Itai Dattner得到了以色列科学基金会387/15号拨款的支持。埃伯哈德·沃伊特获得了NSF-MCB 1615373（PI：戴安娜·唐斯）和NIH 2P30ES019776-05（PI：卡门·马西特）的资助。资助机构不对本文的内容负责。
