2023年6月21日发(作者:)
第十三章曲线拟合以及常用的药代动力学软件
第十三章 药动学数据的曲线拟合以及常用软件
计算药动学参数,是药代动力学进一步应用的基础。如何测定有关的动力学参数呢?常用的方法是:首先在用药后的若干不同时间,采取血样(或尿样),测定其血药浓度值或尿中药量(这些数值称为实测值或观察值,用C i 表示),这样就有了药物浓度经时曲线数据;然后,依据半对数坐标图,选定一种模型方程
(是时间t 的曲线函数)计算理论估算值(用i C )表示),按照观察值和理论估算值之
差的平方和(即残差平方和)或加权残差平方和(均用Re 表示)最小的原则,采用适当的算法,求出有关的动力学参数。这种方法,在数学上称为曲线拟合(fitting a curve)。由于所采用的线性药代动力学的模型方程是多指数项之和的函数形式,并且是所含动力学参数的非线性函数,所以这种曲线拟合方法称为非线性最小二乘法。
一、最小二乘法的一般原理
设y 是变量x 的函数,含有m 个待定参数a 1,a 2,…,a m 。记为
y =f (x ;a 1,a 2,…,a m ) 若对x 和y 作n 次观察,测得观察值(x 1,y 1),(x 2,y 2),…,(x n ,y n )。根据这样一组二维数据,即平面上的若干点,要求确定这个一元函数y =f (x ;a 1,a
2,…,a m );(i=1,2,…,n),即一条曲线,使这些点与曲线总体来说尽量接近。并使y
的观察值y i 与理论估算值=i y )f (x i ;a 1,a 2,…,a m )
;(i=1,2,…,m)的误差平方和,即残差平方和21
()n e i i i y y ==?∑)2121((,,,,)n i i m i y f x a a a ==?…∑R )i 取得最小值,
或者加权残差平方和21()/n e i i i R y y w ==?∑)2121
((,,,,))/n i i m i i y f x a a a w ==?…∑ 取得
最小值。其中w i 称为权重系数,在后面的段落会详细讲解。这时Re 有时候也称为目标函数。这就是数据拟合成曲线的思想,简称为曲线拟合。曲线拟合的目的是根据实验获得的数据去建立因变量和自变量之间有效的函数关系,这个函数关系对于药动学来讲就是通过房室模型推导出来的药时曲线公式。根据观察值求出待 311
定参数(因而也就确定了曲线y =f (x ;a 1,a 2,…,a m ))的问题。我们称f (x ;a 1,a 2,…,a m )为拟合函数。特别地,当f 为x 的线性函数时,则称为直线拟合(或直线回归)。
如下图所示,Re 是点(t i ,C i )与曲线)(x f y =的距离,曲线拟合实际含义就是寻求一个函数,使在某种准则下与所有数据点最为接近,即曲线拟合得最好。最小二乘法准则就是使所有离散点到曲线的距离的平方和最小。例
如对于一房室静注模型函数)(x f y =)(x f i C )=C 0 e -kt i ,就是确定待定系数K 和C 0,使当所有时
间数据点t i ,i =0,1,……,n ,代入函数i C )=C 0 e -kt i 后,按公式计算使求算的结果最小,此时待定系数的值K 和C ∑=?=n i i i
e C C R 1
2
)()0的值就是拟合所求的结果。 度g m l
图13-1 药时曲线拟合示意图
由高等数学中的极值原理知,待定参数a 1,a 2,…,a m 应满足下列方程组(一般称为正规方程)
0Re =??j
a (j=1,2,…,n) (14-1) 这是含有m 个未知数的m 个方程,解这一方程组可求得a 1,a 2,…,a m ,从而确定了拟合函数。当f 是参数的线性函数时,上述正规方程为参数的线性代数方程组,这种情况称为线性最小二乘法。当f 是参数的非线性函数时,上述正规方程则为参数的非线性方程组,这种情况称为非线性最小二乘法。
曲线拟合寻求一个函数,使在某种准则下与所有观察值最为接近,这种搜索最小值的算法目前常采用高斯—牛顿迭代法、单纯形法等。在具体)(x f y =)(x f 312 计算时,有时候采用简化方法—对数回归法或残数法,将非线性方程转化为线性方程,计算虽然比较简单,但计算的误差往往比较大,同时手工计算比较费时费力。目前经常采用高斯—牛顿迭代法、单纯形法等算法编制成计算机程序,当数据比较符合理论情况时,能够比简化计算方法计算出更精确、合理的动力学参数。
二、非线性最小二乘法算法的比较
1、经典高斯—牛顿迭代法以及改进方法—哈特莱方法(Levenberg-Hartley法)、阻
尼最小二乘法
高斯—牛顿迭代法将目标函数(Re)在待定药动学参数初值(a1,a2,…,a m)附近的微分方程用泰勒级数的一次项展开,得正规方程组,采用列主元高斯消元法解出该方程组,就可计算得理论上使得目标函数(Re)最小的最佳药动学参数。
该方法的优点:在某种程度上,按最佳梯度方向搜索,效果较好,虽对初值有一定的依赖性,但依赖程度远远低于其他类型的方法,如单纯形法。所以开始运算时往往收敛较快,运算时间短,这是此方法的突出优点。
缺点:由于泰勒展开中丢弃了高次项等等的原因,使得此法往往不能精确收敛,甚至会引起发散,不能求出解。基于这个缺点,对于经典高斯—牛顿迭代法进行了种种的改进,如哈特莱方法、阻尼最小二乘法等等,可以有效地改进拟合发散的缺点,但还是不能完全避免。
2、单纯形法
本方法是一种多维搜索的直接方法,不需要计算目标函数的导数。通过对n 维空间的n+1个点(它们构成一个初始单纯形)上的函数值进行比较,去掉其中函数值最大的点,代之以新的点,从而构成一个新的单纯形。这样,通过多次迭代逐步逼近极小点。
单纯形法的优点是:原理简单,避免了求导运算以及解正规方程组等步骤,从而不会有经典高斯—牛顿迭代法以及改进方法固有的缺点,基本不会发散。缺点是:收敛速度慢,初值依赖性较大,有陷入局部最小值的弊端,一般不单独使用,有人主张用高斯—牛顿迭代法求出大致的结果,再采用单纯形法作局部区域的精确搜索寻优。
三、估算药代动力学参数中的若干问题
313
1、如何选择模型
在曲线拟合时,需要事先选择合适的药动学模型方程再来进行拟合。也就是说同一组数据可以选择不同的药动学模型方程进行拟合(亦即拟合函数),那么如何选择和确定适当的模型呢?一种简单和直观的方法,就是根据实测的血药浓度的对数值对时间作图(称为对数浓度一时间散点图),作粗略的直观的分析。当散点图的分布比较有规则地反映出某一种单指数或多指数态势时,可确定出相应的房室数。但是,在许多情形下,散点图的分布往往似有多指数态势,或者多指数态势不明显,这时就要从中加以挑选。
选择房室数的方法,最常采用赤池信息判据最小准则AIC(Akaike’s Information Criterion)。
假设观测数据点数为n,拟合函数中所含待定参数的个数为m,拟合所得残差平方和(或加权残差平方和)为Re,
∑=?
=
n
i
i i
e
c c
R
1
2
) ();则
AIC=n ln(Re)+2 m (14-2) 赤池提出按AIC值最小的准则来挑选模型。就是说,在几种预定的模型中,较佳的模型其AIC值最小。根据这一准则,我们就可以对几种不同的房室模型分别进行曲线拟合,算出各自的Re和AIC值,比较AIC值的大小,然后选AIC 值最小的房室,作为优选的模型。另外,根据这一准则,我们对同一模型,在用不同的计算方法作曲线拟合时,AIC值最小的那种方法,算得的参数较为精确;不过这时AIC值最小就相当于Re值最小,因此,只要比较Re值的大小即可。
从公式看,AIC值的大小主要和两个因素有关,一个是Re,一个是待定参数个数。不同的房室模型,其待定参数不一样。例如,一房室静注,待定参数为K 和Vd两个参数,对数浓度一时间散点图为直线;二房室静注,待定参数为A、B、α、β四个参数,对数浓度一时间散点图有一个拐点的两根直线。模型中房室个数越多,待定参数越多,拟合函数越复杂,曲线包含的拐点越多,拟合时Re值就会越低。如果仅仅采用Re作为模型选择依据,非常复杂的多房室模型其Re值往往会更低,产生所谓过拟合现象。AIC值兼顾了模型房室个数和Re两个因素,既使得Re值低,同时又限制了复杂函数的使用,增加了拟合函数的实用
314
价值。从实用角度看,房室模型的房室个数不易超过3个。
在选取模型时要特别注意,房室模型中房室的划分具有相对性。当实验数据比较准确和充分时,可以将药物在体内分布的较小的速度差异区分开来,从而可以将体内分成更多的房室;但是当实验数据比较少或者误差较大时,药物在体内的速度差异就无法区分,只能将机体分成较少的房室或仅单一的中央室。由于上述原因,有可能同一种药物,在不同的文献报道中有不同的房室模型,要理解和容忍这种房室划分的相对性。
2、权重系数
曲线拟合采用的目标函数是Re ,Re 值越小,曲线拟合效果越好。但是,在实验观察值中,高、低血药浓度值相差较大时,Re 值的大小往往过分取决于高浓度
的观察值,而忽视了低浓度值的观测作用。例如,同样是||i i c
c )?=0.5的差值, 对于100ng/ml 相对误差很小,对于1ng/ml 相对误差就很大,所以不能把高低浓度的差值等同看待。这时候就需要采用权重系数的方法,即在加权的情况下,求
加权残差平方和的最小值问题。其中W i n
i i i e w c c R /)(12∑=?=)i 称为权重系数。
权重系数最常采用的是:
W i =1/C 2i
则加权残差平方和为
21Re ∑==n i i i i C C C ) (13-3)
即为相对误差平方和最小。
当W i =1,就是不进行权重。有时候W i =1/C i ,权重介于1和1/C 2之间。至于选择何种权重系数,目前仍然在讨论之中,没有统一的方法,需要根据不同的情况进行具体的分析。需要强调的是,不同的权重系数的Re 值和AIC 值没有可比性。
四、曲线拟合的影响因素
在应用计算机程序进行曲线拟合,估算药动学参数时,影响拟合效果的因素很多,主要有以下几个因素:
1.实验设计问题
药时曲线的采样点数和浓度检测的准确性显著地影响着模型的识别和拟合 315
的质量,因而采样时间点和采样点个数的设计是成功进行拟合的第一步。原则上以较少的采样,尽可能获取血药浓度经时曲线的变化形态信息,如峰时间、峰浓度、拐点、变化趋势等。在曲线变化大的地方,应多采样,在曲线变化小的范围内,则少采样。但是,由于个体差异以及其他因素的存在,可能药时曲线的拐点等位置不一样,所以只能尽量做到符合以上原则。另外经验证明,对于m个参数的模型,采样点个数n至少为2m+1,例如二房室口服有5个待定参数,测11个点以上比较可靠。n≤m是不可行的,否则残差的自由度为n-m-1为-1,不符合统计学要求。事实上,国家药品审评机构的指导原则中明确表明药时曲线采样点个数不得少于11个,一般我们进行药时曲线的点数设计时都会超过11个。
2.样品浓度测定的准确性
目前国家药品审评机构对于体内生物样品检测有比较明确的指导原则,其中明文规定,对于生物样品检测,高、中浓度的检测相对误差必须在±15%以内,低浓度必须在20%以内。由于实验误差的存在,曲线拟合时权重系数的选取必须考虑。如果低浓度样品测试比较准确,就要采取权重拟合,反之普通拟合。对于半衰期等参数,更多地取决于药时曲线末端相的下降速率。样品检测的准确性显著地影响到这类参数的拟合。
3.收敛精度的选取
收敛精度的选取,通常采用由大到小逐步修正的办法,即起初精度可以取得大些,视其计算结果是否满意(残差平方和或加权残差平方和,是否可明显减少;观察值与估算值之间的相关系数,是否可明显增大等),再决定是否提高精度(将精度值取小)。在继续拟合时,应选前一次结果作为初值。
4.初值选择和拟合取值范围问题
目前所采用的非线性最小二乘法的算法都面临一个初值如何选取的问题。初值选得好,迭代次数少,收敛速度快,所得结果也相对理想。初值选取一般首先参考文献上的资料或者以往的经验,也可以应用残数法的计算结果作为初值。但是无论初值如何选取,也还是存在初值依赖性的问题。
初值依赖性是指随着初值选取的不同,拟合所得药动学参数也会随之变化。下图是一个实际的利用单纯形法拟合药动学参数过程中目标函数Re值的分布图。
316
该例为口服一房室模型,有三个待定参数:Ka;Dose*F/V,Ke。在Ka=0.3时,罗列Dose*F/V和Ke两个参数在一定范围内的所有对应组合情况下的Re值;从图中可以看到颜色越深,Re值越小;低于某一限度的Re值区域为一个不规则的长条状。寻优过程中,当参数组合点到达最深的区域的边缘1就可以终止拟合。由图可见,随着起始位置不同,从三个方向开始拟合,所到达的目标值并不是归结为一个点,而是三个点,在Re值最小区域的边缘,这就存在药动学参数拟合的初值依赖性。
改变Ka值,可以绘出一系列图,通过多媒体动画技术,将这些图叠加在一起依次出现,依据Re值最小这个标准,就可以找出一个最佳的三个待定参数的组合范围,使得Re值最小。选取该组合范围的中点,就可以避免初值依赖性。但是随着待定参数数目的增加,待定参数的组合数呈几何级数增长,大大增加计算工作量,获得所需结果就仅仅成为理论上的可能了。
另外,由图1可见,参数Dose*F/V的取值在小于0没有实际意义的情况下,Re值依然很小,说明了数据拟合寻优过程中必须限定取值范围,否则拟合出来的数据没有实际意义。
图13-2 拟合中的初值依赖性
Re值变化率最小,并不是Re值最小。Re变化率最小和Re最小的区
1此处的最深区域为Re值最小,但是真正的拟合收敛条件应该为域虽然不一致,但是也同样是一个不规则区域。
317
初值依赖现象不仅仅在单纯形法模型嵌合中存在,其他任何寻优法中都存在,关键在于药物浓度-时间双指数或多指数曲线本身就存在两个或多个最优解。下面表1中的两组一房室模型参数能产生出如图2所示的完全一样的药时曲线。这两组参数的区别仅仅在于Ka和Ke的大小不同,表观分布容积不同。也就是说,如果Ka和Ke大小不能确定,同一组药时数据会有两组同样最优的药动学参数组合。一般情况下Ka>Ke,所以同一组药时数据最优解只有一组药动学参数组合。但是,对于缓控释制剂、特殊释放制剂,就不可能仅仅通过本身口服后的数据来区分吸收指数项和消除指数项,从而来进行吸收解析。这就是吸收解析的时候最好要同时进行静脉给药实验的原因——能够确切地知道药物的消除相的情况。但是,当药物静脉给药和口服给药消除情况不一致的时候,哪怕进行静脉注射也无法进行口服药物的吸收解析,例如口服给药存在肠肝循环,首过效应强等等的情况。
图13-3 多指数曲线存在两个最优解
100
F(%) 100
口服给药剂量(mg)973 973
0.52
Ka(h-1) 2.5
k(h-1) 0.52
2.5
34
V(l) 163
318
五、目前常用计算机药动学数据拟合程序
目前国内外药动学计算软件比较多。国外比较著名的软件有WinNonlin、Kinetica等,国内比较著名的软件有3p87(3p97)、PKBP—N1、BAPP、DAS等。
(一)WinNonlin 软件
WinNonlin为美国Pharsight公司产品,是目前国外应用最广泛的软件,界面友好,功能强大,与其他软、硬件有很好的兼容性。WinNonlin 有标准版、专业版、企业版等三个版本,标准版中包含了药动学与药效学分析的各种工具,专业版和企业版较标准版增加了几个模块,主要用于商业用途。
目前WinNonlin软件最新版本为4.1,其主要功能如下:1、房室模型分析(PK Fitting,Compartmental Analysis),包括19种模型库,支持用户自定义模型;寻优算法包括Nelder Mead simplex,Gauss-Newton (Levenberg Hartly),Gauss-Newton (Hartly)法,权重方式包括1/obs,1/obs*obs,1/calc,1/calc*calc,初值的设定包括了上下限;模型确认标准包括Akaike,Schwarz,CV,SEM等。2、非房室模型分析(Non Compartmental Analysis),包括各种给药途径,还支持末端相半衰期图解求算,支持参数个体和群体设定,可计算稳态数据参数。3、等效性检验,包括平均等效性、群体等效性、个体等效性(ABE、PBE、IBE)。4、PK-PD数据分析,支持PK、PD分步分析和PK-PD统一嵌合的大模型分析,但是尚不支持群体的PK-PD分析。5、提供了“工具箱”(toolbox)功能及帮助功能:非参数重叠法( nonparametric superposition),用来预测多剂量用药后达到稳态的血药浓度;半房室模型法(semicompartmental modeling),用来估算给定时间和血浆浓度的效应地点浓度。6、数据输入输出支持格式众多,包括ASCII、Excel、SAS?、Transport、Oracle、Watson、ODBC-Treiber等,有非常良好的兼容性;输出的图表能形象化地显示数据,可进行编辑修改。7、提供了一定的数据描述统计功能。
下面两幅图片分别为WinNonlin软件使用过程中模型选择和模型拟合初值输入。
319
图13-4 WinNonlin软件使用过程中模型选择对话框
图13-5 WinNonlin软件使用过程中模型拟合中初值输入对话框
(二)BAPP软件
BAPP软件是由中国药科大学药代中心编写的生物利用度数据处理通用软件(BioAvailability Program Package),专门用于生物利用度研究和生物等效性评价数
320
据处理,最大的特点是在EXCEL的基础上进行二次开发,不但具有EXCEL强大、方便的常用数据处理能力,而且具有强大的专业数据分析能力。目前版本为2.3,其主要功能如下:
图13-6 BAPP软件功能
其中全程自动处理功能提供了全球首创的“傻瓜”操作模式,所有表格、图表、药动学数据处理、等效性检验结果自动转移到WORD中,并且自动排版,可以直接用到用户的实验报告中,并且支持数据缺失,支持低于检测限的数据,支持剂量不一致。
321
图13-7 BAPP软件全程处理功能对话框
另外还支持如下单独功能:
1)在Excel中添加了额外的药动学计算函数,其中包括AUC、MRT、AUMC、t1/2等,支持数据缺失,能将其自动删除后计算。
2)生物等效性检验,支持三交叉试验设计。
3)药代动力学参数的拟合,主要采用改良单纯形法和Gauss-Newton (Levenberg Hartly)法对血药浓度—时间数据进行寻优拟合,找到最小残差平方和,给出血药浓度预报值和一系列药动学参数,并进行作图。
4)权重直线回归,能同时给出不权重(权重因子为1)、权重因子为1/y、权重因子为1/y2三种权重回归结果。
5)自动作图,包括五种生物利用度研究报告中的常用作图模式,并将图形自动拷贝到WORD中。
6)缓释制剂体内外相关性分析。 7)Tmax非参数检验,支持三交叉试验设计。
BAPP软件2002年4月通过了国内著名药动学软件专家组成的评审委员会的审评,目前国内开展生物利用度研究和生物等效性评价的单位中有较广泛的应用。
国内外还有大量的优秀软件,国外如Kinetica,国内如DAS,功能也非常强大,应用也比较广泛,用户选择余地很大。对于药动学软件,通过本章的介绍,我们应该知道即便使用了权威的软件,如果忽视了模型的前提和假设,所得结果也是错误的,所以数据分析处理过程中千万不能迷信软件,轻视对药动学基本概念的掌握。
参考文献:
1.金有豫,药理学(第五版),人民卫生出版社,2001
2.王贤才主译.临床药物大典,青岛出版社,1994
3.吴钟琪、周宏灏、许树梧。全科医学临床药物学,科学出版社,2001
4.戴得哉.临床药理学,中国医药科技出版社,2002
5.Osbourne R,Joel S,Trew D,and Slenvin Morphine and
metabolite behavior after different routes of morphine
tration of the importance of the active
metabolite morphine-6-glucuronide,.,1990,47:12-19
6. 谢海棠黄晓晖孙瑞元中国临床药理学与治疗学 2001,6
(4)289~292
322
2023年6月21日发(作者:)
第十三章曲线拟合以及常用的药代动力学软件
第十三章 药动学数据的曲线拟合以及常用软件
计算药动学参数,是药代动力学进一步应用的基础。如何测定有关的动力学参数呢?常用的方法是:首先在用药后的若干不同时间,采取血样(或尿样),测定其血药浓度值或尿中药量(这些数值称为实测值或观察值,用C i 表示),这样就有了药物浓度经时曲线数据;然后,依据半对数坐标图,选定一种模型方程
(是时间t 的曲线函数)计算理论估算值(用i C )表示),按照观察值和理论估算值之
差的平方和(即残差平方和)或加权残差平方和(均用Re 表示)最小的原则,采用适当的算法,求出有关的动力学参数。这种方法,在数学上称为曲线拟合(fitting a curve)。由于所采用的线性药代动力学的模型方程是多指数项之和的函数形式,并且是所含动力学参数的非线性函数,所以这种曲线拟合方法称为非线性最小二乘法。
一、最小二乘法的一般原理
设y 是变量x 的函数,含有m 个待定参数a 1,a 2,…,a m 。记为
y =f (x ;a 1,a 2,…,a m ) 若对x 和y 作n 次观察,测得观察值(x 1,y 1),(x 2,y 2),…,(x n ,y n )。根据这样一组二维数据,即平面上的若干点,要求确定这个一元函数y =f (x ;a 1,a
2,…,a m );(i=1,2,…,n),即一条曲线,使这些点与曲线总体来说尽量接近。并使y
的观察值y i 与理论估算值=i y )f (x i ;a 1,a 2,…,a m )
;(i=1,2,…,m)的误差平方和,即残差平方和21
()n e i i i y y ==?∑)2121((,,,,)n i i m i y f x a a a ==?…∑R )i 取得最小值,
或者加权残差平方和21()/n e i i i R y y w ==?∑)2121
((,,,,))/n i i m i i y f x a a a w ==?…∑ 取得
最小值。其中w i 称为权重系数,在后面的段落会详细讲解。这时Re 有时候也称为目标函数。这就是数据拟合成曲线的思想,简称为曲线拟合。曲线拟合的目的是根据实验获得的数据去建立因变量和自变量之间有效的函数关系,这个函数关系对于药动学来讲就是通过房室模型推导出来的药时曲线公式。根据观察值求出待 311
定参数(因而也就确定了曲线y =f (x ;a 1,a 2,…,a m ))的问题。我们称f (x ;a 1,a 2,…,a m )为拟合函数。特别地,当f 为x 的线性函数时,则称为直线拟合(或直线回归)。
如下图所示,Re 是点(t i ,C i )与曲线)(x f y =的距离,曲线拟合实际含义就是寻求一个函数,使在某种准则下与所有数据点最为接近,即曲线拟合得最好。最小二乘法准则就是使所有离散点到曲线的距离的平方和最小。例
如对于一房室静注模型函数)(x f y =)(x f i C )=C 0 e -kt i ,就是确定待定系数K 和C 0,使当所有时
间数据点t i ,i =0,1,……,n ,代入函数i C )=C 0 e -kt i 后,按公式计算使求算的结果最小,此时待定系数的值K 和C ∑=?=n i i i
e C C R 1
2
)()0的值就是拟合所求的结果。 度g m l
图13-1 药时曲线拟合示意图
由高等数学中的极值原理知,待定参数a 1,a 2,…,a m 应满足下列方程组(一般称为正规方程)
0Re =??j
a (j=1,2,…,n) (14-1) 这是含有m 个未知数的m 个方程,解这一方程组可求得a 1,a 2,…,a m ,从而确定了拟合函数。当f 是参数的线性函数时,上述正规方程为参数的线性代数方程组,这种情况称为线性最小二乘法。当f 是参数的非线性函数时,上述正规方程则为参数的非线性方程组,这种情况称为非线性最小二乘法。
曲线拟合寻求一个函数,使在某种准则下与所有观察值最为接近,这种搜索最小值的算法目前常采用高斯—牛顿迭代法、单纯形法等。在具体)(x f y =)(x f 312 计算时,有时候采用简化方法—对数回归法或残数法,将非线性方程转化为线性方程,计算虽然比较简单,但计算的误差往往比较大,同时手工计算比较费时费力。目前经常采用高斯—牛顿迭代法、单纯形法等算法编制成计算机程序,当数据比较符合理论情况时,能够比简化计算方法计算出更精确、合理的动力学参数。
二、非线性最小二乘法算法的比较
1、经典高斯—牛顿迭代法以及改进方法—哈特莱方法(Levenberg-Hartley法)、阻
尼最小二乘法
高斯—牛顿迭代法将目标函数(Re)在待定药动学参数初值(a1,a2,…,a m)附近的微分方程用泰勒级数的一次项展开,得正规方程组,采用列主元高斯消元法解出该方程组,就可计算得理论上使得目标函数(Re)最小的最佳药动学参数。
该方法的优点:在某种程度上,按最佳梯度方向搜索,效果较好,虽对初值有一定的依赖性,但依赖程度远远低于其他类型的方法,如单纯形法。所以开始运算时往往收敛较快,运算时间短,这是此方法的突出优点。
缺点:由于泰勒展开中丢弃了高次项等等的原因,使得此法往往不能精确收敛,甚至会引起发散,不能求出解。基于这个缺点,对于经典高斯—牛顿迭代法进行了种种的改进,如哈特莱方法、阻尼最小二乘法等等,可以有效地改进拟合发散的缺点,但还是不能完全避免。
2、单纯形法
本方法是一种多维搜索的直接方法,不需要计算目标函数的导数。通过对n 维空间的n+1个点(它们构成一个初始单纯形)上的函数值进行比较,去掉其中函数值最大的点,代之以新的点,从而构成一个新的单纯形。这样,通过多次迭代逐步逼近极小点。
单纯形法的优点是:原理简单,避免了求导运算以及解正规方程组等步骤,从而不会有经典高斯—牛顿迭代法以及改进方法固有的缺点,基本不会发散。缺点是:收敛速度慢,初值依赖性较大,有陷入局部最小值的弊端,一般不单独使用,有人主张用高斯—牛顿迭代法求出大致的结果,再采用单纯形法作局部区域的精确搜索寻优。
三、估算药代动力学参数中的若干问题
313
1、如何选择模型
在曲线拟合时,需要事先选择合适的药动学模型方程再来进行拟合。也就是说同一组数据可以选择不同的药动学模型方程进行拟合(亦即拟合函数),那么如何选择和确定适当的模型呢?一种简单和直观的方法,就是根据实测的血药浓度的对数值对时间作图(称为对数浓度一时间散点图),作粗略的直观的分析。当散点图的分布比较有规则地反映出某一种单指数或多指数态势时,可确定出相应的房室数。但是,在许多情形下,散点图的分布往往似有多指数态势,或者多指数态势不明显,这时就要从中加以挑选。
选择房室数的方法,最常采用赤池信息判据最小准则AIC(Akaike’s Information Criterion)。
假设观测数据点数为n,拟合函数中所含待定参数的个数为m,拟合所得残差平方和(或加权残差平方和)为Re,
∑=?
=
n
i
i i
e
c c
R
1
2
) ();则
AIC=n ln(Re)+2 m (14-2) 赤池提出按AIC值最小的准则来挑选模型。就是说,在几种预定的模型中,较佳的模型其AIC值最小。根据这一准则,我们就可以对几种不同的房室模型分别进行曲线拟合,算出各自的Re和AIC值,比较AIC值的大小,然后选AIC 值最小的房室,作为优选的模型。另外,根据这一准则,我们对同一模型,在用不同的计算方法作曲线拟合时,AIC值最小的那种方法,算得的参数较为精确;不过这时AIC值最小就相当于Re值最小,因此,只要比较Re值的大小即可。
从公式看,AIC值的大小主要和两个因素有关,一个是Re,一个是待定参数个数。不同的房室模型,其待定参数不一样。例如,一房室静注,待定参数为K 和Vd两个参数,对数浓度一时间散点图为直线;二房室静注,待定参数为A、B、α、β四个参数,对数浓度一时间散点图有一个拐点的两根直线。模型中房室个数越多,待定参数越多,拟合函数越复杂,曲线包含的拐点越多,拟合时Re值就会越低。如果仅仅采用Re作为模型选择依据,非常复杂的多房室模型其Re值往往会更低,产生所谓过拟合现象。AIC值兼顾了模型房室个数和Re两个因素,既使得Re值低,同时又限制了复杂函数的使用,增加了拟合函数的实用
314
价值。从实用角度看,房室模型的房室个数不易超过3个。
在选取模型时要特别注意,房室模型中房室的划分具有相对性。当实验数据比较准确和充分时,可以将药物在体内分布的较小的速度差异区分开来,从而可以将体内分成更多的房室;但是当实验数据比较少或者误差较大时,药物在体内的速度差异就无法区分,只能将机体分成较少的房室或仅单一的中央室。由于上述原因,有可能同一种药物,在不同的文献报道中有不同的房室模型,要理解和容忍这种房室划分的相对性。
2、权重系数
曲线拟合采用的目标函数是Re ,Re 值越小,曲线拟合效果越好。但是,在实验观察值中,高、低血药浓度值相差较大时,Re 值的大小往往过分取决于高浓度
的观察值,而忽视了低浓度值的观测作用。例如,同样是||i i c
c )?=0.5的差值, 对于100ng/ml 相对误差很小,对于1ng/ml 相对误差就很大,所以不能把高低浓度的差值等同看待。这时候就需要采用权重系数的方法,即在加权的情况下,求
加权残差平方和的最小值问题。其中W i n
i i i e w c c R /)(12∑=?=)i 称为权重系数。
权重系数最常采用的是:
W i =1/C 2i
则加权残差平方和为
21Re ∑==n i i i i C C C ) (13-3)
即为相对误差平方和最小。
当W i =1,就是不进行权重。有时候W i =1/C i ,权重介于1和1/C 2之间。至于选择何种权重系数,目前仍然在讨论之中,没有统一的方法,需要根据不同的情况进行具体的分析。需要强调的是,不同的权重系数的Re 值和AIC 值没有可比性。
四、曲线拟合的影响因素
在应用计算机程序进行曲线拟合,估算药动学参数时,影响拟合效果的因素很多,主要有以下几个因素:
1.实验设计问题
药时曲线的采样点数和浓度检测的准确性显著地影响着模型的识别和拟合 315
的质量,因而采样时间点和采样点个数的设计是成功进行拟合的第一步。原则上以较少的采样,尽可能获取血药浓度经时曲线的变化形态信息,如峰时间、峰浓度、拐点、变化趋势等。在曲线变化大的地方,应多采样,在曲线变化小的范围内,则少采样。但是,由于个体差异以及其他因素的存在,可能药时曲线的拐点等位置不一样,所以只能尽量做到符合以上原则。另外经验证明,对于m个参数的模型,采样点个数n至少为2m+1,例如二房室口服有5个待定参数,测11个点以上比较可靠。n≤m是不可行的,否则残差的自由度为n-m-1为-1,不符合统计学要求。事实上,国家药品审评机构的指导原则中明确表明药时曲线采样点个数不得少于11个,一般我们进行药时曲线的点数设计时都会超过11个。
2.样品浓度测定的准确性
目前国家药品审评机构对于体内生物样品检测有比较明确的指导原则,其中明文规定,对于生物样品检测,高、中浓度的检测相对误差必须在±15%以内,低浓度必须在20%以内。由于实验误差的存在,曲线拟合时权重系数的选取必须考虑。如果低浓度样品测试比较准确,就要采取权重拟合,反之普通拟合。对于半衰期等参数,更多地取决于药时曲线末端相的下降速率。样品检测的准确性显著地影响到这类参数的拟合。
3.收敛精度的选取
收敛精度的选取,通常采用由大到小逐步修正的办法,即起初精度可以取得大些,视其计算结果是否满意(残差平方和或加权残差平方和,是否可明显减少;观察值与估算值之间的相关系数,是否可明显增大等),再决定是否提高精度(将精度值取小)。在继续拟合时,应选前一次结果作为初值。
4.初值选择和拟合取值范围问题
目前所采用的非线性最小二乘法的算法都面临一个初值如何选取的问题。初值选得好,迭代次数少,收敛速度快,所得结果也相对理想。初值选取一般首先参考文献上的资料或者以往的经验,也可以应用残数法的计算结果作为初值。但是无论初值如何选取,也还是存在初值依赖性的问题。
初值依赖性是指随着初值选取的不同,拟合所得药动学参数也会随之变化。下图是一个实际的利用单纯形法拟合药动学参数过程中目标函数Re值的分布图。
316
该例为口服一房室模型,有三个待定参数:Ka;Dose*F/V,Ke。在Ka=0.3时,罗列Dose*F/V和Ke两个参数在一定范围内的所有对应组合情况下的Re值;从图中可以看到颜色越深,Re值越小;低于某一限度的Re值区域为一个不规则的长条状。寻优过程中,当参数组合点到达最深的区域的边缘1就可以终止拟合。由图可见,随着起始位置不同,从三个方向开始拟合,所到达的目标值并不是归结为一个点,而是三个点,在Re值最小区域的边缘,这就存在药动学参数拟合的初值依赖性。
改变Ka值,可以绘出一系列图,通过多媒体动画技术,将这些图叠加在一起依次出现,依据Re值最小这个标准,就可以找出一个最佳的三个待定参数的组合范围,使得Re值最小。选取该组合范围的中点,就可以避免初值依赖性。但是随着待定参数数目的增加,待定参数的组合数呈几何级数增长,大大增加计算工作量,获得所需结果就仅仅成为理论上的可能了。
另外,由图1可见,参数Dose*F/V的取值在小于0没有实际意义的情况下,Re值依然很小,说明了数据拟合寻优过程中必须限定取值范围,否则拟合出来的数据没有实际意义。
图13-2 拟合中的初值依赖性
Re值变化率最小,并不是Re值最小。Re变化率最小和Re最小的区
1此处的最深区域为Re值最小,但是真正的拟合收敛条件应该为域虽然不一致,但是也同样是一个不规则区域。
317
初值依赖现象不仅仅在单纯形法模型嵌合中存在,其他任何寻优法中都存在,关键在于药物浓度-时间双指数或多指数曲线本身就存在两个或多个最优解。下面表1中的两组一房室模型参数能产生出如图2所示的完全一样的药时曲线。这两组参数的区别仅仅在于Ka和Ke的大小不同,表观分布容积不同。也就是说,如果Ka和Ke大小不能确定,同一组药时数据会有两组同样最优的药动学参数组合。一般情况下Ka>Ke,所以同一组药时数据最优解只有一组药动学参数组合。但是,对于缓控释制剂、特殊释放制剂,就不可能仅仅通过本身口服后的数据来区分吸收指数项和消除指数项,从而来进行吸收解析。这就是吸收解析的时候最好要同时进行静脉给药实验的原因——能够确切地知道药物的消除相的情况。但是,当药物静脉给药和口服给药消除情况不一致的时候,哪怕进行静脉注射也无法进行口服药物的吸收解析,例如口服给药存在肠肝循环,首过效应强等等的情况。
图13-3 多指数曲线存在两个最优解
100
F(%) 100
口服给药剂量(mg)973 973
0.52
Ka(h-1) 2.5
k(h-1) 0.52
2.5
34
V(l) 163
318
五、目前常用计算机药动学数据拟合程序
目前国内外药动学计算软件比较多。国外比较著名的软件有WinNonlin、Kinetica等,国内比较著名的软件有3p87(3p97)、PKBP—N1、BAPP、DAS等。
(一)WinNonlin 软件
WinNonlin为美国Pharsight公司产品,是目前国外应用最广泛的软件,界面友好,功能强大,与其他软、硬件有很好的兼容性。WinNonlin 有标准版、专业版、企业版等三个版本,标准版中包含了药动学与药效学分析的各种工具,专业版和企业版较标准版增加了几个模块,主要用于商业用途。
目前WinNonlin软件最新版本为4.1,其主要功能如下:1、房室模型分析(PK Fitting,Compartmental Analysis),包括19种模型库,支持用户自定义模型;寻优算法包括Nelder Mead simplex,Gauss-Newton (Levenberg Hartly),Gauss-Newton (Hartly)法,权重方式包括1/obs,1/obs*obs,1/calc,1/calc*calc,初值的设定包括了上下限;模型确认标准包括Akaike,Schwarz,CV,SEM等。2、非房室模型分析(Non Compartmental Analysis),包括各种给药途径,还支持末端相半衰期图解求算,支持参数个体和群体设定,可计算稳态数据参数。3、等效性检验,包括平均等效性、群体等效性、个体等效性(ABE、PBE、IBE)。4、PK-PD数据分析,支持PK、PD分步分析和PK-PD统一嵌合的大模型分析,但是尚不支持群体的PK-PD分析。5、提供了“工具箱”(toolbox)功能及帮助功能:非参数重叠法( nonparametric superposition),用来预测多剂量用药后达到稳态的血药浓度;半房室模型法(semicompartmental modeling),用来估算给定时间和血浆浓度的效应地点浓度。6、数据输入输出支持格式众多,包括ASCII、Excel、SAS?、Transport、Oracle、Watson、ODBC-Treiber等,有非常良好的兼容性;输出的图表能形象化地显示数据,可进行编辑修改。7、提供了一定的数据描述统计功能。
下面两幅图片分别为WinNonlin软件使用过程中模型选择和模型拟合初值输入。
319
图13-4 WinNonlin软件使用过程中模型选择对话框
图13-5 WinNonlin软件使用过程中模型拟合中初值输入对话框
(二)BAPP软件
BAPP软件是由中国药科大学药代中心编写的生物利用度数据处理通用软件(BioAvailability Program Package),专门用于生物利用度研究和生物等效性评价数
320
据处理,最大的特点是在EXCEL的基础上进行二次开发,不但具有EXCEL强大、方便的常用数据处理能力,而且具有强大的专业数据分析能力。目前版本为2.3,其主要功能如下:
图13-6 BAPP软件功能
其中全程自动处理功能提供了全球首创的“傻瓜”操作模式,所有表格、图表、药动学数据处理、等效性检验结果自动转移到WORD中,并且自动排版,可以直接用到用户的实验报告中,并且支持数据缺失,支持低于检测限的数据,支持剂量不一致。
321
图13-7 BAPP软件全程处理功能对话框
另外还支持如下单独功能:
1)在Excel中添加了额外的药动学计算函数,其中包括AUC、MRT、AUMC、t1/2等,支持数据缺失,能将其自动删除后计算。
2)生物等效性检验,支持三交叉试验设计。
3)药代动力学参数的拟合,主要采用改良单纯形法和Gauss-Newton (Levenberg Hartly)法对血药浓度—时间数据进行寻优拟合,找到最小残差平方和,给出血药浓度预报值和一系列药动学参数,并进行作图。
4)权重直线回归,能同时给出不权重(权重因子为1)、权重因子为1/y、权重因子为1/y2三种权重回归结果。
5)自动作图,包括五种生物利用度研究报告中的常用作图模式,并将图形自动拷贝到WORD中。
6)缓释制剂体内外相关性分析。 7)Tmax非参数检验,支持三交叉试验设计。
BAPP软件2002年4月通过了国内著名药动学软件专家组成的评审委员会的审评,目前国内开展生物利用度研究和生物等效性评价的单位中有较广泛的应用。
国内外还有大量的优秀软件,国外如Kinetica,国内如DAS,功能也非常强大,应用也比较广泛,用户选择余地很大。对于药动学软件,通过本章的介绍,我们应该知道即便使用了权威的软件,如果忽视了模型的前提和假设,所得结果也是错误的,所以数据分析处理过程中千万不能迷信软件,轻视对药动学基本概念的掌握。
参考文献:
1.金有豫,药理学(第五版),人民卫生出版社,2001
2.王贤才主译.临床药物大典,青岛出版社,1994
3.吴钟琪、周宏灏、许树梧。全科医学临床药物学,科学出版社,2001
4.戴得哉.临床药理学,中国医药科技出版社,2002
5.Osbourne R,Joel S,Trew D,and Slenvin Morphine and
metabolite behavior after different routes of morphine
tration of the importance of the active
metabolite morphine-6-glucuronide,.,1990,47:12-19
6. 谢海棠黄晓晖孙瑞元中国临床药理学与治疗学 2001,6
(4)289~292
322
发布评论