# # Neuron model # # Matlab code taken from [1] translated into R. # # [1] Izhikevich, E. M. (2007), Dynamical Systems in Neuroscience: # The Geometry of Excitability and Bursting, Chapter 8, p274, MIT # Press, Cambridge, MA # # ODE routine setwd("c:/R/bme_ode/chap4") source("rkf45.R") source("neuron_1.R") # # Model parameters C=100; vr=-60; vt=-40; k=0.7; # parameters used for RS an=0.03; bn=-2; cn=-50; dn=100; # neocortical pyramidal neurons vpeak=35 # spike cutoff # # Pulse of input DC current nout=1000; In=rep(0,nout+1); for(i in 1:nout+1){ if(i< 101)In[i]=0; if(i>=101)In[i]=70;} # # Initial condition v =rep(0,nout+1); w =rep(0,nout+1); v[1]=vr;w[1]=0; t0=seq(from=0,to=nout,by=1); ncall=0; # # Parameters for t integration nt=1;h=1; # # rkf45 integration nint=1; for(i in 1:nout){ yout=rkf45(nt,h,t0[i],c(v[i],w[i])) v[i+1]=yout[1]; w[i+1]=yout[2]; if(v[i+1]>=vpeak){ # a spike is fired v[i]=vpeak; # padding the spike amplitude v[i+1]=cn; # membrane voltage reset w[i+1]=w[i+1]+dn;} # recovery variable reset } # # Selected output to evaluate solution cat(sprintf("\n\n t v w" )); cat(sprintf("\n %5.1f%10.4f%10.4f",t0[1] ,v[1] ,w[1] )); cat(sprintf("\n %5.1f%10.4f%10.4f",t0[251] ,v[251] ,w[251] )); cat(sprintf("\n %5.1f%10.4f%10.4f",t0[501] ,v[501] ,w[501] )); cat(sprintf("\n %5.1f%10.4f%10.4f",t0[751] ,v[751] ,w[751] )); cat(sprintf("\n %5.1f%10.4f%10.4f",t0[1001],v[1001],w[1001])); cat(sprintf("\n\n")); # # Calls to neuron_1 cat(sprintf("\n\n ncall = %5d\n\n",ncall)) # # Three plots for v(t),w(t) vs t, w(t) vs v(t) # # v(t) par(mfrow=c(1,1)) plot(t0,v,xlab="t",ylab="v(t)",xlim=c(0,1000), type="l",lwd=2,main="v(t), rkf45 integration") # # w(t) par(mfrow=c(1,1)) plot(t0,w,xlab="t",ylab="w(t)",xlim=c(0,1000), type="l",lwd=2,main="w(t), rkf45 integration") # # w(t) vs v(t) (phase plane plot) par(mfrow=c(1,1)) plot(v,w,xlab="v(t)",ylab="w(t)", type="l",lwd=2,main="w(t) vs v(t), rkf45 integration")