Spiking neural networks#

Spiking neural networks are nonlinear systems that respond to input by integrating signals over time. Put differently, they react to signals with some delay. As an example, the leaky integrate-and-fire (LIF) neuron model integrates some incoming current until a threshold is reached where it fires a signal (1). Until that point in time, the neuron was silent (0).

Imports#

Before we get started, we need to import some dependencies: numpy and matplotlib.

import numpy as np
import matplotlib.pyplot as plt
import torch
from urllib.request import urlretrieve

Leaky integrator neuron models#

One of the simplest neuron-like models is the leaky integrator. It takes some input current and integrates it over time, but with a leak term trending towards its resting state.

leaky integrator

Fig. 1 A leaky integrator that integrates current over time into its membrane potential. The left panel shows a constant input signal while the right shows a current oscillating between \([-0.1;0.1]A\). Note that the leaky integrator can produce negative membrane voltages.#

This type of neuron provides a continuous voltage, indicating the amount of current it builds up over time. The neuron voltage (\(v\)) evolves as follows:

\[\begin{split} \begin{align*} \dot{v} &= 1/\tau (i - v) \\ \end{align*} \end{split}\]

Task 0: Generate input data#

To begin with, generate at least 1000 consecutive steps of input data. To keep things simple, the data should simply be 1000 ones (1), similar to the left panel above.

input_data = ...

Task 1: Implement a leaky integrator#

Use the equation above to implement a leaky-integrator. To keep things simple, we will use Euler’s method to calculate the evolution of the equation. That is, the value of the next timestep is simply the value of the previous timestep plus the function evaluated at the next timestep

\[ y_{n + 1} = y_n + hf(t_n, y_n) \]

Where \(h\) controls the size of the step. To keep things simple, we’ll assume \(h\) is equivalent to the inverse of our time constant, \(1/\tau\).

The simplest way to implement this is by looping through the timesteps one by one. Note that you must somehow keep track of the current value of \(v\) in the equation above.

def leaky_integrator(input_data, tau):
    # Loop through the input one step at the time
    return ...


# Plot the output
# plt.plot(leaky_integrator(input_data, 1))

Leaky integrate-and-fire neuron models#

The leaky integrate-and-fire model (LIF) is a popular choice for a neuron models because it is relatively simple to simulate, while still exhibiting interesting spiking dynamics.

no-current

Fig. 2 Three examples of how the LIF neuron model responds to three different, but constant, input currents: 0.0, 0.1, and 0.3. At 0.3, we see that the neuron fires a series of spikes, followed by a membrane “reset”. Note that the neuron parameters are non-biological and that the memebrane voltage threshold is 1.#

The LIF model works like the leaky integrator, except when the neuron membrane voltage crosses the threshold \(v_{th}\). If that happens, the membrane voltage resets and a spike (1) is emitted. In all other instances, the neuron outputs (0’s). From a control perspective, that is a quite interesting property because it resembles a binary “gate” that is either open or closed.

The LIF dynamics are implemented as follows (note the similarity with the leaky integrator)

\[\begin{split} \begin{align*} \dot{v} &= 1/\tau (i - v) \\ \end{align*} \end{split}\]

with the jump equation

\[ \begin{align*} v &= (1-z) v + z \end{align*} \]

Where \(z\) represents a spiking tensors with values \(\in \{0, 1\}\).

Task 2: Implement a leaky integrate-and-fire neuron#

Use Euler’s method as you did for the leaky integrator above to implement the LIF dynamics.

def leaky_integrate_and_fire(input_data, tau, threshold):
    # Loop through the input one step at the time
    return output, membrane

# Plot the output
# output, membrane = leaky_integrate_and_fire(input_data)
# plt.plot(membrane)
# plt.scatter(range(len(input_data)), output)

Task 3: High-pass filter#

LIF neurons are essentially high-pass filters that filter out signals below and pass signals above a certain frequency.

Here is a signal at a certain frequency. What happens if you pass that through your LIF model? Can you the neuron spike?

min_x = 0
max_x = np.pi * 5
samples = 100
data_range = np.linspace(min_x, max_x, samples)
data_signal = (np.cos(data_range) + 1) / 2
plt.scatter(data_range, data_signal)
<matplotlib.collections.PathCollection at 0x7fe350313050>
../_images/b7f8e5e14e39ceafd9743d5fc3a5661b782c363e26f544e6d0f7406fae8f6255.png

Now, what if you changed the frequency of the data signal? You can do that by modifying the max_x to change the domain of the signal, for instance. Can you find a frequency that doesn’t spike? Then you have effectively created a band-pass filter!

Task 3.1: Band-pass filter (optional)#

A band-pass filter allows a specific range of frequencies to pass through. Can you combine two high-pass filters to create one single band-pass filter? You would have to create two high-pass filters, (1) one with a high cutoff and (2) one with a low cutoff, and combine them somehow.

Task 3.2: LIF neurons as band-pass filters (optional)#

We wrote that LIF neurons are technically high-pass filters because they filter away frequencies below a certain cut-off. They are actually also low-pass filters in some sense. Can you figure out why?

Task 4: Represent a spiking time-surface in Jupyter#

We are now going to look at events / spikes in time. Events are nothing but discrete values in time. But they represent important information, both in the spatial and temporal dimension. As an example, when sound hits your inner ear, it is decomposed and converted into spikes, roughly representing a discrete spectrogram.

Here’s an example of an audio recording of a digit from the Spiking Heidelberg Digit dataset. Every green dot represents a spike.

filename, _ = urlretrieve("https://github.com/ncskth/phd-course/raw/main/book/module1/sound_data.dat")
data = torch.load(filename)
plt.figure(figsize=(16, 4))
plt.gca().set_xlabel("Time")
plt.gca().set_ylabel("Frequency")
plt.imshow(data, interpolation="none")
<matplotlib.image.AxesImage at 0x7fe2e1d40990>
../_images/624d1354cf09648815f8f9662189eab2ae284dc8550d93b4ef8139bd084598d8.png

Note the vmax=0.1. What happens if you remove it? Why do we need it? Can you explain what that parameter is doing?

Hint: plotting sparse data can hide a tremendous amount of detail if we’re not careful.

What you see is a beautiful visualization of a signal in space and in time. In this kind of representation, time is tremendously important. One way to think about it is to view sound signals as space-time “surfaces”; we need to “see” the whole surface to make sense of it. Put differently, any downstream process that is sent your signal needs to somehow integrate time, otherwise it cannot capture the necessary information.

Time surfaces can be “stretched” and “compacted” in time by squeezing the signals closer together or prying them further apart. This is similar to sharpening and blurring effects in pictures. You’ll see how that can work now.

Task 4.1: Continuous time surface with leaky integrators#

Run the data above through your leaky integrator model and plot the output as a two-dimensional picture, as above. What does the plot represent? What happens when you increase/decrease the time constant \(\tau\)?

Task 4.2: Discrete time surface with leaky integrator-and-fire models#

Run the data above through you leaky integrate-and-fire model and plot the output as a two-dimensional picture, as above. What does the plot represent? What happens when you increase/decrease the time constant \(\tau\)?

Task 4.3: Use notebook sliders to visualize the difference (optional)#

To understand the effects of \(\tau\) for the above time surfaces it would be helpful to have some form of sliders to continuously change the value and visualize the effect. Luckily, this can be done with Jupyter Notebooks! Use the Jupyter Widgets to create a slider that changes the value of the time constant. Every time the time constant changes, you should regenerate and visualize the plot so we can observe the result.import ipywidgets as widgets

import ipywidgets as widgets

# ...
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[7], line 1
----> 1 import ipywidgets as widgets
      3 # ...

ModuleNotFoundError: No module named 'ipywidgets'