褚 娟, 萬建平, 高秀娟, 陳 匯△
1 華中科技大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院,武漢 430074
2 華中科技大學(xué)同濟醫(yī)學(xué)院基礎(chǔ)醫(yī)學(xué)院藥理學(xué)系,武漢 430030
藥物動力學(xué)是一門綜合學(xué)科,由最初的采用數(shù)學(xué)分析手段研究藥物的體內(nèi)過程,現(xiàn)已發(fā)展成為運用數(shù)學(xué)、計算機、生理學(xué)與生命科學(xué)等多學(xué)科知識和技術(shù)來研究藥物體內(nèi)過程的學(xué)科。
在目前的藥物動力學(xué)建模中,較多的是利用常微分方程描述藥物在體內(nèi)的動力學(xué)過程,此方法在描述簡單、線性的藥物體內(nèi)過程時確實具有簡易、合理等優(yōu)勢,但在處理一些更為復(fù)雜的非線性藥物動力學(xué)模型時就遇到了一些問題。
首先,利用常微分方程建模時,極易出現(xiàn)模型過于簡單或過于復(fù)雜的情況。若模型方程過于簡單,不能真實地反映藥物的體內(nèi)過程,就會產(chǎn)生有偏的參數(shù)估計值和相關(guān)的殘差誤差[1];若模型方程過于復(fù)雜,如待估參數(shù)過多,就要求大量數(shù)據(jù)用于參數(shù)估計,而臨床數(shù)據(jù)是有限的。
其次,在常微分方程建模中,常假設(shè)某些重要的藥理參數(shù)(如吸收、消除速率)是常數(shù),而實際上,它們可能是隨時間變化的,這樣也會產(chǎn)生相關(guān)的殘差誤差。
再次,利用常微分方程建立藥物動力學(xué)方程時只考慮了不相關(guān)的殘差,如實驗中產(chǎn)生的誤差,測量誤差等[2],而模型結(jié)構(gòu)設(shè)定的不足,模型中參數(shù)隨時間變化,未知的隨機干擾等都可能引起隨時間變化的相關(guān)殘差,忽略這些相關(guān)的殘差直接影響了參數(shù)的估計,最終導(dǎo)致模型的缺陷。所以,如果只把模型中出現(xiàn)的誤差簡單歸結(jié)為測量誤差顯然是不合理,不科學(xué)的,這是受研究方法所限。
經(jīng)過不斷地探索、研究,在藥物動力學(xué)模型中引入隨機微分方程理論較好地解決了這些問題。
首先,隨機微分方程模型將個體內(nèi)誤差分為兩部分:一部分是不相關(guān)的測量誤差,主要由試驗誤差、測量誤差等造成的;另一部分是系統(tǒng)誤差,由可能存在的模型定義的不合理或體內(nèi)存在的隨機干擾引起的[3]。
其次,對于過于簡單的模型,引入隨機微分方程可以解釋給定數(shù)據(jù)出現(xiàn)的大量變異;對于過于復(fù)雜的模型,在參數(shù)估計時會有很大困難,隨機微分方程可簡化模型并優(yōu)化參數(shù)估計。
再次,隨機微分方程中擴散項參數(shù)的引入,不僅考慮到了模型的不合理性而且可以對其進行定量和定位描述,進而很大程度地改進原有的模型[2]。
因此,本文也是從最初建立的簡單的藥物動力學(xué)的線性模型出發(fā),逐步得到一個既盡可能真實地模擬了藥物在體內(nèi)變化情況又不至太復(fù)雜的藥物動力學(xué)模型,并對建模方法給予了理論上的支持和解釋。新模型在揭示藥物動力學(xué)的非線性等復(fù)雜過程的本質(zhì)時體現(xiàn)了巨大的優(yōu)勢。
這里主要針對常見的幾種給藥方式來建立基于常微分方程的藥動學(xué)模型[4]。目前,大多數(shù)藥物建模時首選的模型通常如表1所示。
表1 藥物建模方式Table 1 The modeling method of drug
以口服給藥為例,建立藥物動力學(xué)模型,通常假設(shè)為一級吸收和一級消除速率過程:
a.一室模型(圖1)
圖1 一室模型結(jié)構(gòu)Fig.1 The structure of one apartment model
吸收部位的藥動學(xué)過程:
體內(nèi)(血液內(nèi))的藥動學(xué)過程:
b.二室模型(圖2)
圖2 二室模型結(jié)構(gòu)Fig.2 The structure of two apartment model
吸收部位的藥動學(xué)過程:
體內(nèi)藥動學(xué)過程:
在實際問題中,許多藥物的吸收、消除并不是簡單的線性變化,而是具有非線性吸收、消除甚至是混合吸收、消除的特征。其中混合吸收、混合消除是指藥物在吸收或消除時不僅具有線性過程而且具有非線性過程。
當線性藥動學(xué)模型不能擬合實驗數(shù)據(jù)時,就可以通過引入隨機微分方程理論,在原有模型基礎(chǔ)上加上一個隨機擾動項來避免由于最初假設(shè)錯誤導(dǎo)致后續(xù)的研究錯誤,以建立更合理的藥物動力學(xué)模型。
利用SDEs建立PK 模型的具體過程如下[1]。
通過探索性數(shù)據(jù)分析和已知的生理學(xué)知識建立關(guān)于狀態(tài)變量的常微分方程,即為式(1.1)~式(1.4)的形式,一般記為:
其中,Xt∈Rn為狀態(tài)向量(如血藥濃度),ut∈Rm為輸入變量(如輸入速率),θ∈Rp為參數(shù)向量,f(·)∈Rn為向量函數(shù)。
其中,f(·)為漂移項,σ(·)為擴散項,σ(·)∈Rn×n為矩陣函數(shù),Bt為n維標準Wiener過程,且Bt1-Bt2~N(0,|t1-t2|Ι),Ι為單位矩陣。
其中,Wiener過程可理解為來源于同分布的隨機事件的和產(chǎn)生的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無關(guān),所以無論是Ito積分還是Stratonovich 積分,最終的結(jié)果均一樣。
其中,yk∈Rl為測量得到的輸出變量(如血液中藥物濃度),h(·)∈Rl為向量函數(shù),tk,k=0…N為抽樣時間,ek為l維白噪聲過程,即ek~N(0,S(θ)),S(θ)為協(xié)方差矩陣。
根據(jù)極大似然估計原理,目前常結(jié)合EKF(Extended Kalman filter)等方法進行參數(shù)估計[2],并通過軟件CTSM、R 等軟件實現(xiàn)。但是,EKF 使用的前提是假設(shè)觀測值的條件密度可以用Gaussian密度近似;當f(·)為關(guān)于Xt的非線性函數(shù)時,EKF是利用f(·)的一階Taylor展式求得Gaussian分布的均值、方差的近似值來計算,基于此得到的是一個近似的極大似然函數(shù)。另外,EKF作為一種迭代算法,其停止規(guī)則未考慮估計的效果和效率,而是完全由數(shù)據(jù)個數(shù)決定的。
目前,國外有一些文獻[6-7]利用由EM 算法衍生出的SAEM-MCMC方法求相關(guān)參數(shù),此算法可依靠Monolix軟件實現(xiàn)。SAEM-MCMC 算法相對于其他的算法(線性化方法,積分近似及普通的EM算法)計算速度更快,精度更高;其對初始值的設(shè)定要求較低,對初始值的選取非常穩(wěn)健,通過若干步迭代之后,基本上收斂到要模擬的平穩(wěn)分布;對觀測值的條件密度沒有特別要求,因此適用范圍更廣。
擴展模型狀態(tài)方程,指出模型的不足。反復(fù)步驟2.4直至模型的系統(tǒng)誤差顯著地為零為止。
若σ顯著為零,則說明新得到的藥物動力學(xué)模型較好地解釋了藥物的體內(nèi)變化情況;若σ不顯著為零,則說明前面建立的ODEs模型不能很好地模擬藥物體內(nèi)變化的實際情況,存在著某些無法給出合理解釋的σ(ut,t,θ)dωt。這部分誤差和干擾,從建模的速率角度分析,可能是有關(guān)吸收、消除速率過程的假設(shè)有不合理之處。因此就面臨著如何求解SDEs的問題。對于求解SDEs,我們可以從以下3方面考慮:
(1)直接利用Ito積分求精確解。雖然直接利用Ito積分求得的解最準確,但只有極少部分特殊的SDEs可以求出精確解,而大部分SDEs不可直接求解。
(2)求SDEs的數(shù)值解。SDEs數(shù)值解只是用隨機過程無數(shù)條可能的軌道中的一個軌道表示過程的特征參數(shù)的性質(zhì)[8],過于特殊而缺乏普遍性;目前,多用隨機微分方程在固定時刻附近的隨機Taylor展開或Runge-Kutta隨機差分方程來求解SDEs的數(shù)值解,這種采用隨機過程的一段相當長時間的數(shù)據(jù)來估值的方法對于反映藥物在各個時刻的濃度情況不是十分理想;另外與解ODEs相比,求隨機微分方程的數(shù)值解仍是較復(fù)雜、費時的。
(3)利用ODEs的解去近似SDEs的解。這種近似不能隨意確定,在一定條件下才能執(zhí)行。
因此,采用擴展模型狀態(tài)方程來合理地修正吸收、消除過程使σ(·)顯著為零,此時就可利用常微分方程的解去近似隨機微分方程的解。
此方法的優(yōu)點在于,ODEs近似SDEs的解是在均方意義下的近似,因而精度可以得到保證;另外,求解常微分方程的理論和方法已十分成熟,為計算帶來很大便利??傊?,此方法在不同程度上彌補了方法(1)、(2)的不足,是十分有意義的。
解常微分方程dXt=f*(Xt,ut,t,θ)dt,得到藥物濃度隨時間變化表達式,用于預(yù)測體內(nèi)藥物濃度變化,從而改進給藥方式、給藥時間間隔,使藥濃度維持在治療窗內(nèi),以優(yōu)化藥效并防止毒副作用。
本文從速率論的角度,在傳統(tǒng)的藥物動力學(xué)線性模型的基礎(chǔ)上,利用隨機微分方程完善復(fù)雜藥物動力學(xué)模型,使其更加接近真實地描述藥物在體內(nèi)的變化情況,由此得到的血藥濃度隨時間變化規(guī)律,對臨床藥物研究更具有積極意義。但本文只是對基于隨機微分方程的藥動學(xué)模型給出了理論上的合理解釋、建模步驟和求解方法,而沒有利用真實數(shù)據(jù)按照隨機微分方程建模的步驟通過統(tǒng)計軟件模擬并建立一個藥物動力學(xué)模型,這些問題有待進一步的研究。
[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] 魏敏吉,趙明.創(chuàng)新藥物藥代動力學(xué)研究與評價[M].北京:北京大學(xué)醫(yī)學(xué)出版社,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] 龔光魯.隨機微分方程及其應(yīng)用概要[M].北京:清華大學(xué)出版社,2008:124-125
[9] Evans L C.An introduction to stochastic differential equations[M].Berkely:Department of Mathematics UC Berkeley,2006:91.