藥物代謝動力學軟件設計（三室模型）

范景行

（甘肅民族師范學院，甘肅 合作 747000）

藥物代謝動力學軟件設計（三室模型）

范景行

（甘肅民族師范學院，甘肅 合作 747000）

本文的研究范疇是三室模型線性藥物動力學，包括一次靜脈推注給藥、一次血管外給藥、靜脈滴注給藥3種模型。首先根據模型列出微分方程組，求出方程組的特解方程，然后將參數整合，根據已知數據求出參數值，也就是擬合出一條曲線。擬合曲線使用最小二乘法和Gauss-Newton迭代法求解，針對奇異矩陣導致Gauss-Newton法迭代結果發散的缺點，引進了Levenberg-Marquardt法，針對殘差平方和偏重高密度數據點的缺點，引進了權重。所有算法都用C語言實現。

藥物代謝動力學；最小二乘法；Gauss-Newton法；Levenberg-Marquardt法

藥物代謝動力學（也稱藥動學、藥物動力學、藥代動力學，PK）是研究藥物及其代謝物在體內吸收、分布、轉化和排泄過程中，隨著時間不同，不斷進行運動變化規律的科學。藥物代謝動力學研究的對象是人體，對患病的機體以及正常人體用藥后，觀察藥物在體內的變化過程。根據該科學可對患者提供安全、有效的治療方案，包括給藥途徑、用藥劑型、用法、用量、給藥間隔等，實現給藥方案個體化；可重新審查給藥計劃，對不良反應做出解釋；可按計劃對正在進行血液、腹膜透析的患者給藥，預防中毒。根據該科學可及時地進行血藥濃度的監測，廣泛地收集藥學情報，為臨床提供科學的給藥方案，進一步提高療效，減少藥物的不良反應。所以該學科所肩負的任務非常重要，其研究是采用數學手段（如用圖像、公式、參數等方法）來描述藥物在體內運動過程中的規律，因而可為用藥方案的設計提供重要的參考依據。

藥物代謝動力學從數學的角度可分為兩類：線性藥物代謝動力學和非線性藥物代謝動力學。在線性藥物代謝動力學中，藥物代謝動力學參數如半衰期與劑量無關；而在非線性藥物代謝動力學中，藥物代謝動力學參數隨劑量（或體內藥物濃度）而變化，如半衰期與劑量有關。本文所討論的內容準確地說屬于線性藥物代謝動力學范疇，并以三室模型為研究對象，對一次靜脈推注給藥、一次血管外給藥、靜脈滴注給藥3種常見的給藥方式做詳細介紹。

1 模型的建立

1.1 一次靜脈推注給藥的三室模型

藥物的三室模型可以這樣設想：假定藥物進入中心室后，逐漸向兩個周邊室轉運，在中央室與周邊室之間藥物進行著可逆的運動，藥物在中央室按一級過程消除[1]，其體內過程模型見圖1。

圖1 一次靜脈推注給藥的三室模型

圖1 顯示：X0為靜脈注射給藥劑量；XC為中央室的藥量；XP和XB為周邊室的藥量；k12和k21為中央室與淺外室之間的運轉速度常數；k13和k31為中央室與深外室之間的運轉速度常數；k10為藥物從中央室消除的一級速度常數；VC、VP、VB分別為中央室、淺外室、深外室的表觀分布容積。

假如藥物的轉運過程均服從一級速度過程，即藥物的轉運速度與該室藥物濃度（或藥量）成正比，則各室藥物的轉運可用下列微分方程組定量描述[2]。

當t=0時，注射的藥物全部存放于中心室內，尚未轉運至其他各室。所以XC（0）=X0，XP（0）=XB（0）=0，這就是以上方程組的初始條件。

因為實際測量的數據普遍都以中央室為主，所以我們只討論中央室的血藥濃度與時間的關系。設C為中央室藥物濃度，通過求解化簡可得需要擬合的數學方程：

1.2 靜脈滴注給藥的三室模型

設劑量為X0的藥物，在τ0這段時間內，以恒速k0=D0/τ0滴入中心室，其體內過程模型[3]見圖2。

圖2 靜脈滴注給藥的三室模型

圖2 顯示：X0為靜脈注射給藥劑量；XC為中央室的藥量；XP和XB為周邊室的藥量；k12和k21為中央室與淺外室之間的運轉速度常數；k13和k31為中央室與深外室之間的運轉速度常數；k10為藥物從中央室消除的一級速度常數；VC、VP、VB分別為中央室、淺外室、深外室的表觀分布容積。

在經過滴注時間t（0≤t≤τ）時，中心室的藥量為XC（t），周邊室的藥量為XP（t）和XB（t）。除滴注是恒速K0之外，如果其余各轉運過程服從一級動力學過程，則各室間的藥物轉運方程為[4]：

當t=0時，XC（0）=XP（0）=XB（0）=0，這就是以上方程組的初始條件。

因為實際測量的數據普遍都以中央室為主，所以我們只討論中央室的血藥濃度與時間的關系。同理可得：

這就是滴注后中心室的藥物濃度與時間關系。其中，是從滴注完成時算起。

1.3 一次血管外給藥的三室模型

有些藥物不能做靜脈推注給藥時，則采用血管外給藥，僅需在靜脈推注給藥的三室模型中心室前增加一個吸收室，讓藥物逐漸從吸收室轉運至中心室，再轉運至兩個周邊室[5]（見圖3）。

圖3 一次血管外給藥的三室模型

圖3 顯示：X0為給藥劑量；Xa為吸收室藥量；XC為中央室的藥量；XP和XB為周邊室的藥量；kа為吸收速度常數；k12和k21為中央室與淺外室之間的運轉速度常數；k13和k31為中央室與深外室之間的運轉速度常數；k10為藥物從中央室消除的一級速度常數；Vа、VC、VP、VB分別為吸收室、中央室、淺外室、深外室的表觀分布容積。

假定吸收過程仍是一級過程，則各室間的藥物運轉可用下列方程表示[6]：

因為實際測量的數據普遍都以中央室為主，所以我們只討論中央室的血藥濃度與時間的關系。同理可得：

2 模型的求解算法

Gauss-Newton 算法步驟歸納如下[7，8]：

（1）取初始近似值A0，允許誤差ε；

（2）計算 F（XY，A0），DF（X，Y，A0），DF（X，Y，A0）T；

（3）解線性方程組：

得△A；

（4）以△ A修正 A0，即計算：A=A0+△ A；

（5）若|△ A|＞ε，則 A0=A，即將 A 作為初值 A0，重復步驟（2）、（3）、（4），直至|△ A|＜ε。

由于Gauss-Newton法采用的是將非線性函數在初值附近做Taylor展開，然后略去二次及二次以上諸項，做簡單線性近似的方法，所以無法避免兩種情況的發生：一是Gauss-Newton法對初值的依賴非常嚴重，如果初值選取不當，往往發散，即越迭代越遠離目標值；二是如果Taylor展開后略去的二次及二次以上諸項的值相對過大，依然會造成Gauss-Newton法不收斂。

克服這個困難的一種途徑是采用Levenberg-Marquardt的方法[9，10]，即加大 DF（X，Y，A0）TDF（X，Y，A0）的主對角元素，而得到改善，將步驟（3）中的方程改寫成：

式中：μ——阻尼因子；

In——n階單位矩陣或DF（X，Y，A0）TDF（X，Y，A0）的主對角元素構成的矩陣（本論文中的算法選擇后者）。相應的算法步驟歸納如下：

（1）取初始近似值A0，計算

（2）取一個適中的μ值，如μ=0.001；

（3）計算 F（X，Y，A0），DF（X，Y，A0），DF（X，Y，A0）T；

（4）解線性方程組[DF（X，Y，A0）TDF（X，Y，A0）+μ In]△A=-DF（X，Y，A0）TF（X，Y，A0），得△ A，并計算 V（X，Y，a1，+V a1a2+V a2，L，an+V an，）；

（5）如果 V （X，Y，a1+V a1，α2+V a2，L，an+V an）≥ V（X，Y，a1，a2，L，an），則將 μ擴大 10 倍（或其他的具有意義的倍數），再回到步驟（4）；

（6）如果V（X，Y，a1+V a1，a2+V a2，L，an+V an）＜V（X，Y，a1，a2，L，an），則將μ縮小10倍（或者其他的具有意義的倍數），將△A+A0代替 A0，再回到步驟（4）；

（7）直至|△ A|＜ε或 V（X，Y，a1，a2，L,an）的變化量＜ε，迭代停止，輸出結果。

[1]劉昌孝，劉定遠.藥物動力學概論[M].北京：中國學術出版社，1984.

[2]梁文權.生物藥劑學與藥物動力學[M].北京：人民衛生出版社，2000.

[3]袁亞湘，孫文瑜.最優化理論與方法[M].北京：科學出版社，2003.

[4]關治，陸金甫.數值分析基礎[M].北京：高等教育出版社，2002.

[5]韋鶴平.最優化技術應用[M].上海：同濟大學出版社，1987.

[6]譚浩強.C程序設計[M].2版.北京：清華大學出版社，2002.

[7]蔣長錦.科學計算和C程序集[M].合肥：中國科學技術大學出版社，1998.

[8]William H P，Saul A T，William T V，et al.C數值算法[M].2 版.北京：電子工業出版社，2004.

[9]張賢達.矩陣分析與應用[M].北京：清華大學出版社，2004.

[10]肖愛玲.Some algorithms of nonlinear least squares[J].數學理論與應用，2004，24（2），86～90.

G431

A

1671-1246（2012）17-0042-03