Spiking threshold as a parameter


#1

Hello,

For a project I need to create LIF ensembles with various spiking threshold. As far as I know, it is not possible to choose this parameter in Nengo. Therefore, I changed the source code to provide the possibility of manually setting the threshold when creating a LIF:

`class LIF(LIFRate):
"""Spiking version of the leaky integrate-and-fire (LIF) neuron model.

Parameters
----------
tau_rc : float
    Membrane RC time constant, in seconds. Affects how quickly the membrane
    voltage decays to zero in the absence of input (larger = slower decay).
tau_ref : float
    Absolute refractory period, in seconds. This is how long the
    membrane voltage is held at zero after a spike.
min_voltage : float
    Minimum value for the membrane voltage. If ``-np.inf``, the voltage
    is never clipped.
amplitude : float
    Scaling factor on the neuron output. Corresponds to the relative
    amplitude of the output spikes of the neuron.
"""

probeable = ('spikes', 'voltage', 'refractory_time')

min_voltage = NumberParam('min_voltage', high=0)

def __init__(self, tau_rc=0.02, tau_ref=0.002, min_voltage=0, amplitude=1, spiking_threshold=1):
    super(LIF, self).__init__(
        tau_rc=tau_rc, tau_ref=tau_ref, amplitude=amplitude)
    self.min_voltage = min_voltage
    self.spiking_threshold=spiking_threshold

def step_math(self, dt, J, spiked, voltage, refractory_time):
    # reduce all refractory times by dt
    refractory_time -= dt

    # compute effective dt for each neuron, based on remaining time.
    # note that refractory times that have completed midway into this
    # timestep will be given a partial timestep, and moreover these will
    # be subtracted to zero at the next timestep (or reset by a spike)
    delta_t = (dt - refractory_time).clip(0, dt)

    # update voltage using discretized lowpass filter
    # since v(t) = v(0) + (J - v(0))*(1 - exp(-t/tau)) assuming
    # J is constant over the interval [t, t + dt)
    voltage -= (J - voltage) * np.expm1(-delta_t / self.tau_rc)

    # determine which neurons spiked (set them to 1/dt, else 0)
    spiked_mask = voltage > self.spiking_threshold
    spiked[:] = spiked_mask * (self.amplitude / dt)

    # set v(0) = 1 and solve for t to compute the spike time
    t_spike = dt + self.tau_rc * np.log1p(
        -(voltage[spiked_mask] - self.spiking_threshold) / (J[spiked_mask] - self.spiking_threshold))

    # set spiked voltages to zero, refractory times to tau_ref, and
    # rectify negative voltages to a floor of min_voltage
    voltage[voltage < self.min_voltage] = self.min_voltage
    voltage[spiked_mask] = 0
    refractory_time[spiked_mask] = self.tau_ref + t_spike` 

I basically replaced every 1 by self.spiking_threshold in step_math, although I’m not sure about what exactly is computed.

Although it successfully sets the threshold, this solution seems to affect the way the ensemble represents values. With the simple code below, the decoded value is inversely scaled by the given spiking_threshold:

`with model:
stim = nengo.Node([0])
a = nengo.Ensemble(n_neurons=50, dimensions=1, 
    neuron_type=nengo.neurons.LIF(spiking_threshold=2.))
nengo.Connection(stim, a)`

Is there a clean way to choose the spiking threshold? Can you give me some hints to implement it otherwise?


#2

Hi @Vaurien96. Thanks for the question and for taking the time to try to figure out these details and get this code working. I think the most important thing to note is that, for the case of the LIF model, the subthreshold dynamics—what happens to the voltage while it’s operating within the open interval (0, threshold)—are completely linear with respect to the current. Within this operating range, the voltage vector is governed by a leaky integrator (a first-order lowpass filter) – a linear dynamical system applied to the input current.

This observation is important because it tells us that scaling the threshold is mathematically equivalent to inversely scaling the input current by the exact same factor. Another way to say this is that: driving the neuron with an input current of $J(t)$ until it reaches $V(t)=1$ is the same as driving an identical neuron with an input current $R \cdot J(t)$ until it reaches a threshold of $V(t)=R$, by linearity. In other words, you can effectively modify the threshold to $R$ (i.e., get the same mathematical outcome) while keeping the actual threshold fixed at $1$ by scaling the input current by $1 / R$.

Since the input current in the NEF is defined by: $$J = \alpha \left \langle {\bf e}\text{,}\, {\bf x}(t) \right \rangle + \beta$$ and we want to scale $J$ by $1 / R$, we can do this by scaling $\alpha$ and $\beta$ by $1 / R$, as in: $$\frac{J}{R} = \left( \frac{\alpha}{R} \right) \left \langle {\bf e}\text{,}\, {\bf x}(t) \right \rangle + \left( \frac{\beta}{R} \right) \text{.}$$

In Nengo, $\alpha$ and $\beta$ are referred to as the gain and bias, respectively. These parameters are by default derived from the max_rates and intercepts of the desired tuning curves. This makes the important point that modifying the threshold is mathematically equivalent to modifying your max_rates and intercepts, i.e., your desired tuning curves. These are all completely configurable in Nengo (e.g., see tuning_curves and nef_summary).

Here is some example code that uses this trick to obtain an effective threshold of $R = 2$, and compares the voltage trace between two otherwise-identical neurons:

import numpy as np
import matplotlib.pyplot as plt

import nengo

threshold = 2.0

with nengo.Network(seed=0) as model1:
    stim1 = nengo.Node(1)
    
    x1 = nengo.Ensemble(n_neurons=1, dimensions=1, encoders=[[1]], max_rates=[100])

    nengo.Connection(stim1, x1)
    p1 = nengo.Probe(x1.neurons, 'voltage')
    
with nengo.Simulator(model1, dt=1e-5) as sim1:
    sim1.run(0.05)
    
gain = sim1.data[x1].gain / threshold
bias = sim1.data[x1].bias / threshold

with nengo.Network(seed=0) as model2:
    stim2 = nengo.Node(1)
    
    x2 = nengo.Ensemble(n_neurons=1, dimensions=1, encoders=[[1]],
                        gain=gain, bias=bias)

    nengo.Connection(stim2, x2)
    p2 = nengo.Probe(x2.neurons, 'voltage')
    
with nengo.Simulator(model2, dt=1e-5) as sim2:
    sim2.run(0.05)
    
plt.figure(figsize=(14, 5))
plt.plot(sim1.trange(), sim1.data[p1],
         label="threshold=1")
plt.plot(sim2.trange(), sim2.data[p2]*threshold,
         label="threshold=%s" % threshold)
plt.xlabel("Time (s)")
plt.ylabel("Voltage")
plt.legend()
plt.show()

To answer your question a bit more directly, though:

The reason you were seeing this behaviour is because: you not only need to change step_math, but you also need to modify the rates method inherited from LIFRate:

        j = J - 1
        output[:] = 0  # faster than output[j <= 0] = 0
        output[j > 0] = self.amplitude / (
            self.tau_ref + self.tau_rc * np.log1p(1. / j[j > 0]))

The decoding weights are learned from the tuning curves produces by rates, and so these must be consistent with step_math. You can use the observation from above to make the appropriate change by rescaling the input current (J) according to the desired threshold.

To summarize, you can either effectively change the LIF threshold by rescaling the input current, or by using your code and making the rates method consistent with your change. But, in either case, this is mathematically equivalent to rescaling the gain and bias parameters of your tuning curves, which you could also achieve by playing with max_rates and intercepts. If your desired tuning curves remain the same (i.e., if you also updated the gain_bias and max_rates_intercepts methods to be consistent with your new rates method), then this change will have no apparent effect, apart from rescaling the voltage trace.

Please feel free to ask more details about the code, math, or how to use this to do something!


#3

Thank you @arvoelke for your detailed answer! Scaling the inputs is indeed a simple yet efficient solution.

EDIT: I’ll create a new thread for my new questions.