余曉美, 張海永，沈永昌

基于貝葉斯推斷估計(jì)藥物動(dòng)力學(xué)模型參數(shù)

余曉美, 張海永，沈永昌

為克服藥理觀測(cè)數(shù)據(jù)的不確定性和未知參數(shù)高維化給非線性藥物動(dòng)力學(xué)模型參數(shù)估計(jì)帶來(lái)的困難，基于貝葉斯定理推導(dǎo)了模型參數(shù)的后驗(yàn)分布。采用馬爾科夫鏈蒙特卡羅(MCMC)方法對(duì)參數(shù)的后驗(yàn)分布進(jìn)行采樣，從而獲得參數(shù)的估計(jì)值。抽樣方法為延遲拒絕自適應(yīng)抽樣算法(DRAM)，實(shí)證結(jié)果表明，DRAM方法能夠?qū)λ幬飫?dòng)力學(xué)模型進(jìn)行良好的參數(shù)估計(jì)。對(duì)比DRAM算法和常用、經(jīng)典的Metropolis Hastings(MH)模擬結(jié)果，可以發(fā)現(xiàn)DRAM算法在求解藥物動(dòng)力學(xué)模型參數(shù)估計(jì)問(wèn)題具有更小的誤差，同時(shí)具有穩(wěn)定、收斂速度更快，高效的特點(diǎn)。

貝葉斯推理；藥物動(dòng)力學(xué)模型；DRAM；MCMC

1 前言

藥物動(dòng)力學(xué)主要用于建立監(jiān)測(cè)個(gè)體的體內(nèi)藥量或藥物濃度隨時(shí)間變化的數(shù)學(xué)表達(dá)式，并求算出有關(guān)藥動(dòng)學(xué)參數(shù)，從而研究藥物、毒物及其代謝在體內(nèi)的吸收、分布、代謝和排泄的動(dòng)態(tài)過(guò)程及這些過(guò)程與藥理反應(yīng)間的定量規(guī)律。結(jié)合藥動(dòng)學(xué)模型和藥動(dòng)學(xué)參數(shù)，調(diào)整和確定用藥方案，保證藥物治療的有效性和安全性[1]。目前用于刻畫(huà)藥物在體內(nèi)的動(dòng)力學(xué)規(guī)律的動(dòng)力學(xué)模型類(lèi)別多樣，從常微分方程到隨機(jī)微分方程，從房室模型到非房室模型，從線性到非線性模型[2]。非線性藥動(dòng)學(xué)參數(shù)的計(jì)算方法主要有加權(quán)非線性最小二乘法進(jìn)行曲線擬合，Gauss-Newton、Hartley等求導(dǎo)類(lèi)算法和單純形算法、遺傳算法等非求導(dǎo)類(lèi)算法[3]。對(duì)現(xiàn)實(shí)研究而言，經(jīng)常面臨臨床實(shí)驗(yàn)數(shù)據(jù)相對(duì)比較少，或者是實(shí)驗(yàn)條件昂貴，甚至是實(shí)驗(yàn)情況不允許大規(guī)模進(jìn)行的情況。由于貝葉斯估計(jì)方法不僅能夠充分結(jié)合模型、數(shù)據(jù)的信息和參數(shù)的先驗(yàn)信息，而且還能對(duì)缺失數(shù)據(jù)進(jìn)行簡(jiǎn)明化處理，因此，貝葉斯估計(jì)方法在藥物動(dòng)力學(xué)研究中具有無(wú)語(yǔ)倫比的優(yōu)勢(shì)，成為藥物動(dòng)力學(xué)參數(shù)估計(jì)的新興方法[4]。

馬爾科夫鏈蒙特卡羅(MCMC)方法的提出解決了貝葉斯估計(jì)中復(fù)雜的后驗(yàn)概率分布積分的數(shù)值求解問(wèn)題，從而將貝葉斯統(tǒng)計(jì)方法從理論拓展到了藥物動(dòng)力學(xué)、水動(dòng)力學(xué)等非線性高維度模型的應(yīng)用。為適應(yīng)更加復(fù)雜的現(xiàn)實(shí)領(lǐng)域問(wèn)題的模型參數(shù)求解，MCMC方法在基本的Gibbs抽樣和Metropolis-Hastings(MH)抽樣算法的基礎(chǔ)上，發(fā)展演化出切變抽樣(Slice Sampling)算法、適應(yīng)性Metropolis(AM)算法、延遲拒絕(DR)算法、延遲拒絕適應(yīng)性Metropolis(DRAM)算法等。DRAM耦合了DR和AM的優(yōu)點(diǎn)，被認(rèn)為在處理非線性、參數(shù)高維、參數(shù)高度相關(guān)的模型，或者是在目標(biāo)分布遠(yuǎn)遠(yuǎn)偏離于高斯分布時(shí)有突出的優(yōu)勢(shì)[5]。目前，DRAM在國(guó)外和國(guó)內(nèi)的藥物動(dòng)力學(xué)等領(lǐng)域應(yīng)用研究都比較少[6]。

本文以非線性藥物代謝動(dòng)力學(xué)模型未知參數(shù)的貝葉斯估計(jì)為目標(biāo)，結(jié)合DRAM算法實(shí)現(xiàn)模型反演并獲得良好的藥物代謝規(guī)律的模擬效果。同時(shí)，對(duì)比DRAM這一新的抽樣方法與MH進(jìn)行對(duì)比，為高度非線性、高維模型參數(shù)估計(jì)提供參考標(biāo)準(zhǔn)。

2 藥物代謝動(dòng)力學(xué)模型的貝葉斯參數(shù)估計(jì)

2.1 藥物動(dòng)力學(xué)模型

假設(shè)藥物的藥代動(dòng)力學(xué)滿足口服單隔室分布，且模型是非線性吸收動(dòng)力學(xué)及一級(jí)消除動(dòng)力學(xué)過(guò)程：

(1)

(2)

其中，Q(mg)代表口服藥物的藥量，C(mg/l)為血漿中的藥物濃度，Vmax(mg/min)為體內(nèi)藥物最大轉(zhuǎn)運(yùn)速率常數(shù)，Km稱(chēng)為米氏常數(shù)，CL(l/min)為清除速率常數(shù)，Vd(l)為表觀分布容積，它是聯(lián)系著體內(nèi)藥量Q與血藥濃度C的橋梁?；颊叻盟幬锖竽骋粫r(shí)刻的體內(nèi)藥量無(wú)法知道，但是血藥濃度可以根據(jù)即時(shí)抽取的血樣測(cè)定出來(lái)，兩者聯(lián)系的數(shù)學(xué)公式表示為Vd=Q/C。

藥物代謝模型的形式多樣。研究表明，許多新的類(lèi)別的模型如隨機(jī)微分模型，分?jǐn)?shù)階模型等的改進(jìn)和確定都可以基于以上典型常微分藥物代謝動(dòng)力學(xué)模型進(jìn)行。因此，不失一般性，我們的計(jì)算模型將以此為例，利用實(shí)測(cè)的時(shí)間-血藥濃度的數(shù)據(jù)構(gòu)建有效合理的貝葉斯算法，以獲取藥物代謝動(dòng)力學(xué)特征參數(shù)θ=(Vmax,Km,V,CL)，為新藥的研制、劑量的確定和給藥方案的設(shè)計(jì)提供理論依據(jù)。在計(jì)算中，記初始藥量為Q0=Q(0)，初始血藥濃度為C0=C(0)。

2.2 模型參數(shù)估計(jì)的貝葉斯算法

為獲取模型特征參數(shù)θ=(Vmax,Km,V,CL)的信息，根據(jù)貝葉斯原理，可以知道參數(shù)的后驗(yàn)分布πpost(θ|Y)滿足以下計(jì)算公式：

(3)

其中π(θ)為參數(shù)θ的先驗(yàn)分布，L(Y|θ)為給定樣本觀察數(shù)據(jù)Y的似然函數(shù)。由貝葉斯算法的基本原理，參數(shù)θ=(Vmax,Km,V,CL)的識(shí)別問(wèn)題將轉(zhuǎn)化為下述給定觀測(cè)數(shù)據(jù)Y的統(tǒng)計(jì)反問(wèn)題。

記前向藥動(dòng)力學(xué)模型式(1)的解為f(θ)，假設(shè)觀測(cè)導(dǎo)致的誤差模型表示為

(4)

(5)

根據(jù)模型(4)，給定n個(gè)獨(dú)立同分布樣本觀察數(shù)據(jù)Y=(y1,y2,…,yn)'的似然函數(shù)L(Y|θ,σ2)的似然函數(shù)可以表示為：

(6)

由貝葉斯定理(3)，可以推導(dǎo)出所有參數(shù)的后驗(yàn)分布為

2.3 DRAM抽樣算法

一旦獲得了所有參數(shù)的后驗(yàn)分布，就可以獲得參數(shù)的后驗(yàn)特征。然而一般情況下是無(wú)法得到其相關(guān)的解析形式，因此我們通常借助MCMC抽樣算法來(lái)實(shí)現(xiàn)復(fù)雜積分的數(shù)值求解。經(jīng)典的MH算法計(jì)算步驟是[5,7]：

1)給出參數(shù)初始值θ(0)和建議分布q；

2)從建議分布q(θ(i+1)|θ(i))中抽取變量θ(i+1)；

3)以概率α(θ(i+1),θ(i))接受θ(i+1)，其中

(7)

4)重復(fù)步驟2，直至獲取足夠的抽樣。

顯然，M-H算法中建議分布一旦設(shè)定變保持不變，研究證明建議分布會(huì)影響馬爾科夫鏈的收斂速度，建議分布越接近目標(biāo)分布，則馬爾科夫鏈?zhǔn)諗吭娇靃1]。馬爾科夫鏈的收斂速度在高非線性和高維參數(shù)估計(jì)模型中顯得尤為重要。在AR算法中采取不斷更新建議分布的協(xié)方差的方式，使其不斷接近目標(biāo)分布。DR算法采取延遲拒絕抽樣值的方式來(lái)降低拒絕概率，為了提高效率，Haario等吸收AR和DR的抽樣思想，設(shè)計(jì)提出了DRAM算法[1]。以高斯型建議分布為例，基于藥物動(dòng)力學(xué)模型(1)(2)的貝葉斯參數(shù)估計(jì)的DRAM抽樣算法具體步驟如下：

2)DR循環(huán)抽樣。從建議分布qi(θ*|θ(i))中抽取新的參數(shù)θ*，以第1層的接受概率來(lái)接受θ*；

3)如果接受θ*，令θ(i+1)=θ*；如果拒絕θ*，令θ(i+1)=θ(i)；

4)如果循環(huán)次數(shù)i>n0，提出新的建議分布的協(xié)方差C(i)，從而得到更高層的建議分布q，并根據(jù)新的建議分布。

5)重復(fù)步驟2，直至獲取足夠的抽樣。

3 模型的貝葉斯模擬

在藥物動(dòng)力學(xué)模型的參數(shù)估計(jì)和DRAM抽樣算法的模擬中，選取文獻(xiàn)中的臨床試驗(yàn)數(shù)據(jù)[8]取參數(shù)的初始值為

CL=0.05(l/min)，Q0=5(mg)，C0=0(mg/l)

樣本的時(shí)間點(diǎn)取值為1,2,4,6,8,10,15,20,30,40,50,60,70,80,90, 100,150,200,250,300；相應(yīng)的血藥濃度為0.018081，0.035952，0.071065，0.10535，0.1388，0.17145，0.24949，0.32249，0.45305，0.56133，0.64181，0.68151，0.6663，0.61523，0.55859，0.50569，0.3067，0.186，0.1128，0.068409。根據(jù)上述貝葉斯算法，進(jìn)行了10000次迭代，為減少統(tǒng)計(jì)誤差，舍棄前1000次迭代結(jié)果。圖1和圖2分別給出了MH抽樣算法和DRAM抽樣算法的迭代過(guò)程，顯然，DRAM算法相對(duì)于MH算法收斂更快，具有更好的混合性。如圖3所示，DRAM算法取得良好的參數(shù)估計(jì)效果。為更細(xì)致考察兩種抽樣算法生成的馬爾科夫鏈的收斂性，表1和表2給出了兩種算法相關(guān)統(tǒng)計(jì)指標(biāo)的結(jié)果。從參數(shù)的均值、標(biāo)準(zhǔn)差和估計(jì)誤差來(lái)看， MH算法和DRAM算法都取得了比較好的估計(jì)效果。從Geweke檢驗(yàn)方法結(jié)果來(lái)看，兩種方法的馬爾科夫鏈均已收斂，收斂程度有微小的差異。從序列的自相關(guān)時(shí)間間隔(τ)來(lái)看，DRAM方法顯然更有效率。以參數(shù)V的τ值來(lái)看，大約每間隔15個(gè)取樣點(diǎn)能取得獨(dú)立同分布的樣本，而MH算法大約需要每隔20個(gè)點(diǎn)采樣才能取得獨(dú)立同分布的樣本。

圖1 MH抽樣算法參數(shù)估計(jì)的迭代過(guò)程和誤差

圖2 DRAM抽樣算法參數(shù)估計(jì)的迭代過(guò)程和誤差

表1 MH算法模擬的統(tǒng)計(jì)數(shù)據(jù)

參數(shù)均值標(biāo)準(zhǔn)差估計(jì)誤差τGewekeVmax0．110610．0039250．0001272215．0410．99615Km0．764050．137410．004920716．5880．97035V4．94680．0795320．0031120．6570．99522CL0．0502060．000727093．0009e?00515．9720．99817

表2 DRAM算法模擬的統(tǒng)計(jì)數(shù)據(jù)

圖3 DRAM抽樣算法模型參數(shù)擬合效果

4 結(jié)論

本文針對(duì)非線性藥物動(dòng)力學(xué)模型的參數(shù)估計(jì)問(wèn)題，給出了MCMC的DRAM算法，仿真結(jié)果表明：MCMC算法能夠很好解決藥物動(dòng)力學(xué)模型的參數(shù)估計(jì)問(wèn)題，為藥理分析和藥物臨床試驗(yàn)提供較好的參考，藥物動(dòng)力學(xué)模型擬合結(jié)果和實(shí)測(cè)結(jié)果擬合效果良好，說(shuō)明了用MCMC及其相關(guān)方法應(yīng)用于藥物動(dòng)力學(xué)模型的有效性；結(jié)果顯示DRAM算法具有更好的靈活性和高效性，在收斂速度和混合性方面都具有明顯的優(yōu)勢(shì)，適合用在高度非線性和高維參數(shù)估計(jì)模型中。因此，基于DRAM

算法的非線性藥物動(dòng)力學(xué)模型反演結(jié)果能在現(xiàn)實(shí)應(yīng)用中具有重要的價(jià)值。

[1] 張萍, 陸濤, 言方榮等. 基于隨機(jī)微分方程的群體藥物代謝動(dòng)力學(xué)參數(shù)的極大似然估計(jì) [J]. 藥學(xué)研究, 2013, 26(3), 304-309.

[2] 曾文藝,孫曉穎. 藥物動(dòng)力學(xué)房室模型的改進(jìn)[J]. 北京師范大學(xué)學(xué)報(bào)(自然科學(xué)版), 2010, 46(5): 541-545.

[3] 蔡曄,孫春萌,沈雁. 非線性藥物動(dòng)力學(xué)參數(shù)的計(jì)算方法研究進(jìn)展[J]. 藥學(xué)研究, 2014, 33(7), 401-404.

[4] Fanni Natanegara,Beat Neuenschwander,John W, Seamen et al. The current state of Bayesian methods in medical product development: survey results and recommendations from the DIA Bayesian scientific working group[J]. Pharmaceutical Statistics, 2014, 13：3-12.

[5] Marko Lain. Adaptive MCMC methods with applications in environmental and geophysical models[D]. Finland, 2008,12-14.

[6] 張慶慶, 許月萍, 張徐杰等. 基于DRAM的水質(zhì)模擬不確定分析和風(fēng)險(xiǎn)決策[J]. 浙江大學(xué)學(xué)報(bào)(工學(xué)版), 2012, 46(12): 2231-2236,2242.

[7] Zuev K M, Katafygiotis I S. Modified Metropolis-Hastings algorithm with delayed rejection[J]. Probability Engineering Mechanics,2011, 26(3): 405-412.

[8] 褚娟. 利用隨機(jī)微分方程建立藥物代謝動(dòng)力學(xué)模型及其應(yīng)用[D]. 華中科技大學(xué), 2012.

責(zé)任編輯：劉海濤

Application of Bayesian Inference to Estimate the Parametersin Pharmacokinetics Model

Yu Xiaomei, Zhang Haiyong, Shen Yongchang

To overcome the parameter estimation's difficulty from the measurement uncertainty and high dimensions of unknown parameters, mathematical model of posterior probability distribution based on Bayesian inference was set up to estimate the parameters in pharmacokinetics model. Markov chain Monte Carlo(MCMC) method is used to explore the posterior distribution to get the estimation. Delay Rejection Adaptive MCMC is used during sampling. The results show that DRAM method could give a better estimate results comparison to traditional sampling MCMC method such as Metropolis Hastings(MH). DRAM method has been proved to be more stable and effective than MH for pharmacokinetics parameters calculation.

Bayesian inference; Pharmacokinetics model; DRAM; MCMC

O212.8;R969.1

A

1673-1794(2017)02-0004-04

余曉美，滁州學(xué)院數(shù)學(xué)與金融學(xué)院講師，博士，研究方向：應(yīng)用統(tǒng)計(jì)；張海永，沈永昌，滁州學(xué)院數(shù)學(xué)與金融學(xué)院(安徽 滁州 239000)。

安徽省人才基金項(xiàng)目(1408085QA06，KJ2014A180)；滁州學(xué)院科研啟動(dòng)項(xiàng)目(2012qd04)；滁州學(xué)院科學(xué)研究項(xiàng)目(2015GH34)

2016-12-07