library(multicore) # Für multicore Rechnungen library(deSolve) library(ggplot2) ggplot <- function(...) ggplot2::ggplot(...) + theme_bw() ## Definition der Funktionen --------------------------- ## Hodgkin-Huxley HH <- function(time, state, pars, alpha, beta) { with(as.list(c(state, pars)), { gK <- gKbar * n^4 dV = - gK/Cm * (V-VK) dn <- alpha(V)*(1-n) - beta(V)*n return(list(c(dV, dn), alpha=alpha(V), beta=beta(V), gK = gK)) }) } ## Ergebnis von ode() in data.frame umschreiben out2config <- function(out, parameter="none") { ncols <- dim(out)[2] nrows <- dim(out)[1] data.frame(time = out[,1], value = as.vector(out[,2:ncols]), species = rep(colnames(out)[-1], each=nrows), parameter = parameter) } ## Parameter festlegen und integrieren ---------------------------- pars <- c( Cm = 10^(-1.4), VK = 20, gKbar = 10^2) times <- seq(0, 40, len=200) alpha <- function(V) 0.01*(10-V)/(exp((10-V)/10)-1) beta <- function(V) 0.125*exp(-V/80) n0 <- c(0.095, 0.1, 0.105) # Für jeden Wert von n0 eine Integration durchführen # gibt eine Liste zurück out <- lapply(n0, function(myn) { yini <- c( V = -50, n = myn ) result <- ode(yini, times, HH, pars, alpha=alpha, beta=beta) return(result) }) # out2config auf jeden Listeneintrag anwenden outconfig <- lapply(1:3, function(i) out2config(out[[i]], paste0("n0=", n0[i]))) # Die Liste mit rbind zu einem data.frame() kombinieren outconfig <- do.call(rbind, outconfig) # Aus jedem Listeneintrag eine data.frame für die Phasenraumdarstellung machen outphase <- lapply(1:3, function(i) data.frame(out[[i]][,c("V", "gK")], parameter=paste0("n0=", n0[i]))) # Die Liste mit rbind zu einem data.frame() kombinieren outphase <- do.call(rbind, outphase) ## Ergebnisse darstellen --------------------------------- P1 <- ggplot(outconfig, aes(x=time, y=value, group=parameter, color=parameter)) + facet_wrap(~species, scales="free") + geom_line() print(P1) out2 <- subset(outconfig, species%in%c("V", "gK")) P2 <- ggplot(out2, aes(x=time, y=value, group=species, color=species)) + facet_wrap(~parameter, ncol=1) + geom_line() print(P2) P3 <- ggplot(outphase, aes(x=V, y=gK, group=parameter, color=parameter)) + geom_line() + scale_y_log10() print(P3) out3 <- subset(outconfig, species%in%c("alpha", "beta")) P4 <- ggplot(out3, aes(x=time, y=value, group=species, color=species)) + facet_wrap(~parameter, ncol=1) + geom_line() print(P4)