脳内のレセプターと薬物の結合のシミュレーション
odesolveで血中濃度が時間依存的に減少するとき(時間の関数になっているとき)どういう風に書けばいいのか悩んだが
Cp <- dose/mwt*(50*exp(-0.1*t) + 125*exp(-0.025*t))
とtそのまま使えることがわかったので解決
library(odesolve)
times=c(0,(1:120))
params = c(fb=0.2,mwt=250, Kp=8.0, kon=2,koff=0.2,Rt=1,Vt=1,Q=1,dose=0.2)
dydt <- function(t,y,p){
fb <- p['fb']
Kp <- p['Kp']
kon <- p['kon']
koff <- p['koff']
Rt <- p['Rt']
Vt <- p['Vt']
Q <- p['Q']
dose <- p['dose']
mwt <- p['mwt']
Cp <- dose/mwt*(50*exp(-0.1*t) + 125*exp(-0.025*t))
Ct <- (Q*(Cp-y[1]/Kp)-kon*fb*y[1]/Kp*Rt*(1-y[2])+koff*Rt*y[2])/Vt
PHAI <- kon*fb*y[1]/Kp*(1-y[2])-koff*y[2]
return(list(c(Ct,PHAI)))
}
y <- lsoda(c(Ct=0.0,PHAI=0.0),times,dydt,params)
y02 <- y
params = c(fb=0.2,mwt=250, Kp=8.0, kon=2,koff=0.2,Rt=1,Vt=1,Q=1,dose=1)
y <- lsoda(c(Ct=0.0,PHAI=0.0),times,dydt,params)
y1 <- y
params = c(fb=0.2,mwt=250, Kp=8.0, kon=2,koff=0.2,Rt=1,Vt=1,Q=1,dose=5)
y <- lsoda(c(Ct=0.0,PHAI=0.0),times,dydt,params)
y5 <- y
plot(y02[,1],y02[,3],ylim=c(0,1.0),xlab="time(min)",ylab="Phai",type="l",col="red")
par(new=T)
plot(y1[,1],y1[,3],ylim=c(0,1.0),xlab="",ylab="",axes=F,type="l",col="blue")
par(new=T)
plot(y5[,1],y5[,3],ylim=c(0,1.0),xlab="",ylab="",axes=F,type="l",col="green")
plot(y2[,2],y5[,3],ylim=c(0,1.0),xlab="Cb",ylab="Phai",type="l",col="red")
ヒステリシスがみられる