Code: Select all
// $Id: myizh.hoc,v 1.4 2010/01/26 19:20:38 billl Exp $
// load_file("izh.hoc")
load_file("stdgui.hoc")
nrnmainmenu()
objref g[10]
tstop=500
//* setup the cell
create acell
access acell
objref izh,stim,nc,fih[2]
cvode_local(1)
acell izh=new IZH(0.5)
//* parameters for different cell types
objref avec,bvec,cvec,dvec,vviv,tstv,butnl,bub
avec=new Vector()
{bvec=avec.c cvec=avec.c dvec=avec.c vviv=avec.c tstv=avec.c}
// parameters a,b,c,d for different models
avec.append(0.02,0.02,0.02,0.02,0.02,0.01,0.02,0.2,0.02,0.05,0.1,0.02,0.03,0.03,0.03,0.1,1,0.02,-0.02,-0.026)
bvec.append(0.2,0.25,0.2,0.25,0.2,0.2,-0.1,0.26,0.2,0.26,0.26,-0.1,0.25,0.25,0.25,0.26,0.2,1,-1,-1)
cvec.append(-65,-65,-50,-55,-55,-65,-55,-65,-65,-60,-60,-55,-60,-52,-60,-60,-60,-55,-60,-45)
dvec.append(6,6,2,0.05,4,8,6,0,6,0,-1,6,4,0,4,0,-21,4,8,-2)
objref playvec,playtvec
{playvec=new Vector() playtvec=new Vector()}
vviv.append(-70,-64,-70,-64,-70,-70,-60,-64,-70,-62,-62,-60,-64,-64,-64,-61,-70,-65,-63.8,-63.8)
tstv.append(100,200,220,200,160,85,300,300,100,200,400,100,200,200,100,300,50,400,350,350)
butnl=new List()
for ii=0,19 butnl.append(new String())
{butnl.object(0).s="tonic spiking" butnl.object(1).s="phasic spiking" butnl.object(2).s="tonic bursting" butnl.object(3).s="phasic bursting" butnl.object(4).s="mixed mode" butnl.object(5).s="spike frequency adaptation" butnl.object(6).s="Class 1" butnl.object(7).s="Class 2" butnl.object(8).s="spike latency"}
{butnl.object(9).s="subthreshold oscillations" butnl.object(10).s="resonator" butnl.object(11).s="integrator" butnl.object(12).s="rebound spike" butnl.object(13).s="rebound burst" butnl.object(14).s="threshold variability" butnl.object(15).s="bistability" butnl.object(16).s="Depolarizing afterpotential" butnl.object(17).s="accomodation" butnl.object(18).s="inhibition-induced spiking" butnl.object(19).s="inhibition-induced bursting"}
proc prpars () { local ix
if (numarg()==1) ix=$1 else ix=rnum
printf("\t%s\n",butnl.object(ix).s)
printf("a:%g\tb:%g\tc:%g\td:%g\tvv0:%g\n",avec.x[ix],bvec.x[ix],cvec.x[ix],dvec.x[ix],vviv.x[ix])
}
//* run routine
proc p () { local ix
ix=rnum=$1 // global rnum
izh.a=avec.x[ix] izh.b=bvec.x[ix] izh.c=cvec.x[ix] izh.d=dvec.x[ix]
tstop=tstv.x[ix]
cvode_simgraph()
g.size(0,tstop,-100,50)
playinit()
run()
}
proc playinit () { local T1,ix
ix=rnum
izh.f=5 izh.g=140 // standard params: vv'=0.04*vv^2+5*vv+140-u+I
sprint(bub.label,"%s (%c) #%d",butnl.object(ix).s,65+ix,ix)
if (ix==16) sprint(bub.label,"%s -- REPEATED SPIKING",bub.label)
if (ix==17) sprint(bub.label,"%s -- NOT IMPLEMENTED (different functional form;see izh.mod)",bub.label)
if (ix==19) sprint(bub.label,"%s -- NOT IMPLEMENTED (convergence problems)",bub.label)
g.erase_all
g.addvar("izh.vv",2,2)
g.label(0.1,0.9,bub.label)
playvec.play_remove()
playtvec.resize(0) playvec.resize(0)
if (ix==6) {
T1=30
playtvec.append(0,T1,tstop)
playvec.append(0,0,0.075*(tstop-T1))
} else if (ix==7) { // (H) Class 2 exc.
T1=30
playtvec.append(0,T1,tstop)
playvec.append(-0.5, -0.5,-0.05+0.015*(tstop-T1))
} else if (ix==17) { // (R) accomodation
playtvec.append(0, 200, 200.001, 300, 312.5, 312.501, tstop)
playvec.append( 0, 200/25, 0 , 0 , 4 , 0 , 0)
}
if (ix==6 || ix==7 || ix==17) playvec.play(&izh.I,playtvec,1)
if (ix==6 || ix==11) { izh.f=4.1 izh.g=108 }
}
//* box of buttons
begintemplate bubox
public label
external p
objref hbox,vbox
strdef tstr,label
func transpose () { return int($1/rows) + $1%rows*cols }
endtemplate bubox
//* plotting & printing
// newPlotV()
// g=graphItem
// g.erase_all
// g.addvar("izh.vv",2,2)
// g.addvar("izh.u",3,1)
// g.addvar("izh.I",4,2) // will show false ramps
// nrnpointmenu(izh)
// bub = new bubox(butnl)
//* initialization
fih=new FInitializeHandler("uvvset()")
fih[1]=new FInitializeHandler(0,"Isend()")
// initialization routines
proc uvvset () {
izh.vv=vviv.x[rnum]
izh.u=izh.vv*izh.b
if (rnum==17) izh.u=-16 // example 17 also requires different mod file
}
// current injections for specific models
proc Isend () { local T1,T2,T3,T4
T1=tstop/10
izh.I=0
if (rnum==0) { // (A) tonic spiking
Isend1(T1,14)
} else if (rnum==1) { // (B) phasic spiking
T1=20
Isend1(T1,0.5)
} else if (rnum==2) { // (C) tonic bursting
T1=22
Isend1(T1,15)
} else if (rnum==3) { // (D) phasic bursting
T1=20
Isend1(T1,0.6)
} else if (rnum==4) { // (E) mixed mode
Isend1(T1,10)
} else if (rnum==5) { // (F) spike freq. adapt
Isend1(T1,30)
} else if (rnum==6) { // (G) Class 1 exc. -- playvec
} else if (rnum==7) { // (H) Class 2 exc. -- playvec
} else if (rnum==8) { // (izh.I) spike latency
Isend1(T1,7.04)
Isend1(T1+3,0.0)
} else if (rnum==9) { // (J) subthresh. osc.
Isend1(T1,2)
Isend1(T1+5,0)
} else if (rnum==10) { // (K) resonator
T2=T1+20 T3 = 0.7*tstop T4 = T3+40
Isend1(T1,0.65) Isend1(T2,0.65) Isend1(T3,0.65) Isend1(T4,0.65)
Isend1(T1+4,0.) Isend1(T2+4,0.) Isend1(T3+4,0.) Isend1(T4+4,0.)
} else if (rnum==11) { // (L) integrator
T1=tstop/11 T2=T1+5 T3 = 0.7*tstop T4 = T3+10
Isend1(T1,9) Isend1(T2,9) Isend1(T3,9) Isend1(T4,9)
Isend1(T1+2,0.) Isend1(T2+2,0.) Isend1(T3+2,0.) Isend1(T4+4,0.)
} else if (rnum==12) { // (M) rebound spike
T1=20
Isend1(T1,-15)
Isend1(T1+5,0)
} else if (rnum==13) { // (N) rebound burst
T1=20
Isend1(T1,-15)
Isend1(T1+5,0)
} else if (rnum==14) { // (O) thresh. variability
T1=10 T2=70 T3=80
Isend1(T1,1) Isend1(T2,-6) Isend1(T3,1)
Isend1(T1+5,0.) Isend1(T2+5,0.) Isend1(T3+5,0.)
} else if (rnum==15) { // (P) bistability
T1=tstop/8 T2=216
izh.I=0.24
Isend1(T1,1.24) Isend1(T2,1.24)
Isend1(T1+5,0.24) Isend1(T2+5,0.24)
} else if (rnum==16) { // (Q) DAP depolarizing afterpotential
T1 = 10
Isend1(T1-1,20)
Isend1(T1+1,0)
} else if (rnum==17) { // (R) accomodation -- playvec
} else if (rnum==18) { // (S) inhibition induced spiking
izh.I=80
Isend1(50,75)
Isend1(250,80)
} else if (rnum==19) { // (T) inhibition induced bursting
izh.I=80
Isend1(50,80) // Isend1(50,75) -- will crash simulator
Isend1(250,80)
}
}
proc Isend1 () {
sprint(tstr,"izh.I=%g cvode.re_init()",$2)
cvode.event($1,tstr)
}
// izhstim() sets up a single stim into izh cell
// effect easily seen by running "Class 1" -- p(6)
proc izhstim () {
stim=new NetStim(0.5)
stim.number = stim.start = 1
nc = new NetCon(stim,izh)
nc.delay = 2
nc.weight = 0.1
izh.erev = -5
}
//* ==== net.hoc starts ==== Declarations
objref g,rdm[2],ind,XO,YO,animv[2],tmpobj,tvec,vec[10],scr[2],tmpvec,veclist,ind
objref cells, nil
objref nclist,netcon
objref sl,gx[2],drv
{cells = new List() nclist = new List()}
for ii=0,1 rdm[ii] = new Random()
for ii=0,9 vec[ii]=new Vector()
for ii=0,1 scr[ii]=new Vector()
for ii=0,1 animv[ii]=new Vector()
tvec=new Vector()
ind=new Vector()
veclist=new List()
sl = new SectionList() // empty list
drv = new Vector() // drawing vector
Gshp=0
tveclen=0
//* Cell template
begintemplate Cell
public M
public pp, connect2target, x, y, z, position, is_art
external acell
objref pp
proc init () {
acell pp = new IZH(0.5)
// pp.tau = 10
// pp.invl = 20
}
func is_art () { return 1 }
proc connect2target () { $o2 = new NetCon(pp, $o1) }
proc position (){ x=$1 y=$2 z=$3}
func M () {
// return pp.M ()
}
endtemplate Cell
//* Network specification
rcell = 0
func cell_append() {cells.append($o1) $o1.position($2,$3,$4) return cells.count-1 }
proc cell_append() {cells.append($o1) $o1.position($2,$3,$4) }
func nc_append () { //srcindex, tarcelindex, synindex
if ($3 >= 0) {
cells.object($1).connect2target(cells.object($2).synlist.object($3),netcon)
netcon.weight = $4 netcon.delay = $5
} else {
cells.object($1).connect2target(cells.object($2).pp,netcon)
netcon.weight = $4 netcon.delay = $5
}
nclist.append(netcon)
return nclist.count - 1
}
//** createnet()
proc createnet () { local i, j
netcon = nil
nclist.remove_all()
cells.remove_all()
for i=0, ncell-1 cell_append(new Cell(), i, 0, 0)
wire()
}
//** wire():: full non-self connectivity
// artificial cell templates have obj.pp
// params: ncell
// creates nclist: list of NetCons
proc wire () {
nclist.remove_all()
for i=0,ncell-1 for j=0,ncell-1 if (i!=j) {
netcon = new NetCon(cells.object(i).pp,\
cells.object(j).pp)
nclist.append(netcon)
}
}
//** ranwire()
proc ranwire () { local proj,beg,end
nclist.remove_all()
rdm.discunif(0,ncell-1)
for i=0,ncell-1 {
proj=rdm.repick()
if (proj<ncell/2) { beg=proj end=proj+int(ncell/2)-1
} else { beg=proj-int(ncell/2)+1 end=proj }
for j=beg,end if (i!=j) {
netcon = new NetCon(cells.object(i).pp,\
cells.object(j).pp)
nclist.append(netcon)
}
}
}
//* Parameters
//** parameters
tstop = 500
ncell = 10
ta=10
w=-0.01
del=4
low=10
high=11
//** routines: weight(),delay(),tau(),interval()
proc weight () { local i localobj vv
vv=new Vector()
if (numarg()==0) {
for i=0, nclist.count-1 vv.append(nclist.object(i).weight)
print vv.min,vv.max,vv.mean,vv.stdev
} else {
w = $1
for i=0, nclist.count-1 nclist.object(i).weight = w
}
}
//** weight2(WT,EXCLUDE_VEC) :: set weight to WT
// unless in EXCLUDE_VEC then set wt. to 0
proc weight2 () { local i,ww
w = $1
for i=0,nclist.count-1 {
if ($o2.contains(i)) ww=0 else ww=w
nclist.object(i).weight = ww
}
}
proc delay () { local i
del = $1
for i=0, nclist.count-1 {
nclist.object(i).delay = $1
}
}
proc tau () { local i
ta = $1
for i=0, cells.count-1 {
// cells.object(i).pp.tau = $1
}
}
//** interval(low,high) randomly sets cells to have periods between low and high
proc interval () { local i, x, dx // args are low and high
low = $1
high = $2
rdm.uniform(low,high)
vec.resize(ncell)
vec.setrand(rdm)
for i=0, cells.count-1 {
cells.object(i).pp.I = vec.x[i]
}
}
//** setparams() sets weight, delay, tau and intervals
proc setparams () { local ii
if (numarg()>=3) {w=$1 del=$2 ta=$3}
if (numarg()>=5) {low=$4 high=$5}
weight(w)
delay(del)
// tau(ta)
interval(low, high)
}
//* Run code
//** savspks() -- save to a vector
proc savspks () { local ii
for ii=0,1 animv[ii].resize(0)
for ii=0,cells.count-1 {
XO=cells.object(ii)
if (cvode.netconlist(XO, "", "").count==0) {
tmpobj=new NetCon(XO.pp, nil)
} else {
tmpobj=cvode.netconlist(XO, "", "").object(0)
}
tmpobj.record(tvec,ind)
animv[0].append(objnum(tmpobj))
animv[1].append(objnum(tmpobj.pre)) // in this case are all in a row anyway
}
}
//** showspks() -- show spikes on graph g
proc markspks () { local ii,x
g.erase
ind.mark(g,tvec,"O",4,2,1)
g.flush()
}
proc showspks () { local ii,x
// ind.mark(g,tvec,"O",8,2,1)
for ii=0,tvec.size-1 {
x=tvec.x[ii]
g.beginline(3,1) g.line(x,ind.min) g.line(x,ind.max)
}
g.flush()
}
//** syncer() :: returns sync measure 0 to <1
// measures how well spikes "fill up" the time
// assumes spike times in tvec, tstop
// param: width
// syncer doesn't take account of prob of overlaps
// due to too many spikes stuffed into too little time
func syncer () { local t0,tt,cnt,width
t0=-1 width=1 cnt=0
for ii=0,tvec.size-1 {
tt=tvec.x[ii]
if (tt>=t0+width) {t0=tt cnt+=1}
}
return 1-cnt/(tstop/width)
}
//** autorun() changes weights
proc autorun () { local pij,x
veclist.remove_all()
rvc() // clear the storage vectors
for case(&x,-0.5,-0.3,-0.1,-0.01,-0.001) {
weight(x)
run()
savevec(ind,tvec)
print w,syncer()
vec[1].append(w) vec[2].append(syncer())
}
}
//** autorun1() changes connections density
proc autorun1 () { local ii,S,pij_inc,C
veclist.remove_all()
rvc()
max=-0.1 pij_inc=0.1
S=ncell*ncell
vec[3].resize(0)
for ii=0,9 {
C = (1-ii*pij_inc) // percent convergence
w=max/C // scale weight up as convergence goes down
setparams()
weight2(w,vec[3])
run()
savevec(ind,tvec)
print S-vec[3].size,syncer()
vec[1].append((S-vec[3].size)/S) vec[2].append(syncer())
rdm.discunif(0,S-1)
rdmunq(vec[3],0.1*S,rdm) // increase those set to 0
}
}
//* Utility functions
//** savevec(vec1,vec2,...) add vectors onto veclist
proc savevec () { local i
for i=1, numarg() {
tmpvec = new Vector($oi.size)
tmpvec.copy($oi)
veclist.append(tmpvec)
tmpvec = nil
}
}
//** rvc() clears vec[0..9]
proc rvc () { for ii=0,9 vec[ii].resize(0) }
//** setdensity(pij) sets connection density to 0<pij<1
proc setdensity () { local ii,S,pij,C,max
max=-0.1 pij=$1
S=nclist.count
rdm.discunif(0,S-1)
vec[3].resize(0)
rdmunq(vec[3],(1-pij)*S,rdm) // number to set to 0
C = (1-pij) // percent convergence
w=max/C // scale weight up as convergence goes down
setparams()
weight2(w,vec[3])
}
//** rdmunq(vec,n,rdm) -- augment vec1 by n unique vals from rdm
proc rdmunq () { local n,num,flag,xx,loop
n=$2 num=0 flag=1 loop=0
scr.resize(n*4) // hopefully will get what we want
while (flag) {
scr.setrand($o3)
for ii=0,scr.size-1 {
xx=scr.x[ii]
if (! $o1.contains(xx)) { $o1.append(xx) num+=1 }
if (num==n) { flag=0 break }
}
loop+=1
if (loop==10) { print "rdmunq ERR; inf loop" flag=0 break }
}
}
//** rdmord (vec,n) randomly ordered numbers 0->n-1 in vec
proc rdmord () { local n
n=$2
rdm.uniform(0,100)
scr.resize(n)
scr.setrand(rdm)
scr.sortindex($o1)
}
// vcount (num,vec)
func vcount () {
scr.where($o2,"==",$1)
return scr.size
}
//* Mapping functions
//** objnum(OBJ) -- find object number
func objnum () {
sprint(temp_string_,"%s",$o1)
if (sscanf(temp_string_,"%*[^[][%d]",&x) != 1) x=-1
return x
}
//** getcnum(CELL_OBJ) return index given cell object
func getcnum () {
sprint(tstr,"%s",$o1)
if (sscanf(tstr,"IntervalFire[%d]",&x) != 1) x=-1
return x
}
//** fconn(PREVEC,POSTVEC) places values of
// pre- and post-syn cells in parallel vectors
// only lists pairs with non-zero connections
// getcnum() returns index of cell obj
proc fconn () {
$o1.resize(0) $o2.resize(0)
for ii=0,nclist.count-1 {
XO=nclist.object(ii)
if (XO.weight!=0) {
$o1.append(getcnum(XO.pre))
$o2.append(getcnum(XO.syn))
}
}
}
//** showconns() -- show all the connections as line segments
proc showconns () { local pr,po
g.erase
fconn(scr[0],scr[1])
for ii=0,scr.size-1 {
pr=scr.x[ii] po=scr[1].x[ii]
drawline(pr,po,10)
}
g.flush()
}
//** showconv1() -- show convergence to one cell as line seg
proc showconv1 () { local pr,po,id,colr
id=$1
if (numarg()==2) colr=$2 else colr=2
fconn(scr[0],scr[1])
for ii=0,scr.size-1 {
pr=scr.x[ii] po=scr[1].x[ii]
if (po==id) { drawline(pr,po,10,colr,3) printf("%d ",pr)}
}
print ""
g.flush()
}
//** showdiv1() -- show divergence from one cell as line seg
proc showdiv1 () { local pr,po,id,colr
id=$1
if (numarg()==2) colr=$2 else colr=2
fconn(scr[0],scr[1])
for ii=0,scr.size-1 {
pr=scr.x[ii] po=scr[1].x[ii]
if (pr==id) { drawline(pr,po,10,colr,8) printf("%d ",po) }
}
print ""
g.flush()
}
//*** pos (cell_ind,&xi,&yi,cols) -- sets xi, yi to location on the grid
proc pos () { local ci,cols
ci=$1 cols=$4
$&2 = ci%cols $&3=int(ci/cols)
print $&2,$&3
}
func xpos () { return $1%$2 }
func ypos () { return int($1/$2) }
//*** func distn () calc distance
func distn () { local c1,c2,xd,yd,cols
c1=$1 c2=$2 cols=$3
xd=xpos(c1,cols)-xpos(c2,cols)
yd=ypos(c1,cols)-ypos(c2,cols)
return sqrt(xd*xd+yd*yd)
}
//*** distwire(pij)
proc distwire () { local syn,prob,maxdist,cnt,total,C,loop
pij=$1
allsyns = nclist.count
total=pij*allsyns // how many syns to set
rdm[1].uniform(0,1) // for flipping coin
// maxdist==12.728 for 10x10; mindist=1 (neighbors)
maxdist=0.33*distn(0,ncell-1,sqrt(ncell)) // the full dist from lower left to upper right
maxwt=-0.9/pij // norm wt by convergence
loop=cnt=0 // counters
for i=0,allsyns-1 nclist.object(i).weight = 0 // clear weights
while (cnt<total && loop<4) {
rdmord(vec[3],allsyns) // test each synapse in random order
for ii=0,vec[3].size-1 {
XO=nclist.object(vec[3].x[ii]) // pick a synapse
// max prob of connection is 0.8*(1-mindist/maxdist)~74% for 10x10
// zero prob of diag connection from corner to corner
prob = 1.0*(1-(distn(getcnum(XO.pre),getcnum(XO.syn),sqrt(ncell))/maxdist))
if (rdm[1].repick<prob) {
XO.weight=maxwt
cnt+=1
}
if (cnt>=total) break // finished
}
}
print cnt,total
if (cnt<total) printf("distwire ERR: target %d, set %d\n",total,cnt)
}
//*** drawline(beg,end,columns[,color,line_width])
proc drawline () { local beg,end,cols,clr,lwid
beg=$1 end=$2 cols=$3
if (numarg()>=4) clr=$4 else clr=4
if (numarg()>=5) lwid=$5 else lwid=1
g.beginline(clr,lwid)
g.line(xpos(beg,cols),ypos(beg,cols))
g.line(xpos(end,cols),ypos(end,cols))
}
//* Animation
//** animplot() put up the shape plot for hinton diagram
proc animplot () {
gx[Gshp] = new PlotShape(sl,0)
flush_list.append(gx[Gshp])
ctern(gx[Gshp],0,2)
gx[Gshp].view(1,0,1.1,1.1,500,200,100,100)
drawcells()
}
//* Ternary color map
proc ctern () {
$o1.colormap(3)
$o1.colormap(0, 255, 0, 0)
$o1.colormap(1, 255, 255, 0)
$o1.colormap(2, 0, 0, 255)
$o1.scale($2, $3)
}
//** drawcells() draw squares of hinton diagram
proc drawcells () { local i,j,ny,nx
if (ncell!=100) { print "ERROR: drawcells() currently written for ncell=100" return }
gx[Gshp].erase_all
drv.resize(ncell)
xoff=1.1 yoff=0.1 wdt=0.1
nx=ny=10
for i=0,nx-1 for j=0,ny-1 gx[Gshp].hinton(&drv.x[j*nx+i],(i+.5)*wdt+xoff,(j+.5)*wdt+yoff,wdt)
gx[Gshp].size(1, 2.2, 0, 1.2)
gx[Gshp].exec_menu("Shape Plot")
}
//** chkhint(cell#) light up a location for a single cell
proc chkhint () { drv.fill(0) drv.x[$1]=2 gx.flush() }
//** anim() animates sim stored in tvec,ind
proc anim () { local tt,tstep,ii,jj,sz
tstep=0.1
sz = ind.size-1
gx[Gshp].exec_menu("Shape Plot")
scr.copy(tvec)
scr.add(500*tstep) // how many steps to keep it illuminated
drv.fill(0)
ii=jj=0
for (tt=0;tt<=tstop;tt+=tstep) {
for (;ii<sz && tt>tvec.x[ii];ii+=1) drv.x[animv.indwhere("==",ind.x[ii])]=2
for (;jj<sz && tt> scr.x[jj];jj+=1) drv.x[animv.indwhere("==",ind.x[jj])]=0
gx[Gshp].flush
doEvents()
}
}
//* Run sequences
cvode_active(1)
proc runseq () {
createnet(ncell)
setparams()
savspks()
run()
}
// ncell=10
// w=-1e-6
// runseq()
// g=new Graph()
// g.erase_all
// markspks()
// showspks()
// load_file("net.hoc")
ncell=100
createnet()
savspks()
w=-1e-5
setparams()
// run() // don't repeate runseq() unless want to change # of cells
// animplot()
// anim()