PK4-1

PBPKモデリング

ProductName ファーマコキネティクス―演習による理解
杉山 雄一,山下 伸二,加藤 基浩
南山堂 / ¥ 6,300 ()
在庫あり。

i.v.

library(odesolve)
params = c(Vh=44,Vr=8,Qh=58,Qr=45.6,Vsys=1000,Kph=1,Kpr=1,RB=1.5,fp=0.1,Fa=1,Vmax=1000,Km=0.1,GFR=6)
times=c(0,(1:240))

dydt <- function(t,y,p){
Vh <- p['Vh']
Vr <- p['Vr']
Qh <- p['Qh']
Qr <- p['Qr']
Vsys <- p['Vsys']
Kph <- p['Kph']
Kpr <- p['Kpr']
RB <- p['RB']
fp <- p['fp']
Fa <- p['Fa']
Vmax <- p['Vmax']
Km <- p['Km']
GFR <- p['GFR']
fb <- fp/RB
Cb <- (Qh*y[2]/Kph + Qr*y[3]/Kpr - Qh*y[1] - Qr*y[1])/Vsys
Ch <- (Qh*y[1] - (Vmax*fb*y[2]/Kph)/(Km+fb*y[2]/Kph) - Qh*y[2]/Kph)/Vh
Cr <- (Qr*y[1] - fb*GFR*y[1] - Qr*y[3]/Kpr)/Vr
Xu <- fb*GFR*y[1]
return(list(c(Cb,Ch,Cr,Xu)))
}

y <- lsoda(c(Cb=0.1,Ch=0.0,Cr=0.0,Xu=0.0),times,dydt,params)

par(mfrow=c(2,2))
plot(y[,1],y[,2],ylim=c(0,0.1),xlab="time(min)",ylab="Cb",type="l",col="red")
plot(y[,1],y[,3],ylim=c(0,0.01),xlab="time(min)",ylab="Ch",type="l",col="black")
plot(y[,1],y[,4],ylim=c(0,0.1),xlab="time(min)",ylab="Cr",type="l",col="green")
plot(y[,1],y[,5],ylim=c(0,0.8),xlab="time(min)",ylab="Xu",type="l",col="blue")

pk41_iv

p.o.

params = c(Vh=44,Vr=8,Qh=58,Qr=45.6,Vsys=1000,Kph=1,Kpr=1,RB=1.5,fp=0.1,Fa=1,Vmax=1000,Km=0.1,GFR=6,Dose=0.1,ka=0.1)
times=c(0,(1:240))

dydt <- function(t,y,p){
Vh <- p['Vh']
Vr <- p['Vr']
Qh <- p['Qh']
Qr <- p['Qr']
Vsys <- p['Vsys']
Kph <- p['Kph']
Kpr <- p['Kpr']
RB <- p['RB']
fp <- p['fp']
Fa <- p['Fa']
Vmax <- p['Vmax']
Km <- p['Km']
GFR <- p['GFR']
Dose <- p['Dose']
ka <- p['ka']
fb <- fp/RB
Cb <- (Qh*y[2]/Kph + Qr*y[3]/Kpr - Qh*y[1] - Qr*y[1])/Vsys
Ch <- (Qh*y[1] - (Vmax*fb*y[2]/Kph)/(Km+fb*y[2]/Kph) - Qh*y[2]/Kph + Dose*ka*Fa*exp(-ka*t))/Vh
Cr <- (Qr*y[1] - fb*GFR*y[1] - Qr*y[3]/Kpr)/Vr
Xu <- fb*GFR*y[1]
return(list(c(Cb,Ch,Cr,Xu)))
}

y <- lsoda(c(Cb=0.0,Ch=0.0,Cr=0.0,Xu=0.0),times,dydt,params)

par(mfrow=c(2,2))
plot(y[,1],y[,2],xlab="time(min)",ylab="Cb",type="l",col="red")
plot(y[,1],y[,3],xlab="time(min)",ylab="Ch",type="l",col="black")
plot(y[,1],y[,4],xlab="time(min)",ylab="Cr",type="l",col="green")
plot(y[,1],y[,5],xlab="time(min)",ylab="Xu",type="l",col="blue")

pk41_po

母集団薬物データの解析

薬物動態解析と統計の基礎から始まる概論的な読み物。

ProductName 母集団薬物データの解析 (統計科学選書)
矢船 明史,石黒 真木夫,赤池 弘次
朝倉書店 / ¥ 3,045 ()


  1. 序論
  2. 薬物動態解析の概略
  3. 統計理論の概略
  4. 情報量規準
  5. 母集団薬物動態解析
  6. 血液中薬物濃度推移のベイズ推定
  7. おわりに

1-4の概略はまとまっていてわかりやすいんだけど、5,6の実例がいきなり難しくなるという。MCMCとかPPKがどういった場面で使われるかとかそういうことが知りたい場合には役立つが、実際どうやって使うかとかは他の本のほうがよいかな。

  • Kullback-Leibler情報量の解釈
  • AICはその値に意味があるのではなく比較するモデルのAICの差に意味がある
  • ベイズの定理を利用して推定されるものは薬物動態パラメータの事後分布

PK3-5

クリアランス

ProductName ファーマコキネティクス―演習による理解
杉山 雄一,山下 伸二,加藤 基浩
南山堂 / ¥ 6,300 ()
在庫あり。

設問1

CLr = Dose/AUC * fe = Dose/AUC * [排泄量]/Dose = [排泄量]/AUC

> data <- read.csv("/Users/kzfm/PK/pk35a.csv")
> attach(data)
> data
  Dose  AUC     e
a  200  7.3  71.9
b  400 15.3 128.0
c  800 34.3 217.0

> CLr <- e/AUC
> CLr
[1] 9.849315 8.366013 6.326531

> F <- 0.85
> QH <- 1.5
> QR <- 1.2
> GFR <- 0.12
> fb <- 0.2
> Vd <- 100

設問2

> CLtot <- Dose*0.85/AUC
> CLtot
[1] 23.28767 22.22222 19.82507
> CLh <- CLtot - CLr
> CLh
[1] 13.43836 13.85621 13.49854

解答では E = 1-F = CLh/Qhから求めてた。 BAと消化管の吸収性から肝クリアランスの見積もりも可能。

設問3

> data2 <- read.csv("/Users/kzfm/PK/pk35b.csv")
> data2
   ve   conc
a   5  0.194
b  10  0.404
c  20  0.863
d  50  2.480
e 100  5.540
f 200 12.000

> CLtot <- data2[,"ve"]/data2[,"conc"]
> CLtot
[1] 25.77320 24.75248 23.17497 20.16129 18.05054 16.66667

> CLr <- CLtot -13.5
> CLr
[1] 12.273196 11.252475  9.674971  6.661290  4.550542  3.166667

設問4

CLrはGFRと腎尿細管分泌クリアランスの和 CLr = (1-FR)(fb*GFR + CLrs) CLrs = CLr - fb*GFR

> CLrs <- CLr -fb*GFR*60
> CLrs
[1] 10.833196  9.812475  8.234971  5.221290  3.110542  1.726667

設問5

CLrs_uint = CLrs*QR/fb*(QR-CLrs)

> CLrs_uint <- CLrs*QR*60/(fb*(QR*60-CLrs))
> CLrs_uint
[1] 63.75927 56.80385 46.49241 28.14766 16.25495  8.84546

> e <- nls(CLrs_uint ~ Vmax/(Km+data2[,"conc"]*fb),start=list(Vmax=10,Km=10))
> e
Nonlinear regression model
  model:  CLrs_uint ~ Vmax/(Km + data2[, "conc"] * fb) 
   data:  parent.frame() 
   Vmax      Km 
23.2331  0.3267 
 residual sum-of-squares: 0.2010

Number of iterations to convergence: 12 
Achieved convergence tolerance: 2.395e-07 

設問6

省略というかわからん

設問7

省略というかわからん

PK3-4

クリアランス

ProductName ファーマコキネティクス―演習による理解
杉山 雄一,山下 伸二,加藤 基浩
南山堂 / ¥ 6,300 ()
在庫あり。

設問1

定常状態ではe*(-kelt) ~ Css_min/Css_maxであることを利用して

Css_avg = Css_max-Css_min/(ln(Css_max) - ln(Css_min))

> data <- read.csv("/Users/kzfm/PK/pk34.csv")
> data
     D Css_max Css_min    fe
a   60    14.1    9.55 0.082
b  120    22.2   14.80 0.098
c  300    45.0   30.50 0.076
d  600    90.8   62.30 0.090
e 1200   114.0   81.30 0.110
f 2400   146.0  105.00 0.091

> attach(data)
> Css_avg <- (Css_max - Css_min)/log(Css_max/Css_min)
> Css_avg
[1]  11.67764  18.25065  37.28122  75.65745  96.73056 124.37574

> plot(D,Css_avg,ylim=c(0,150))
> par(new=T)
> plot(D,Css_max,ylim=c(0,150),xlab="",ylab="",pch=2)
> par(new=T)
> plot(D,Css_min,ylim=c(0,150),xlab="",ylab="",pch=3)
> par(new=T)

plot

DoseとCpssに線形性がない。

設問2

定常状態の尿中排泄速度を定常状態の平均血中濃度で割る

> CLr <- (D/12*fe/60*1000)/Css_avg
> CLr
[1] 0.5851641 0.8949455 0.8494000 0.9913102 1.8952991 2.4388464

設問3

設問2からCLrが増加しており、これはfbの上昇、つまり血中蛋白結合の飽和が示唆される。

PK3-3

クリアランスとメカニズム

ProductName ファーマコキネティクス―演習による理解
杉山 雄一,山下 伸二,加藤 基浩
南山堂 / ¥ 6,300 ()
在庫あり。

設問1

> data <- read.csv("/Users/kzfm/PK/pk33.csv")
> data
  drag  dose      AUC   fe
a    A 1e+00    169.0 0.13
b    A 1e+01   1720.0 0.14
c    A 1e+02  22400.0 0.18
d    A 1e+03 530000.0 0.42
e    B 1e-02      2.4 0.17
f    B 1e-01     23.4 0.17
g    B 1e+00    167.0 0.19
h    B 1e+01    658.0 0.25

> attach(data)
> CLtot <- dose/AUC * 1000
> CLtot
[1]  5.917160  5.813953  4.464286  1.886792  4.166667  4.273504  5.988024 15.197568

> CLr <- CLtot * fe
> CLr
[1] 0.7692308 0.8139535 0.8035714 0.7924528 0.7083333 0.7264957 1.1377246 3.7993921

> CLh <- CLtot - CLr
> CLh
[1]  5.147929  5.000000  3.660714  1.094340  3.458333  3.547009  4.850299 11.398176

設問2

CLr = fb * GFR つまりfb = VLr/GFR

> GFR <- 8
> QH <- 60

> fb_a <- CLr[1]/8
> fb_b <- CLr[5]/8
> fb_a
[1] 0.09615385
> fb_b
[1] 0.08854167

CLh = (QH*fb*CLh_uint)/(QH+fb*CLh_uint)つまりCLh_uint=CLh*QH/(fb*(QH-CLh))

> CLh_uint_a <- CLh[1]*QH/(fb_a*(QH-CLh[1]))
> CLh_uint_a
[1] 58.56311
> CLh_uint_b <- CLh[5]*QH/(fb_b*(QH-CLh[5]))
> CLh_uint_b
[1] 41.44783

設問3

DrugA

CLrの投与量変化が見られないことからfbの変化はない。つまり投与量増加に伴うCLhの減少は CLh_uintの低下、つまり代謝過程の飽和

DrugB

投与量増加に伴い、CLh,CLrが増加。CLrの増加はfbの増加によるもので、これはつまり血中蛋白結合の 飽和をあらわしている。 肝固有クリアランスとの関係はなんとなく。

PK3-2(設問4,5)

繰り返し経口投与の続き

ProductName ファーマコキネティクス―演習による理解
杉山 雄一,山下 伸二,加藤 基浩
南山堂 / ¥ 6,300 ()
在庫あり。

設問 4

肝クリアランスと肝血流の大きさから律速を判断。

v = Vmax*[S] / (Km + [S])

肝クリアランスに飽和が見られるから

CLh_int = Vmax/Km * 1/Cu

ここでCu = Cpss * fp

vh = Cb*CLh vr = Cpss*Rb*CLr

設問 5

v=D/t*(F-fe),C=Cpss*fp でMicaehlis-Menten式の形でフィッティング

> data <- read.csv("/Users/kzfm/PK/pk32b.csv")
> data
  Dose      C     v
a   10  0.212 0.308
b   30  0.675 0.850
c  100  2.400 2.420
d  300  7.850 5.750
e  500 13.800 7.920
> attach(data)
> e <- nls(v ~ Vmax*C/(Km+C),start=list(Vmax=10,Km=10))
> e
Nonlinear regression model
  model:  v ~ Vmax * C/(Km + C) 
   data:  parent.frame() 
 Vmax    Km 
15.09 12.59 
 residual sum-of-squares: 0.01301

PK3-2(設問1-3)

繰り返し経口投与

ProductName ファーマコキネティクス―演習による理解
杉山 雄一,山下 伸二,加藤 基浩
南山堂 / ¥ 6,300 ()
在庫あり。

設問1

省略

設問 2

CL = (F * (Dose/t)) / Cp,ss 10mgのデータからFを決定する

CLb=CLp/RbからRbをつかって全血基準のクリアランスを求める

t1/2をもとめる

一般に、薬物を持続投与する場合、消失半減期の5倍の時間が経っている時には定常状態に達しているとみてよいので これをもって定常状態の判断に利用してよい

設問 3

> data <- read.csv("/Users/kzfm/PK/pk3-2.csv")
> data
  dose     cpss   fe
A   10 0.000424 0.63
B   30 0.001350 0.66
C  100 0.004800 0.71
D  300 0.015700 0.77
E  500 0.027500 0.81

> attach(data)
> fp <- 0.5
> Rb <- 0.7
> Vd <- 30
> ka <- 1.5

> CLb <- (dose/12)/(cpss*1000*Rb)
> CLb
[1] 2.807727 2.645503 2.480159 2.274795 2.164502

> CLr <- CLb*fe
> CLr
[1] 1.768868 1.746032 1.760913 1.751592 1.753247

> CLh <- CLb - CLr
> CLh
[1] 1.0388589 0.8994709 0.7192460 0.5232029 0.4112554

PK3-1

設問1はクリアランスのパラメータを求める

設問2,3はi.v., p.o.それぞれの状況でパラメータを変えた場合どう変化するかを推定する

ProductName ファーマコキネティクス―演習による理解
杉山 雄一,山下 伸二,加藤 基浩
南山堂 / ¥ 6,300 ()
在庫あり。

設問1

> data <- read.csv("/Users/kzfm/PK/pk31.csv")
> data
  CLtot   Vd   fe  fb  Ka
A   0.2  0.1 1.00 0.1 1.5
B   0.2 10.0 1.00 0.1 1.5
C   1.3  0.1 0.02 0.1 1.5
D   1.3 10.0 0.02 0.1 1.5

> Qh <- 1.4
> Qr <- 1.1
> GFR <- 0.17
> Vb <- 0.077

> CLh <- (1-fe)*CLtot
> CLh
[1] 0.000 0.000 1.274 1.274

> CLr <- fe*CLtot
> CLr
[1] 0.200 0.200 0.026 0.026

> CLh_uint <- Qh*CLh/(fb*(Qh-CLh))
> CLh_uint
[1]   0.0000   0.0000 141.5556 141.5556

> CLr_uint <- Qr*(CLr -fb*GFR)/((Qr-CLr+fb*GFR)*fb)
> CLr_uint
[1] 2.19520174 2.19520174 0.09074244 0.09074244

設問2

AUCiv = D/CLtot

  1. fbの変動によって影響を受けるパラメータはVd,CLtot
  2. Qhの変動により影響を受けるパラメータはCLh
    • CLh = Qh*fb*CLh_uint/(Qh + fb*CLh_uint)なので分母の項の大小で律速が決まる
  3. CLr_uintの変動により影響をうけるのはCLr
  4. CLh_uintの変動により影響をうけるのはCLh

設問3

flip-flopを考える

ファーマコキネティクス 3章メモ

クリアランスの章

  • 固有クリアランスは血流や蛋白結合性などの影響を受けない本質的なパラメータ
  • 組織クリアランスと固有クリアランスの関係は、流入、流出する薬物濃度と細胞内薬物濃度の関係
  • 全身クリアランスは静注速度を定常状態の血中濃度で除する
  • 糸球体での薬物濾過速度は薬物によらずほぼ一定
  • 血液血漿濃度比

ProductName ファーマコキネティクス―演習による理解
杉山 雄一,山下 伸二,加藤 基浩
南山堂 / ¥ 6,300 ()
在庫あり。

pythonでdeconvolution

演習本にデコンボリューション法が載っていたのだけど、実装までは載ってなかったので、調べたら、論文が見つかった

ProductName ファーマコキネティクス―演習による理解
杉山 雄一,山下 伸二,加藤 基浩
南山堂 / ¥ 6,300 ()
在庫あり。

GOTOとIFの組み合わせがwhileなのかforなのか悩むとこが幾つかあったが、pythonで書きなおした。結構、fortranっぽさが残っている気もするが。

#!/usr/bin/env python

def calc_eps():
    eps = 1.0
    while eps + 1.0 > 1.0: eps = eps/2.0
    return eps

def sign(a1,a2):
    if a2 >= 0: return abs(a1)
    else:       return -abs(a1)

def deconv(a,div,b,dose):
    n = len(a); m = len(b); l = n + m -1
    eps = calc_eps()
    u = [0.0 for i in range(l)]
    v = [0.0 for i in range(l)]
    gamma = [0.0 for i in range(20)]
    g = [0.0 for i in range(20)]

    a.sort(lambda x,y:cmp(x[1],y[1]))

    for II in range(n-1):
        x1 = a[II][1] * -1
        x2 = a[II+1][1] * -1
        qx1 = 0.0
        qx2 = 0.0

        for i in range(n):
            s1 = 1.0
            s2 = 1.0
            for j in range(n):
                if j != i:
                    s1 = s1*(a[j][1] + x1)
                    s2 = s2*(a[j][1] + x2)
            qx1 = qx1 + a[i][0] * s1
            qx2 = qx2 + a[i][0] * s2

        x3 = x1
        qx3 = qx1
        d = x2 -x1
        e = d

        while(1):

            if abs(qx3) < abs(qx2):
                x1  = x2; x2  = x3; x3  = x1
                qx1 = qx2; qx2 = qx3; qx3 = qx1

            tol = 4.0*eps*abs(x2)
            xm  = 0.5*(x3-x2)

            if abs(xm) <= tol or qx2 == 0.0: break

            if abs(e) >= tol and abs(qx1) > abs(qx2):
                if x1 == x3:
                    s = qx2/qx1; p = 2.0*xm*s; q = 1.0-s
                else:
                    q = qx1/qx3; r = qx2/qx3; s = qx2/qx1
                    p = s*(2.0*xm*q*(q-r) - (x2-x1)*(r-1.0))
                    q = (q-1.0)*(r-1.0)*(s-1.0)

                if p > 0.0: q = -q
                p = abs(p)

                if 2.0*p < 3.0*xm*q - abs(tol*q) and p < abs(0.5*e*q):
                    e = d
                    d = p/q
                else:
                    d = xm
                    e = d
            else:
                d = xm
                e = d

            x1 = x2
            qx1 =qx2
            x2 = x2 + d
            if abs(d) <= tol: x2 = x2 - d + sign(tol,xm)
            qx2 = 0.0
            for i in range(n):
                s2 = 1.0
                for j in range(n):
                    if j != i: s2 = s2*(a[j][1] +x2)
                qx2 = qx2 + a[i][0] * s2
            if qx2*(qx3/abs(qx3)) > 0.0:
                x3 = x1
                qx3 = qx1
                d = x2 - x1
                e = d

        gamma[II] = x2
        s2 = 0.0
        for j in range(n):
            s1 = 0.0
            for k in range(n):
                if k != j:
                    s1 = s1 + 1.0/(gamma[II] + a[k][1])
            s2 = s2 + a[j][0]*s1/(gamma[II] + a[j][1])
        g[II] = 1.0/s2

    ak1 = 0.0
    ak2 = 0.0
    ak3 = 0.0
    for i in range(n):
        ak1 = ak1 + a[i][0];
        ak2 = ak2 + a[i][0]*a[i][1]
        if n > 1 and i != n-1: ak3 = ak3 + g[i]/gamma[i]
    ak1 = 1.0/ak1
    ak2 = ak1*ak1*ak2
    ak3 = ak2 - ak3
    ak4 = 100.0*div/dose
    uz = 0.0

    for i in range(l):
        s = 0.0
        if i > m-1:
            v[i] = -gamma[i-m]
            for j in range(m): 
                s = s + b[j][0]/(gamma[i-m] + b[j][1])
            u[i] = ak4*g[i-m]*s/gamma[i-m]
        else:
            v[i] = b[i][1]
            u[i] = ak4*b[i][0]*(ak1 - ak3/b[i][1])
            if n != 1:
                for j in range(n-1):
                    s = s + g[j]/(gamma[j]*(gamma[j] + b[i][1]))
                u[i] = u[i] - ak4*b[i][0]*s
        uz = uz - u[i]

    return uz,u,v

if __name__ =='__main__':
    print "=== Data set 1 ==="
    a     = [[1.0427,5.0526],[0.97617,0.97567]]
    div   = 1.0
    b     = [[-0.80547,3.6373],[0.80547,0.83033]]
    dose  = 0.6

    print deconv(a,div,b,dose)

    print "=== Data set 2 ==="
    a     = [[1.3377,3.1851],[0.52381,0.61065]]
    div   = 1.0
    b     = [[-2.2637,2.0800],[2.2637,1.2513]]
    dose  = 0.6

    print deconv(a,div,b,dose)

    print "=== Data set 3 ==="
    a     = [[1.0308,5.9131],[1.0487,1.0262]]
    div   = 1.0
    b     = [[-1.6103,3.3563],[1.6103,1.3917]]
    dose  = 0.6

    print deconv(a,div,b,dose)

    print "=== Data set 4 ==="
    a     = [[1.3377,3.1851],[0.52381,0.61065]]
    div   = 1.0
    b     = [[-5.5440,2.2488],[5.5440,1.7521]]
    dose  = 0.6

    print deconv(a,div,b,dose)

Table.1 のテスト計算を実行。論文のデータセット4のa2はタイポでしょう。

=== Data set 1 ===
(103.38116460095375, 
[99.756988719460651, -23.218031509886707, -179.92012181052769], 
[3.6373000000000002, 0.83033000000000001, 2.9469592648362699])
=== Data set 2 ===
(94.012343334076377, 
[212.39909284817011, 2395.320119577972, -2701.7315557602187], 
[2.0800000000000001, 1.2513000000000001, 1.3350740721242431])
=== Data set 3 ===
(94.362883301636657, 
[-1704.7606134762168, 73.013373899587776, 1537.3843562749923], 
[3.3563000000000001, 1.3916999999999999, 3.4906828227939406])
=== Data set 4 ===
(91.15932635164495, 
[370.51731307178676, -1111.1884882257646, 649.51184880233291], 
[2.2488000000000001, 1.7521, 1.3350740721242431])