基于随机微分方程的药物动力学建模新方法＊

褚 娟， 万建平， 高秀娟， 陈 汇△

1 华中科技大学数学与统计学院，武汉 430074

2 华中科技大学同济医学院基础医学院药理学系，武汉 430030

药物动力学是一门综合学科，由最初的采用数学分析手段研究药物的体内过程，现已发展成为运用数学、计算机、生理学与生命科学等多学科知识和技术来研究药物体内过程的学科。

在目前的药物动力学建模中，较多的是利用常微分方程描述药物在体内的动力学过程，此方法在描述简单、线性的药物体内过程时确实具有简易、合理等优势，但在处理一些更为复杂的非线性药物动力学模型时就遇到了一些问题。

首先，利用常微分方程建模时，极易出现模型过于简单或过于复杂的情况。若模型方程过于简单，不能真实地反映药物的体内过程，就会产生有偏的参数估计值和相关的残差误差［1］；若模型方程过于复杂，如待估参数过多，就要求大量数据用于参数估计，而临床数据是有限的。

其次，在常微分方程建模中，常假设某些重要的药理参数（如吸收、消除速率）是常数，而实际上，它们可能是随时间变化的，这样也会产生相关的残差误差。

再次，利用常微分方程建立药物动力学方程时只考虑了不相关的残差，如实验中产生的误差，测量误差等［2］，而模型结构设定的不足，模型中参数随时间变化，未知的随机干扰等都可能引起随时间变化的相关残差，忽略这些相关的残差直接影响了参数的估计，最终导致模型的缺陷。所以，如果只把模型中出现的误差简单归结为测量误差显然是不合理，不科学的，这是受研究方法所限。

经过不断地探索、研究，在药物动力学模型中引入随机微分方程理论较好地解决了这些问题。

首先，随机微分方程模型将个体内误差分为两部分：一部分是不相关的测量误差，主要由试验误差、测量误差等造成的；另一部分是系统误差，由可能存在的模型定义的不合理或体内存在的随机干扰引起的［3］。

其次，对于过于简单的模型，引入随机微分方程可以解释给定数据出现的大量变异；对于过于复杂的模型，在参数估计时会有很大困难，随机微分方程可简化模型并优化参数估计。

再次，随机微分方程中扩散项参数的引入，不仅考虑到了模型的不合理性而且可以对其进行定量和定位描述，进而很大程度地改进原有的模型［2］。

因此，本文也是从最初建立的简单的药物动力学的线性模型出发，逐步得到一个既尽可能真实地模拟了药物在体内变化情况又不至太复杂的药物动力学模型，并对建模方法给予了理论上的支持和解释。新模型在揭示药物动力学的非线性等复杂过程的本质时体现了巨大的优势。

1 基于常微分方程（ODEs）的药动学线性模型

这里主要针对常见的几种给药方式来建立基于常微分方程的药动学模型［4］。目前，大多数药物建模时首选的模型通常如表1所示。

表1 药物建模方式Table 1 The modeling method of drug

以口服给药为例，建立药物动力学模型，通常假设为一级吸收和一级消除速率过程：

a.一室模型（图1）

图1 一室模型结构Fig.1 The structure of one apartment model

吸收部位的药动学过程：

体内（血液内）的药动学过程：

b.二室模型（图2）

图2 二室模型结构Fig.2 The structure of two apartment model

吸收部位的药动学过程：

体内药动学过程：

在实际问题中，许多药物的吸收、消除并不是简单的线性变化，而是具有非线性吸收、消除甚至是混合吸收、消除的特征。其中混合吸收、混合消除是指药物在吸收或消除时不仅具有线性过程而且具有非线性过程。

2 基于随机微分方程（SDEs）的药动学模型

当线性药动学模型不能拟合实验数据时，就可以通过引入随机微分方程理论，在原有模型基础上加上一个随机扰动项来避免由于最初假设错误导致后续的研究错误，以建立更合理的药物动力学模型。

利用SDEs建立PK 模型的具体过程如下［1］。

2.1 建立传统的线性PK 模型

通过探索性数据分析和已知的生理学知识建立关于状态变量的常微分方程，即为式（1.1）～式（1.4）的形式，一般记为：

其中，Xt∈Rn为状态向量（如血药浓度），ut∈Rm为输入变量（如输入速率），θ∈Rp为参数向量，f（·）∈Rn为向量函数。

2.2 引入干扰建立SDEs模型

其中，f（·）为漂移项，σ（·）为扩散项，σ（·）∈Rn×n为矩阵函数，Bt为n维标准Wiener过程，且Bt1－Bt2～N（0，｜t1－t2｜Ι），Ι为单位矩阵。

其中，Wiener过程可理解为来源于同分布的随机事件的和产生的Xt的真实变化与f（·）描述的Xt的变化之差。换言之，扩散项σ（·）可解释模型的不确定性［1］。

由于dB（t）dt＝0，（dt）2＝0，（dB（t））2＝d［B，B］（t）＝dt［5］，所以，式（2.2）简化为：

在Ito积分意义下：

在Stratonovich积分意义下：

显然，本文中σ（ut，t，θ）与Xt无关，所以无论是Ito积分还是Stratonovich 积分，最终的结果均一样。

2.3 建立测量方程

其中，yk∈Rl为测量得到的输出变量（如血液中药物浓度），h（·）∈Rl为向量函数，tk，k＝0…N为抽样时间，ek为l维白噪声过程，即ek～N（0，S（θ）），S（θ）为协方差矩阵。

2.4 参数估计判别扩散项的显著性

根据极大似然估计原理，目前常结合EKF（Extended Kalman filter）等方法进行参数估计［2］，并通过软件CTSM、R 等软件实现。但是，EKF 使用的前提是假设观测值的条件密度可以用Gaussian密度近似；当f（·）为关于Xt的非线性函数时，EKF是利用f（·）的一阶Taylor展式求得Gaussian分布的均值、方差的近似值来计算，基于此得到的是一个近似的极大似然函数。另外，EKF作为一种迭代算法，其停止规则未考虑估计的效果和效率，而是完全由数据个数决定的。

目前，国外有一些文献［6－7］利用由EM 算法衍生出的SAEM－MCMC方法求相关参数，此算法可依靠Monolix软件实现。SAEM－MCMC 算法相对于其他的算法（线性化方法，积分近似及普通的EM算法）计算速度更快，精度更高；其对初始值的设定要求较低，对初始值的选取非常稳健，通过若干步迭代之后，基本上收敛到要模拟的平稳分布；对观测值的条件密度没有特别要求，因此适用范围更广。

2.5 求解随机微分方程

扩展模型状态方程，指出模型的不足。反复步骤2.4直至模型的系统误差显著地为零为止。

若σ显著为零，则说明新得到的药物动力学模型较好地解释了药物的体内变化情况；若σ不显著为零，则说明前面建立的ODEs模型不能很好地模拟药物体内变化的实际情况，存在着某些无法给出合理解释的σ（ut，t，θ）dωt。这部分误差和干扰，从建模的速率角度分析，可能是有关吸收、消除速率过程的假设有不合理之处。因此就面临着如何求解SDEs的问题。对于求解SDEs，我们可以从以下3方面考虑：

（1）直接利用Ito积分求精确解。虽然直接利用Ito积分求得的解最准确，但只有极少部分特殊的SDEs可以求出精确解，而大部分SDEs不可直接求解。

（2）求SDEs的数值解。SDEs数值解只是用随机过程无数条可能的轨道中的一个轨道表示过程的特征参数的性质［8］，过于特殊而缺乏普遍性；目前，多用随机微分方程在固定时刻附近的随机Taylor展开或Runge－Kutta随机差分方程来求解SDEs的数值解，这种采用随机过程的一段相当长时间的数据来估值的方法对于反映药物在各个时刻的浓度情况不是十分理想；另外与解ODEs相比，求随机微分方程的数值解仍是较复杂、费时的。

（3）利用ODEs的解去近似SDEs的解。这种近似不能随意确定，在一定条件下才能执行。

因此，采用扩展模型状态方程来合理地修正吸收、消除过程使σ（·）显著为零，此时就可利用常微分方程的解去近似随机微分方程的解。

此方法的优点在于，ODEs近似SDEs的解是在均方意义下的近似，因而精度可以得到保证；另外，求解常微分方程的理论和方法已十分成熟，为计算带来很大便利。总之，此方法在不同程度上弥补了方法（1）、（2）的不足，是十分有意义的。

2.6 用微分方程模型解释药动学过程

解常微分方程dXt＝f＊（Xt，ut，t，θ）dt，得到药物浓度随时间变化表达式，用于预测体内药物浓度变化，从而改进给药方式、给药时间间隔，使药浓度维持在治疗窗内，以优化药效并防止毒副作用。

3 结语

本文从速率论的角度，在传统的药物动力学线性模型的基础上，利用随机微分方程完善复杂药物动力学模型，使其更加接近真实地描述药物在体内的变化情况，由此得到的血药浓度随时间变化规律，对临床药物研究更具有积极意义。但本文只是对基于随机微分方程的药动学模型给出了理论上的合理解释、建模步骤和求解方法，而没有利用真实数据按照随机微分方程建模的步骤通过统计软件模拟并建立一个药物动力学模型，这些问题有待进一步的研究。

［1］ Kristensen N R，Madsen H，Ingwersen S H.Using stochastic differential equations for PK／PD model development［J］.J Pharmacokinet Pharmacodyn，2005，32（1）：109－113.

［2］ Tornøe C W，Overgaard R V，AgersøH，et al.Stochastic differential equations in NONMEM ®：Implementation application and comparison with ordinary differential equations［J］.Pharm Res，2005，22（8）：1248－1256.

［3］ Overgaard R V，Jonsson N，Torne C W，et al.Non－linear mixed－effects models with stochastic differential equations：implementation of an estimation algorithm［J］.J Pharmacokinet Pharmacodyn，2005，32（1）：85－107.

［4］ 魏敏吉，赵明.创新药物药代动力学研究与评价［M］.北京：北京大学医学出版社，2008：9－52.

［5］ Klebaner F C.Introduction to stochastic calculus with applications［M］.2版.北京：人民邮电出版社，2008：110.

［6］ Marc Lavielle.Monolix Version 3.1 User Guide［EB／OL］.［2011－9－29］.http：／／www.monolix.org／2010.

［7］ Panhard X，Samson A.Extension of the SAEM algorithm for nonlinear mixed models with 2levels of random effects［J］.Biostatistics，2009，10（1）：121－135.

［8］ 龚光鲁.随机微分方程及其应用概要［M］.北京：清华大学出版社，2008：124－125

［9］ Evans L C.An introduction to stochastic differential equations［M］.Berkely：Department of Mathematics UC Berkeley，2006：91.