フェンタニルの持続静注をすることで術後鎮痛を行うことがあるが、フェンタニルを低用量で開始するだけでは目標とする鎮痛濃度に達するにはものすごい時間がかかる。
例えば、身長160cm、体重50kg、年齢40歳の男性に20ug/hr でフェンタニルの持続静注をすると、効果部位の平衡濃度（の0.05 ng/ml 誤差範囲）に達するのに27時間程度かかる。
術後鎮痛に必要なフェンタニルの部位効果濃度は1-3 ng/ml 程度らしいので20ug/hr の投与はそもそも足りてないのだが、こんなものである。
というわけで、投与前もしくは投与してから早めの段階で、けっこうな用量をボーラス投与して、部位効果濃度を高くする。
体重と同用量の投与（この場合は50ug のフェンタニル）では、24時間
x3 （150ug）では、15時間
x5 （250ug）では、4時間程度になる。
library(shiny)
library(shinyTime)
library(googleVis)
weight <- 50
height <- 160
age <- 40
hour <- 60
Cont_conc <- 20
starttime <- strptime("09:00:00", "%T")
gender <- "M"
iv1 <- 3*weight
comp <- c(1:3, "e")
cname <- c("Bolus", "Cont", sprintf("A%s", comp), mapply(function(z) sprintf("dA%s_%d", comp, z), 1:4), sprintf("C%s", comp))
X <- matrix(0, hour*60, length(cname), dimnames=list(NULL, cname))
X <- as.data.frame(X)
timestep <- 60
k_1N <- c(0.0827, 0.471, 0.225)
k_N1 <- c(0.102, 0.006)
k_e0 <- 0.114
k <- c(k_1N, k_N1, k_e0)/(60/timestep)
names(k) <- c(sprintf("k1%d", c(0, 2, 3)), sprintf("k%d1", 2:3), "ke0")
LBM <- ifelse(gender == "M", 1.1*weight-128*(weight/height)^2,1.07*weight-148*(weight/height)^2)
k1 <- k["k10"] + k["k12"] + k["k13"]
V1 <- 0.105*weight
CLs <- k[c("k10", "k12", "k13")]*(60/timestep)*V1
V2 <- CLs[2]/(k["k21"]*(60/timestep))
V3 <- CLs[3]/(k["k31"]*(60/timestep))
t0 <- starttime - 61*0
tx <- t0 + (seq(nrow(X))-1)*60
X[, "Cont"] <- Cont_conc/60*1000
X <- na.omit(X)
X <- rbind(0, X)
X[2, "Bolus"] <- iv1*1000*(60/timestep)
for(i in 1:(nrow(X)-1)){
X$dA1_1[i+1] <- k["k21"]*X$A2[i]+k["k31"]*X$A3[i]-k1*X$A1[i]+X$Cont[i+1]/(60/timestep)
X$dA2_1[i+1] <- k["k12"]*X$A1[i] - k["k21"]*X$A2[i]
X$dA3_1[i+1] <- k["k13"]*X$A1[i] - k["k31"]*X$A3[i]
X$dAe_1[i+1] <- k["ke0"]*(X$A1[i] - X$Ae[i])
X$dA1_2[i+1] <- k["k21"]*(X$A2[i]+X$dA2_1[i+1]/2)+k["k31"]*(X$A3[i]+X$dA3_1[i+1]/2)-k1*(X$A1[i]+X$dA1_1[i+1]/2)+X$Cont[i+1]/(60/timestep)
X$dA2_2[i+1] <- k["k12"]*(X$A1[i]+X$dA1_1[i+1]/2) - k["k21"]*(X$A2[i]+X$dA2_1[i+1]/2)
X$dA3_2[i+1] <- k["k13"]*(X$A1[i]+X$dA1_1[i+1]/2) - k["k31"]*(X$A3[i]+X$dA3_1[i+1]/2)
X$dAe_2[i+1] <- k["ke0"]*(X$A1[i]+X$dA1_1[i+1]/2 - (X$Ae[i]+X$dAe_1[i+1]/2))
X$dA1_3[i+1] <- k["k21"]*(X$A2[i]+X$dA2_2[i+1]/2)+k["k31"]*(X$A3[i]+X$dA3_2[i+1]/2)-k1*(X$A1[i]+X$dA1_2[i+1]/2)+X$Cont[i+1]/(60/timestep)
X$dA2_3[i+1] <- k["k12"]*(X$A1[i]+X$dA1_2[i+1]/2) - k["k21"]*(X$A2[i]+X$dA2_2[i+1]/2)
X$dA3_3[i+1] <- k["k13"]*(X$A1[i]+X$dA1_2[i+1]/2) - k["k31"]*(X$A3[i]+X$dA3_2[i+1]/2)
X$dAe_3[i+1] <- k["ke0"]*(X$A1[i]+X$dA1_2[i+1]/2 - (X$Ae[i]+X$dAe_2[i+1]/2))
X$dA1_4[i+1] <- k["k21"]*(X$A2[i]+X$dA2_3[i+1])+k["k31"]*(X$A3[i]+X$dA3_3[i+1])-k1*(X$A1[i]+X$dA1_3[i+1])+X$Cont[i+1]/(60/timestep)
X$dA2_4[i+1] <- k["k12"]*(X$A1[i]+X$dA1_3[i+1]) - k["k21"]*(X$A2[i]+X$dA2_3[i+1])
X$dA3_4[i+1] <- k["k13"]*(X$A1[i]+X$dA1_3[i+1]) - k["k31"]*(X$A3[i]+X$dA3_3[i+1])
X$dAe_4[i+1] <- k["ke0"]*(X$A1[i]+X$dA1_3[i+1] - (X$Ae[i]+X$dAe_3[i+1]))
X$A1[i+1] <- X$A1[i]+(X$dA1_1[i+1]+2*X$dA1_2[i+1]+2*X$dA1_3[i+1]+X$dA1_4[i+1])/6+X$Bolus[i+1]/(60/timestep)
X$A2[i+1] <- X$A2[i]+(X$dA2_1[i+1]+2*X$dA2_2[i+1]+2*X$dA2_3[i+1]+X$dA2_4[i+1])/6
X$A3[i+1] <- X$A3[i]+(X$dA3_1[i+1]+2*X$dA3_2[i+1]+2*X$dA3_3[i+1]+X$dA3_4[i+1])/6
X$Ae[i+1] <- X$Ae[i]+(X$dAe_1[i+1]+2*X$dAe_2[i+1]+2*X$dAe_3[i+1]+X$dAe_4[i+1])/6
X$Ce[i+1] <- X$Ae[i+1]/V1/1000
X$C1[i+1] <- X$A1[i+1]/V1/1000
X$C2[i+1] <- X$A2[i+1]/V2/1000
X$C3[i+1] <- X$A3[i+1]/V3/1000
}
X <- X[-1,]
yl <- c(0, 2)
par(mar=c(5, 5, 2, 1), las=1, cex.lab=1.5)
plot(tx, X$Ce, type="n", ylim=yl, xlab="Time", ylab="Effect concentratoin [ng/ml]")
lines(tx, X$Ce, lwd=7)
abline(h=1.0, lty=3)
abline(h=2.0, lty=3, col="red", lwd=2)
abline(v=seq.POSIXt(tx[1], tail(tx, 1)+100, by=15*60), lty=3)
X <- data.frame(time=tx, X)
plot(X$time, X$Ce, type="l")
Cs <- grep("C[123e]", colnames(X), value=TRUE)
matplot(tx, X[, Cs], type="l")
cols <- c("yellow", "orange", "hotpink", "blue")
h24 <- seq.POSIXt(tx[1], tail(tx, 1)+100, by=60*60*24)
par(mar=c(5, 5, 2, 1), las=1, cex.lab=1.5)
plot(tx, X$Ce, type="n", ylim=yl, xlab="Time", ylab="Effect concentratoin [ng/ml]", xaxt="n")
for(i in seq(Cs)) lines(tx, X[, Cs[i]], lwd=7, col=cols[i])
axis(1, at=h24, labels=24*(seq_along(h24)-1))
abline(h=1.0, lty=3)
abline(h=2.0, lty=3, col="red", lwd=2)
abline(v=seq.POSIXt(tx[1], tail(tx, 1)+100, by=120*60), lty=3)
abline(v=h24[-1], lty=1, lwd=3)
legend("bottomright", legend=Cs, col=cols, pch=15, cex=2, bg="white")
xd <- tail(which(abs(X$Ce-tail(X$Ce, 1)) > 0.05), 1)
abline(v=tx[xd], col=2)
eqt <- difftime(tx[xd], tx[1], units="hours")
text(tx[xd], par()$usr[4], sprintf("%0.2f hrs", eqt), pos=3, xpd=TRUE)
