基于生理学动力学模型的R软件pksensi：全局灵敏度分析和校准工具

软件X 12（2020）100609原始软件出版物pksensi：一个R软件包，用于在基于生理学的动力学建模放大图片作者：谢南鸿a，Brad Reisfeldb，Weihsueh A.邱亚a兽医综合生物科学，兽医和生物医学科学学院，德克萨斯州农工大学，学院站，德克萨斯州，美国b美国科罗拉多州柯林斯堡科罗拉多州立大学化学与生物工程学院ar t i cl e i nf o文章历史记录：收到2020年收到修订版2020年7月8日接受2020年保留字：基于生理学的动力学建模R软件包a b st ra ct灵敏度分析是建模者了解模型参数影响的重要工具。设置模型输出。它也越来越多地用于开发和评估基于生理学的动力学（PBK）模型。例如，一些研究已经应用全局SA来减少基于贝叶斯马尔可夫链蒙特卡罗的校准过程PBK模型中的计算负担。 虽然有几种SA算法和软件包可用，但没有全面的软件包存在 它允许用户无缝地求解PBK模型中的微分方程，进行SA结果并将其可视化，并区分可以固定的非影响性模型参数和需要校准的参数。因此，我们开发了一个R包，名为pksensi，使全局SA更容易在PBK建模。 该软件包可以研究PBK模型中的不确定性和灵敏度，包括具有多变量模型输出的模型。 它还包括检查全局SA结果的收敛性的功能。总体而言，pksensi改善了执行全局SA的用户体验，并可以为PBK模型校准中的决策制定创建稳健且可重复的结果©2020作者由爱思唯尔公司出版这是CC BY-NC-ND下的开放获取文章许可证（http://creativecommons.org/licenses/by-nc-nd/4.0/）。代码元数据当前代码版本1.2.0用于此代码版本的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX_2020_29法律代码许可证LPGL-3.0使用git的代码版本控制系统软件代码语言、工具和服务使用R编译要求，操作环境依赖性ggplot2，data.table，deSolve，doppler，getPass，magrittr，foreach，parallel，doParallel，重塑如果可用，链接到开发人员文档/手册https://nanhung.rbind.io/pksensi/问题支持电子邮件d99622005@ntu.edu.tw软件元数据当前软件版本1.2.0此版本可执行文件的永久链接https://cran.r-project.org/web/packages/pksensi法律软件许可证GPL-3计算平台/操作系统Linux，MacOS，Microsoft Windows安装要求依赖项ggplot2、data.table、deSolve、dapr、getPass、magrittr、foreach、parallel、doParallel、reshape如果可用，请链接到用户手册-如果正式出版，请在参考列表中引用该出版物https://cran.r-project.org/web/packages/pksensi/pksensi.pdf问题支持电子邮件d99622005@ntu.edu.tw*通讯作者。电子邮件地址：wchiu@cvm.tamu.edu（W.A.Chiu）。https://doi.org/10.1016/j.softx.2020.1006091. 动机和意义灵敏度分析（SA）是一种研究模型参数变化如何影响模型输出的数学技术2352-7110/©2020作者。由爱思唯尔公司出版。这是一篇开放获取的文章，使用CC BY-NC-ND许可证（http://creativecommons.org/licenses/by-nc-nd/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softxN.- H. 谢湾，巴西-地 Reisfeld和W.A.Chiu软件X 12（2020）1006092=越来越多的研究使用SA来确定哪些模型参数导致模型预测的高度变化[1]。该技术也已应用于药理学和毒理学研究[2，3]。药代动力学模型描述了体内物质的浓度或量随时间的变化（动力学模型的各种术语，包括“药代动力学”、“毒代动力学”和“生物动力学”模型，在此互换使用）。生理动力学（PBK）研究中SA的目标是检查输出变量（如血液或组织中的化学浓度）相对于输入参数（如解剖学、生理学和动力学常数）的灵敏度[4]。 此外，SA可以指导实验设计和参数估计过程[5]。例如，SA可以识别 在模型标定之前，可以进一步确定参数的优先级和参数的固定.在我们早期的工作中[7]，我们开发了一个更新的工作-流程应用全局SA以减少基于生理学的药代动力学（PBPK）模型的基于贝叶斯、马尔可夫链蒙特卡罗（MCMC）的我们使用GNU MCSim [8]，一个用于贝叶斯群体PBPK建模的有效模拟我们发现，扩展傅立叶振幅敏感性检验（eFAST）是一种基于方差的全局SA算法，对于复杂的多室、多数据集和多代谢物PBPK模型，具有效率和准确性的最佳平衡这种方法最早是在1973年引入的[9]，用于研究由具有众多速率系数的微分方程构成的耦合反应系统的灵敏度（经典的FAST方法仅具有一阶效应）。然后，它演变为估计SobolFAST算法通过“搜索曲线”创建输入参数的多维空间。它在计算模型参数的影响方面比其他全局SA方法（如基于Monte Carlo的方法）更有效。2. 软件描述我们的软件包可以研究具有多变量模型输出的PBK模型中的参数不确定性和灵敏度。此外，它与R或外部C代码中可用的通用ODE求解器以及用于进行复杂PBK建模的GNU MCSim软件（包括MCMC）无缝集成pksensi的分析框架不仅包括全局SA，还包括支持建模过程的不确定性分析和诊断工具2.1. 安装该软件包可以与RdeSolve软件包连接，后者包含一个综合积分器来解决ODE问题。GNU MCSim可以按照www.example.com上的GNU MCSim手册中的说明进行安装，或者使用pksensi中的内置函数进行安装https://www.gnu.org/software/mcsim/mcsim.html2.2. 参数矩阵生成我们采用了eFAST，广泛使用的方差分析（ANOVA）类似的全球SA方法，在这个包中的数值模型中的多维参数的方差分解为部分方差。大多数可用的eFAST功能，例如生成参数空间和设置采样频率的算法，都来自R灵敏度软件包[17]。eFAST方法通过对函数进行类似方差分析的分解来计算灵敏度指数（SI）。这种类型的全局SA方法和相应的SI也被称为Sobol由于使用较小样本量时估计的SI不稳定（高变化），因此我们采用随机相移方法从参数空间的随机起点进行重复采样，考虑在y f（xi）下建立的模型。抽样方案可以描述为，我们以前的工作也发现了一些有效的可视化AP-X1 1可以用来区分“有影响力的”和“非影响性”参数。这种方法解决了PBPK模型参数的常见可识别性问题，这些问题增加了贝叶斯分析的计算负担[13]，从而降低了计算效率。我们还开发了在决策中传达参数敏感性大多数SA工具都是使用科学计算软件构建的，如Matlab[14]，Python [15]和R [16]。R软件被广泛使用，包括许多开发的软件包，如灵敏度[17]，其中包括进行本地和全球SA的工具。灵敏度软件包包括一些函数，通过使用不同的计算算法来生成参数序列。同时，为集成外部建模结果提供了一种可行的方法。这种计算方法可以有效地解决纯R环境下数值解的计算负担其 他 几 个 R 软 件 包 包 括 SA 工 具 ， 可 以 在 Comprehensive RArchive Network（CRAN）中免费获得。尽管这些工具提供了执行SA的各种方法，但它们各自个别地具有应用于PBK建模的限制为i=2+πarcsin（sin（ωis+πi））其中xi是第i个参数的标称值，ωi是给出频率集合的向量，每个参数一个频率，并且π是范围从0到2π的随机相移系数。对于给定的样本量，SA应在具有随机相移的每次重复中给出类似的SI，以便被认为是稳健的或“收敛的”。对于大多数然而，由于其构造方式，eFAST估计值始终为正为了研究SI的收敛性，我们采用了以前提出的该方法通过计算所有时间点和输出变量的所有参数的计算SI的95%置信区间宽度来定量评估收敛性。通过随机相移过程，重复的SI值可以独立地重复采样。因此，它们可以用来研究SI的收敛性，CIi，t=max（SIub−SIlb）例如，SA中的估计索引在使用i，ti，t样本量不足然而，没有软件包中SIub和SIub是SI的上限和下限，i， t i， t提供评估趋同性的工具。此外，现有的eFAST功能不包括生成评估收敛所需的随机采样曲线的功能。我们在这里提出了一个R包，称为pksensi，其目的是使SA更容易获得和重现的药理学和毒理学研究。在时间（t）的特定模型输出（i）。CIi，t是相应的收敛指数。在pksensi中，我们将随机相移与eFAST集成对整个建模过程进行灵敏度测试。函数可以根据给定的样本大小和概率N.- H. 谢湾，巴西-地 Reisfeld和W.A.Chiu软件X 12（2020）1006093++Fig. 1. 在贝叶斯MCMC模型校准过程中执行不确定性和灵敏度分析的拟议工作流程分布（例如，平均值/标准差对数正态分布的平均对数/标准对数，以及均匀分布的最小值/最大值）。模型评估的数量等于样本量乘以模型参数的数量2.3. 建模pksensi还提供了一个有用的函数，可以使用分析和数值方法求解PBK模型中的ODE。通过与R deSolve包[19]链接，可以在纯R编程环境下计算然而，实现高计算速度的最因此，pksensi包含一个函数，用于在Windows上编译和创建动态链接库（.dll），在类Unix系统上编译和创建共享对象（.so）（例如，Linux和MacOS）。这个编译后的程序可以通过R中的内置函数加载和执行。但是，Windows用户需要安装Rtools才能使用GNU GCC编译器编译源代码。示例模型的C实现可以在GNU MCSim [8]中找到pksensi还可以与GNU MCSim链接以创建模型程序，用于求解每个方程组。2.4. 不确定性分析pksensi中不确定性和SA的一般工作流程如图所示。1.一、不确定性分析是SA中一个重要的有时，很难有信息数据来建立参数分布。参数分布的范围过宽可能会导致极值的数值误差。另一方面，过窄的范围可能导致预测无法覆盖校准数据中观察到的范围。在数据驱动的建模过程中，使用蒙特卡洛进行不确定性分析是解决这一问题的适当方法。在将数据拟合到模型之前，可以将Monte Carlo应用于给定的模型参数分布，以检查实验数据是否在模型输出分布的“合理”范围内pksensi还有一个函数，可以链接到GNU MCSim进行Monte Carlo分析，该分析内置于GNU MCSim中，然后检查模拟结果。此外，用户还可以检查由eFAST的参数矩阵生成的模型输出。2.5. 敏感性分析通过pksensi中的绘图功能，可以同时绘制所有参数的第一效应和总效应的时间相关灵敏度测量值。如果计算的敏感性指数在相同样本量下的重复试验中显示不同的值，则表明样本数量不足，会产生不可靠的结果必须使用足够的样本量，以防止参数固定中的错误判断。该软件包包括检查模型参数的收敛性和灵敏度的功能，提供了一种评估灵敏度测量的鲁棒性的方法。我们还在这个包中开发了一个基于“截止”的方法来区分“有影响”和“无影响”参数。最后，该包还提供了可视化工具，用于有效的调查和决策支持结果的交流3. 说明性实例我们展示了两个示例模型来说明在PBK建模中使用pksensi。在第一个例子中，我们使用了一个通用的，一室模型从httk包[20]来解释整个工作流程和细节在模拟过程中。用户可以参考CRAN中 的 在 线 小 插 曲 和 用 户 手 册 （ https ： //cran.r-project.org/web/packages/pksensi/index.html ） 或 pksensi 网 站（https://nanhung.rbind.io/pksensi/）获取详细的使用示例。在附录和图中给出了R码的部分例子及其结果。2，其在该过程中包括以下基本步骤。1. 构建模型2. 定义初始条件、输出时间步长和变量3. 集参数分布4. 进行模拟5. 不确定性分析6. 检查并可视化SA结果第一步是准备和构建模型。模型代码可以基于RdeSolve和GNUMCSim格式编写，为用户提供灵活性。 然后，必须设置模型变量及其初始条件。此外，本部分还必须提供模拟时间点下一步是识别将在全局SA中检查的模型参数并创建参数矩阵。在参数选择之后，我们必须为每个参数分配一个概率分布，以便进行评估，例如具有指定的最小值和最大值的均匀分布。然后，用户可以使用rfast99（）进行上述赋值（例如，参数名称、样品编号、参数分布及其相应的性质）来创建参数矩阵。最 后 一 步 是 使 用 以 下 公 式 求 解 模 型 中 的 常 微 分 方 程 ：solve_fun（）函数并可视化结果。一室模型的结果表明，体重（BW）模拟中对血药浓度无影响的参数（见附录）。此外，使用来自 GNU MCSim 的模型代码和求解引擎（solve_mcsim（））可以提供最快的计算速度。第二个实施例使用对乙酰氨基酚（APAP）-PBPK模型，一种具有多个隔室的复杂模型，以证明我们先前出版物[7]的重现性和pksensi的能力。APAP-PBPK模型描述了母体APAP及其两种代谢产物（葡萄糖醛酸AG和硫酸盐AS）的药代动力学。我们将全面的全局SA工作流程应用于具有21个模型参数的原始发布模型[21]。图3示出了针对全局SA方法计算的时间相关的总SI。热图提供了一种可视化估计结果的替代方法。一些参数（例如，在整个调查时间间隔内，所有三个输出变量均无影响N.- H. 谢湾，巴西-地 Reisfeld和W.A.Chiu软件X 12（2020）1006094图二. 实施例R脚本及其在单室药代动力学建模中应用pksensi的 可以通过check（）函数直接在R控制台查看结果，通过plot（）函数输出可视化结果图3.第三章。 根据A P A P - P B P K 模 型 中 的 21个 物质特异性参数计 算 的 时 间 依 赖 性总SI的热图表示。4. 影响和结论计算建模已成为现代科学研究中支持假设、实验和发现的基本技术。不确定性和SA是建模者进行模型开发和评估的重要工具。CRAN存储库和其他软件平台上的几个软件包已经被构建为使用各种算法执行全局SA，但没有一个包括与基于ODE的模型（如PBK模型）的集成。我们发起pksensi是为了在PBK建模中促进全面高效的全球SA工作流程，填补开源建模社区的关键空白。它提供了一个简单的应用程序来研究参数对模拟时间点和变量的模型输出的影响。我们采用了eFAST方法，这是一种常见的基于方差的SA，我们发现它在效率和准确性方面取得了最佳平衡。此外，该软件包还包括评估SA结果的收敛性的能力，这在大多数全球SA研究和软件包中很少涉及。我们还开发函数来可视化输出结果，并帮助区分“有影响”的参数和“无影响”的参数，这些参数可以在模型校准中固定。我们的软件包使全局SA更容易在PBK建模。pksensi的当前版本包括GNU MCSim和deSolve包作为主要的ODE求解器。我们选择将SA与GNU MCSim集成，因为后者是一个强大的开源软件包，用于基于生物模型的动态模拟，主要用于PBK研究和概率方法[8]。与deSolve软件包相比，GNU MCSim可以提供更好的速度和效率来求解PBK模型中的常微分方程。虽然pksensi提供了安装和链接GNU MCSim的基本特性和功能，在PBK建模中进行不确定性和SA，但对于不熟悉GNU MCSim的用户来说，从模型开发，调试和测试，此工作流程的整体学习曲线可能很陡峭。因此，通用包，如deSolve，可以成为R用户的替代选择。我们的软件包的灵活性可以帮助建模者量化模型参数对模型输出的影响此外，pksensi还拥有N.- H. 谢湾，巴西-地 Reisfeld和W.A.Chiu软件X 12（2020）1006095library（pksensi）library（httpk）#Example：One-compartmentpbtkmdel------------------------------------1. Construct1-cptpbtkmododelforordeSolve包pbtk1cpt<-function（t，state，parameters）{with（as. list（c（state，parameters）），{dAgutlumment=-kgutabs*AgutlummentdAcompartment=kgutabs*Agutlumment-ke*一个合作伙伴dAmeabolized=ke*AcmparmetCcmparmet=Acmparmet/vdist*BW;list（c（dAgutlumment，dAcompartment，dAmetabolized），“Ccompartment“=Ccompartment）}）}#Step2. 定义定义、输出时间表和变量有可能应用于最近研究中提到的交互式PBPK建模平台[22]。值得注意的是，当前的全局SA结果取决于特定模型和给定的模型参数范围。其他因素，如检查的参数数量、参数分布/范围以及PBK模型中的给定具有高参数维数的复杂模型对计算机的要求也较高随着解决这些复杂性的方法的发展，我们希望将它们纳入pksensi。总体而言，pksensi提供了一套全面的特性和功能，用于执行PBK模型的全局SA，并可以为模型开发和校准中的决策制定创建强大且可重复的结果。虽然这里使用的pksensi主要用于PBK建模，但它也可以应用于其他基于ODE的动态模型，以研究模型输出对输入参数的灵敏度。CRediT作者贡献声明谢南宏：概念化，软件，写作-初稿。Brad Reisfeld：资金获取，项目管理。魏雪湾邱：督导，方法论，写作-评论&编辑.竞合利益作者声明，他们没有已知的竞争性财务利益或个人关系，可能会影响本文报告的工作致谢这项工作的部分资金来自美国食品药品监督管理局（FDA）的1U01FD005838 号 拨 款 和 美 国 国 家 环 境 健 康 科 学 研 究 所 的 P42ES027704号拨款。这篇文章反映了作者的观点， 不代表FDA的观点或 政策 。 我们 感 谢德 克 萨斯 A M 大 学的 Yasha Hartberg 博 士和Barbara Gastel 博 士 、 Certara 的 Frederic Bois 博 士 和 美 国 FDA 的Eleftheria Tsakalozou博士审查手稿和提供咨询。阑尾 示例代码其中，State-c（Agutlummet=10，Acmpartmet=0，Ametabolized=0）t<-seq（from=0. 01，to=24。01，by=1）输出s-c（“Ccompartment“）第三章. Gennerrateparametermatr a rix##3号。1 .一、 Extractparamete r amevalefmtkpackagepars1cmp<-（parameterize_1cmp（chem.name=“acetaminophen“））##3号。二、SETPARAMERDISTRIBUT IONSq-c（“qunif“，“qunif“，“qunif“，“qnorm“）Q. arg-list（list（ min=pars1comp$Vdist/2，max=pars1comp$Vdist*2），list（ min=pars1comp$kelim/2，max=pars1comp$kelim*2），list（ min=pars1comp$kgutabs/2，max=pars1comp$kgutabs*2），list（ mean=pars1comp$BW，sd=5））##3号。3 .第三章。Createparametermatixset.seed（1234）params-c（“vdist“，“ke“，“kgutabs“，“BW“）X<-射线衍射（XRD）分析表明，X-射线衍射（XRD）分析表明，arg=q。arg，重复= 10）第四章.条件是imulat ion（将kefewmintes）out<-solve_fun（x，t ime=t，fu nc=pbtk1cpt，initState=initState，outnames=outputts）第5页。你知道吗？pksim（outt）#Usetocomparewith“real“data（ifany）第6页。敏感性分析结果的确认和验证plot（out） #Visualizeresultcheck（out）#Printresultoconsolee引用[1]Ferretti F，Saltelli A，Tarantola S.近十年来敏感性分析实践的发展趋势科学总环境2016;568：666-70。http://dx.doi.org/10的网站。1016/j.scitotenv.2016.02.133.[2] Loizou GD，McNally K，Jones K，Cocker J.整体敏感性分析在开发人体间二甲苯和乙醇联合暴露的生理药代动力学模型中的应用Front Pharmacol2015;6：135.http://dx.doi.org/10.3389/fphar.2015.00135网站。[3]McNally K、Cotton R、Cocker J、Jones K、Bartels M、Rick D等人，使用PBPK建模、贝叶斯推断和马尔可夫链蒙特卡罗模拟从人体生物监测数据中重建间二甲苯暴露。J Toxicol 2012. http://dx.doi.org/10.1155/2012/760281网站。[4]McNally K，Cotton R，Loizou GD. PBPK模型全局敏感性分析的工作流程。前沿药理学2011;2：31。http://dx.doi.org/10.3389/fphar.2011.00031。[5]张晓宇，特莱姆M，莱斯科L，施密特S. Sobol敏感性分析：指导系统药理学模型开发和评价的工具。CPT Pharmacomet Syst Pharmacol 2015;4（2）：69-79。http://dx.doi.org/10的网站。1002/psp4.6。[6] Chu Y，Hahn J.通过准线性化使用全局灵敏度分析的定量最优实验设计。IndEngChemRes2010;49（17）：7782-94.http://dx.doi.org/10.1021/ie9009827网站。[7] 谢N-H，Reisfeld B，Bois FY，Chiu WA.应用全局灵敏度分析工作流程提高基于生理学的药代动力学建模的计算效率前沿药理学2018;9：588。得双曲正切值. doi.org/10.3389/fphar.2018.00588网站。[8] 博伊斯·菲。GNU MCSim：SBML编码系统生物学模型的贝叶斯统计推断 。 Bioinformatics 2009;25 （ 11 ） ： 1453-4. http://dx.doi 的 网 站 。org/10.1093/bioinformatics/btp162.[9]Cukier R，Fortuin C，Shenyke，Petschek A，Schaibly J.耦合反应系统对速率系数不确定性的敏感性研究。我的理论。J Chem Phys 1973;59（8）：3873-8. http://dx.doi.org/10.1063/1的网站。1680571。[10] Saltelli A，Tarantola S，Chan K-S.模型输出全局灵敏度分析的一种定量模型无关方法。Technometrics1999;41（1）：39-56.http://dx.doi.org/10.2307/1270993网站。[11] Jansen MJ 模型输出的方差分析设计。 COM-put PhysComm 1999;117（1-2）：35-43. http://dx.doi.org/10.1016/S0010-4655（98）00154-4.N.- H. 谢湾，巴西-地 Reisfeld和W.A.Chiu软件X 12（2020）1006096[12]欧文AB。更好地估计小sobol'敏感性指数。ACM跨模型计算模拟2013;23（2）：11. http://dx.doi.org/10.1145/2457459的网站。2457460。[13]Garcia RI，Ibrahim JG，Wambaugh JF，Kenyon EM，Setzer RW. PBPK模型的可识别性及其在二甲基胂酸暴露中的应用。 药物动力学杂志，2015;42（6）：591-609。http://dx.doi.org/10的网站。1007/s10928-015-9424-2。[14]Pianosi F，Sarrazin F，Wagener T. Matlab工具箱，用于全局灵敏度分析。环境 模 型 软 件 2015;70 ： 80-5.http://dx.doi.org/10.1016/j 的 网 站 。envsoft.2015.04.009网站。[15] 赫尔曼J，亚瑟W。SALib：用于敏感性分析的开源Python库。J OpenSource Softw 2017;2 （ 9 ） . http://dx.doi.org/10.21105/joss 的 网 站 。00097。[16]Looss B，Saltelli A.敏感性分析导论。不确定性量化手册。Springer; 2017，p.1103-22. http://dx.doi.org/10.1007/978-3-319-12385-1_31。[17]Looss B，Janon A，Pujol G，Boumhaout K，DaVeiga S，Delage T，et al.Sensitivity ： Global sensitivity analysis of model outputs. 2018 ， R 包 版 本1.15.1，URLhttps://CRAN.R-project.org/package=sensitivity。[18] Sarrazin F，Pianosi F，Wagener T.环境模型的全局敏感性分析：收敛性和有效性。EnvironModelSoftw2016;79：135-52.http://dx.doi.org/10.1016/j.envsoft.2016.02.005网站。[19] Soetaert K ， Petzoldt T ， Setzer RW. 在 R 中 求 解 微 分 方 程 ： PackagedeSolve。统计软件杂志2010;33（9）：1-25.http://dx.doi.org/10.18637/jss.v033.i09.[20]Pearce R，Setzer RW，Strope C，Sipes N，Wambaugh J. Httk：R packagefor high-throughput toxicoketics. J Stat Softw 2017;79（4）：1-26. 网址：//dx.doi.org/10.18637/jss.v079.i04网站。[21] 作者：J.对乙酰氨基酚及其主要代谢物在人体内的药理学建模，使用贝叶斯群 体 方 法 。 Eur J Drug Metab Pharmacokinet 2016;41 ： 267-80 。http://dx.doi.org/10.1007/s13318-015-0253-x网站。[22]Li M，Cheng Y-H，Chittenden JT，Baynes RE，Tell LA，Davis JL，et al.Integration of food animal residue avoidance databank（FARAD）empiricalmethods for drug withdrawal interval determination with a mechanism- ticpopulation-based interactive physiological based pharmacokinetics（iPBK）modeling platform：example for flunixin meglumine administration.毒理学文献2019;1-16。http://dx.doi.org/10.1007/s00204-019-02464-z.