Fano factor

The basics of how to develop, test, and use models.
Post Reply
oren
Posts: 55
Joined: Fri Mar 22, 2013 1:03 am

Fano factor

Post by oren »

Hello, I am trying to see if the synaptic events that are created from netstim with noise=1 have a Fano Factor http://en.wikipedia.org/wiki/Fano_factor of 1 (as should be for a Poisson distribution)
This is my code

Code: Select all

// Test Fano factor

objref syn, sti , con, Vec, Gra, Vec1
create soma
insert pas
Vec = new Vector()
syn = new Exp2Syn(0.5)
sti = new NetStim(0.5)
sti.interval = 100
sti.start 	 = 10
sti.number   = 500000
sti.noise     = 1
con  = new NetCon(sti,syn)
con.weight=1
con.record(Vec)  // record all events of synaptic release 

tstop = 55000
run()

Vec1 = new Vector()
for i=1,(Vec.size()-1){  // getting only the difference between each synaptic event
	Vec1.append(Vec.x[i] - Vec.x[i-1])
}

Vec1.var()/Vec1.mean()

The result that I should get from Vec1.var()/Vec1.mean() should be ~1
as We can see from the simple matlab code

Code: Select all

rand_sm = poissrnd(1650,1,70);
var(rand_sm)/mean(rand_sm)
But I get form Vec1.var()/Vec1.mean() = 91.550064

What am I doing wrong?

Thank You.
ted
Site Admin
Posts: 5786
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Fano factor

Post by ted »

What am I doing wrong?
Misinterpreting the definition of Fano factor. It applies to discrete distributions, e.g. counts of events that occur within a particular time interval, such as spike counts. It does not apply to continuous random variables, such as interspike intervals. A NetStim with random=1 generates events at intervals that are described by the negexp (negative exponential) distribution, which is a continuous distribution. Over a particular time interval, the number of spikes should be found to follow the Poisson distribution.

"But I have seen/heard neuroscientists describe spike times as being Poisson!"

As Orwell wrote, "the slovenliness of our language makes it easier for us to have foolish thoughts." He was writing about politics, but the principle applies much more generally.
oren
Posts: 55
Joined: Fri Mar 22, 2013 1:03 am

Re: Fano factor

Post by oren »

Thank You Ted.

I have fixed the code to get so now it give the right results.

Code: Select all

// Test Fano factor

objref syn, sti , con, Vec, Vec1, hist
create soma
insert pas
Vec = new Vector()
syn = new Exp2Syn(0.5)
sti = new NetStim(0.5)
sti.interval = 100
sti.start 	 = 10
sti.number   = 500000
sti.noise     = 1

con  = new NetCon(sti,syn)
con.weight=1
con.record(Vec)


tstop = 35000
timewindow = 250 // Should work for any time window!
run()

Vec1 = new Vector()
for i=1,(Vec.size()-1){
	Vec1.append(Vec.x[i] - Vec.x[i-1])
}

hist = Vec.histogram(0,tstop,20)
hist.var()/hist.mean()  //This is the Fano factor --> Should be 1 for possion...

Vec1.stdev()/Vec1.mean() //Should be 1 for possion...

ted
Site Admin
Posts: 5786
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Fano factor

Post by ted »

Good work. Here's a further revision that runs faster yet produces the same results.

Code: Select all

objref sti
sti = new NetStim()
sti.interval = 100
sti.start = 10
sti.number = 1e9 // no need to limit this
sti.noise = 1

objref con, nil, vec
con  = new NetCon(sti,nil)
vec = new Vector()
con.record(vec)

tstop = 35000
cvode_active(1)
run()

objref isivec, hist

proc analyze() {
  isivec = new Vector()
  isivec.deriv(vec,1,1) // Euler deriv with dx == 1 yields vector of ISIs
  print "mean ISI ", isivec.mean(), "based on ", isivec.size(), "samples"
  print "ISI stdev ", isivec.stdev()
  print "stdev/mean ", isivec.stdev()/isivec.mean()
  print " "
  hist = vec.histogram(0,tstop,20)
  print "estimate from spike counts over 20 ms intervals: ", hist.var()/hist.mean()
}

analyze()
Comments:
1. The biophysical model cell requires numerical integration but adds nothing to the statistical evaluation. Just use the NetStim as an artificial spiking cell.
2. After the biophysical model cell has been eliminated, adaptive integration will execute extremely fast because the simulation is now an event-driven simulation (no numerical integration involved).
3. There's no need to limit the maximum number of spikes generated by the NetStim. Specifying a large number allows you to explore a much wider range of simulation times and total spike counts.
4. Vector class's deriv method using Euler method and "dx" parameter == 1 transforms recorded spike times to interspike intervals at machine language speeds.

You might try a run with tstop 1e6 ms and NetStim interval 1 ms, just to see how fast event-driven simulation can be.
Just execute
tstop = 1e6
sti.interval = 1
run()
analyze()
in that sequence at the interpreter's oc> prompt.
ted
Site Admin
Posts: 5786
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Fano factor

Post by ted »

Two more items:
1. I'd insert
load_file("stdgui.hoc")
as the first line in the file. This ensures that the standard run library is loaded. Otherwise MSWin users may find that double clicking on the hoc file will generate an error message about cvode_active()
2. cvode_active() is a procedure that is included in the standard run library (i.e. not built into hoc). It calls cvode.active and does several other things that ensure the standard run system works properly when adaptive integration is used.

"What other things?"

Read lib/hoc/stdrun.hoc to find out.
Post Reply