基于模拟退火算法的适当集总法简化茶碱生理药动学模型*

梁露花,黄叶真,普思静,郑浩东,刘欣然,项荣武

(1.沈阳药科大学医疗器械学院生物医药信息教研室,沈阳 110016;2.辽宁省医药大数据与人工智能工程技术研究中心,沈阳 110016;3.沈阳药科大学制药工程学院,沈阳 110016;4.沈阳药科大学药学院,沈阳 110016)

全身生理药动学模型 (whole-body pharmacokinetics model,WBPBPK)是高维复杂的,通常是非线性的,表示为一个常微分方程(ordinary differential equation,ODE)系统,其中每个常微分方程表示一个组织或器官腔室的状态[1]。在临床研究中,通常希望建立的WBPBPK能吻合动物或临床数据,进而改善药物相关的特异性参数,但是,由于伦理限制,通常不能直接测定组织浓度,许多参数(例如组织-血液分配系数)不能仅仅通过血浆数据得到[2]。因此,根据模型参数先验的信息,在高维系统中,统计模型可能无法被识别,将WBPBPK拟合到临床数据通常会产生数值上不稳定的分析。故一个常见的策略是,将全身生理药动学系统进行降阶,把高维复杂的模型简化为更少的状态,从而减少生理药动学模型中ODE和参数的数量,实现对WBPBPK的简化,以稳定临床数据[3]。目前常用简化生理药动学模型的4种方法来搜索最优集总模型,包括完全枚举法(参考方法)、非自适应随机搜索(non-adaptive random search,NARS)、碎石图加非自适应随机搜索(scree plot plus NARS)、模拟退火算法(simulated annealing algorithm,SA)[4]。笔者通过已有的茶碱成人WBPBPK,运用SA对全身模型进行适当集总简化。对简化模型和原始模型预测的血药浓度-时间分布进行比较,以保留对系统的机制解释。

1 材料与方法

1.1数据来源 模型数据来源于KRAUSS等[5]提供的非线性混合效应模型用户指南的第五部分12例30岁的健康志愿者的身高、体质量等基本参数,将这些基本参数输入到PK-sim数据库,得到志愿者的系统特定的解剖学参数。在开放系统药理学(www.open-systems-pharmacology.org)中找到茶碱的吸收、分布、代谢和排泄(absorption distribution, metabolism, and excretion,ADME)参数,放入PK-sim数据库(版本10.0,是Open Systems药理学套件的一部分)中查看,得到药物特定的参数。按照12例志愿者不同的体质量对其分别给予不同剂量,为了计算机模型的简单性,文献[5]给出计算其给药剂量的平均值(单位:mg),这并不会影响最终集总模型的通用性。

1.2茶碱的原始模型 从BJÖRKMAN 等[6]已发表的文献中确定茶碱的WBPBPK。茶碱的PBPK模型在MATLAB®(版本 R2016)中被编码为常微分方程,使用MATLAB矩阵指数求解其常微分方程。茶碱的WBPBPK是由14个隔室组成的PBPK模型[7]。为了方便将茶碱的原始PBPK模型转化为矩阵的形式,首先就是列出原始微速率常数矩阵K。把每个隔室进行单独的编号,其中静脉血为输出状态,做为非集总状态,不对其进行编号,剩余的13个隔室对其进行编号,肺为1号、动脉血为2号、大脑为3号、心脏为4号、肾脏为5号、肝脏为6号、肌肉为7号、皮肤为8号、脂肪为9号、酮体为10号、肠道(小肠)为11号、脾脏和胰脏合并为12号、胃为13号。

1.3模拟退火算法适当集总以简化模型 茶碱的WBPBPK为14个隔室的高维复杂模型,为了保证原始系统动力学不会改变,选用一种对于线性或线性化模型来说是一种潜在的更灵活的技术——集总,其中状态的合并是根据模型的性能特征来确定的,集总技术参考文献[7]。在本文中SA被用来在寻找集总矩阵时最小化相对差异的绝对值(absolute value of the relative difference,ARD%),使用SA对M矩阵进行最优搜索[8]。SA的变量设置为:初始温度为104℃,负责在搜索过程中监管接受新的能量状态的概率;温度下降速率为每周期0.999,即第一次循环之后温度从104℃下降到104×0.999 ℃,线性冷却方案继续每个循环;最大迭代次数为2×104。直到寻找到最优M矩阵停止。

1.4WBPBPK和集总的生理药动学模型的评估 对集总前茶碱的WBPBPK和利用SA集总得到的适当集总简化的茶碱生理药动模型进行评估。简化模型和原始模型之间茶碱浓度-时间曲线下面积的ARD%<0.002%和模拟动脉浓度-时间的分布图中原始模型和集总模型的曲线拟合程度,若二者比较良好则认为简化模型可以较好地替代原始PBPK模型,并且可以应用到临床[8]。

1.5软件 PK-sim(版本10.0,是Open Systems药理学套件的一部分,www.open-systems-pharmacology.org)用于查找系统特定解剖学参数和药物的特异性参数。MATLAB®(版本 R2016)用于SA进行茶碱PBPK模型的简化,寻找最优集总矩阵M和进行生理药动学建模。

2 结果

2.1SA适当集总算法简化的生理药动学模型 将SA适当集总算法代码,分别保存为m文件,输入12例志愿者生理解剖学参数和茶碱的相关参数进行运行,最终得到最优集总M矩阵为集总成7态的全身药动学模型,最终简化模型的集总矩阵具体如下:

将最优集总M矩阵转化为具体的隔室,得到茶碱集总后的简化PBPK模型集总成7态时,大多数集总状态下的预测与原始系统预测相当一致,动脉血液单独为一态;第2个伪状态是肌肉、酮体、肠道;大脑、合并的脾脏和胰脏、胃合并为一个伪状态;肾脏程单独的一态;肝脏为单独的一态;肺为单独的一态;第7个伪状态是心脏、皮肤、脂肪合并为一态。第二伪状态是肌肉、酮体、肠道的集总,第三伪状态是大脑和内脏器官的集总。然而,预测这2个集总隔室中的药物量对茶碱没有特别的临床意义。

2.2适当集总简化PBPK模型的评估 运用MATLAB进行SA搜索集总得到一个7态的M最优集总矩阵,此矩阵结果输出的血药浓度-时间曲线下面积为332.47 ng·h·mL-1,ARD%为0.000 190 82%。在对茶碱PBPK模型简化之前人为规定集总模型和原始模型之间茶碱浓度-时间曲线下面积的最大差异为0.002%,即ARD%<0.002%,表示原始模型和集总模型血药浓度-时间偏差不大,此时认为简化后的茶碱模型有意义。额外运行不同初始温度(即分别为103和105℃)在集总数量方面给出了相同的结果状态。图1温度变化趋势和经过迭代后对最优ARD%解的收敛显示了模拟退火算法对数尺度的温度变化的趋势,和经过迭代后对最优ARD%解的收敛。随着温度指数的下降即退火,在连续迭代中形成ARD%的卷积曲面,在经过20 000次数的迭代后,ARD%达到了稳定,表明迭代收敛到最佳状态。上坡跳跃随着温度的降低而减小,最终使搜索收敛到最优解。

图1 温度变化趋势(左)和经过迭代后对最优ARD%解的收敛(右)

2.3WBPBPK和集总的生理药动学模型的评估 在MATLAB中使用模拟退火算法得到最优M矩阵之后,继续在MATLAB中调用茶碱的集总模型参数,模拟动脉血药浓度-时间分布图(图2),并与原始的茶碱PBPK模型的模拟动脉血药浓度-时间分布图进行比较,可以看出集总后隔室与原始系统的偏差几乎不会影响血浆中药物浓度的时间过程,认为集总模型与原始模型相关生理意义继续保持,此次集总的茶碱简化成人模型有生理学意义。

图2 原始模型和7态集总模型模拟动脉浓度-时间分布

3 讨论

本研究通过查阅茶碱的相关文献得到关于该药物的相关参数、临床实验数据以及茶碱的WBPBPK模型,运用SA适当集总算法对PBPK模型进行简化,对简化后的模型进行评估。尽管简化后模型与原始模型比较模拟浓度-时间曲线偏差较小。但研究中仍然存在局限性,本研究所使用的参数为12例30岁男性欧洲人生理参数,只考虑了肝、肾的清除率忽略了肺、皮肤等其他组织器官的清除率。总体看来,本文利用SA适当集总算法对茶碱的成人模型的简化相对有意义,经评估与原始模型预测血药浓度的差异性较小,可以用于表征潜在生理和ADME相关参数以及药物特定参数的个体间变异性。后续还需要更多的药物例子来证明这种简化方法在药物开发中的优势。