欢迎来到网络研讨会。我的名字叫Sietse Braakman。我是一个在MathWorks应用工程师。今天,我将谈论与SimBiology全局灵敏度分析。
今天我想讨论的话题是第一个带你通过全局灵敏度分析的一些概念,以确保每个人都理解并与这些概念在相同的页面上。然后,我将解释如何执行在SimBiology全局灵敏度分析。我们将使用两种方法,Sobol法和multiparametric方法。我将花时间在如何解释这些情节,以及如何你们的样品大小,这样你得到可靠的结果。
我们会有时间的问答。有一个问答窗口网讯的一部分,你可以输入你的问题在会议期间。我的一些同事,Fulden Buyukozturk和杰里米•Huard能够回答。和其他人,我们可能会保持到最后的问答环节。
,让我们开始吧。所以我们先的一些概念。你们中的大多数可能熟悉这一事实有本地和全局灵敏度分析。当我谈论当地的敏感性分析,讨论分析在一个操作点的参数空间。所以你可能会有多个参数,实际上,我们只是执行此灵敏度分析为一组参数值。
全局灵敏度分析,另一方面,跨域参数空间中执行。所以对于每一个参数,现在有一个上限和下限。之间的界限,我们探讨模型输出与输入参数很敏感。
敏感性分析是样本中一个重要的概念。为当地的敏感性分析,这是一次完成的。我们在一个操作点。我们扰乱一个参数,看看影响模型的输出。我们把参数改回到原来的值,我们会扰乱下一个,等等。这是一次。因此,你无法从单一吸收天文台参数之间的相互作用分析。
全局灵敏度分析,另一方面,都是一次。所以你带参数的随机样本空间来计算敏感性指数。因此,你可以观察到参数之间的相互作用。然后也有多种方式,你可以样本参数空间。一般来说,这些树苗是统一的,但也有所谓的低偏差抽样方法,如Sobol、拉丁超立方体,哈尔顿序列,您可以使用它来执行抽样。
最后,有多种方法可以计算敏感性指标。对于当地的敏感性,我们使用导数或比率,我们如何看模型输出的变化对模型输入的变化。Sobol和eFAST, Fourier-based方法,你用方差。所以你想属性为每个参数方差。然后有发布出去是依靠方法,如multiparametric全局灵敏度分析和correlation-based部分等级相关系数等方法。
接下来要讲的是为什么我们要使用本地或全局灵敏度分析。这里我有一个非常简单的例子,某个部位模型与吸收ka,分布体积Vd,和消除——参数化酶消除Vm和K马克斯,Vm和公里。现在,如果我每一个参数设置为1,这个模型模拟了10小时,我得到以下当地灵敏度值。然而,如果我ka设为0.1,所以我k除以10,我做同样的分析,我得到这些结果。
所以你可以看到,来自当地的敏感性分析的结果强烈依赖于参数空间操作点。这就是为什么建议使用全局灵敏度分析操作点的如果你不知道多少信心。所以在这种情况下,全局灵敏度分析最合适当你探索敏感性参数域。
不过,有一个原因你可能想要使用当地的敏感性分析,例如,对于目标识别。如果你有一个好主意的校准,你有,例如,一个校准,代表一个kidney-impaired耐心,你可能想要使用局部灵敏度分析对目标识别,如Birgit Schoeberl和我们的同事在下面的纸。而且,SimBiology MATLAB使用当地的敏感性计算梯度优化算法。
现在,另一方面,全球灵敏度分析更适合当你想要了解哪些模型输入驱动模型响应或基于模型的决策指标,如果你不知道操作点。然后你可以使用全局灵敏度分析的结果也通知参数估计策略。模型的参数非常敏感,可以最有可能,你想估计,确保你有一个好的理解它们的价值是什么,你的模型是正确的校准。
,我想继续Sobol全局灵敏度分析。这是其中一个在SimBiology实现的两种方法,我想解释第一Sobol指数是如何计算出来的。的想法是分摊方差模型的输出输入X Y的模型。
和有多个指标可以计算。一阶灵敏度指数代表了个人的贡献,因此个人可以方差,方差可以分配到每个参数。通过这个表达式的计算。
分母,你看到的无条件方差。我们假设模拟单舱模型,我们有,我们采样两个参数,吸收系数ka和间隙。现在,如果我这些参数样本,每个样本模拟,得到模拟的合奏。在每一个时间点,我可以计算总体的方差,这就是发生在分母。
现在,在分子上,你看到的条件方差。方差不是因为ξ。习说,卡是我们的参数。那么我们要做的就是解决ka在一个点。我们仍然允许间隙变化,得到方差。
现在,我们能做的所有值的ka,取平均值,然后我们可以说是由于方差不是ka。所以1 -整个的价值比率给出了方差ka所致。这是第一个订单的方式Sobol指数计算。所以它的条件方差的比值无条件方差,然后1 -。
其他敏感性指数可以计算的总效应,和显示交互,因为它是所有的低阶效应的总和,包括你感兴趣的参数。如果你的参数或兴趣是第一或ka,一阶,二阶与每个其他的参数,每个组合的三阶参数,等等。和这个数字基本上是有一个更简单的方法来计算,是通过这个表达式,我不会讲太多细节。但计算的概念类似于一阶灵敏度指数。
好,这给了你一个想法Sobol指数的计算,以及不同于当地的敏感性分析,在那里你的比率变化模型输出变化模型的输入。让我们——之前在SimBiology显示如何做到这一点,我只是想带你通过工作流,我们要遵循当我们搬到SimBiology。
所以我们要做的第一件事是我们要在参数空间中定义感兴趣的领域。所以我想说,这些都是我感兴趣的参数将在我的全局敏感性分析,对于它们中的每一个,我需要决定低,上限是什么。然后定义我的参数空间,我可以样品。我可以模拟每个样本。一旦完成我所有的模拟,我可以计算灵敏度的措施。
我要证明这一点由Sergey Aksenov和他的同事使用一个模型。这是一个模型,描述了lesinurad febuxostat,这是两个批准的药物来治疗痛风。和模型是这样的。我们有一个两舱制模型lesinurad顶部和底部的febuxostat两舱制模式。febuxostat可以看到,中央集中,会影响血清尿酸的生产,而lesinurad增加肾小球滤过,从而基本上增加尿酸的清除,以避免积累。
今天我要做什么,我要从这一模型,选择四个参数尤其是PD模型的一部分,以确保我们——我们探索这四个参数,看他们在多大程度上影响输出,我选择血清尿酸。所以这四个参数是k1 fc50 fmax, e0。当然,我可能需要更多的参数,然后我们要花大量的时间模拟——今天的目的,它不是有用的。
好了,有两种方法,你可以执行这个分析。你可以执行它完全以编程方式,所以您可以编写自己的脚本,或者您可以使用一个应用程序来执行分析。应用程序,您可以下载,如果你去MATLAB,家庭选项卡,并单击插件浏览器,和你搜索全球SimBiology灵敏度分析,可以直接下载并安装这个应用程序到你的机器。
如果你不能访问这个出于某种原因,注意,函数本身是建立在和安装在2020年——2020年SimBiology。所以这两个函数是sbiosobol sbiompgsa,应用基本上是一个用户界面,这两个功能。所以我们要做的是,今天我们将使用应用程序,这样更容易让你跟着我做什么。但是最终的演示,我将与你分享我的代码和模型,这样你也可以试试使用代码。
为了运行这个应用程序,我需要做的就是我需要这个模型转移到MATLAB。所以我可以将此模型导出到工作区,MATLAB的工作区。例如,称之为m1。然后在MATLAB,我可以退出剂量。所以在这种情况下,我想用剂量1和2。我已经做了所有这些,开始应用,但你这样的调用应用程序,startGlobalSensi tivityAnalysisApp,然后通过模型和剂量。启动应用程序,它看起来像这样。
让我们走过这一步一步。首先,正如我所说的,我们需要做的就是我们需要定义参数我们感兴趣。所以我可以去这里。基本上打开另一个窗口的所有对象在我的模型中,隔间,参数、初始条件的物种,等等,我可以使用作为输入来全局灵敏度分析。我可以排序。
所以我现在感兴趣的是选择这四个参数感兴趣的我,e0 fc50, fmax, k1。现在我可以选择值,上界和下界值,每一个。我已经做了,最后我基本上把模型值,标准模型中的参数值,然后除以1又1/2乘以1又1/2,主要因为它更容易总是绑定1。所以我想确保所有的参数也有类似的宽度数量级的上限和下限之间的传播。
你可以编辑和写作0.5,然后改变参数值。然后当你完成,你选择所有参数和输入所有的上界和下界,我们继续前进。在这里,您可以选择样本的数量。默认情况下,这是1000年,我要保持这种方式。
然后你需要定义输出的时候,这基本上是MATLAB代码。所以你可以使用linspace,但您也可以编写自己的时间向量。所以我要用向量从0到180年的0.5小时的步骤。
然后这里有几个下拉菜单。你可以选择不同的抽样方法。就像我说的,Sobol,哈尔顿和拉丁超立方体,那些低抽样方法,差异——我建议使用这三个,因为它的一个更有效的比仅使用一个标准的随机均匀分布。
报告输出次仿真完成之后,他们可能不配合,ODE求解器的步骤。这就是为什么你可以选择一个插入器。我只会推荐使用默认值。
然后有不同的方法可以加快这个全局灵敏度分析。首先,你可以并行化模拟。我将使用一个平行四个工人在这里,在我的酷睿i7与四核笔记本电脑。我们也将加速模型,所以,编译模型寻求代码为了加快模拟。然后我们可以以后再讨论。
这就是输入全局灵敏度分析的一部分,我们已经建立了。接下来我们能做的就是我们可以定义为我们感兴趣的输出是什么。好,这个模型痛风,兴趣是血清尿酸水平的输出。在这里,你可以选择你感兴趣的任何输出。如果你只是感兴趣看lesinurad pKa,您可以使用这个——中央lesinurad浓度。
但在这种情况下,我们非常感兴趣的临床相关的输出模型,即血清尿酸。所以我就选择并点击完成。然后我们可以开始模拟模型。和我将讨论这一点。但基本上,有四个参数的组合和绘画1000个样本意味着我们需要做1000 * 4 + 2,所以6000模拟。
所以我要继续,现在开始,与此同时,我会回到SimBiology的幻灯片来讨论另一个新特性是与这个特定的相关情况。这些被称为可见。所以可见背后的想法是,他们取代和消耗计算统计功能。
所以你可能熟悉计算统计功能,你在任务编辑器中模拟模型在2019年之前,你可以计算,例如,cmax什么的。它一直是一个标量值。可见,建立在最重要的是,他们不一定导致标量,但他们也可以基于时间的。所以基于时间相关,例如,你可以有不同的物种,你想加起来,例如,总肿瘤体积,或全部药物浓度,之类的。
和这个想法是可见的,可以适用于你的仿真数据。所以说,100年你做了模拟,你这100的AUC从每个模拟。所以你有一个100年1 SimData数组,你可以添加一个可观测SimData数组,就说,好吧,给我的AUC等于trapz时间,中央。drug_central,然后它就会给你100 auc,基本上。这将使你的生活更容易。你可以用一行代码。
现在,有些人可能意识到,这也可以通过使用重复作业。所以重复赋值的区别和可见,可见计算常微分方程后解决。所以他们不是方程组的一部分。所以可见不能系统依赖的变量。如果你的模型或系统取决于肿瘤体积,例如,改变隔间肿瘤体积,然后它需要重复任务。但如果没有在模型中取决于它,最好使用一个可观测的,因为它是更少的计算昂贵。
可见可以应用到你的模型和SimData对象。如果你添加一个可观测模型和模拟模型中,可观察到的自动计算。SimData,如果您添加可观测到SimData,然后计算并返回,例如,AUC。
好的,所以我希望有意义可见。让我们回到SimBiology。它看起来像它现在所做的,我们有我们的结果。上次,我这今天早些时候。6000年花了大约90秒来模拟这些模拟。所以,好的。
我们这里有结果。我们有蓝色的一阶敏感性,Sobol敏感性,在红总敏感性。所以你可以看到e0显然是最重要的参数模型中,可能随后k1, fmax,然后fc50。
现在,如果没有相互作用,这些都应该加起来1在任何时刻或大约1。你可以看到,并不完全如此。所以——这是无法解释的方差不是0的一小部分。这给你一个指标之间有相互作用参数。
还有其他的原因可能是零,这是,你可能会有一些在你的模拟数值漂移什么的。但最常见的解释是,有相互作用。你可以看到,有相互作用通过比较总订单值,例如,在这里,一阶值。
如果我们比较这两个,你可以看到,这个比这个,这是你所希望的。你期望总要相同或高于第一个订单,因为总订单是一阶的总和,所有的第二,第三,等等,订单。因此,通过比较这两个,你可以看到哪些参数显示,多数交互。
好吧,这为我们提供了一个全球敏感性的概念。是血清尿酸的输出,这是一个连续变量中,随时间变化的模型。和海浪,你在这里看到的是不同剂量事件。所以我们可以看到的是这里e0看起来更重要的早期比后期的仿真。
和其他的事情,当然,我们可以做我们可以使用一个实际的观测,标量观测,如AUC或血清尿酸的最小值作为我们的模型输出。当然,我们不会有时间,但是我们会有一个数字对于每一个标量值,为每个第一和总指数。
你可以画出的另一件事,通常是不错的,这种完整性检查,你可以画出数据。如果你画出数据中,您可以看到在这里所有的模拟结果与90年百分比地区蓝色,和一些个人的痕迹点缀。在红色的,你看到的意思是模拟值的所有样本,我们的1000个样本。
好,那么,我们将进入multiparametric全局灵敏度分析。和multiparametric全局灵敏度分析背后的想法是,你使用一个标识符来分析参数的相对重要性,分类器是一个基于模型的决策指标基于模型输出,但却导致一个正确或错误的结果。因此,必须有一个不等式。所以效果——例如,分类器可以,我的药效学效应大于70%或最终的浓度大于最终浓度的均值。
所以,你基本上可以分类之间的模拟,他们满足这一分类器吗?他们是事实上超过70%吗?他们还是失败了吗?他们拒绝了吗?
您可以使用这些字符的组合。,您可以创建一个单一,效果更大,最终的浓度比。或者你也可以执行multiparametric全球敏感性与多个输出,与多个分类器,看看哪些是相关的。
这是怎样的一个例子——它是什么样子。说你有一个模型有两个参数,你不同,凯尔和IC50。这两个值,然后模拟模型,你看看效果大于70。如果是这样,那么你说,是的,否则,你说不。因此,我们得到了一组模拟,要么接受或拒绝。
然后,我们可以计算一个经验累积分布函数。如果,我们这样做样品的接受和拒绝。所以你可以看到一个例子。k1,是这种情况。所以在k1值低,看起来更多的样本是接受,而在更高的价值,更多的样品被拒绝。
我们可以做我们可以计算最大距离。这是垂直两个函数之间的最大距离,基本上。我们可以使用一个Kolmogorov-Smirnov测试是否这些分布在统计上显著不同。所以这样做的好处在Sobol实际上是,你得到一个答案,有显著不同的结果。当然,你可以选择指标,分类器为你的案子有关。
现在,有一件事我还没接触过的阈值。阈值,例如,70%,应该是在你的模拟。大约1/2的模拟应该通过和1/2的模拟应该会失败,为了能够构建这两个累积分布函数。
因为如果比通过更多的人失败了,那么你会得到更少的传递。和你会得到一个非常参差不齐的运作,然后你Kolmogorov-Smirnov测试是不可靠的。所以你可以做的是,你可以使用这里的模拟,情节的模拟,例如,用红线想出一个阈值。
好吧,那么,我们如何在实践中把这个?好吧,其实很简单。我们可以回到这里,我们可以重用的模拟。我们已经做了蒙特卡罗模拟。所以我们可以重用这些模拟。我唯一需要定义是我的分类器。
所以我可以输入分类器。我可以,我不需要运行模拟了。我可以计算出multiparametric全局灵敏度分析。
所以在这里,你看到multiparametric全局灵敏度分析的结果。现在你看到k1和fmax——两者之间有显著差异。e0,几英里远。fc50,他们非常接近。
你可以比较这直方图。累积分布函数和直方图是相关的。所以你可以看到,例如,e0,大多数接受样品发生在较低的e0的价值观。这就是为什么你看到这个蓝色线上升。而拒绝样品只发生在更高的值,这就是为什么红线的只有0,如高于0,e0的这些更高的值。这就是直方图和eCDF是相关的。
我们看这个Kolmogorov-Smirnov测试是如何不同的这些直方图或不同的是——真的,Kolmogorov-Smirnov测试不同如何这些累积分布函数。和我们如何做,我们可以画出——被称为酒吧。这情节——计算ks统计,你看到蓝色的,无论是统计上显著不同的p值。
然后你看到fmax, e0, k1,效果是显著的。所以这两个分布明显不同。而对于fc50,它不会达到p小于0.05的阈值。
好的,这总结了multiparametric全局灵敏度分析。还有一个主题我想讨论迭代的数量。为了执行这个全局灵敏度分析,计算,你需要进行大量的模拟。这是计算昂贵。但同时,为了得到可靠的结果,你需要——你不能undersample。所以你知道,什么是好的样品?
如果我们假设N是样品的数量,我们从参数输入空间和P是输入参数的数量,这是参数输入空间的维数,然后我们能说哪个更大。所以N应大于2 p的力量是绝对的最低限度。
如果我有两个参数,然后我带四个样品。只允许我介绍我的输入空间的角落。如果我有三个参数,这是八个,所以这是立方体的角落。所以,你想要的是一个较高的基数。所以,3 P的力量或4 P的力量。
然后你可以看到,你知道如上P增加,15,你看大量的模拟。一般而言,——如果你有——如果你执行这不到五个或更少的参数,我建议1000 - 5000模拟,像我所做的。我做1000这四个参数的模拟。
所以尝试。有几个方法,你会发现你实际上欠采样,这是如果你Sobol结果低于0。这是略高于0,但如果这是消极的,那么我将会担心我,欠采样——或者如果他们以上的1。
因为Sobol指数总是试图看看比例的方差可以总方差可以归因于单一参数。它应该是在0和1之间。这无疑是症状的东西是错误的,你应该增加样本的数量。
当然,另一个是,如果你重新运行分析,你要有不同的结果,但是你可能没有计算资源尝试多次。MPGSA最后,如果你有,,如果你有一个锯齿状或阶梯形状,我之前提到了这一点,那么你可能有一个低的阈值,阈值的选择可能不是这样,大约1/2是拒绝,1/2被接受,但它也可以是由于样本量不足。你需要知道,如果你有锯齿状的楼梯,你不应该相信从你MPGSA结果。
最后,因为我们所做的,我们要做很多模拟和因为我们有模拟所有的样品在我们可以计算灵敏度的措施之前,那么所有的仿真结果需要在内存中以计算敏感性指标。所以如果你可以减少内存占用你的模拟,您可以执行更多的样品。你可以模拟更多的样品。
所以你怎么能最小化内存占用?嗯,首先,你可以减少输出。如果你只有一个输出,最终的浓度或最终效果,将大大减少输出的数量。并限制登录状态的数量。所以确保没有——你不记录你所有的物种,等等。然后最后,当然,你可以使用一个可观测的为了减少足迹,因为整个模拟基本上可以归结为一个标量使用,可见,这将允许您采取更多,运行更大的全局敏感性分析。
好,总之,我们谈论我们需要做什么来执行这个敏感性分析。我们需要定义该域,每个参数的上界和下界。然后我们可以样本域。对于每一个样品,我们可以模拟模型。一旦我们模拟所有的样本,我们得到的模拟。从那里,我们可以计算灵敏度的措施。
所以总的来说,你应该考虑的一些事情,你想要包括哪些参数,这些参数的范围。生理上和范围可以定义,也可以在文学,发现不同的值取最小值你发现你发现和最高价值。然后你要样品。你需要选择一个抽样法。所有这些样品必须统一为基础的假设Sobol指数计算不违反。
所以拉丁超立方体,这Sobol序列,序列和哈尔顿,他们都是统一的取样方法。但是他们有一个优势的随机均匀采样,他们需要更少的样本,同时还覆盖相同的输入空间。所以我建议使用拉丁超立方体、Sobol或哈尔顿的抽样。然后你还需要选择样本的数量。
然后你要模拟模型,所以你需要确保你添加剂量。如果你不添加剂量,你不扰动模型,你不会看到任何结果。变异需要是正确的。您需要定义是什么时间点,感兴趣的模型的输出是什么,你的分类器是什么,等等,这里思考你的内存占用。
最后,我们可以谈论,sbiosobol sbiompgsa,基本上工作引擎盖下面的两个函数的应用程序,我给你,他们将计算一阶和全序Sobol指数和eCDFs,并执行钴测试所有的模拟,然后运行。
