## model dM1 <- function(dt, M1, M2, e01, k12, k21, k10) dt*(e01 - k12*M1 + k21*M2 - k10*M1) dM2 <- function(dt, M1, M2, k12, k21) dt*( + k12*M1 - k21*M2) ## `sol2C` implements the Euler method of ## numerical integration, where h is the step size ## in time, and f(t, y) is the differential equation dy/dt: ## y_n = y_n-1 + h*f(t_n-1, y_n-1) sol2C <- function(e01, k12=2, k21=2, k10=0.7, dt=1/60) { M1 <- M2 <- time <- rep(0, length(e01)) for(t in 1:(length(e01)-1)) { d1 <- dM1(dt, M1[t], M2[t], e01[t], k12, k21, k10) d2 <- dM2(dt, M1[t], M2[t], k12, k21) M1[t+1] <- M1[t] + d1 M2[t+1] <- M2[t] + d2 time[t+1] <- time[t] + dt } ret <- data.frame(M1=M1, M2=M2, e01=e01, time=time) class(ret) <- c("sol2C", class(ret)) return(ret) } ## plot the results of `sol2C` plot.sol2C <- function(obj) { par(mfrow=c(2,1), mai=rep(0,4), omi=rep(0.75,4)) plot(obj\$time, obj\$M1, type="l", ylim=range(obj\$M1, obj\$M2), ylab="Mass", xlab="Time", xaxt="n", yaxt="n") mtext("Mass", 2, 1, xpd=NA) lines(obj\$time, obj\$M2, lty=2) abline(h=0, lty=2, col="gray") legend("topright", legend=c("M1", "M2"), lty=c(1,2), cex=0.7, bty="n") plot(obj\$time, obj\$e01, type="l", ylab="e01", xlab="Time", xaxt="n", yaxt="n") mtext("e01", 2, 1, xpd=NA) abline(h=0, lty=2, col="gray") } ## stimulus functions dt <- 1/60 time <- seq(0,8,dt) ## single IV infusion e01_1iv <- ifelse(time <= 30/60, 1, 0) ## double IV infusion e01_2iv <- ifelse(time <= 30/60 | (time > 60/60 & time <= 90/60), 1, 0) ## single IV followed by removal e01_1iv1rm <- ifelse(time <= 30/60, 1, 0) e01_1iv1rm <- ifelse(time > 60/60 & time <= 90/60, -0.45, e01_1iv1rm) ## single IV followed by removal then additional IV e01_2iv1rm <- ifelse(time <= 30/60, 1, 0) e01_2iv1rm <- ifelse(time > 60/60 & time <= 90/60, -0.45, e01_2iv1rm) e01_2iv1rm <- ifelse(time > 120/60 & time <= 150/60, 1.0, e01_2iv1rm) ## sawtooth e01_saw <- ifelse(time <= 30/60, time, 0) e01_saw <- ifelse(time > 60/60 & time <= 90/60, time-60/60, e01_saw) ## compute the Jacobian using a numerical derivative sen2C <- function(e01, k12=2, k21=2, k10=0.7, dt=1/60) { sol0 <- sol2C(e01, k12, k21, k10, dt) delk12 <- sqrt(.Machine\$double.eps)*k12 solk12 <- sol2C(e01, k12+delk12, k21, k10, dt) delk21 <- sqrt(.Machine\$double.eps)*k21 solk21 <- sol2C(e01, k12, k21+delk21, k10, dt) delk10 <- sqrt(.Machine\$double.eps)*k10 solk10 <- sol2C(e01, k12, k21, k10+delk10, dt) ret <- data.frame(M1=sol0\$M1, M2=sol0\$M2, k12=(solk12\$M1-sol0\$M1)/delk12, k21=(solk21\$M1-sol0\$M1)/delk21, k10=(solk10\$M1-sol0\$M1)/delk10, e01=e01, time=sol0\$time) class(ret) <- c("sen2C", class(ret)) return(ret) } ## compute the parameter correlations cor2C <- function(sen, times=c(seq(0.25, 3, 0.25),4,6,8)) { senm <- as.matrix(sen[,c("k12","k21","k10")]) covm <- t(senm)%*%senm return(cov2cor(covm)) } ## plot the model solutions and Jacobian plot.sen2C <- function(obj) { par(mfrow=c(5,1), mai=rep(0,4), omi=rep(0.5,4)) plot(obj\$time, obj\$M1, type="l", ylim=range(obj\$M1, obj\$M2), ylab="Mass", xaxt="n", yaxt="n") mtext("Mass", 2, 1, xpd=NA, cex=0.8) lines(obj\$time, obj\$M2, lty=2) abline(h=0, lty=2, col="gray") legend("topright", legend=c("M1", "M2"), lty=c(1,2), bty="n") plot(obj\$time, obj\$e01, type="l", ylab="e01", xlab="Time", xaxt="n", yaxt="n") mtext("e01", 2, 1, xpd=NA, cex=0.8) abline(h=0, lty=2, col="gray") plot(obj\$time, obj\$k12, type="l", ylab="dM1/dk12", xlab="Time", xaxt="n", yaxt="n") mtext("dM1/dk12", 2, 1, xpd=NA, cex=0.8) abline(h=0, lty=2, col="gray") plot(obj\$time, obj\$k21, type="l", ylab="dM1/dk21", xlab="Time", xaxt="n", yaxt="n") mtext("dM1/dk21", 2, 1, xpd=NA, cex=0.8) abline(h=0, lty=2, col="gray") plot(obj\$time, obj\$k10, type="l", ylab="dM1/dk10", xlab="Time", xaxt="n", yaxt="n") mtext("dM1/dk10", 2, 1, xpd=NA, cex=0.8) mtext("Time", 1, 1, xpd=NA, cex=0.8) abline(h=0, lty=2, col="gray") } ## plot the solutions plot(sol2C(e01_1iv)) plot(sol2C(e01_2iv)) plot(sol2C(e01_saw)) plot(sol2C(e01_1iv1rm)) plot(sol2C(e01_2iv1rm)) ## plot the solutions with Jacobian plot(sen2C(e01_1iv)) plot(sen2C(e01_2iv)) plot(sen2C(e01_saw)) plot(sen2C(e01_1iv1rm)) plot(sen2C(e01_2iv1rm)) ## show the correlations cor2C(sen2C(e01_1iv)) cor2C(sen2C(e01_2iv)) cor2C(sen2C(e01_saw)) cor2C(sen2C(e01_1iv1rm)) cor2C(sen2C(e01_2iv1rm))