基于隨機(jī)微分方程的藥物動(dòng)力學(xué)建模新方法＊

褚 娟， 萬建平， 高秀娟， 陳 匯△

1 華中科技大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，武漢 430074

2 華中科技大學(xué)同濟(jì)醫(yī)學(xué)院基礎(chǔ)醫(yī)學(xué)院藥理學(xué)系，武漢 430030

藥物動(dòng)力學(xué)是一門綜合學(xué)科，由最初的采用數(shù)學(xué)分析手段研究藥物的體內(nèi)過程，現(xiàn)已發(fā)展成為運(yùn)用數(shù)學(xué)、計(jì)算機(jī)、生理學(xué)與生命科學(xué)等多學(xué)科知識(shí)和技術(shù)來研究藥物體內(nèi)過程的學(xué)科。

在目前的藥物動(dòng)力學(xué)建模中，較多的是利用常微分方程描述藥物在體內(nèi)的動(dòng)力學(xué)過程，此方法在描述簡單、線性的藥物體內(nèi)過程時(shí)確實(shí)具有簡易、合理等優(yōu)勢(shì)，但在處理一些更為復(fù)雜的非線性藥物動(dòng)力學(xué)模型時(shí)就遇到了一些問題。

首先，利用常微分方程建模時(shí)，極易出現(xiàn)模型過于簡單或過于復(fù)雜的情況。若模型方程過于簡單，不能真實(shí)地反映藥物的體內(nèi)過程，就會(huì)產(chǎn)生有偏的參數(shù)估計(jì)值和相關(guān)的殘差誤差［1］；若模型方程過于復(fù)雜，如待估參數(shù)過多，就要求大量數(shù)據(jù)用于參數(shù)估計(jì)，而臨床數(shù)據(jù)是有限的。

其次，在常微分方程建模中，常假設(shè)某些重要的藥理參數(shù)（如吸收、消除速率）是常數(shù)，而實(shí)際上，它們可能是隨時(shí)間變化的，這樣也會(huì)產(chǎn)生相關(guān)的殘差誤差。

再次，利用常微分方程建立藥物動(dòng)力學(xué)方程時(shí)只考慮了不相關(guān)的殘差，如實(shí)驗(yàn)中產(chǎn)生的誤差，測(cè)量誤差等［2］，而模型結(jié)構(gòu)設(shè)定的不足，模型中參數(shù)隨時(shí)間變化，未知的隨機(jī)干擾等都可能引起隨時(shí)間變化的相關(guān)殘差，忽略這些相關(guān)的殘差直接影響了參數(shù)的估計(jì)，最終導(dǎo)致模型的缺陷。所以，如果只把模型中出現(xiàn)的誤差簡單歸結(jié)為測(cè)量誤差顯然是不合理，不科學(xué)的，這是受研究方法所限。

經(jīng)過不斷地探索、研究，在藥物動(dòng)力學(xué)模型中引入隨機(jī)微分方程理論較好地解決了這些問題。

首先，隨機(jī)微分方程模型將個(gè)體內(nèi)誤差分為兩部分：一部分是不相關(guān)的測(cè)量誤差，主要由試驗(yàn)誤差、測(cè)量誤差等造成的；另一部分是系統(tǒng)誤差，由可能存在的模型定義的不合理或體內(nèi)存在的隨機(jī)干擾引起的［3］。

其次，對(duì)于過于簡單的模型，引入隨機(jī)微分方程可以解釋給定數(shù)據(jù)出現(xiàn)的大量變異；對(duì)于過于復(fù)雜的模型，在參數(shù)估計(jì)時(shí)會(huì)有很大困難，隨機(jī)微分方程可簡化模型并優(yōu)化參數(shù)估計(jì)。

再次，隨機(jī)微分方程中擴(kuò)散項(xiàng)參數(shù)的引入，不僅考慮到了模型的不合理性而且可以對(duì)其進(jìn)行定量和定位描述，進(jìn)而很大程度地改進(jìn)原有的模型［2］。

因此，本文也是從最初建立的簡單的藥物動(dòng)力學(xué)的線性模型出發(fā)，逐步得到一個(gè)既盡可能真實(shí)地模擬了藥物在體內(nèi)變化情況又不至太復(fù)雜的藥物動(dòng)力學(xué)模型，并對(duì)建模方法給予了理論上的支持和解釋。新模型在揭示藥物動(dòng)力學(xué)的非線性等復(fù)雜過程的本質(zhì)時(shí)體現(xiàn)了巨大的優(yōu)勢(shì)。

1 基于常微分方程（ODEs）的藥動(dòng)學(xué)線性模型

這里主要針對(duì)常見的幾種給藥方式來建立基于常微分方程的藥動(dòng)學(xué)模型［4］。目前，大多數(shù)藥物建模時(shí)首選的模型通常如表1所示。

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

以口服給藥為例，建立藥物動(dòng)力學(xué)模型，通常假設(shè)為一級(jí)吸收和一級(jí)消除速率過程：

a.一室模型（圖1）

圖1 一室模型結(jié)構(gòu)Fig.1 The structure of one apartment model

吸收部位的藥動(dòng)學(xué)過程：

體內(nèi)（血液內(nèi)）的藥動(dòng)學(xué)過程：

b.二室模型（圖2）

圖2 二室模型結(jié)構(gòu)Fig.2 The structure of two apartment model

吸收部位的藥動(dòng)學(xué)過程：

體內(nèi)藥動(dòng)學(xué)過程：

在實(shí)際問題中，許多藥物的吸收、消除并不是簡單的線性變化，而是具有非線性吸收、消除甚至是混合吸收、消除的特征。其中混合吸收、混合消除是指藥物在吸收或消除時(shí)不僅具有線性過程而且具有非線性過程。

2 基于隨機(jī)微分方程（SDEs）的藥動(dòng)學(xué)模型

當(dāng)線性藥動(dòng)學(xué)模型不能擬合實(shí)驗(yàn)數(shù)據(jù)時(shí)，就可以通過引入隨機(jī)微分方程理論，在原有模型基礎(chǔ)上加上一個(gè)隨機(jī)擾動(dòng)項(xiàng)來避免由于最初假設(shè)錯(cuò)誤導(dǎo)致后續(xù)的研究錯(cuò)誤，以建立更合理的藥物動(dòng)力學(xué)模型。

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

2.1 建立傳統(tǒng)的線性PK 模型

通過探索性數(shù)據(jù)分析和已知的生理學(xué)知識(shí)建立關(guān)于狀態(tài)變量的常微分方程，即為式（1.1）～式（1.4）的形式，一般記為：

其中，Xt∈Rn為狀態(tài)向量（如血藥濃度），ut∈Rm為輸入變量（如輸入速率），θ∈Rp為參數(shù)向量，f（·）∈Rn為向量函數(shù)。

2.2 引入干擾建立SDEs模型

其中，f（·）為漂移項(xiàng)，σ（·）為擴(kuò)散項(xiàng)，σ（·）∈Rn×n為矩陣函數(shù)，Bt為n維標(biāo)準(zhǔn)Wiener過程，且Bt1－Bt2～N（0，｜t1－t2｜Ι），Ι為單位矩陣。

其中，Wiener過程可理解為來源于同分布的隨機(jī)事件的和產(chǎn)生的Xt的真實(shí)變化與f（·）描述的Xt的變化之差。換言之，擴(kuò)散項(xiàng)σ（·）可解釋模型的不確定性［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é)果均一樣。

2.3 建立測(cè)量方程

其中，yk∈Rl為測(cè)量得到的輸出變量（如血液中藥物濃度），h（·）∈Rl為向量函數(shù)，tk，k＝0…N為抽樣時(shí)間，ek為l維白噪聲過程，即ek～N（0，S（θ）），S（θ）為協(xié)方差矩陣。

2.4 參數(shù)估計(jì)判別擴(kuò)散項(xiàng)的顯著性

根據(jù)極大似然估計(jì)原理，目前常結(jié)合EKF（Extended Kalman filter）等方法進(jìn)行參數(shù)估計(jì)［2］，并通過軟件CTSM、R 等軟件實(shí)現(xiàn)。但是，EKF 使用的前提是假設(shè)觀測(cè)值的條件密度可以用Gaussian密度近似；當(dāng)f（·）為關(guān)于Xt的非線性函數(shù)時(shí)，EKF是利用f（·）的一階Taylor展式求得Gaussian分布的均值、方差的近似值來計(jì)算，基于此得到的是一個(gè)近似的極大似然函數(shù)。另外，EKF作為一種迭代算法，其停止規(guī)則未考慮估計(jì)的效果和效率，而是完全由數(shù)據(jù)個(gè)數(shù)決定的。

目前，國外有一些文獻(xiàn)［6－7］利用由EM 算法衍生出的SAEM－MCMC方法求相關(guān)參數(shù)，此算法可依靠Monolix軟件實(shí)現(xiàn)。SAEM－MCMC 算法相對(duì)于其他的算法（線性化方法，積分近似及普通的EM算法）計(jì)算速度更快，精度更高；其對(duì)初始值的設(shè)定要求較低，對(duì)初始值的選取非常穩(wěn)健，通過若干步迭代之后，基本上收斂到要模擬的平穩(wěn)分布；對(duì)觀測(cè)值的條件密度沒有特別要求，因此適用范圍更廣。

2.5 求解隨機(jī)微分方程

擴(kuò)展模型狀態(tài)方程，指出模型的不足。反復(fù)步驟2.4直至模型的系統(tǒng)誤差顯著地為零為止。

若σ顯著為零，則說明新得到的藥物動(dòng)力學(xué)模型較好地解釋了藥物的體內(nèi)變化情況；若σ不顯著為零，則說明前面建立的ODEs模型不能很好地模擬藥物體內(nèi)變化的實(shí)際情況，存在著某些無法給出合理解釋的σ（ut，t，θ）dωt。這部分誤差和干擾，從建模的速率角度分析，可能是有關(guān)吸收、消除速率過程的假設(shè)有不合理之處。因此就面臨著如何求解SDEs的問題。對(duì)于求解SDEs，我們可以從以下3方面考慮：

（1）直接利用Ito積分求精確解。雖然直接利用Ito積分求得的解最準(zhǔn)確，但只有極少部分特殊的SDEs可以求出精確解，而大部分SDEs不可直接求解。

（2）求SDEs的數(shù)值解。SDEs數(shù)值解只是用隨機(jī)過程無數(shù)條可能的軌道中的一個(gè)軌道表示過程的特征參數(shù)的性質(zhì)［8］，過于特殊而缺乏普遍性；目前，多用隨機(jī)微分方程在固定時(shí)刻附近的隨機(jī)Taylor展開或Runge－Kutta隨機(jī)差分方程來求解SDEs的數(shù)值解，這種采用隨機(jī)過程的一段相當(dāng)長時(shí)間的數(shù)據(jù)來估值的方法對(duì)于反映藥物在各個(gè)時(shí)刻的濃度情況不是十分理想；另外與解ODEs相比，求隨機(jī)微分方程的數(shù)值解仍是較復(fù)雜、費(fèi)時(shí)的。

（3）利用ODEs的解去近似SDEs的解。這種近似不能隨意確定，在一定條件下才能執(zhí)行。

因此，采用擴(kuò)展模型狀態(tài)方程來合理地修正吸收、消除過程使σ（·）顯著為零，此時(shí)就可利用常微分方程的解去近似隨機(jī)微分方程的解。

此方法的優(yōu)點(diǎn)在于，ODEs近似SDEs的解是在均方意義下的近似，因而精度可以得到保證；另外，求解常微分方程的理論和方法已十分成熟，為計(jì)算帶來很大便利?？傊?，此方法在不同程度上彌補(bǔ)了方法（1）、（2）的不足，是十分有意義的。

2.6 用微分方程模型解釋藥動(dòng)學(xué)過程

解常微分方程dXt＝f＊（Xt，ut，t，θ）dt，得到藥物濃度隨時(shí)間變化表達(dá)式，用于預(yù)測(cè)體內(nèi)藥物濃度變化，從而改進(jìn)給藥方式、給藥時(shí)間間隔，使藥濃度維持在治療窗內(nèi)，以優(yōu)化藥效并防止毒副作用。

3 結(jié)語

本文從速率論的角度，在傳統(tǒng)的藥物動(dòng)力學(xué)線性模型的基礎(chǔ)上，利用隨機(jī)微分方程完善復(fù)雜藥物動(dòng)力學(xué)模型，使其更加接近真實(shí)地描述藥物在體內(nèi)的變化情況，由此得到的血藥濃度隨時(shí)間變化規(guī)律，對(duì)臨床藥物研究更具有積極意義。但本文只是對(duì)基于隨機(jī)微分方程的藥動(dòng)學(xué)模型給出了理論上的合理解釋、建模步驟和求解方法，而沒有利用真實(shí)數(shù)據(jù)按照隨機(jī)微分方程建模的步驟通過統(tǒng)計(jì)軟件模擬并建立一個(gè)藥物動(dòng)力學(xué)模型，這些問題有待進(jìn)一步的研究。

［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)新藥物藥代動(dòng)力學(xué)研究與評(píng)價(jià)［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］ 龔光魯.隨機(jī)微分方程及其應(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.