Izhikevich Model: Connecting Two Cells

[charlie]

I am looking to connect two Izhikevich model cells with a synapse. I can already connect a NetStim input to an Izhikevich cell using NetCon, but I can't access the izh.vv variable to use an Izhikevich cell as the presynaptic cell. If I try to use the izh.vv variable as the first argument of a NetCon I get an error. Are there any examples or models that show how to connect two Izhikevich models, using one as a presynaptic cell?

[wwlytton]

It's often easier to work within the context of a cell template. This is the way that the CellBuilder and NetworkBuilder sets things up so particularly helpful if you want to end up connecting this cell to some other cell type (eg a compartment cell or an event-driven ArtCell)

Anyway, that's how I've set it up in the following example file (myizh.hoc) which was an attempt to get an IZH model to synchronize via the Hopfield-Brody mechanism. This did not succeed with the parameters I used but I did not do much parameter search. If someone gets a chance to do that and finds a successful IZH parameterization please post.

See next post for code

[wwlytton]

reposting code as <code>

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()

[charlie]

I noticed in your forum code post that you used a NetCon for izhikevich cells that simply put the izh object as the first argument for the postsynaptic cell without using &izh or referring to the izh.vv variable. When I set it up as:

NetCon(izh,izh2,0,0,1)

it connects izh as the presynaptic cell to izh2 as the postsynaptic cell and seems to work fine, using vv as the variable to monitor for the threshold. My question is, how does it know to monitor the vv variable from the above message command. If vv is not specified in the NetCon does it simply go to the first STATE variable of izh and use that, which happens to be vv? Or is it monitoring all STATE variables? Im just a little unclear of how that works.

[ted]

In NEURON-implemented network models that use spike-triggered synaptic communication, presynaptic cells ("spike sources") are connected to postsynaptic targets by NetCons. Spike sources fall into three broad categories.

1. Biophysical model cells. These are mechanistic models of biological neurons, which are represented by sections which have membrane currents, and membrane potential which is computed by numerical integration of the cable equation. To use a biophysical model cell as a spike source, one must attach a NetCon to it and, at a minimum, tell that NetCon the name of the variable that it is to monitor for threshold crossing.

2. ARTIFICIAL_CELLs. An ARTIFICIAL_CELL has two chief distinguishing features. First, it is not associated with a section location or numerical integrator; this means that it doesn't have a BREAKPOINT block, refer to v or any ions, or have a POINTER. Second, the only way it can affect, or be affected by, anything outside itself is by sending or receiving events; this means it has a NET_RECEIVE block. For an ARTIFICIAL_CELL to serve as a spike source, its NET_RECEIVE block must contain a net_event() statement. Execution of a net_event() statement will cause events to be delivered, with appropriate latencies and weights, to all targets of the NetCons for which the ARTIFICIAL_CELL is a spike source. Since the spike event is generated directly by calling net_event(), the NetCons don't have to monitor a variable in the ARTIFICIAL_CELL for threshold crossing.

3. POINT_PROCESSes. A POINT_PROCESS may refer to v or ions, have a POINTER coupled to some external variable, or have one or more state variables that require numerical integration--any of which is enough to require it to have a BREAKPOINT block. A POINT_PROCESS may also have a NET_RECEIVE block that allows it to use events to communicate with ARTIFICIAL_CELLs or POINT_PROCESSes. That's what is going on with instances of the IZH class (examine izh.mod and you'll see that its NET_RECEIVE block contains a net_event() statement that is executed when vv crosses a threshold). Since an IZH cell generates spike events by calling net_event(), NetCons that use an IZH cell as their spike source do not have to be told to monitor a variable in the IZH mechanism for threshold crossing.

The "acell" section in this model exists only because POINT_PROCESSes have to be attached to a section. It is merely a "host" for all IZH cell instances. That's just an accident of NEURON's developmental history (POINT_PROCESSes were originally used to represent "point sources" of current like a microelectrode for a current or voltage clamp, or the postsynaptic membrane at a site of synaptic attachment).

[wwlytton]

In the specific case of izh.mod the threshold monitoring is internal to the "cell" and the NetCon does no monitoring

the thresholding is done here:
WATCH (vv>thresh) 2

so the threshold is a RANGE variable that can be set different for each cell -- the default value is
thresh=30

any post-connected NetCon is notified via the net_event() call -- this call is present in any mechanism that directly generates synaptic events

[charlie]

Thank you we have been able to successfully connect two Izh cells. We are now interested in connecting multiple cells to one Izh cell. We are also interested in learning through synaptic plasticity for these connections. Is there anyway to have multiple separate synaptic connections from several cells to one Izh cell so that the weights of those connections can be varied independently?

[wwlytton]

I don't offhand see why there would be any problem with convergence? What difficulty did you run into?

[charlie]

One example is that we want to use both inhibitory connections and excitatory connections. If all other cells are connected to a single Izhi cell through multiple NetCon statements that all have our izhi cell as the target, how will the Izhi cell know in the Net_recieve code block if it should update an inhibitory or excitatory synapse? The synapses do have different time constants of synaptic activation.

[charlie]

Would it be possible to use a nonspecific_current to connect the izhi equation to the cell voltage. If the cell voltage equation has a leak conductance of zero and a capactiance of 1, with an area of 1, then the dV/dt equation would be exactly the izhikevich equation. In this way the izh.mod could just use the equation as the current without a derivative state for vv. vv would instead be a nonspecific_current so that the cell would actually be executing the equation. In this way, synaptic currents could be added to the cell for different synapse types and they would effect the voltage separately. Does this seem to be a possible solution to the problem?

[wwlytton]

in the mod file the weight is handled by summation
so the sign of the weight determines whether the input is excitatory or inhibitory

} else { : synaptic activation
gsyn = gsyn+w
}

I may be missing the point of this thread -- what code are you trying to run that does not run?

[charlie]

We have different dynamics for different types of connections to the cell. Excitatory and inhibitory connections have different time constants of activation and decay. We also have synaptic learning for each of the synaptic connections that differs depending on where the connection is coming from. We are looking for a way to add different kinds of alpha synapses that have varying dynamics of their dual exponential onset and decay, as well as being able to keep track of synaptic plasticity separately for each synapse. Is this possible with the current model setup?

[wwlytton]

take a look at my intf.mod in various models in modeldb
this uses 4 different synapse types for a single net_receive() block