Regression in NEST 2.8.0 for 'aeif_cond_exp' with Delta_T = 0

[apdavison]

When Delta_T = 0, the membrane voltage should go to infinity as soon as "V_th" is reached (Scholarpedia). This is handled correctly in NEST 2.6.0, but in NEST 2.8.0 there is no spike until the voltage reaches "V_peak".

Example:

import matplotlib.pyplot as plt
import nest

neuron = nest.Create('aeif_cond_exp')
nest.SetStatus(neuron,
               dict(Delta_T=0.0, I_e=1000.0, V_th=-50.0, V_peak=-45.0))

recorder = nest.Create('multimeter')
nest.SetStatus(recorder, {'record_from': ['V_m'],
                          'interval': nest.GetKernelStatus('resolution')})
nest.Connect(recorder, neuron)

nest.Simulate(100.0)

data = nest.GetStatus(recorder, 'events')[0]
t = data['times']
vm = data['V_m']

plt.plot(t, vm)
plt.xlabel("Time (ms)")
plt.ylabel("V_m (mV)")
plt.title(nest.version())
plt.ylim(-70, -45)
plt.savefig("aeif_conf_exp_{}.png".format(nest.version().replace(" ", "_")))

[heplesser]

@apdavison Thanks for spotting and reporting this. I can confirm the observation.

I have traced it back to a change in how we handle excessively large exponents (svn r11833).

In 2.6.0, we had in aeif_cond_exp_dynamics()

static const double_t largest_exp=std::exp(10.);
const double_t exp_arg=(V - node.P_.V_th) / node.P_.Delta_T;
const double_t I_spike = (exp_arg>10.)? largest_exp : node.P_.Delta_T * std::exp(exp_arg);

Note that that code was not entirely correct, since we, for exp_arg > 0 did not multiply by Delta_T, we simply set I_spike to e^10 == 22026.5. For Delta_T==0, the division by zero gives exp_arg==inf,I_spike becomes 22026.5, and the neuron spikes in the next step.

In 2.8.0, we have instead

  const double_t exp_arg = ( V - node.P_.V_th ) / node.P_.Delta_T;
  const double_t MAX_EXP_ARG = 10.;
  const double_t I_spike = node.P_.Delta_T * std::exp( std::min( exp_arg, MAX_EXP_ARG ) );

so we now (correctly) multiply with Delta_T under all conditions, but for Delta_T==0 this leads to I_spike==0 under all conditions.

This affects all aeif_cond_* models.

[heplesser]

@gewaltig Could you take a look at this. I am not entirely certain what would be the best way to handle this case without adding if ( Delta_T == 0 ) in several places, particularly since the adaptive stepsize algorithm needs to be able to handle whatever we do in the dynamics() function. Maybe we can solve this in conjunction with #35?

[gewaltig]

Thanks for pointing this out. I’ll look at the problem.

Best Marc-Oliver

On 29 Sep 2015, at 16:50, Hans Ekkehard Plesser notifications@github.com<mailto:notifications@github.com> wrote:

@gewaltighttps://github.com/gewaltig Could you take a look at this. I am not entirely certain what would be the best way to handle this case without adding if ( Delta_T == 0 ) in several places, particularly since the adaptive stepsize algorithm needs to be able to handle whatever we do in the dynamics() function. Maybe we can solve this in conjunction with #35https://github.com/nest/nest-simulator/pull/35?

— Reply to this email directly or view it on GitHubhttps://github.com/nest/nest-simulator/issues/114#issuecomment-144084129.

[heplesser]

@gewaltig: I talked some more about this with @apdavison. From a PyNN perspective, it would be ok to not support Delta_T==0. NEURON does not support that edge case either, PyNN then substitutes the corresponding I&F model. So throwing an exception when the user tries to set Delta_T==0 and recommending iaf_psc_* is probably the easiest solution.

Another interesting question is what happens if Delta_T is small but not strictly zero.The default value is 2 mV, so

const double_t I_spike = node.P_.Delta_T * std::exp( std::min( exp_arg, MAX_EXP_ARG ) );

is larger than the exponential. But if we had Delta_T == 0.0001 mV, then for the limiting case, we have I_spike = 0.0001 * std::exp(10) = 2.2, which probably will not take membrane potential to V_peak within a single time step.

[apdavison]

@heplesser I think you misunderstood me this morning. NEURON didn't support this edge case, but it does now (see https://github.com/NeuralEnsemble/PyNN/issues/367). It was running the PyNN regression tests with NEST 2.8.0 which uncovered the issue. Nevertheless, if it is too difficult to support this in NEST going forward, I am happy with the Exception approach.

[apdavison]

Looking at this again, it is not true that when Delta_T==0 the model is the same as iaf_psc_*, because there is still the adaptation term - w

[gewaltig]

I am somewhat hesitant here. As a fix, I would just go back to the old code, as it behaved correctly. At the time I deliberately "precomputed" the result to achieve the desired threshold crossing. On a more general note, I would like to remove the old GSL based adex versions in favor of the faster and more accurate RK5 versions.

[heplesser]

@apdavison I have tried for some time now to find a numerically stable solution supporting Delta_T==0 that carries minimal overhead, but have not been entirely successful. So I fear we won't be able to fix this for the next NEST release. It will have to wait for a new implementation of aeif-numerics in context of #35.

[jougs]

@heplesser, @apdavison, @gewaltig please also take into account https://github.com/nest/nest-simulator/pull/35#issuecomment-186309378

[tammoippen]

See also #229

[heplesser]

In the dynamics function computing the derivative, we must use min( y [ S::V_m ], node.P_.V_peak ) wherever V_m appears in expressions. I will fix this in the standard AEIF-model shortly. This ensures that w will not be "polluted" by super-peak values of V_m and keep numerics stable.

[heplesser]

This have been fixed by #474 .