Open zbarni opened 2 years ago
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
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?
Issue automatically marked stale!
@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).
Issue automatically marked stale!
@zbarni After some more analysis of this problem (to be provided later via #2755), I have come to the conclusion that the most pragmatic solution is to limit V_peak
. Indeed, in what I would consider reference implemenations of the model (by Romain Brette in Brian2, Brette & Gerstner (2005), Touboul & Brette (2008)) they use V_peak = V_T + 5 Delta_T
(they call it v_cut
). My experiments indicate that this is rather a bit low. Given that dV/dt in real neurons seems to go up to about 1200 V/s (Borst & Sakmann, 1998), I would suggest to use
$$
V_{\text{peak}} = V_T + \Delta_T \ln \frac{C_m \times 1200 mV/ms}{g_L \Delta_T}
$$
i.e., set V_peak
to the voltage at which the exponential current drives the membrane potential at the maximum biologically plausible rate.
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:
nest.ResetKernel() nest.set_verbosity('M_ALL') neuron_pars = { 'E1': { # IT cell firing profile
"model": "aeif_psc_alpha",
}
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
]
syn_specs = [
module 1 internal
]
connections = [
module 1 internal
]
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")