Irregular spiking threshold when using `aeif_cond` models

[francesshei]

Adaptive exponential integrate-and-fire neurons seem to have an irregular spiking threshold with a certain set of parameters and a relatively strong input current.

I am trying to port a spiking neural network and its corresponding single neuron models implemented using the Brian simulator to NEST, but I am not able to reproduce the same firing patterns.

Minimal reproducible example:

import nest
from matplotlib import pyplot as plt
import numpy as np

I_AMPLITUDE  = 5000.0

# Create the AdEX model
neuron = nest.Create("aeif_cond_exp")

dc_current = nest.Create("dc_generator")
dc_current.set({"start": 1200.0, "stop": 1600.0, "amplitude": I_AMPLITUDE})

multimeter = nest.Create("multimeter")
multimeter.set(record_from=["V_m"])

w_multimeter = nest.Create("multimeter")
w_multimeter.set(record_from=["w"])

spikerecorder = nest.Create("spike_recorder")

nest.Connect(multimeter, neuron)
nest.Connect(w_multimeter, neuron)
nest.Connect(neuron, spikerecorder)
nest.Connect(dc_current, neuron)

# Run the simulation for 1000.0 ms (1 sec)
nest.Simulate(2500.0)

# Access the multimeter quantities (dictionary)
dmm = multimeter.get()
# Specifically, retrieve the voltage values and the timestamps of recording
Vms = dmm["events"]["V_m"]
ts = dmm["events"]["times"]

# Retrieve the adapting variable w changes
wmm = w_multimeter.get()
ws = wmm["events"]["w"]

# Plot the result
fig = plt.figure()

ax1 = fig.add_subplot(211)
ax1.plot(ts, Vms, "#1e2b61")
spike_times = spikerecorder.get()["events"]["times"]
ax1.scatter(spike_times, np.array([-30] * len(spike_times)), s=3)
ax1.set_xlim(1000, 2500)
ax1.set_ylim(-100, 0)

ax2 = fig.add_subplot(212)
ax2.plot(ts, ws, "#1e2b61")
ax2.set_xlim(1000, 2500)
plt.show()

The parameters used are as follows:

{
    'C_m': 1000.0,
    'g_L': 50.0,
    'E_L': -60.0,
    'Delta_T': 2.5,
    'V_th': -50.0,
    'V_reset': -60.0, 
    't_ref': 2.5,
    'I_e': 0.0,
    'tau_w': 600.0,
    'a': 400.0,
    'b': 20.0 
}

Simulating the same neuron in Brian and NEST yields different results: Brain (expected behaviour) [note that the V threshold is equivalent, but the potential at each spike time is plotted at +20 mV]

Nest

Not only the neuron does not seem to describe the same level of spike frequency adaptation, but the spike times (reported as blue dots superimposed over the voltage trace) appear to follow an irregular threshold, or the reset condition is irregularly computed after surpassing the threshold.

This seems to compromise the ability of the aeif_models to capture the same spiking dynamics that can be described with Brian. Indeed, I wasn't able to reproduce the same figure reported in the reference paper ( doi:10.1371/journal.pone.0161934 , fig. 1B)

Reference (NOTE: the legend might contain a typo, the stimulating current is 1-5 nA, as confirmed by running simulations with Brain able to replicate the figure)

NEST implementation

As it can be seen from the picture, the ISI are too regular in NEST and not adapting their frequency as much as those in Brian, resulting in a flatter histogram.

• OS: MacOS 13.1, Apple M1 Pro
• Shell: zsh
• Python-Version: Python 3.10.6
• NEST-Version: 3.3

[clinssen]

Hi, thanks for writing in!

It looks like, given the large current injection, the neuron is firing at a very high rate (around 400 Hz). The default resolution in NEST is not adequate to capture this fast behaviour. Please change the resolution on both the kernel

nest.resolution = 0.01    # [ms]

as well as on the multimeters:

multimeter.interval = nest.resolution
w_multimeter.interval = nest.resolution

This should hopefully allow a more close match between Brian and NEST. Note that there might still be small numerical discrepancies due to use of different numerical solvers, the exact time that the threshold condition is checked, etc.

Please let us know if you run into any further trouble, or otherwise please close this issue if everything is resolved. Cheers!

[francesshei]

Thanks @clinssen for the answer!

The original resolution value was chosen to match the one used in Brian. I changed both the kernel and multimeters resolutions and indeed the voltage trace looks better now (although the reset condition seems to be still irregular):

The ISI histogram still looks flatter though, I assumed it because of some intrinsic solver differences? The default solver in Brian is the Euler solver. Is there a way to force NEST to use a similar implementation?

I also have another question: I am attempting to replicate the two neurons network proposed in the study, which has synaptic weights values in the order of uS. Which unit of measurement do connection weights have in NEST?

[clinssen]

If you would like to try with the forward Euler solver, I have added some support to NESTML for this, please see https://github.com/nest/nestml/pull/940. Please note that the threshold condition check could still be different, so you might still see small numerical differences.

The default synaptic weight unit in NEST is nS or pA, depending on which model you are using.

[heplesser]

@francesshei The reason for the irregular looking spike threshold in NEST is that the trace from NEST shows the last sub-threshold voltage value before a spike. Since the membrane potential is driven by an exponentially growing current at the aeif neuron approaches threshold, a rather wide range of membrane potential values are possible at the last time step before a spike. So the variation you observe in the plot is expected behavior.

The difference between NEST and Brian is most likely due to the fact that Brain uses a fixed time-step forward Euler integration scheme for its aeif examples. Since the model with its extremely stiff differential equation (due to the exponentially rising current), forward Euler is not really a good fit for the aeif models. NEST, on the other hand, uses an adaptive step-size Runge-Kutta-Fehlberg integrator. Thus, the results obtained with NEST should be closer to the mathematically correct solution of the aeif equations than the results produced by Brian. For details on NEST's aeif implementation, see https://nest-simulator.readthedocs.io/en/stable/model_details/aeif_models_implementation.html.

We are currently working on a version of the aeif models that would allow the user to set an upper limit for the exponentially growing spike current (https://github.com/nest/nest-simulator/pull/2755). According to my calculations (see in the PR), in the original Brette & Gerstner implementation of the model (from their 2005 paper), this current could never exceed about 9 nA.

On a deeper level, this raises the question of how the aeif model is actually defined: In terms of the equations, assumed solved as exactly as possible, or in terms of the equations plus the original implementation, i.e., forward Euler with a fixed (but in the paper not specified) time step?

[github-actions[bot]]

Issue automatically marked stale!