Improving Accuracy

Improving Accuracy#

Our implementation of LIF neurons is good enough for many applications but there are lots of ways in which our implementation is imprecise. Namely:

If the refractory period ends in the middle of a step (if 0 < self.refractory_time < t_step) then we pretend that the refractory time ended at the beginning of the step. In reality, the voltage shouldn’t have updated for the part of the step where we were in refractory time.
If a neuron fires in the middle of a timestep—between \((t - T_{step})\) and \(t\)—then we are effectively “rounding” the firing time to the nearest timestep, since we assign self.refractory_time = self.tau_ref at the end of the step (whereas some of the refractory time should have passed before the step end).
We are using Euler’s method as a first-order approximation of self.v, which has some error (which is magnified for larger values of \(T_{step}\)).

We can actually “correct” all of these inaccuracies at little computational cost. In fact, by adding accuracy, we can reduce computational cost by running our simulation in larger time steps. You can look at the optional descriptions below to see how we can rectify all of these:

1. Accurately Setting the Refractory Time

Accurately Setting the Refractory Time (Optional)

Recall that:

\[v(t + \Delta{}t) \approx v(t)(1 - \frac{\Delta{}t}{\tau_{rc}}) + \frac{\Delta{}t}{\tau_{rc}}I(t)\]

Now, consider the image below, which represents a neuron crossing the threshold voltage (\(v_{th}\) in the diagram) at some time between \((t - \Delta{}t)\) and \(t\). We will call this time \(t_{th}\).

As this image shows, our neuron should exit the refractory period at time \((t - \Delta{}t + t_{th} + \tau_{ref})\). However, it is currently incorrectly exiting the refractory period at time \((t + \tau_{ref})\). Let’s rectify this. First, we need to figure out \(t_{th}\).

We know that \(v(t - \Delta{}t + t_{th}) = v_{th}\) . Using Euler’s above, we get:

\[v(t - \Delta{}t + t_{th}) = v_{th} = (1 - \frac{t_{th}}{\tau_{rc}})v(t-\Delta{}t) + \frac{I(t)t_{th}}{\tau_{rc}}\]

Moving things around, we get:

\[t_{th} = \tau_{rc}\frac{v_{th} - v(t - \Delta{}t)}{I(t) - v(t-\Delta{}t)}\]

It can be a little more convenient to rewrite this using \(v(t)\) (the potential now) rather than \(v(t - \Delta{}t)\) (the potential in the last timestep). We know that they are related by \(v(t) = v(t-\Delta{}t)(1 - \frac{\Delta{}t}{\tau_{rc}}) + \frac{I(t)\Delta{}t}{\tau_{rc}}\). Or converseley, \(v(t - \Delta{}t) =\frac{\tau_{rc}v(t) - I(t)\Delta{}t}{\tau_{rc}-\Delta{}t}\). Doing the math, you will find that we can solve for \(t_{th}\):

\[t_{th} = \frac{\tau_{rc}(v_{th} - v(t)) + \Delta{}t(I(t) - v_{th})}{I(t) - v(t)}\]

Now, rather than assigning self.refractory_time = self.tau_ref, we can be more precise and take into account the actual spike time:

spike_time = (self.tau_rc * (self.v_th - self.v) + t_step * (I - self.v_th)) / (I - self.v)
self.refractory_time = self.tau_ref + spike_time # Start refractory countdown

Note that we need to do this before setting self.v = 0. Let’s see what this code actually looks like:

class FirstOrderLIF1: # First Order Leaky Integrate and Fire
    def __init__(self, tau_rc=0.2, tau_ref=0.002, v_init=0, v_th=1): # Default values for tau_rc and v_init
        self.tau_rc  = tau_rc  # Set instance variables
        self.tau_ref = tau_ref
        self.v       = v_init
        self.v_th    = v_th

        self.output = 0
        self.refractory_time = 0
        

    def step(self, I, t_step): # Advance one time step (input I and time step size t_step)
        self.refractory_time -= t_step # Subtract the amount of time that passed from our refractory time

        if self.refractory_time < 0: # If we aren't in our refractory period, update the voltage
            self.v = self.v * (1 - t_step / self.tau_rc) + I * t_step / self.tau_rc

        if self.v >= self.v_th:      # Voltage is above the threshold
            # vvv CHANGED vvv #
            spike_time = (self.tau_rc * (self.v_th - self.v) + t_step * (I - self.v_th)) / (I - self.v)
            self.refractory_time = self.tau_ref + spike_time # Start refractory countdown
            # ^^^ CHANGED ^^^ #
        
            self.output = 1 / t_step                         # Fire
            self.v = 0                                       # Reset potential
        else:
            self.output = 0          # Don't fire

        return self.output

2. Accurately Updating Voltage

Accurately Updating Voltage (Optional)

If the refractory period ends in the middle of a timestep, our voltage should only be updated for the part of the timestep that was not in the refractory period. So we need to update our step() function again so that we take this into account. This means we can no longer use t_step as the change in time (since the actual time where I was applied might be smaller).

Instead, we’re going to create a variable delta_t that represents the actual amount of time that our neuron’s value updated. If 0 < self.refractory_time < t_step, we can assign delta_t to t_step - self.refractory_time to represent the actual amount of time that I was applied for:

delta_t = t_step - self.refractory_time # How much time passed where we can actually accept input?

However, if this.refractory_time > t_step then delta_t would be negative. We don’t want this so we add:

if delta_t < 0:      delta_t =  0       # Make sure it's not negative

The way that we defined self.refractory_time, it can be negative when our neuron finishes its refractory period. In that case, we want delta_t = t_step. If self.refractory_time is negative then we would end up with delta_t > t_step so we can add:

if delta_t > t_step: delta_t = t_step   # Make sure it's not greater than our time step

Then, we need to replace our use of t_step in our equations to instead use delta_t:

self.v = self.v * (1 - delta_t / self.tau_rc) + I * delta_t / self.tau_rc # CHANGED: Note that we use delta_t rather than t_step

class FirstOrderLIF2: # First Order Leaky Integrate and Fire
    def __init__(self, tau_rc=0.2, tau_ref=0.002, v_init=0, v_th=1): # Default values for tau_rc and v_init
        self.tau_rc  = tau_rc  # Set instance variables
        self.tau_ref = tau_ref
        self.v       = v_init
        self.v_th    = v_th

        self.output = 0
        self.refractory_time = 0
        

    def step(self, I, t_step): # Advance one time step (input I and time step size t_step)
        self.refractory_time -= t_step # Subtract the amount of time that passed from our refractory time

        # vvv CHANGED vvv
        delta_t = t_step - self.refractory_time # How much time passed where we can actually accept input?
        if delta_t < 0:      delta_t =  0       # Make sure it's not negative
        if delta_t > t_step: delta_t = t_step   # Make sure it's not greater than our time step
        # ^^^ CHANGED ^^^

        self.v = self.v * (1 - delta_t / self.tau_rc) + I * delta_t / self.tau_rc # CHANGED: Note that we use delta_t rather than t_step

        if self.v >= self.v_th:      # Voltage is above the threshold
            spike_time = (self.tau_rc * (self.v_th - self.v) + t_step * (I - self.v_th)) / (I - self.v)
            self.refractory_time = self.tau_ref + spike_time # Start refractory countdown
        
            self.output = 1 / t_step                         # Fire
            self.v = 0                                       # Reset potential
        else:
            self.output = 0          # Don't fire

        return self.output

After these improvements, we end up with the following code:

import numpy as np

class AnalyticalLIF: # First Order Leaky Integrate and Fire
    def __init__(self, tau_rc=0.2, tau_ref=0.002, v_init=0, v_th=1): # Default values for tau_rc and v_init
        self.tau_rc  = tau_rc  # Set instance variables
        self.tau_ref = tau_ref
        self.v       = v_init
        self.v_th    = v_th

        self.output = 0
        self.refractory_time = 0
        

    def step(self, I, t_step): # Advance one time step (input I and time step size t_step)
        self.refractory_time -= t_step # Subtract the amount of time that passed from our refractory time

        delta_t = t_step - self.refractory_time # How much time passed where we can actually accept input?
        if delta_t < 0:      delta_t =  0       # Make sure it's not negative
        if delta_t > t_step: delta_t = t_step   # Make sure it's not greater than our time step
        # ^---CHANGED to get a more accurate estimate how how much time was not part of the refractory time

        self.v = I + (self.v - I) * np.exp(-delta_t / self.tau_rc) # <-CHANGED to more accurately update the voltage

        if self.v >= self.v_th:      # Voltage is above the threshold
            spike_time = delta_t + self.tau_rc * np.log((self.v - I) / (self.v_th - I)) #<-CHANGED to get a more precise spike time
            
            self.refractory_time = self.tau_ref + spike_time # Start refractory countdown
        
            self.output = 1 / t_step                         # Fire
            self.v = 0                                       # Reset potential
        else:
            self.output = 0          # Don't fire

        return self.output

Measuring the Benefit#

So far, we’ve made some changes to get a more accurate model…but is there really a benefit? In short, that depends. Let’s measure. We’re going to take a few random time step options (defined using TIME_STEP_SIZES below):

…for each of these time steps, we are going to evaluate a rough idea of “accuracy” by measuring the number of times our neuron spikes given a few inputs. We will measure “ground truth” using a very small time step size, count how many times our neuron fires in the “ground truth” case, and compare that to how many times it fires for larger time steps. We will measure the error as the difference in the number of spikes (divided by the total number of spikes to scale it accurately)

We will do this for:

The first order approximation (in orange below), which was our original code
The first order approximation (in green below) that allows for refractory times that are fractions of the time step
Our analytical solution (in blue below), which is the AnalyticalLIF class defined above

Show code cell source Hide code cell source

import math
import matplotlib.pyplot as plt
import numpy as np

class FirstOrderLIF: # First Order Leaky Integrate and Fire
    def __init__(self, tau_rc=0.2, tau_ref=0.002, v_init=0, v_th=1): # Default values for tau_rc and v_init
        self.tau_rc  = tau_rc  # Set instance variables
        self.tau_ref = tau_ref
        self.v       = v_init
        self.v_th    = v_th

        self.output = 0
        self.refractory_time = 0
        

    def step(self, I, t_step): # Advance one time step (input I and time step size t_step)
        self.refractory_time -= t_step # Subtract the amount of time that passed from our refractory time

        if self.refractory_time < 0: # If we aren't in our refractory period, update the voltage
            self.v = self.v * (1 - t_step / self.tau_rc) + I * t_step / self.tau_rc

        if self.v >= self.v_th:      # Voltage is above the threshold
            self.refractory_time = self.tau_ref # Start refractory countdown
            self.output = 1 / t_step            # Fire
            self.v = 0                          # Reset potential
        else:
            self.output = 0          # Don't fire

        return self.output

class FirstOrderLIF2: # First Order Leaky Integrate and Fire
    def __init__(self, tau_rc=0.2, tau_ref=0.002, v_init=0, v_th=1): # Default values for tau_rc and v_init
        self.tau_rc  = tau_rc  # Set instance variables
        self.tau_ref = tau_ref
        self.v       = v_init
        self.v_th    = v_th

        self.output = 0
        self.refractory_time = 0
        

    def step(self, I, t_step): # Advance one time step (input I and time step size t_step)
        self.refractory_time -= t_step # Subtract the amount of time that passed from our refractory time

        # vvv CHANGED vvv
        delta_t = t_step - self.refractory_time # How much time passed where we can actually accept input?
        if delta_t < 0:      delta_t =  0       # Make sure it's not negative
        if delta_t > t_step: delta_t = t_step   # Make sure it's not greater than our time step
        # ^^^ CHANGED ^^^

        self.v = self.v * (1 - delta_t / self.tau_rc) + I * delta_t / self.tau_rc # CHANGED: Note that we use delta_t rather than t_step

        if self.v >= self.v_th:      # Voltage is above the threshold
            spike_time = (self.tau_rc * (self.v_th - self.v) + t_step * (I - self.v_th)) / (I - self.v)
            self.refractory_time = self.tau_ref + spike_time # Start refractory countdown
        
            self.output = 1 / t_step                         # Fire
            self.v = 0                                       # Reset potential
        else:
            self.output = 0          # Don't fire

        return self.output

NUM_SECONDS = 2
TIME_STEP_SIZES = []

currSize = 0.0001
while currSize < 0.002:
    TIME_STEP_SIZES.append(currSize)
    currSize *= 1.2

def getSpikeCount(neuron, I, numSeconds, timeStepSize):
    numSteps = int(numSeconds / timeStepSize)
    count = 0
    for i in range(numSteps):
        output = neuron.step(I(i * timeStepSize), timeStepSize)
        if output > 0:
            count += 1
    return count

neurons = {
    'analytical': lambda: AnalyticalLIF(),
    'firstOrder': lambda: FirstOrderLIF(),
    'firstOrderFractional': lambda: FirstOrderLIF2()
}

avgErrors = { k: [] for k in neurons.keys() }
errors = { k: {} for k in neurons.keys() }

for timeStepSize in TIME_STEP_SIZES:
    for k in errors:
        errors[k][timeStepSize] = []

for i in [1.2, 2, 4, 8, 16, 32, 64, 128]:
    I = lambda t: i
    groundTruthCount = getSpikeCount(AnalyticalLIF(), I, NUM_SECONDS, 1e-6)

    for neuronType in neurons:
        for timeStepSize in TIME_STEP_SIZES:
            neuron = neurons[neuronType]()
            err = abs(groundTruthCount - getSpikeCount(neuron, I, NUM_SECONDS, timeStepSize))
            errors[neuronType][timeStepSize].append(err)
            
        
for neuronType in neurons:
    for timeStepSize in TIME_STEP_SIZES:
        avgErrors[neuronType].append(np.mean(errors[neuronType][timeStepSize]))

# Plotting the error vs time step size
plt.figure()

for neuronType in neurons:
    plt.plot(TIME_STEP_SIZES, avgErrors[neuronType], label=neuronType)

plt.xscale('log')
plt.xlabel('Time Step Size')
plt.ylabel('Average Error')
plt.title('Error vs Time Step Size')
plt.legend()
plt.grid(True)
plt.show()

../_images/0db92fcb3e660083a4a9b1dbdcdb987ea4d703b06780ccae6b6f3af984a3e7a0.png

In the graph below, a lower error rate is better.

As we can see, our analytical solution (blue line) performs better than our original solution (orange line). This tells us that with the analytical model of LIFs, we might be able to use slightly larger time steps while maintaining accuracy.