# # 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 # # Library of R ODE solvers library("deSolve") # # ODE routine setwd("c:/R/bme_ode/chap4") source("euler.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; # # Euler integration for(i in 1:nout){ yout=euler(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)) # # Analysis of derivatives dvdt=rep(0,nout+1); dwdt=rep(0,nout+1); out=rep(0,2); for(i in 1:nout+1){ out=neuron_1(t0[i],c(v[i],w[i])); dvdt[i]=out[1]; dwdt[i]=out[2]; } # # 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), Euler 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), Euler 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), Euler integration") # # Two plots for dv(t)/dt, dw(t)/dt vs t # # dv/dt(t) vs v(t) par(mfrow=c(1,1)) plot(t0,dvdt,xlab="t",ylab="dv(t)/dt", type="l",lwd=2,main="dv(t)/dt") # # dw/dt(t) vs v(t) par(mfrow=c(1,1)) plot(t0,dwdt,xlab="t",ylab="dw(t)/dt", type="l",lwd=2,main="dw(t)/dt")