Import Packages

from elephant.spike_train_generation import StationaryPoissonProcess
import matplotlib.pyplot as plt
from neo.core import SpikeTrain
import numpy as np
import quantities as pq

Utility Functions

These are some helper functions we’ll use in the exercises. Please run the cell to get the functions:

from numpy.typing import NDArray

class utils:

    @staticmethod
    def generate_spiketrain(rate: float, t_stop: float) -> NDArray:
        neuron = StationaryPoissonProcess(rate=rate*pq.Hz, t_stop=t_stop*pq.s)
        st = neuron.generate_spiketrain()
        st_times = st.times
        return st_times
    
    @staticmethod
    def plot_poisson_distribution(mu: float) -> None:
        import numpy as np
        from scipy.stats import poisson, expon
        from scipy.interpolate import interp1d
        
        low = int(poisson.ppf(.001, mu))
        high = int(poisson.ppf(.999, mu))
        
        x = np.arange(low, high + 1, dtype=int)
        prob = poisson.pmf(k=x, mu=mu)
        
        plt.plot(x, prob, color="black")
        

In this first section, you will generate individual spike trains from a stationary Poisson process. This is a simple model where spikes occur randomly over time, but with a constant average firing rate. For example, a neuron firing at 5 Hz is expected to produce about 5 spikes per second on average, but that does not mean the spikes will be evenly spaced or that every simulation will produce exactly the same number of spikes. The exercises here are meant to help you connect the firing rate, the observation time, the number of spikes you expect, and the spike times you actually observe.

Code Reference

Exercises

Please fill in the blank cells with working code, according to each exercise below.

Note: “Example” exercises already have the solutions, and will help show how the code is used for the exercises that follow.

Example: Generate a spikeetrain from a cell spiking an average rate of 5 times per second (“Hz”) over a 1000-second period, and count the number of spikes.

• How many spikes would you expect to see?
• Is the measured amount close to your expectation?

st = utils.generate_spiketrain(rate=5, t_stop=1000)
len(st)

4972

Exercise: Generate a spiketrain from a cell spiking an average rate of 10 times per second (“Hz”) over a 5-second period, and count the number of spikes.

• How many spikes would you expect to see?
• Is the measured amount close to your expectation?

st = utils.generate_spiketrain(rate=10, t_stop=5)
len(st)

52

Exercise: Generate a spiketrain from a cell spiking an average rate of 7 times per second (“Hz”) over a 10-second period, then use the plt.eventplot() function to make a rasterplot and view each spike in the spiketrain over time.

• Are the spikes evenly-spaced in time?
• If you re-run the code, do the spikes always come at the same times?

Solution

st = utils.generate_spiketrain(rate=7, t_stop=5)
plt.eventplot(st)

Exercise: Generate a spiketrain from a cell spiking an average rate of 20 times per second (“Hz”) over a 100-second period, then use the plt.eventplot() function to make a rasterplot and view each spike in the spiketrain over time.

Now, it’s too difficult to see the spikes because of a problem called “over-plotting”! Fix the plot: add the following parameters to the plt.eventplot() function to make each event more visible:

• linewidths=0.5: Makes the lines narrower (so less overlap)
• colors='black': Makes the lines darker (so stronger contrast with background)
• alpha=0.5: Makes the lines translucent (so overlap increases contrast)

Solution

st = utils.generate_spiketrain(rate=5, t_stop=100)
plt.eventplot(st, colors='black', linewidths=.5, alpha=.5)

One spike train can be quite noisy, so it is hard to learn much from a single observation. In this section, we will generate many spike trains from the same underlying process and look at how the spike counts vary across repetitions. This is closer to the way we often reason about neural data: individual trials may differ, but repeated trials can reveal the statistical structure of the process that generated them. We will compare the simulated spike-count distribution to the distribution expected from a Poisson model, where the expected number of spikes is given by the firing rate multiplied by the recording duration.

Code Reference

Exercises

Example: Generate 3 spike trains, each with an average firing rate of 6 Hz and a measurement duration of 9 seconds. Then plot them as a raster plot.

sts = [utils.generate_spiketrain(rate=6, t_stop=9) for _ in range(3)]
plt.eventplot(sts, linelengths=.95, linewidths=.5, colors='black');

Exercise: Generate 10 spike trains, each with an average firing rate of 4 Hz and a measurement duration of 15 seconds. Then plot them as a raster plot.

Solution

sts = [utils.generate_spiketrain(rate=9, t_stop=15) for _ in range(10)]
plt.eventplot(sts, linelengths=.95, linewidths=.5, colors='black');

Example: For the following spike trains, use np.mean() to calculate the average number of spikes per spike train.

sts = [utils.generate_spiketrain(10, 10) for _ in range(100)]

counts = [len(st) for st in sts]
np.mean(counts)

np.float64(100.91)

Exercise: For the following spike trains, use np.std() to calculate the standard deviation of spikes observed between each spike train.

sts = [utils.generate_spiketrain(10, 10) for _ in range(100)]

counts = [len(st) for st in sts]
np.std(counts)

np.float64(9.897570408943803)

Exercise: For the following spike trains, use plt.hist() to visualize the distribution of spike counts measured in each spike train.

sts = [utils.generate_spiketrain(10, 10) for _ in range(100)]

Solution

counts = [len(st) for st in sts]
plt.hist(counts);

Example: What is the range of spike counts we would expect to see from a stationary Poisson-like neuron firing at an average rate of 4 Hz over a 3-second period? Use utils.plot_poisson_distribution() to plot the expected distribution.

• Run the code multiple times. Does this distribution change?

utils.plot_poisson_distribution(mu=12)  # 4Hz x 3secs

Exercise: What is the range of spike counts we would expect to see from a stationary Poisson-like neuron firing at an average rate of 8 Hz over a 10-second period? Use utils.plot_poisson_distribution() to plot the expected distribution.

• Run the code multiple times. Does this distribution change?

utils.plot_poisson_distribution(mu=80) 

Exercise: Does the distribution of observed spike counts match what would be expected from a Poisson model? For the spike counts calculated below:

• Plot the proportional distribution of spike counts with plt.hist(density=True).
• Plot the expected distribution with utils.plot_poisson_distribution(). This function needs a parameter called mu; this is the expected number of spikes (i.e. average spike rate, multiplied by time duration)
• Does the shape of the expected distrubtion roughly follow the measured spike distribution?

sts = [utils.generate_spiketrain(4, 2) for _ in range(500)]
counts = [len(st) for st in sts]

Solution

plt.hist(counts, bins=25, density=True);
utils.plot_poisson_distribution(mu = 8)

This section introduces Neo, a package that provides standardized data objects for electrophysiology. So far, we have mostly treated spike times as simple arrays or lists. That is fine for small examples, but real analysis often needs more context: when the recording started, when it stopped, what units the times are measured in, and how the data should be passed between tools. A Neo SpikeTrain object stores spike times together with this kind of metadata. In this section, you will create SpikeTrain objects manually, inspect their contents, plot their spike times, and then use Elephant to generate Neo-compatible spike trains automatically.

Code Reference

Exercises

from neo.core import SpikeTrain
import quantities as pq

Example: Construct a Neo SpikeTrain object that describes the following spike train:

• spike times were observed at 0.2, 0.5, 0.9, and 1.3 seconds
• measurements started at 0 seconds and ended at 2 seconds

st = SpikeTrain(
    times=[.2, .5, .9, 1.3] * pq.s, 
    t_start=0 * pq.s, 
    t_stop=2 * pq.s
)
st

SpikeTrain containing 4 spikes; units s; datatype float64  time: 0.0 s to 2.0 s

Exercise: Construct a Neo SpikeTrain object that describes the following spike train:

• spike times were observed at 0.2, 0.5, 0.9, and 1.3 seconds
• measurements started at 0 seconds and ended at 2 seconds

st = SpikeTrain(
    times=[.2, .5, .9, 1.3] * pq.s, 
    t_start=0 * pq.s, 
    t_stop=2 * pq.s
)
st

SpikeTrain containing 4 spikes; units s; datatype float64  time: 0.0 s to 2.0 s

Example: Take the following SpikeTrain object and display the starting time (e.g. t_start) contained inside it.

st = SpikeTrain(times=[.2, .5, .9, 1.3] * pq.s, t_start=0 * pq.s, t_stop=2 * pq.s)
st

st.t_start

array(0.) * s

Exercise: Take the following SpikeTrain object and display the ending time (e.g. t_stop) contained inside it.

st = SpikeTrain(times=[.2, .5, .9, 1.3] * pq.s, t_start=0 * pq.s, t_stop=2 * pq.s)
st

st.t_stop

array(2.) * s

Exercise: Take the following SpikeTrain object and display the spike times (e.g. times) contained inside it.

st = SpikeTrain(times=[.2, .5, .9, 1.3] * pq.s, t_start=0 * pq.s, t_stop=2 * pq.s)
st

st.times

array([0.2, 0.5, 0.9, 1.3]) * s

Exercise: Using plt.eventplot(), make a raster plot of the spike times from the following SpikeTrain object:

st = SpikeTrain(times=[.2, .28, .5, .9, .97, 1.3, 1.9] * pq.s, t_start=0 * pq.s, t_stop=2 * pq.s)
st

SpikeTrain containing 7 spikes; units s; datatype float64  time: 0.0 s to 2.0 s

Solution

plt.eventplot(st.times);

The elephant package contains many useful utility classes and functions that work with neo objects.

Example: Use the StationaryPoissonProcess class to generate a spike train with an average firing rate of 5 Hz and ending at 3 seconds of measurement time. Display a raster plot of the generated spike times.

from elephant.spike_train_generation import StationaryPoissonProcess

process = StationaryPoissonProcess(rate=5*pq.Hz, t_stop=3 * pq.s)
st = process.generate_spiketrain()
plt.eventplot(st.times);

Exercise: Use the StationaryPoissonProcess class to generate a spike train with an average firing rate of 8 Hz and ending at 4 seconds of measurement time. Display a raster plot of the generated spike times.

Solution

process = StationaryPoissonProcess(rate=8*pq.Hz, t_stop=4 * pq.s)
st = process.generate_spiketrain()
plt.eventplot(st.times);

Exercise: Instead of generating a single spike train with proces.generate_spiketrain(), use process.generate_n_spiketrains(as_array=True) to generate 50 spike trains’ times, then display them as a raster plot.

Solution

process = StationaryPoissonProcess(rate=5*pq.Hz, t_stop=3 * pq.s)
st_times = process.generate_n_spiketrains(50, as_array=True)
plt.eventplot(st_times);

Neuroscience data almost always come with physical units: time might be measured in seconds or milliseconds, firing rate in Hertz, and voltage in volts or millivolts. It is easy to make mistakes when these units are tracked only in variable names or comments. The quantities package helps by attaching units directly to values. This means Python can convert compatible units automatically and raise an error when you try to combine quantities that do not make physical sense. In this section, you will practice creating unit-aware variables, inspecting their units and magnitudes, and using them to compute quantities such as expected spike count.

Code Reference

Exercises

Example: Create a variable called firing_rate and assign it the value 5 * pq.Hz.

firing_rate = 5 * pq.Hz

Exercise: Create a variable called t_start and assign it the value 100*pq.ms.

t_start = 100 * pq.ms

Exercise: Create a variable called t_stop and assign it the value 3 * pq.s

t_stop = 3 * pq.s

Exercise: Compute the difference between t_stop and t_start and assign the result to a new variable duration. Then, print that variable.

duration = t_stop - t_start
duration

array(2.9) * s

Exercise: Get the .magnitude and the .units of duration.

duration.units, duration.magnitude

(array(1.) * s, array(2.9))

Exercise: Try adding the variables firing_rate and duration. What error message do you observe?

firing_rate+duration

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
File ~/repositories/new-learning-platform/.pixi/envs/lfp/lib/python3.12/site-packages/quantities/quantity.py:234, in Quantity.rescale(self, units, dtype)
    233 try:
--> 234     cf = get_conversion_factor(from_u, to_u)
    235 except AssertionError:

File ~/repositories/new-learning-platform/.pixi/envs/lfp/lib/python3.12/site-packages/quantities/quantity.py:56, in get_conversion_factor(from_u, to_u)
     55 to_u = to_u._reference
---> 56 assert from_u.dimensionality == to_u.dimensionality
     57 return from_u.magnitude / to_u.magnitude

AssertionError: 

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
Cell In[8], line 1
----> 1 firing_rate+duration

File ~/repositories/new-learning-platform/.pixi/envs/lfp/lib/python3.12/site-packages/quantities/quantity.py:66, in scale_other_units.<locals>.g(self, other, *args)
     64     other = other.view(type=Quantity)
     65 if other._dimensionality != self._dimensionality:
---> 66     other = other.rescale(self.units, dtype=np.result_type(self.dtype, other.dtype))
     67 return f(self, other, *args)

File ~/repositories/new-learning-platform/.pixi/envs/lfp/lib/python3.12/site-packages/quantities/quantity.py:236, in Quantity.rescale(self, units, dtype)
    234     cf = get_conversion_factor(from_u, to_u)
    235 except AssertionError:
--> 236     raise ValueError(
    237         'Unable to convert between units of "%s" and "%s"'
    238         %(from_u._dimensionality, to_u._dimensionality)
    239     )
    240 if np.dtype(dtype).kind in 'fc':
    241     cf = np.array(cf, dtype=dtype)

ValueError: Unable to convert between units of "s" and "Hz"

Exercise: Multiply firing_rate by duration and assign the result to a new variable n_spikes. Print that variable.

n_spikes = firing_rate * duration
n_spikes

array(14.5) * s*Hz