library(deSolve) library(ggplot2) ## Definition der Funktionen --------------------------- ## Massenwirkungsgesetz ODE enzMassAct <- function(time, state, pars) { with(as.list(c(state, pars)), { ds <- k_m1*c - k_1*s*e de <- (k_m1 + k_2)*c - k_1*s*e dc <- k_1*s*e - (k_m1 + k_2)*c dp <- k_2*c return(list(c(ds, de, dc, dp))) }) } ## Gleichgewichtsapproximation ODE enzGleichGewApprox <- function(time, state, pars) { with(as.list(c(state, pars)), { ds <- -k_2*e_t*s / (k_m1/k_1 + s) dp <- -ds return(list(c(ds, dp))) }) } ## Ergebnis von ode() in data.frame umschreiben out2config <- function(out, ini) { data.frame(time = out[,1], value = as.vector(out[,2:5]), species = rep(c("s", "e", "c", "p"), each=dim(out)[1]), ini = ini) } ## Ergebnis von ode() in data.frame umschreiben out2phase <- function(out, ini) { data.frame(out[,2:5], ini=as.factor(ini)) } ## Ergebnis von ode() in data.frame umschreiben out2config2 <- function(out, ini) { data.frame(time = out[,1], value = as.vector(out[,2:4]), species = rep(c("s", "c", "p"), each=dim(out)[1]), ini = ini) } ## Ergebnis von ode() in data.frame umschreiben out2phase2 <- function(out, ini) { data.frame(out[,2:4], ini=as.factor(ini)) } ## Eigentliches Skript ------------------------------- times <- seq(0, 100, len=1000) sini <- seq(10,100,by=10) outConfig <- data.frame() outPhase <- data.frame() outConfig2 <- data.frame() outPhase2 <- data.frame() outConfig3 <- data.frame() outPhase3 <- data.frame() for(i in 1:length(sini)){ ini <- c(s = sini[i], e = 10, c = 0, p = 0) pars <- c(k_1 = 0.1, k_m1 = 0.1, k_2 = 0.1) out <- ode(ini, times, enzMassAct, pars) outConfig <- rbind(outConfig, out2config(out, ini[1])) outPhase <- rbind(outPhase, out2phase(out, ini[1])) ini2 <- c(s = sini[i], p = 0) e_t <- ini["e"]+ini["c"] pars2 <- c(k_1 = 0.1, k_m1 = 0.1, k_2 = 0.1, e_t = e_t) out2 <- ode(ini2, times, enzGleichGewApprox, pars2) out2 <- cbind(out2[,1:2],c = e_t*out2[,2]/(pars2["k_m1"]/pars2["k_1"] + out2[,2]), p = out2[,3]) outConfig2 <- rbind(outConfig2, out2config2(out2, ini2[1])) outPhase2 <- rbind(outPhase2, out2phase2(out2, ini2[1])) out3 <- cbind(out2[,1:2],c = e_t*out2[,2]/((pars2["k_m1"]+pars2["k_2"])/pars2["k_1"] + out2[,2]), p = out2[,3]) outConfig3 <- rbind(outConfig3, out2config2(out3, ini2[1])) outPhase3 <- rbind(outPhase3, out2phase2(out3, ini2[1])) } # Konfigurationsraumplot P1 <- ggplot(outConfig, aes(x=time, y=value, group=paste0(species,ini))) + geom_line() + facet_wrap(~species, scales="free") print(P1 + theme_bw()) # Phasenraumplot P2 <- ggplot(outPhase, aes(x=s, y=c, group=ini)) + geom_path() print(P2 + theme_bw()) P3 <- ggplot(outConfig2, aes(x=time, y=value, group=paste0(species,ini))) + geom_line() + facet_wrap(~species, scales="free") print(P3 + theme_bw()) P4 <- ggplot(outPhase2, aes(x=s, y=c, group=ini)) + geom_path(size=3, alpha=0.1) print(P4 + theme_bw()) P5 <- ggplot(outConfig3, aes(x=time, y=value, group=paste0(species,ini))) + geom_line() + facet_wrap(~species, scales="free") print(P5 + theme_bw()) P6 <- ggplot(outPhase3, aes(x=s, y=c, group=ini)) + geom_path() print(P6 + theme_bw()) print(P4 + geom_path(data=outPhase3, color="red") + geom_path(data=outPhase, color="blue"))