PK-PDモデリング
ファーマコキネティクス―演習による理解設問1
> data <- read.csv("/Users/kzfm/PK/pk7-1.csv")
> data
time Cp E
k 0.0 0.000 0.0
a 0.2 0.309 NA
b 0.5 0.560 17.4
c 1.5 0.678 35.4
d 2.0 0.556 50.1
e 3.0 0.386 NA
f 4.0 0.259 50.4
g 6.0 0.115 40.3
h 8.0 0.051 27.8
i 9.0 NA 22.0
j 12.0 NA 9.3
> plot(data$Cp,data$E)
半時計周りのヒステリシス
設問2
> nls(Cp ~ alpha*(exp(-beta*time)-exp(-gamma*time)),data[c("time","Cp")],start=list(alpha=100,beta=10,gamma=0.1))
Nonlinear regression model
model: Cp ~ alpha * (exp(-beta * time) - exp(-gamma * time))
data: data[c("time", "Cp")]
alpha beta gamma
-1.3322 1.8587 0.4085
residual sum-of-squares: 2.501e-07
Number of iterations to convergence: 8
Achieved convergence tolerance: 8.835e-07
ka = 1.8587
ke = 0.4085
つまりCp = -1.3322 (exp(-1.8587t)-exp(-0.4085*t))で
Ce = 1.3322*ke0/(ke0-0.4085)*exp(-0.4085*t) - \
1.3322*ke0/(ke0-1.8587)*exp(-1.8587*t) + 1.931666*ke0/((ke0-0.4085)*\
(ke0-1.8587))*exp(-ke0*t)
これをE=Emax*Ce/(EC50+Ce)=100*Ce/(EC50+Ce)に代入してフィッティング
> nls(E ~ 100*(1.3322*ke0/(ke0-0.4085)*exp(-0.4085*time) - \
1.3322*ke0/(ke0-1.8587)*exp(-1.8587*time) + 1.931666*ke0/((ke0-0.4085)* \
(ke0-1.8587))*exp(-ke0*time))/(EC50 + (1.3322*ke0/(ke0-0.4085)*exp(-0.4085*time) \
- 1.3322*ke0/(ke0-1.8587)*exp(-1.8587*time) + 1.931666*ke0/((ke0-0.4085)*\
(ke0-1.8587))*exp(-ke0*time))),data[c("time","E")],start=list(ke0=0.2,EC50=0.1))
Nonlinear regression model
model: E ~ 100 * (1.3322 * ke0/(ke0 - 0.4085) * exp(-0.4085 * time) \
- 1.3322 * ke0/(ke0 - 1.8587) * exp(-1.8587 * time) + 1.931666 * \
ke0/((ke0 - 0.4085) * (ke0 - 1.8587)) * exp(-ke0 * time))/ \
(EC50 + (1.3322 * ke0/(ke0 - 0.4085) * exp(-0.4085 * time) - 1.3322 *
ke0/(ke0 - 1.8587) * exp(-1.8587 * time) + 1.931666 * ke0/((ke0 - 0.4085) * (ke0 - 1.8587)) * exp(-ke0 * time)))
data: data[c("time", "E")]
ke0 EC50
0.4807 0.3595
residual sum-of-squares: 0.01642
Number of iterations to convergence: 6
Achieved convergence tolerance: 2.781e-07
結果
Ce = 8.869647*exp(-0.4085*t) + 0.4647232*exp(-1.8587*t) -9.332967*exp(-0.4807*t)
E = 100*(8.869647*exp(-0.4085*t) + 0.4647232*exp(-1.8587*t) - \
9.332967*exp(-0.4807*t))/(0.3595 + 8.869647*exp(-0.4085*t) + \
0.4647232*exp(-1.8587*t) - 9.332967*exp(-0.4807*t))
> plot(f,xlim=c(0,12),ylim=c(0,60),xlab="Time(hour)")
> plot(f,xlim=c(0,12),ylim=c(0,60),xlab="Time(hour)",ylab="E")
> par(new=T)
> plot(data$time,data$E,xlim=c(0,12),ylim=c(0,60),xlab="",ylab="",axes=F,type="p",col="blue")
設問3
CpとEの計算値も一致していることを確かめる
設問2のCeを求める際のラプラス変換があやふやだったのと、別解として用意されてあった連立微分方程式をRでフィッティングさせる方法が分からなかった。