Spikes analysis

Semantic Pointer Architecture
1

Hi,
I am trying to build a model such that it can mimic experimental results of Rat spiking data where subject is performing different tasks. The SPA question example in the provided hbb ch5 tutorials seems relevant at this point. I have following questions:

  1. To probe spiking data from let’s say Motor state, we need to put model.motor.state_ensembles.add_neuron_output(), right?
  2. After probing, I am having spikes at every instant of time (dt=0.001) which will not be useful, as inter-spiking time would be just 0.001 for all neurons at all instances.
    Can anyone help in this matter, how can I make this model more realistic to match-up with experimental data, where inter-spiking-times are usually exponential or heavy tail distributed.
    Thanks
    Gaurav
2

Hi Gaurav,

To answer your questions:

  1. To probe spiking data from let’s say Motor state, we need to put model.motor.state_ensembles.add_neuron_output(), right?
    For the HBB CH5 question-control example, if you want to probe the spikes from the motor state, you’ll want to do:
motor.add_neuron_output()
motor_spike_probe = nengo.Probe(motor.neuron_output)
  1. After probing, I am having spikes at every instant of time (dt=0.001) which will not be useful, as inter-spiking time would be just 0.001 for all neurons at all instances.

This is not an easy question to answer since the inter-spike interval measurement can depend on the specifics of the behaviour you want to replicate. There are a few things you can try to achieve your goal.

=========================

  1. The interspike interval is typically measured for a single neuron. However the spike data gathered using the spike probe above is for all of the neurons in the ensemble. You can isolate the spikes from a single neuron like so:
neuron_0_spikes = motor_spike_probe[:, 0]
neuron_1_spikes = motor_spike_probe[:, 1]

The data returned like this will look something like:
[0, 0, 0, 0, 1000.0, 0, 0, 0, .... 0, 0]

Each element corresponds to whether or not the specific neuron spiked for a specific timestep, where 0 indicates no spikes occurring during that timestep, and a value of 1 / sim.dt when a spike did occur during that timestep. To convert this data into an array of timestamps where a spike happened, you’ll need to multiply by sim.trange() and sim.dt, and then strip out all of the zeros:

neuron_0_spike_times = neuron_0_spikes * sim.trange() * sim.dt  # Mulitply by trange and dt
neuron_0_spike_times = neuron_0_spike_times[neuron_0_spike_times > 0]  # Strip out zeros

Once you have the spike data for individual neurons, you can then proceed to compute the distribution of interspike intervals for each neuron, and then across the entire model to get the distribution for all of the neurons.

=========================
2. The default neuron firing rates in Nengo are pretty high (200-400Hz), while biological neurons fire at a much lower frequency. This can be done by passing in a max_rates parameter to the initialization of the Ensemble or EnsembleArray:

from nengo import Uniform

# Example EnsembleArray with neurons firing from 10-50 Hz
motor = nengo.networks.EnsembleArray(n_neurons=N, n_ensembles=dim, max_rates=Uniform(10, 50),  label='Motor')

# Example ensemble with neurons firing from 5-20 Hz
ens = nengo.Ensemble(n_neurons=N, max_rates=Uniform(5, 20))

=======================
3. Changing the neuron type from the default LIF neuron to the adaptive LIF neuron (nengo.neurons.AdaptiveLIF) type may also help achieve the desired ISI distribution. Unlike the LIF neuron which has a fixed firing rate for a constant input, for a constant input, the adaptive LIF neuron starts a high firing rate, but the rate falls as it adapts the the constant input. The details of the adaptive LIF neuron can be found here.

To use the adaptive LIF neuron, specify the neuron_type parameter when creating your ensemble or ensemble array:

from nengo import Uniform

# Example EnsembleArray with ALIF neurons
motor = nengo.networks.EnsembleArray(n_neurons=N, n_ensembles=dim, neuron_type=nengo.AdaptiveLIF(...),  label='Motor')

# Example ensemble with ALIF neurons
ens = nengo.Ensemble(n_neurons=N, , neuron_type=nengo.AdaptiveLIF(...))

=========================

I hope these suggestions help with your query! :slight_smile:

1 Like
3

Thanks a lot @xchoo for the suggestions. I can at least see exponential distr for the inter spiking times now.

4

That’s great to hear! :smiley:
If you have any other questions, please don’t hesitate to post here.

5

I need to simulate for longer times like ~6500sec. At this point I exhaust all of memory 46 GB, so is there any way such that I just probe only the spiking times, instead of whole neuron spikes output which has lots of zeros.
Thanks

6

This is a long-standing feature request.

Fortunately there are some quick-and-dirty work-arounds you could take. Here is one that I quickly hacked together. You may wish to modify how the buffer data is structured according to what is convenient or efficient for your use-case (currently it stores a list of spike times for each neuron separately).

import nengo
import numpy as np


class SpikeTimeRecorder:
    """Records a list of spike times for each neuron.
    
    Intended for use with ``nengo==3.0.0``
    
    Source: https://forum.nengo.ai/t/spikes-analysis/1348/5
    """
    
    def __init__(self, n_neurons):
        self.n_neurons = n_neurons
        self.buffer = [[] for _ in range(n_neurons)]

    def __call__(self, t, x):
        if x.shape != (self.n_neurons,):
            raise RuntimeError(
                "Expected n_neurons=%d; got shape: %s"
                % (self.n_neurons, x.shape)
            )
        # Note: if multiple spikes happen in a single time-step
        # then they will be recorded only once!
        for i in np.where(x > 0)[0]:
            self.buffer[i].append(t)


def probe_spike_times(x):
    """Helper method for probing ``x.neurons`` times using a node."""
    obj = SpikeTimeRecorder(x.n_neurons)
    output = nengo.Node(obj, size_in=x.n_neurons)
    nengo.Connection(x.neurons, output, synapse=None)
    return obj.buffer

Example usage:

with nengo.Network(seed=0) as model:
    x = nengo.Ensemble(n_neurons=10, dimensions=1)

    # Get list of spike times for each neuron
    buffer = probe_spike_times(x)


with nengo.Simulator(model) as sim:
    sim.run(1)


for i, neuron_spike_times in enumerate(buffer):
    # neuron_spike_times is a list of spike times (in seconds)
    print("Neuron %d spiked %d times" % (i, len(neuron_spike_times)))

Example output:

Neuron 0 spiked 0 times
Neuron 1 spiked 0 times
Neuron 2 spiked 99 times
Neuron 3 spiked 0 times
Neuron 4 spiked 0 times
Neuron 5 spiked 99 times
Neuron 6 spiked 163 times
Neuron 7 spiked 0 times
Neuron 8 spiked 0 times
Neuron 9 spiked 10 times

Let me know if you have any questions about how to use this snippet. :slight_smile:

3 Likes
7

Thanks @arvoelke, this works like expected :slightly_smiling_face:

8

I had this same problem. For what it’s worth, here is what I have been using. It’s very inelegant compared to @arvoelke’s solution, but it handles multiple channels. It’s also pretty fast – sub millisecond per time step for compression of 1 million channels. That’s probably overkill unless you are using Nengo OCL or Nengo DL.

You can test it out with this notebook. The notebook will work up to 1M neurons, but just make sure to comment out plotting if you want to do that.

event_based.py (8.9 KB)
event-based-representation.ipynb (24.6 KB)