Numerical instability in aeif_psc_alpha neuron model

[zbarni]

Describe the bug
While running some experiments together with @TeEsTeBe, we found the implementation of the aeif_psc_alpha neuron model appears to be numerically unstable for certain RNG seeds, with two different behaviors as described below:

To Reproduce
Steps to reproduce the behavior:

1. Consider the following example script

import nest
import numpy as np

nest.ResetKernel()
nest.set_verbosity('M_ALL')
neuron_pars = {
            'E1': {  # IT cell firing profile
                # "model": "aeif_psc_alpha",
                "C_m": 200.,  # capacity of the membrane (pF)
                "g_L": 12.,  # leak conductance (nS)
                "E_L": -70.,  # leak reversal potential (mV) [in other models: resting potential]
                "V_th": -50.,  # spike initiation threshold (mV)
                "Delta_T": 2.,  # slope factor (mV)
                "a": 2.,  # sub-threshold adaptation (ns)
                "tau_w": 300.,  # adaptation time constant (ms) [large]
                "b": 60.,  # spike-triggered adaptaton (pA)
                "V_reset": -58.,  # reset value for Vm after a spike (mV)
                "I_e": 0.,  # constant external input current (pA)
                'tau_syn_ex': 5.,
                'tau_syn_in': 10.,
                "t_ref": 0.1
            },
            'E2': {  # PT cell firing profile
                # "model": "aeif_psc_alpha",
                "C_m": 200.,  # capacity of the membrane (pF)
                "g_L": 10.,  # leak conductance (nS)
                "E_L": -58.,  # leak reversal potential (mV) [in other models: resting potential]
                "V_th": -50.,  # spike initiation threshold (mV)
                "Delta_T": 2.,  # slope factor (mV)
                "a": 2.,  # sub-threshold adaptation (ns)
                "tau_w": 120.,  # adaptation time constant (ms) [large]
                "b": 100.,  # spike-triggered adaptaton (pA)
                "V_reset": -46.,  # reset value for Vm after a spike (mV)
                "I_e": 0.,  # constant external input current (pA)
                'tau_syn_ex': 5.,
                'tau_syn_in': 10.,
                "t_ref": 0.1
            },
            'I1': {
                # 'model': 'iaf_psc_exp',
                'C_m': 250.,
                'E_L': -70.0,
                'I_e': 0.,
                'V_m': -70.0,
                'V_th': -55.0,
                'V_reset': -70.0,
                't_ref': 2.,
                'tau_m': 20.,
                'tau_syn_ex': 5.,
                'tau_syn_in': 10.
            },
            'I2': {
                # 'model': 'iaf_psc_exp',
                'C_m': 250.,
                'E_L': -70.0,
                'I_e': 0.,
                'V_m': -70.0,
                'V_th': -55.0,
                'V_reset': -70.0,
                't_ref': 2.,
                'tau_m': 20.,
                'tau_syn_ex': 5.,
                'tau_syn_in': 10.
            }
}

kernel_pars = {'resolution': 0.1, 'data_prefix': 'asd',
               'data_path': 'asd',
               'overwrite_files': True, 'print_time': True, 'total_num_virtual_procs': 6,
               'local_num_threads': 6,
               'rng_seed': 130030284    # !!! -> NEST hangs
               # 'rng_seed': 1234       # !!! -> Error message
               }
nest.SetKernelStatus(kernel_pars)
np_rng = np.random.default_rng(kernel_pars['rng_seed'])
print(f"NEST RNG: {nest.GetKernelStatus('rng_seed')}")

NE = 2000
NI = 500

for pop_name, pars in neuron_pars.items():
    model = 'aeif_psc_alpha' if pop_name[0] == 'E' else 'iaf_psc_exp'
    nest.CopyModel(model, pop_name)
    accepted_keys = list(nest.GetDefaults(model).keys())
    accepted_keys.remove('model')
    nest_dict = {k: v for k, v in pars.items() if k in accepted_keys}
    nest.SetDefaults(pop_name, pars)

# create populations
E1 = nest.Create('E1', NE)
I1 = nest.Create('I1', NI)
E2 = nest.Create('E2', NE)
I2 = nest.Create('I2', NI)
populations = [E1, I1, E2, I2]

# randomize Vm
for pop in populations:
    pop.V_m = list(np_rng.uniform(low=pop[0].E_L, high=pop[0].V_th, size=len(pop)))

# create input
pg = nest.Create('poisson_generator', 1, {'rate': 2400.})

# connect input
for pop, w in zip(populations, [10., 40., 10., 20.]):
    nest.Connect(pg, pop, 'all_to_all', syn_spec={'weight': w})

epsilon = 0.1
delay = 1.5
wscale = 1.
wE_IT = 17. * wscale  
wE_PT = 32. * wscale 
wE_IT_to_PT = 14.7 * wscale  
gamma = 10.  # scaling factor for the inhibitory synapses.
nu_x = 12.  # external background input intensity.
wIx = 4.960628287400623

conn_specs = [
    # module 1 internal
    {'allow_autapses': False, 'allow_multapses': False, 'rule': 'pairwise_bernoulli', 'p': epsilon / 2.},
    {'allow_autapses': False, 'allow_multapses': False, 'rule': 'pairwise_bernoulli', 'p': epsilon},
    {'allow_autapses': False, 'allow_multapses': False, 'rule': 'pairwise_bernoulli', 'p': epsilon},
    {'allow_autapses': False, 'allow_multapses': False, 'rule': 'pairwise_bernoulli', 'p': epsilon},
    # module 2 internal
    {'allow_autapses': False, 'allow_multapses': False, 'rule': 'pairwise_bernoulli', 'p': epsilon / 2.},
    {'allow_autapses': False, 'allow_multapses': False, 'rule': 'pairwise_bernoulli', 'p': epsilon},
    {'allow_autapses': False, 'allow_multapses': False, 'rule': 'pairwise_bernoulli', 'p': epsilon},
    {'allow_autapses': False, 'allow_multapses': False, 'rule': 'pairwise_bernoulli', 'p': epsilon},
    # feed-forward excitation
    {'allow_autapses': False, 'allow_multapses': False, 'rule': 'pairwise_bernoulli', 'p': epsilon / 2.},
    # feedback excitation
    {'allow_autapses': False, 'allow_multapses': False, 'rule': 'pairwise_bernoulli', 'p': epsilon / 2.},
]

syn_specs = [
    # module 1 internal
    {'synapse_model': 'static_synapse', 'delay': delay, 'weight': wE_IT},
    {'synapse_model': 'static_synapse', 'delay': delay, 'weight': wE_IT * 2},
    {'synapse_model': 'static_synapse', 'delay': delay, 'weight': -gamma * wE_IT},
    {'synapse_model': 'static_synapse', 'delay': delay, 'weight': -gamma * wIx},
    # module 2 internal
    {'synapse_model': 'static_synapse', 'delay': delay, 'weight': wE_PT},
    {'synapse_model': 'static_synapse', 'delay': delay, 'weight': wE_PT},
    {'synapse_model': 'static_synapse', 'delay': delay, 'weight': -gamma * wE_PT},
    {'synapse_model': 'static_synapse', 'delay': delay, 'weight': -gamma * wIx},

    # # feed-forward
    {'synapse_model': 'static_synapse', 'delay': delay, 'weight': wE_IT_to_PT},
    # feedback
    {'synapse_model': 'static_synapse', 'delay': delay, 'weight': wE_PT},  # no data for this connection, assumed
]

connections = [
    # module 1 internal
    (E1, E1),  # E1 ---> E1 (recurrent E)
    (I1, E1),  # E1 ---> I1
    (E1, I1),  # I1 ---> E1
    (I1, I1),  # I1 ---> I1
    # # module 2 internal
    (E2, E2),
    (I2, E2),
    (E2, I2),
    (I2, I2),
]

for idx, conn_tuple in enumerate(connections):
    nest.Connect(conn_tuple[1], conn_tuple[0], conn_spec=conn_specs[idx], syn_spec=syn_specs[idx])

spike_recorders = nest.Create("spike_recorder", 4)
for idx in range(4):
    nest.Connect([E1, I1, E2, I2][idx], spike_recorders[idx])

nest.Simulate(1.)
nest.Simulate(35.)

# print output
sr = spike_recorders[0]
print("Spikes:", sr.events['times'][:10])
print("Spikes:", sr.events['senders'][:10])

from matplotlib import pyplot as plt

fig, axes = plt.subplots(1, 4, figsize=(20, 6))
cnt = 0
for idx, sr in enumerate(spike_recorders):
    axes[idx].scatter(sr.events['times'], sr.events['senders'], s=1)
    axes[cnt].set_xlim((-5, 40))
    cnt += 1
fig.tight_layout()
fig.savefig(f"raster_merged_test.jpg")

2. Used parameters
   See script. We tried two different RNG seeds, which lead to two different behaviors:

   • for { ... 'rng_seed': 130030284}, the simulation hangs at 31 ms, likely due to the GSL solver not converging.

   • for { ... 'rng_seed': 1234}, the simulation exits with the error message

     nest.lib.hl_api_exceptions.NumericalInstability: NumericalInstability in SLI function Simulate_d: NEST detected a numerical instability while updating aeif_psc_alpha.

3. Command to run
   Execute above script using Python 3.8+.

4. How to see error
   See point 2., for the two different RNG seeds.

Expected behavior
The simulation should exit after 35.0 ms without errors.

Desktop/Environment (please complete the following information):
We tested the script in two different environments on different PCs, also using multiple NEST versions 3.x

• OS: Linux Mint 19.3, kernel 5.4.0-050400-generic / Ubuntu 20.4
• Python-Version: 3.9.0 / 3.8.8
• NEST-Version: 3.0, 3.3, master@e0f8991aa
• Installation: manual install inside Miniconda env, see attached file.

Additional context
Add any other context about the problem here.

[clinssen]

I was able to reproduce this. For 1 thread instead of 6, I get a NumericalInstability exception, which turns out to be due to V_m going below -1000. Dumping the values, it looks like the inhibitory current is indeed quite strong. The result of an earlier integration failure? Setting the timestep to 0.01 doesn't help either.

NumericalInstability: t = 30
         State_::V_M = -1003.7
         State_::I_EXC = 2867.42
         State_::DI_EXC = 329.099
         State_::I_INH = 42880.1
         State_::DI_INH = 5319.17
         State_::W = -53.0438

[heplesser]

Due to the exponentially growing current in the AEIF models, it is fiendishly difficult to integrate numerically. Some other implementations of the model simply side-step this by using forward Euler for integration, but NEST strives to be exact and uses an adaptive-stepsize RKF45 solver. To keep numerics somewhat under control, we limit V_m to V_peak when evaluating the right-hand side of the ODEs, but for your parameters (V_th=-50 mV, V_peak=0 mV, Delta_T=2 mV, g_L=12 mV), this means a maximum current of about 1.7 A, while a single synaptic input current has a peak of 320 pA. The exponential current thus is some nine orders of magnitude larger than the synaptic input currents, and it grows very fast.

The current NEST implementation of AEIF models is discussed in https://nest-simulator.readthedocs.io/en/latest/neurons/model_details/aeif_models_implementation.html.

One easy way of stabilizing the model somewhat would be to reduce V_peak, but that has a side effect, because V_peak not only limits the membrane voltages upwards upon rhs-evaluation, it is also the voltage at which a spike is triggered.

A better approach, which is also fully biologically justified, would be to limit the exponential current. In a biological neuron, this current is obviously bounded (limited number of ion channels and ions, inactivation of Na channels at high voltage), so setting some upper bound is justified. To get a ballpark number for this limit, one could try to estimate the current if all Na-channels on the soma of a reasonably representative neuron were fully open simultaneously while the neuron is still well polarised. This could enter as a new parameter I_exp_max. Could you give it a try?

[github-actions]

Issue automatically marked stale!

[heplesser]

@zbarni What do you think about the idea to limit the exponentially growing current? And what upper limit to the current value would you consider biologically plausible?

Looking at the parameters of the hh_cond_exp_traub, an upper limit seems to be around 20000 nS x 100 mV = 2.000.000 pA = 2 µA, i.e., a maximum current about 10^6 times smaller than the one in the model at present. This would certainly contribute to reducing numerical instability, but also delay spike initiation (in a biologically plausible way).

[github-actions]

Issue automatically marked stale!