HHsim <- function(
	samplitude=c(6.0, -3.0, 2.0), #amplitude of sample command stimulations

	###simulation parameters
	tmax=100.0,     #end time
	dt=0.025,       #delta time
	sdelay=25.0,    #delay to command stimulation onset
	sduration=50.0, #duration of command stimulation
	maxline=10,     #maximum number of the lines

	###physiological parameters (set identical to that of original H-H simulation)
	Cm=1.0,        #membrane capacitance
	Gk0=36.0,      #potassium conductance
	Gna0=120.0,    #sodium conductance
	Gl0=0.3,       #chloride conductance

	###equivalent potentials (set to satisfy the original Eeq relationships)
	Vrest=-65.0,   #resting membrane potential
	Ek=-77.0,      #potassium equivalent potential
	Ena=50.0,      #sodium equivalent potential
	El=-54.387     #chloride equivalent potential
){

cat("Choose a graphic mode.\n\n")
cat("\t(1)6-split window ... presents Vm, CmdCurrent, Ina, Ik, Gna, and Gk.\n")
cat("\t(2)wide window ...... presents only Vm and CmdCurrent, in a time sequence.\n\n")
gmode <- as.numeric( readline("put the NUMBER of the graphic mode : ") )
cat("\n")

cond <- length(samplitude) #number of sample amplitudes
if(gmode==1) {
	col1 <- colorRampPalette(c("skyblue","royalblue4"))(cond)
	col2 <- colorRampPalette(c("orange","red4"))(maxline-cond)
	col  <- c(col1, col2)
	ts <- list()
	length(ts) <- maxline
	ts[] <- list( seq(0, tmax-dt, dt) )
	nsplit <- split.screen(c(3,2))
} else if(gmode==2) {
	col <- colorRampPalette(c("skyblue","royalblue4"))(maxline)
	ts <- list()
	length(ts) <- maxline
	ts[] <- list( seq(0, tmax-dt, dt) )
	for(i in 1:maxline) { ts[[i]] <- ts[[i]] + (i-1)*tmax }
	nsplit <- split.screen(c(2,1))
} else {
	stop(paste("There is no graphic mode numbered ", gmode, ".", sep=""))
}

#data matrix generation
Vms    <- list()
Istims <- list()
Inas   <- list()
Iks    <- list()
Gnas   <- list()
Gks    <- list()


HHsimloop <- function(sampl) { #ti=list(tmax=tmax, dt=dt, sdelay=sdelay, sduration=sduration), phy=list(Vrest=Vrest, Ek=Ek, Ena=Ena, El=El, Cm=Cm, Gk0=Gk0, Gna0=Gna0, Gl0=Gl0)

	Vms <- Istims <- Inas <- Iks <- Gnas <- Gks <- rep(0, tmax/dt)

	###reset parameters
	Vm <- 0.0
	An <- 0.01*(10.0-Vm)/(exp((10.0-Vm)/10.0)-1.0)
	Bn <- 0.125*exp(-1*Vm/80.0)
	Am <- 0.1*(25.0-Vm)/(exp((25.0-Vm)/10.0)-1.0)
	Bm <- 4.0*exp(-1*Vm/18.0)
	Ah <- 0.07*exp(-1*Vm/20.0)
	Bh <- 1.0*exp((30.0-Vm)/10.0+1.0)
	n <- An/(An+Bn)
	m <- Am/(Am+Bm)
	h <- Ah/(Ah+Bh)

	for(j in 1:(tmax/dt) ) { #time sequence loop
		if( ((sdelay/dt)<=j) && (j<=((sdelay+sduration) / dt)) ) { Istim <- sampl }
		else { Istim <- 0.0 }

		An <- 0.01  * (10.0-Vm) / (exp((10.0-Vm)/10.0)-1.0)
		Bn <- 0.125 * exp(-1*Vm/80.0)
		Am <- 0.1   * (25.0-Vm) / (exp((25.0-Vm)/10.0)-1.0)
		Bm <- 4.0   * exp(-1*Vm/18.0)
		Ah <- 0.07  * exp(-1*Vm/20.0)
		Bh <- 1.0 / (exp((30.0-Vm)/10.0)+1.0)
		n  <- n + (An*(1.0-n)-Bn*n) * dt
		m  <- m + (Am*(1.0-m)-Bm*m) * dt
		h  <- h + (Ah*(1.0-h)-Bh*h) * dt

		Gk  <- Gk0  * n^4
		Gna <- Gna0 * m^3 * h
		Gl  <- Gl0

		Ik  <- Gk  * (Vm + Vrest - Ek)
		Ina <- Gna * (Vm + Vrest - Ena)
		Il  <- Gl  * (Vm + Vrest - El)
		Ichannel <- Ik+Ina+Il

		Im <- Istim - Ichannel
		Vm <- Vm + (Im / Cm)*dt

		Vms[j]    <- Vm + Vrest
		Istims[j] <- Istim
		Inas[j]   <- Ina
		Iks[j]    <- Ik
		Gnas[j]   <- Gna
		Gks[j]    <- Gk
	} #time sequence loop

return(list(Vms=Vms, Istims=Istims, Inas=Inas, Iks=Iks, Gnas=Gnas, Gks=Gks))
}


########## START SIMULATION ##########
###
for(i in 1:cond) {

	rslt <- HHsimloop(samplitude[i])

	Vms[[i]]    <- rslt$Vms
	Istims[[i]] <- rslt$Istims
	Inas[[i]]   <- rslt$Inas
	Iks[[i]]    <- rslt$Iks
	Gnas[[i]]   <- rslt$Gnas
	Gks[[i]]    <- rslt$Gks

}
###
########## END OF SIMULATION ##########


########## GRAPHICAL OUTPUT ##########
###

lines1 <- function(x, y, lcol, llty) { lines(x, y, lty=llty, col=lcol, lwd=2) }

if(gmode==1) {
	plot1 <- function(xf, Vmsf,    colf, lt) { par(new=TRUE,mar=c(3,3,2,1)) ; plot(0, 0, xlim=c(0,tmax), ylim=c(-100,50),    xlab="time", ylab="Vm",    main="voltage",               type="n", col=colf) ; mapply(FUN=lines1, xf, Vmsf,    colf, lt) }
	plot2 <- function(xf, Istimsf, colf, lt) { par(new=TRUE,mar=c(3,3,2,1)) ; plot(0, 0, xlim=c(0,tmax), ylim=c(-10,50),     xlab="time", ylab="Istim", main="command current",       type="n", col=colf) ; mapply(FUN=lines1, xf, Istimsf, colf, lt) }
	plot3 <- function(xf, Inasf,   colf, lt) { par(new=TRUE,mar=c(3,3,2,1)) ; plot(0, 0, xlim=c(0,tmax), ylim=c(-1000,1000), xlab="time", ylab="Ina",   main="sodium current",        type="n", col=colf) ; mapply(FUN=lines1, xf, Inasf,   colf, lt) }
	plot4 <- function(xf, Iksf,    colf, lt) { par(new=TRUE,mar=c(3,3,2,1)) ; plot(0, 0, xlim=c(0,tmax), ylim=c(-1000,1000), xlab="time", ylab="Ik",    main="potassium current",     type="n", col=colf) ; mapply(FUN=lines1, xf, Iksf,    colf, lt) }
	plot5 <- function(xf, Gnasf,   colf, lt) { par(new=TRUE,mar=c(3,3,2,1)) ; plot(0, 0, xlim=c(0,tmax), ylim=c(0,35),       xlab="time", ylab="Gna",   main="sodium conductance",    type="n", col=colf) ; mapply(FUN=lines1, xf, Gnasf,   colf, lt) }
	plot6 <- function(xf, Gksf,    colf, lt) { par(new=TRUE,mar=c(3,3,2,1)) ; plot(0, 0, xlim=c(0,tmax), ylim=c(0,35),       xlab="time", ylab="Gk",    main="potassium conductance", type="n", col=colf) ; mapply(FUN=lines1, xf, Gksf,    colf, lt) }

	screen(nsplit[1]) ; plot1(ts[1:cond], Vms,    col[1:cond], "11")
	screen(nsplit[2]) ; plot2(ts[1:cond], Istims, col[1:cond], "11")
	screen(nsplit[3]) ; plot3(ts[1:cond], Inas,   col[1:cond], "11")
	screen(nsplit[4]) ; plot4(ts[1:cond], Iks,    col[1:cond], "11")
	screen(nsplit[5]) ; plot5(ts[1:cond], Gnas,   col[1:cond], "11")
	screen(nsplit[6]) ; plot6(ts[1:cond], Gks,    col[1:cond], "11")
} else if(gmode==2) {
	plot1 <- function(xf, Vmsf,    colf, lt) { par(new=TRUE,mar=c(3,3,2,1)) ; plot(0, 0, xlim=c(0,tmax*maxline), ylim=c(-100,50),    xlab="time", ylab="Vm",    main="voltage",               type="n", col=colf) ; mapply(FUN=lines1, xf, Vmsf,    colf, lt) }
	plot2 <- function(xf, Istimsf, colf, lt) { par(new=TRUE,mar=c(3,3,2,1)) ; plot(0, 0, xlim=c(0,tmax*maxline), ylim=c(-10,50),     xlab="time", ylab="Istim", main="command current",       type="n", col=colf) ; mapply(FUN=lines1, xf, Istimsf, colf, lt) }

	screen(nsplit[1]) ; plot1(ts[1:cond], Vms,    col[1:cond], "11")
	screen(nsplit[2]) ; plot2(ts[1:cond], Istims, col[1:cond], "11")
}

###
########## END OF GPAPH. OUTPUT ##########


########## OPTIONAL GRAPH DRAWING ##########
###

for(i in (cond+1):maxline){

	choice=readline("put a new command current : ")
	if(choice=="") { break }

	rslt <- HHsimloop(as.numeric(choice))

	Vms[[i]]    <- rslt$Vms
	Istims[[i]] <- rslt$Istims
	Inas[[i]]   <- rslt$Inas
	Iks[[i]]    <- rslt$Iks
	Gnas[[i]]   <- rslt$Gnas
	Gks[[i]]    <- rslt$Gks

	if(gmode==1) {
		screen(nsplit[1]) ; plot1(ts[i], Vms[i],    col[i], "11")
		screen(nsplit[2]) ; plot2(ts[i], Istims[i], col[i], "11")
		screen(nsplit[3]) ; plot3(ts[i], Inas[i],   col[i], "11")
		screen(nsplit[4]) ; plot4(ts[i], Iks[i],    col[i], "11")
		screen(nsplit[5]) ; plot5(ts[i], Gnas[i],   col[i], "11")
		screen(nsplit[6]) ; plot6(ts[i], Gks[i],    col[i], "11")
	} else if(gmode==2) {
		screen(nsplit[1]) ; plot1(ts[i], Vms[i],    col[i], "11")
		screen(nsplit[2]) ; plot2(ts[i], Istims[i], col[i], "11")
	}

}

###
########## END OF OPTIONAL GRAPHS ##########

invisible(NULL)
}

comment(HHsim) <- c(
"HHsim - Draw a predicted electrophysiological response of a neuron by Hodgkin-Huxley model.",
"",
"SYNOPSIS",
"    HHsim(samplitude)",
"INPUT",
"    samplitude : a numeric vector which contains the default amplitudes of command current.",
"                 Set sapmlitude=0 if you don't need any default drawing of command.",
"    other variables such as the cell's resting potential or conductances, as well as the stimulus parameters",
"    can be changed, although it isn't necessary in general case.",
"OUTPUT",
"    NULL.",
"",
"See http://www7b.biglobe.ne.jp/~homunculus/r/HHsim.html for detailed instruction in Japanese.",
"Reference",
"    Miyakawa, H., & Inoue, M. 2003. Biophysics of the neuron. (in Japanese) Maruzen: Tokyo.",
"ver.1.10, written by MOCHI, 2009."
)