# first passage set up to compute the firing times # -1 = f(V) T'(V) + sigma^2/2 T''(V) -oo < V < Vspike # we desire T(Vreset) but we cannot get this directly # from AUTO, so we will instead define two variables u,w # defined on [0,1] and scaled so that u(1)=w(0) = T(Vreset) # # # this is a bit weird, but here goes: # u is defined on -a < V < Vreset # a is sort of like infinity. b=Vreset+a is the length of this interval # w is defined on Vreset < V < VSpike c = (Vspike-Vreset) is length of this # u(s=1)=w(s=0) # note d/dV = d/d(bs) for b=(Vreset+a) # d/dV = d/d(cs) for c=(Vspike-Vreset) # c u'(s=1)=b w'(s=0) # now what about the BC at V=-a. # easiest is to make T'(-a)=0 # or u'(s=0)=0 # # for better accuracy, T'(V) f(V) ~ -1 for V large # so for QIF T'(V)= -1/a^2 # and thus u'(s=0) = -b/a^2 # # clearly # w(s=1)=0 since T(Vspike)=0 # # this is set up for sigma=1, I=50 easy to solve! # for I large and negative or sigma small, XPP cannot solve it # using shooting, but AUTO can. Below it is set up for AUTO # with struggle AUTO will take it down to almost -2 # # in XPP: # Click on Bndry Val Show # File AUTO # in AUTO # Run BoundryValue # youll get the FI curve # QIF: f(x)=x^2+I init u=1,w=1,wp=-2,freq=1 par I=50,sig=1,vreset=-10,vspike=10,a=20 !b=(vreset+a) !c=(vspike-vreset) du/dt=up dup/dt=-b*b/sig-f(-a+b*s)*up*b/sig dw/dt=wp dwp/dt=-c*c/sig-f(vreset+c*s)*wp*c/sig ds/dt=1 dfreq/dt=0 b up b w' b c*up'-b*wp b u'-w b s b freq-1/w # I add the actual values of V here to plot the whole thing # on the real V scale aux vleft=-a+b*s aux vright=vreset+c*s # numerics stuff @ total=1,dt=.005,meth=cvode,tol=1e-8,atol=1e-8 # plot both segments @ nplot=2 @ yp=u,yp2=w,xp=vleft,xp2=vright,xlo=-10,xhi=10,ylo=0,yhi=10 # AUTO parameters @ NTST=60,DSMIN=1e-4,NPR=200,NMAX=1000,DS=-.01 @ PARMIN=-2,PARMAX=60 # Assumed here that I will be the parameter to get "f-I" curve # freq=1/w(0) will be plotted - this is 1/T(Vreset) @ AUTOXMIN=-2,AUTOXMAX=50,AUTOYMAX=10,AUTOVAR=freq,AUTOYMIN=0 done