第十三章 曲线拟合以及常用的药代动力学软件

第十三章 药动学数据的曲线拟合以及常用软件计算药动学参数，是药代动力学进一步应用的基础。

如何测定有关的动力学参数呢？常用的方法是：首先在用药后的若干不同时间，采取血样（或尿样），测定其血药浓度值或尿中药量（这些数值称为实测值或观察值，用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 12)()0的值就是拟合所求的结果。

度g m l图13-1 药时曲线拟合示意图由高等数学中的极值原理知，待定参数a 1，a 2，…，a m 应满足下列方程组（一般称为正规方程）0Re =∂∂ja (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个点(它们构成一个初始单纯形)上的函数值进行比较，去掉其中函数值最大的点，代之以新的点，从而构成一个新的单纯形。

这样，通过多次迭代逐步逼近极小点。

单纯形法的优点是：原理简单，避免了求导运算以及解正规方程组等步骤，从而不会有经典高斯—牛顿迭代法以及改进方法固有的缺点，基本不会发散。

缺点是：收敛速度慢，初值依赖性较大，有陷入局部最小值的弊端，一般不单独使用，有人主张用高斯—牛顿迭代法求出大致的结果，再采用单纯形法作局部区域的精确搜索寻优。

三、估算药代动力学参数中的若干问题3131、如何选择模型在曲线拟合时，需要事先选择合适的药动学模型方程再来进行拟合。

也就是说同一组数据可以选择不同的药动学模型方程进行拟合（亦即拟合函数），那么如何选择和确定适当的模型呢？一种简单和直观的方法，就是根据实测的血药浓度的对数值对时间作图（称为对数浓度一时间散点图），作粗略的直观的分析。

当散点图的分布比较有规则地反映出某一种单指数或多指数态势时，可确定出相应的房室数。

但是，在许多情形下，散点图的分布往往似有多指数态势，或者多指数态势不明显，这时就要从中加以挑选。

选择房室数的方法，最常采用赤池信息判据最小准则AIC(Akaike’s Information Criterion)。

假设观测数据点数为n，拟合函数中所含待定参数的个数为m，拟合所得残差平方和（或加权残差平方和）为Re，∑=−=nii iec cR12) ()；则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 ni 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，不符合统计学要求。