Plotting Percent Activity of Artificial Cell Network

[tseredynski]

I have established working hoc code for a network model based on some published connectivity data. I am interested in not just generating the rasterplot, but would like to see a plot of the % of currently firing neurons. I thought I had figured out a way of doing so by manipulating the rasterplot code in Chapter 11 of the NEURON book, by using .histogram to attempt to bin each vector in the variable spikes - followed by a plotactivity() protocol that attempts to sum the bins together and normalize the values. The modified rig.hoc file is shown here:

Code: Select all

sample_rate = 1 //ms	// NEW

objref netcon, vec, spikes, nil, graster
objref spikeBins, activityPlot, vec2	// NEW

proc preprasterplot() {
	spikes = new List()
	spikeBins = new List()	// NEW
	for i=0, cells.count()-1 {
		vec = new Vector()
		vec2 = new Vector()	// NEW
		netcon = new NetCon(cells.object(i).pp,nil)
		netcon.record(vec)
		spikes.append(vec)
		vec2 = spikes.object(i).histogram(0, tstop, sample_rate)	// NEW
		spikeBins.append(vec2)	// NEW
	}
	objref netcon, vec, vec2	// NEW
	
	graster = new Graph(0)
	graster.view(0,0, tstop, cells.count(), 300, 105, 300.48, 500)
	
	activityPlot = new Graph(0)	// NEW
	activityPlot.view(0,0, tstop, 100,500, 105, 300.48, 700)	// NEW
}

preprasterplot()

objref spikey

proc showraster() {
	graster.erase_all()
	for i=0, cells.count()-1 {
		spikey = spikes.object(i).c
		spikey.fill(i+1)
		spikey.mark(graster, spikes.object(i), "|", 6)
	}
	objref spikey
}

	// BELOW ALL NEW

objref percent, x_values

proc showactivity(){
	activityPlot.erase_all()
	percent = spikeBins.object(0).c
	percent.fill(0)
	for i=0, cells.count()-1{
		percent.add(spikeBins.object(i))
	}
	percent.div(cells.count())
	x_values = spikeBins.object(0).c
	x_values.indgen(sample_rate)
	percent.line(activityPlot, x_values)
}

	// END NEW

proc run() {
	stdinit()
	continuerun(tstop)
	showraster()
	showactivity()	// NEW
}

I tried to highlight the parts I added with //NEW.

Unfortunately nothing plots on activityPlot, and when further investigated I find that vec2 in prepraster() never contains values other than 0 (tested with a for-loop):

Code: Select all

proc preprasterplot() {
	spikes = new List()
	spikeBins = new List()
	for i=0, cells.count()-1 {
		vec = new Vector()
		vec2 = new Vector()
		netcon = new NetCon(cells.object(i).pp,nil)
		netcon.record(vec)
		spikes.append(vec)
		vec2 = spikes.object(i).histogram(0, tstop, sample_rate)
		spikeBins.append(vec2)
		for j = 0, vec2.size()-1 {                                    //TEST
			if(vec2.x(j)!=0) print j
		}
	}
	objref netcon, vec, vec2
	
	graster = new Graph(0)
	graster.view(0,0, tstop, cells.count(), 300, 105, 300.48, 500)
	
	activityPlot = new Graph(0)
	activityPlot.view(0,0, tstop, cells.count(),500, 105, 300.48, 700)
}

I'm at a real loss as to what I'm doing wrong here with trying to use .histogram. I attempted running the same code just in nrniv and the concept seems like it should work... are there any glaring mistakes you can see?

Thank you!

[ted]

An aside: presumably your network is built from aritificial spiking cells.

To produce a nonzero result, histogram() must be executed on Vectors that have nonzero elements, i.e. after a run has generated some spike events.

A comment: if your networks get large, you're going to have a lot of Vectors. Your code will scale better if a unique integer (identifier) is assigned to each cell, so that all spike times and associated cell ids can be recorded to a single pair of vectors tvec and idvec. You'll find examples in our paper
Hines, M.L. and Carnevale, N.T.
Translating network models to parallel hardware in NEURON.
J. Neurosci. Methods 169:425-455, 2008
which is available from http://www.neuron.yale.edu/neuron/nrnpubs.
This style lends itself nicely to raster plots (for each element in tvec plot a mark at coordinates (tvec.x, idvec.x+1). Furthermore, the fact that all event times are already aggregated in tvec eliminates some of the steps involved in generating a "percent activity" graph. You might also consider convolving the aggregated tvec with a window function (boxcar for moving average, tophat for weighted average, gaussian etc.).

[tseredynski]

Thank you for the suggestions. I'll give that a shot and see how it goes.