应用数学和力学，第35卷第9期Applied Mathematics and Mechanics2014年9月15日出版Vol.35,No.9,Sep.15,2014文章编号：1000-0887（2014)09-1033-13应用数学和力学编委会，ISSN1000-0887血管外给药的非线性房室模型解的逼近唐三一胡晓虎，（陕西师范大学数学与信息科学学院，西安710062）（我刊编委唐三一来稿）摘要：药物动力学模型的解析求解公式在新药设计特别是药物动力学参数确定等方面具有非常重要的意义.近年来，由于非线性米氏消除速率过程确定的药物动力学模型解析求解公式的获得，使得大多数单房室模型的解析解基本确定.但是，由于刻画血管外给药的非线性米氏药物动力学模型是一个非自治系统，进而不可能寻求其解析求解公式.该文的目的是讨论一次性血管外给药和周期血管外给药下非线性药物动力学模型解的逼近问题.采用微分方程和脉冲微分方程的比较定理并借助LambertW函数的定义以及相关性质给出模型的不同上下界，估计模型解的逼近程度，并通过数值模拟进行验证关键词：房室模型：米氏消除速率；血管外给药：LambertW函数：药物动力学中图分类号：0241.8：0242文献标志码：Adoi:10.3879/j.issn.1000-0887.2014.09.009引言药物的主动转运，肝脏的代谢，胆汁的分泌以及肾小管的分泌都需要酶或载体参与.这些酶或载体参与的过程不但特异性强，而且受酶或载体本身数量的限制，也就是说其转运或消除速度有一定的限制.上述特性决定了一级消除、零级消除速率过程不能准确刻画药物在体内的消除过程，需要采用基于酶催化反应而提出的非线性米氏消除速率过程-来刻画，即有下面的方程：-VC(t)dc(t)D(1), C(0) = C = ,=- + C()dt其中C（t）表示血药浓度，Vmx表示最大反应速率，K是Michaelis-Menten常数，D是一次静脉注射给药的药物剂量，V是中心室表观分布容积.来氏消除速率过程在药物动力学中应用非常广泛，尤其是在酶动力学中-Lundquist和Wolthers在1958年将米氏消除速率过程应用在恒定静脉注射和周期静脉注射中，提出了模型（1）并得到了该模型如下积分形式的解析解：1[DD- C(t) +Kmln((2)=t-loVlv(Vc())收稿日期：2014-03-18；修订日期：2014-06-16基金项目：国家自然科学基金（11171199)作者简介：胡晓虎（1989一），男，湖北宜城人，硕士：唐三一（1970一），男，教授，博士生导师（通讯作者，E-mail：sytang@snnu.edu.cn）1033?1994-2014 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net
书 文章编号:1000-0887(2014)09-1033-13 应用数学和力学编委会，ISSN 1000-0887 血管外给药的非线性房室模型解的逼近* 胡晓虎， 唐三一 (陕西师范大学 数学与信息科学学院，西安 710062) (我刊编委唐三一来稿) 摘要: 药物动力学模型的解析求解公式在新药设计特别是药物动力学参数确定等方面具有非常 重要的意义． 近年来，由于非线性米氏消除速率过程确定的药物动力学模型解析求解公式的获得， 使得大多数单房室模型的解析解基本确定． 但是，由于刻画血管外给药的非线性米氏药物动力学 模型是一个非自治系统，进而不可能寻求其解析求解公式． 该文的目的是讨论一次性血管外给药 和周期血管外给药下非线性药物动力学模型解的逼近问题． 采用微分方程和脉冲微分方程的比较 定理并借助 Lambert W 函数的定义以及相关性质给出模型的不同上下界，估计模型解的逼近程度， 并通过数值模拟进行验证． 关 键 词: 房室模型; 米氏消除速率; 血管外给药; Lambert W 函数; 药物动力学 中图分类号: O241． 8; O242 文献标志码: A doi: 10． 3879 /j． issn． 1000-0887． 2014． 09． 009 引 言 药物的主动转运，肝脏的代谢，胆汁的分泌以及肾小管的分泌都需要酶或载体参与． 这些 酶或载体参与的过程不但特异性强，而且受酶或载体本身数量的限制，也就是说其转运或消除 速度有一定的限制． 上述特性决定了一级消除、零级消除速率过程不能准确刻画药物在体内的 消除过程，需要采用基于酶催化反应而提出的非线性米氏消除速率过程［1-4］来刻画，即有下面 的方程: dC(t) dt = － VmaxC(t) KM + C(t) ，C(0) = C0 = D V ， (1) 其中 C(t) 表示血药浓度，Vmax 表示最大反应速率，KM 是 Michaelis-Menten 常数，D 是一次静脉 注射给药的药物剂量，V 是中心室表观分布容积． 米氏消除速率过程在药物动力学中应用非常 广泛，尤其是在酶动力学中［5-6］． Lundquist 和 Wolthers［7］在 1958 年将米氏消除速率过程应用在恒定静脉注射和周期静脉 注射中，提出了模型(1)并得到了该模型如下积分形式的解析解: 1 Vmax D V － C(t) + KM ln D VC(t [ ( ) ] ) = t － t0 ． (2) 3301 应用数学和力学，第 35 卷 第 9 期 2014 年 9 月 15 日出版 Applied Mathematics and Mechanics Vol． 35，No． 9，Sep． 15，2014 * 收稿日期: 2014-03-18; 修订日期: 2014-06-16 基金项目: 国家自然科学基金(11171199) 作者简介: 胡晓虎(1989—)，男，湖北宜城人，硕士; 唐三一(1970—)，男，教授，博士生导师(通讯作者． E-mail: sytang@ snnu． edu． cn)．
1034胡晓虎唐三一利用上述积分形式的解，作者结合试验数据进行了药物动力学模型的数据拟合和相应药物动力学参数估计的研究。后来，Beal,Godfrey和Fitch[i系统地研究了积分方程（2）关于药物浓度C(t）的可解性。进一步，Goudar等，Schnell和Mendozai2分别给出了模型具有一次静脉推注时如下的解析求解公式： exp(Co - V (t -t0))C(t) = KmW(e(3)KMKM其中，W（）表示LambertW函数，下同.有关LambertW函数定义和有关性质在第1节给出，Tang和Xiaoi首先给出具有周期静脉推注给药方式下模型（1）的解析求解公式，同时建立了具有治疗窗口的药物动力学模型并给出了模型的解析求解公式.特别地，文献首次对如下恒定给药下的非线性药物动力学模型的解析求解公式进行了系统研究：VmaxC(t)dc(t)D(4)ddo =r-k+Co). C(c) = Co = D,其中r代表常数输入率。Tang和Xiao同样借助LambertW函数的定义和性质给出了其相应的解析求解公式：C(t) =TKubK (-1, - exp(D(_rKm + bCo +b (t - t0)b>0,KuVaasKyVux-Km+/(Km+C)*+2Vm（t-t),b = 0,rKM-bKumwfo, -kthcexp(-k+Cet(-)b<0,KVmaxKuVaux其中b=r-Vmx近些年，针对上述解析求解公式的数值计算和近似计算，也得到深入研究[44S对于血管外给药（比如口服给药），由于药物从用药部位进入血液循环要经过吸收过程，因此，通过口服剂量为D的药物首先进入吸收室，在室内逐渐被吸收，同时将吸收到的药物逐渐向中心室输送，根据这一过程可以建立相应的血管外给药的单房室模型.如果记C，为吸收室药物的浓度，V，为吸收室表观分布容积，则药物浓度在吸收室的变化规律为dC, (t)D= - K,C,(t), C,(0) = C =dtV.其中，K，为药物在吸收室的消除速率常数，求解上述方程得到一次性血管外给药吸收室的浓度为C, (t) = Cge-ky.如果假设从吸收室吸收的药物以生物利用度F的比率输送到中心室，因此中心室的血药浓度C（t）的变化规律可以描述为?1994-2014 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net
利用上述积分形式的解，作者结合试验数据进行了药物动力学模型的数据拟合和相应药物动 力学参数估计的研究． 后来，Beal［8-9］，Godfrey 和 Fitch［10］系统地研究了积分方程(2)关于药物浓度 C(t) 的可解 性． 进一步，Goudar 等［11］，Schnell 和 Mendoza［12］分别给出了模型具有一次静脉推注时如下的 解析求解公式: C(t) = KM W C0 KM exp C0 － Vmax(t － t0 ) K ( ( ) ) M ， (3) 其中，W() 表示 Lambert W 函数，下同． 有关 Lambert W 函数定义和有关性质在第 1 节给出． Tang 和 Xiao［13］首先给出具有周期静脉推注给药方式下模型(1)的解析求解公式，同时建 立了具有治疗窗口的药物动力学模型并给出了模型的解析求解公式． 特别地，文献［7］首次对 如下恒定给药下的非线性药物动力学模型的解析求解公式进行了系统研究: dC(t) dt = r － VmaxC(t) KM + C(t) ，C(t + 0 ) = C0 = D V ， (4) 其中 r 代表常数输入率． Tang 和 Xiao［13］同样借助 Lambert W 函数的定义和性质给出了其相应 的解析求解公式: C(t) = － rKM b + KM VmaxW － 1，－ rKM + bC0 KM Vmax exp － rKM + bC0 + b2 (t － t0 ) KM V ( ( ) ) max ， b ＞ 0， － KM + (KM + C0 )2 + 2Vmax(t － t 槡 0 )， b = 0， － rKM b + KM VmaxW 0，－ rKM + bC0 KM Vmax exp － rKM + bC0 + b2 (t － t0 ) KM V ( ( ) ) max ， b ＜ 0 ， 其中 b = r － Vmax ． 近些年，针对上述解析求解公式的数值计算和近似计算，也得到深入研究［14-15］． 对于血管外给药(比如口服给药)，由于药物从用药部位进入血液循环要经过吸收过程． 因此，通过口服剂量为 D 的药物首先进入吸收室，在室内逐渐被吸收，同时将吸收到的药物逐 渐向中心室输送，根据这一过程可以建立相应的血管外给药的单房室模型． 如果记 Ca 为吸收 室药物的浓度，Va 为吸收室表观分布容积，则药物浓度在吸收室的变化规律为 dCa(t) dt = － KaCa(t)，Ca(0) = Ca 0 = D Va ， 其中，Ka 为药物在吸收室的消除速率常数． 求解上述方程得到一次性血管外给药吸收室的浓 度为 Ca(t) = Ca 0 e －Kat ． 如果假设从吸收室吸收的药物以生物利用度 F 的比率输送到中心室，因此中心室的血药 浓度 C(t) 的变化规律可以描述为 4301 胡 晓 虎 唐 三 一
1035血管外给药的非线性房室模型解的逼近V.C(t)dc(t)= FK,C, (t)KM+C()=dtV.C(t)(5)FCgK,e-ka --K + C(t)明显地，上述方程是一个非自治的一阶微分方程，由于变量C和t是不可分离的，因此无法求出该方程的解析解，由于实际问题和药物动力学研究的需要，有必要知道上述血管外给药下药物在体内的浓度如何变化？因此采用微分方程和脉冲微分方程的比较定理并借助LambertW函数的定义以及相关性质给出模型（5）的不同上下界，估计模型解的逼近程度，并通过数值模拟分析解的可行性，1LambertW函数和特殊情形在介绍具体的方法之前，首先给出Lambert W函数的定义和一些基本性质[定义1LambertW函数是函数f(z）=ze的反函数，且满足：W(z)exp(W(z)) = z,为了方便使用，本文将LambertW函数记做W.注意函数f(z）=ze在z>-1时有正导数f（z）=z+1)e.在区间[-1,αo）上定义f(z）=ze的反函数为W(0,z)=W(z)：易得当z>0时有W（0,z）>0，且当z>-1时W(0，z）为增函数.类似地，可以在区间（-，-1）上定义ze的反函数为W（-1，z).通过计算容易得到其导数满足如下关系式：W(2)W (z) =z(1 + W(z))为了研究模型的解的性质，首先考虑以下特殊情形情形1米氏常数远远大于血药浓度时的情形，即K》C.假设K>C，则方程（5）可以简化成Vc(t).dc()= FCK,e-k .(6)dtKu上述方程是一个简单的一阶线性非齐次方程，求解得到血药浓度C（t）满足：, VmaxFCK,KM(e-Vmat/KMK≠e-kt),K,K- VKMC(t):Va[FCK,te-*,K, =KM由此可见，当K》C时方程（5）刻画了血管外给药时体内药物浓度与时间的关系，即药物浓度在时刻KMIn(Vmxtmax = - 7KKM-V.KKM时达到峰值，血药浓度在时间区间（0，t）内由小变大，然后从t开始逐渐变小，其数值结果如图1所示，情形2米氏常数远远小于血药浓度时，即有KC假设K≤C，则方程（5）可以简化成dc(t)= - V + FCgK,e-Ka.dt?1994-2014 China Academic Journal Electronic Publishing House.All rights reserved.http:/www.cnki.net
dC(t) dt = FKaCa(t) － VmaxC(t) KM + C(t) = FCa 0Ka e －Kat － VmaxC(t) KM + C(t) ． (5) 明显地，上述方程是一个非自治的一阶微分方程，由于变量 C 和 t 是不可分离的，因此无法求 出该方程的解析解． 由于实际问题和药物动力学研究的需要，有必要知道上述血管外给药下药 物在体内的浓度如何变化? 因此采用微分方程和脉冲微分方程的比较定理并借助 Lambert W 函数的定义以及相关性质给出模型(5)的不同上下界，估计模型解的逼近程度，并通过数值模 拟分析解的可行性． 1 Lambert W 函数和特殊情形 在介绍具体的方法之前，首先给出 Lambert W 函数的定义和一些基本性质［16］． 定义 1 Lambert W 函数是函数 f(z) = zez 的反函数，且满足: W(z)exp(W(z)) = z， 为了方便使用，本文将 Lambert W 函数记做 W． 注意函数 f(z) = zez 在 z ＞ － 1 时有正导数 f '(z) = (z + 1)ez ． 在区间［－ 1，∞ ) 上定义 f(z) = zez 的反函数为 W(0，z) = W(z)． 易得当 z ＞ 0 时有 W(0，z) ＞ 0，且当 z ＞ － 1 时 W(0，z) 为增函数． 类似地，可以在区间( － ∞ ，－ 1］上定义 zez 的反函数为 W( － 1，z)． 通过计算容易得到其导数满足如下关系式: W'(z) = W(z) z(1 + W(z))． 为了研究模型的解的性质，首先考虑以下特殊情形． 情形 1 米氏常数远远大于血药浓度时的情形，即 KM C． 假设 KM C，则方程(5)可以简化成 dC(t) dt = FCa 0Ka e －Kat － Vmax KM C(t)． (6) 上述方程是一个简单的一阶线性非齐次方程，求解得到血药浓度 C(t) 满足: C(t) = － FCa 0KaKM KaKM － Vmax e －Vmaxt/KM － e ( －Ka ) t ， Ka ≠ Vmax KM ， FCa 0Ka te －Kat ， Ka = Vmax KM { ． 由此可见，当 KM C 时方程(5)刻画了血管外给药时体内药物浓度与时间的关系，即药物浓 度在时刻 tmax = － KM KaKM － Vmax ln Vmax KaK ( ) M 时达到峰值，血药浓度在时间区间 (0，tmax) 内由小变大，然后从 tmax 开始逐渐变小，其数值结果 如图 1 所示． 情形 2 米氏常数远远小于血药浓度时，即有 KM C． 假设 KM C，则方程(5)可以简化成 dC(t) dt = － Vmax + FCa 0Ka e －Kat ． 血管外给药的非线性房室模型解的逼近 5301
1036胡晓虎唐三求解上述微分方程得到血药浓度C(t）为C(t) = C, - Vt - FCfe-kt.容易看出，若Vm≥FCK，此时C（t）≤0，即在此时血药浓度始终是降低的8(ru/Bu)/)5410203040501/hCo=15mg/mL,Vmx=1 mg/(h-mL),Km =10mg/mL,K, =0.4mg/(hamL),F =0.8图1当KC时方程（5)刻画的血药浓度与时间的关系Fig. 1Plasncording to model (5) when Ky Ctratinnthiy2.01.5C-00.5C(0)C,(0)~0.51051520253040503545t/hCg=2mg/mL,Vmx=1 mg/(hmL),Km=1.2mg/mL,K,=0.1 mg/(hamL),F=1图2模型（5)血药浓度的一个粗略上下界Fig. 2A rough upper and lower bounds of model (5)情形3方程（5）刻画的血药浓度的一个粗略上下界情形1和情形2在两种非常特殊的情况下得到了模型（5）的解析求解公式.正如引言中所强调的那样，一般形式的模型（5）已经不可能得到其解析公式了.因此，笔者希望借助引言中讨论的两类可求解非线性米氏消除速率过程刻画的药物动力学模型，即模型（1）和（4），通过微分方程的比较定理得到模型（5）解的一个粗略下界和上界为此记C,（）和C,（）分别表示下述两个模型?1994-2014 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net
求解上述微分方程得到血药浓度 C(t) 为 C(t) = C0 － Vmax t － FCa 0 e －Kat ． 容易看出，若 Vmax ≥ FCa 0Ka，此时 C'(t) ≤ 0，即在此时血药浓度始终是降低的． Ca 0 = 15 mg /mL，Vmax = 1 mg /(h·mL)，KM = 10 mg /mL，Ka = 0． 4 mg /(h·mL)，F = 0． 8 图 1 当 KM C 时方程(5)刻画的血药浓度与时间的关系 Fig． 1 Plasma concentration vs． time according to model(5) when KM C Ca 0 = 2 mg /mL，Vmax = 1 mg /(h·mL)，KM = 1． 2 mg /mL，Ka = 0． 1 mg /(h·mL)，F = 1 图 2 模型(5)血药浓度的一个粗略上下界 Fig． 2 A rough upper and lower bounds of model(5) 情形 3 方程(5)刻画的血药浓度的一个粗略上下界． 情形 1 和情形 2 在两种非常特殊的情况下得到了模型(5)的解析求解公式． 正如引言中 所强调的那样，一般形式的模型(5)已经不可能得到其解析公式了． 因此，笔者希望借助引言 中讨论的两类可求解非线性米氏消除速率过程刻画的药物动力学模型，即模型(1)和(4)，通 过微分方程的比较定理得到模型(5)解的一个粗略下界和上界． 为此记 C1 (t) 和 C2 (t) 分别表示下述两个模型 6301 胡 晓 虎 唐 三 一
1037血管外给药的非线性房室模型解的逼近dC,(t)VC,()dC, (t)Vm C, (t)=FCoK,-dtK+CdtKM+C,()的解，其解析求解公式引言中已经给出.由于VmC, (t)VC, (t)VC(t)< FC,K,e-Ky ≤FCK.-+C()Kμ + C()K+C,(t)根据微分方程的比较定理，模型（5）的解C(t）满足关系C,（t）<C(t）≤C,（t)，这样就得到了模型（5）解的一个粗略下界和上界，如图2所示由于在寻找模型（5）解的一个上下界时放缩幅度过大，得到了如图2所示的精确解与上下界相距甚远.这样不能用这两个上下界确定或判别模型（5）解的一些行为特点.为此，在下面将采用更加技巧性的方法确定模型（5）的一个非常精确的上下界2一次性血管外给药下模型（5）解的精确上下界为了寻求模型（5）的一个非常精确的上下界，采用如下的分析办法：在每一个整数区间[n，n+1]（n=0，1，2，）上考虑模型（5），寻找该模型在区间[n，n+1]上的一个上下界，然后推广到整个区间上，实现寻求模型（5）解的一个精确上下界的目标.在区间[n，n+1]上首先考虑如下微分方程：dC,()V.C, (t)= FCK,e-K,-(7)ten,n+i],dtK+C()其中表示取整，即在区间[n，n+1]上=n或=n+1.为此先考虑=n时的情形，即在区间[n，n+1]上先考虑如下的微分方程：dC (t)VC, (t)= FCgK,e-Ka --(8)te[n,n+l)dtK+C,()将上述方程重新整理得dC,(t)u+g,C, ()(9)te[n,n+1),dt=Km+ C,(t)其中un=FCK,Kue-ka,ga=FCgK,e-k-Vmx下面根据g，的符号考虑几种不同情形.首先考虑g，=0的情形，此时关于n求解g，=0得到一个根，记做N，，且In(FCK.)1(10)N, =K因此，如果(11)FCK,>V'则N，是正的.当g,=0时容易得到方程(9）在区间[n,n+1）上的解为C,(o) =-K + /(C,(n) + K)? +2u. (t -n).若g。+0时，模型（9）在区间[n，n+1）上的解可以转化为下述等价形式：(1 + Kug,)dc,() = g.d..u,+g,C,())将上述方程右端从n到t开始积分得+Kug. -"s dc,(t) = g, J.dt1+u.+g,C,(o))?1994-2014 China Academic Journal Electronic Publishing House.All rights reserved.http:/www.cnki.net
dC1 (t) dt = － VmaxC1 (t) KM + C1 (t) ，dC2 (t) dt = FCa 0Ka － VmaxC2 (t) KM + C2 (t) 的解，其解析求解公式引言中已经给出． 由于 － VmaxC1 (t) KM + C1 (t) ＜ FCa 0Ka e －Kat － VmaxC(t) KM + C(t) ≤ FCa 0Ka － VmaxC2 (t) KM + C2 (t) ． 根据微分方程的比较定理，模型(5)的解 C(t) 满足关系 C1 (t) ＜ C(t) ≤ C2 (t)，这样就得到 了模型(5)解的一个粗略下界和上界，如图 2 所示． 由于在寻找模型(5)解的一个上下界时放缩幅度过大，得到了如图 2 所示的精确解与上 下界相距甚远． 这样不能用这两个上下界确定或判别模型(5)解的一些行为特点． 为此，在下 面将采用更加技巧性的方法确定模型(5)的一个非常精确的上下界． 2 一次性血管外给药下模型(5)解的精确上下界 为了寻求模型(5)的一个非常精确的上下界，采用如下的分析办法:在每一个整数区间 ［n，n + 1］(n = 0，1，2，.) 上考虑模型(5)，寻找该模型在区间［n，n + 1］上的一个上下界， 然后推广到整个区间上，实现寻求模型(5)解的一个精确上下界的目标． 在区间［n，n + 1］上首先考虑如下微分方程: dC1 (t) dt = FCa 0Ka e －Ka ［t］ － VmaxC1 (t) KM + C1 (t) ， t ∈［n，n + 1］， (7) 其中［t］表示取整，即在区间［n，n + 1］上［t］ = n 或［t］ = n + 1． 为此先考虑［t］ = n 时的 情形，即在区间［n，n + 1］上先考虑如下的微分方程: dC1 (t) dt = FCa 0Ka e －Kan － VmaxC1 (t) KM + C1 (t) ， t ∈［n，n + 1)． (8) 将上述方程重新整理得 dC1 (t) dt = un + gnC1 (t) KM + C1 (t) ， t ∈［n，n + 1)， (9) 其中 un = FCa 0KaKM e －Kan ，gn = FCa 0Ka e －Kan － Vmax ． 下面根据 gn 的符号考虑几种不同情形． 首先考虑 gn = 0 的情形，此时关于 n 求解 gn = 0 得到一个根，记做 N1，且 N1 = 1 Ka ln FCa 0Ka V ( ) max ． (10) 因此，如果 FCa 0Ka ＞ Vmax， (11) 则 N1 是正的． 当 gn = 0 时容易得到方程(9) 在区间［n，n + 1) 上的解为 C1 (t) = － KM + (C1 (n) + KM )2 槡 + 2un ·(t － n)． 若 gn ≠ 0 时，模型(9) 在区间［n，n + 1) 上的解可以转化为下述等价形式: 1 + KM gn － un un + gnC1 (t ( ) ) dC1 (t) = gn dt． 将上述方程右端从 n 到 t 开始积分得 ∫ C1(t) C1(n) 1 + KM gn － un un + gnC1 (t [ ] ) dC1 (t) = gn ∫ t n dt， 血管外给药的非线性房室模型解的逼近 7301