Multiple events for a single spike in hh_cond_exp_traub when increasing time resolution

[astoeckel]

When running a simulation with the hh_cond_exp_traub neuron model, an increased time resolution (a time step smaller than the default 0.1ms), causes single output spikes in the HH model to be issued multiple times. The cause of issue should be

void
nest::hh_cond_exp_traub::calibrate()
{
  B_.logger_.init(); // ensures initialization in case mm connected after Simulate
  V_.RefractoryCounts_ = 20;
  V_.U_old_ = S_.y_[ State_::V_M ];
}

in hh_cond_exp_traub.cpp where the number of ticks in V_.RefractoryCounts_ is initialised with the constant 20 instead of calculating the number of ticks depending on the time step.

However, I'm not convinced that introducing an artificial, time-dependent refractoriness to the Hodgkin-Huxley model is actually the right way to go (actually, this causes bug #330). In my implementation of the HH model [1] I simply set a boolean in_refrac to true once a spike is issued (at the maximum membrane potential above a certain threshold) and reset this variable once the membrane potential drops below the threshold potential. This way no assumption about timings is made.

[1] https://github.com/hbp-sanncs/cinder/blob/ede53f1070b2b0a4eeed90278b61eea6db4c775d/cinder/models/neurons/hodgkin_huxley.hpp#L197

[Silmathoron]

I agree with @astoeckel that it seems the way refractoriness is considered in the hh_cond_exp model is rather weird: on hybrid models such as the iaf or aeif, the voltage is clamped during the whole refractory period so there is no risk of missing spike. Such a feature can be necessary to get closer to real neural behaviour because these simple models have no intrinsic equivalent of refractoriness.

On the contrary, the hh model has a current that induces to a "real" refractory period, so implementing an "artificial" refractory period which is not related to the intrinsic parameters and does not clamp the potential seems a bit strange... and indeed seems to lead to problems.

Wouldn't it be better to propose either no such refractory period or one that (as in the other models) clamps the potential?

EDIT: just for information, here is the implementation for hh in PyNN for Brian and I don't see any refractory periods in the parameters (regardless of the simulator).

[heplesser]

@astoeckel Thank you for reporting this and #330. I suggest we discuss both issues here. @Silmathoron Thanks for the reference to the PyNN for Brian implementation.

I agree that the current implementation of the hh_cond_exp_traub model badly abuses the refractoriness mechanism. @astoeckel's solution to emit a spike when a threshold crossing is detected for the first time and then to suppress all further spikes until V_m < V_th, make sense. But before committing to a solution, we should take a look at the original definition of the model. Does anyone have the original paper or book chapter available?

[ynodem]

Hi @heplesser , here is a good reference about this model.

Biol Cybern. 2008 Nov;99(4-5):427-41. doi: 10.1007/s00422-008-0263-8. Epub 2008 Nov 15. Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Pospischil M1, Toledo-Rodriguez M, Monier C, Piwkowska Z, Bal T, Frégnac Y, Markram H, Destexhe A.

[heplesser]

@yan1982 Thank you, I will take a look at it shortly!

[heplesser]

See also #473 .

[JonThom]

Hi,

I was wondering what the status is for hh_cond_exp_traub refractoriness.

Issue: neurons fire 'bursts' of 3-4 spikes with 0.1ms inter-spike-intervals

Setup:

• Python 3.6
• NEST 2.12

Example:

A balanced network, in part based on Brunel's network. The bursts have occurred with a range of network scales and parameters. The attached script should reproduce the problem and print a spiketrain.

As far as I can see hh_cond_exp_traub is the only conductance based HH model. Any comments or suggestions?

Many thanks in advance!

Jonatan

hh_cond_exp_traub_test.py.tar.gz

[Silmathoron]

Hi Jonatan, the hh_cond_exp_traub's refractoriness was updated for 2.12, is that the version you are using ?

As a general rule, it is best to always provide:

1. the Python and NEST versions you are using
2. a small example that reproduces the problem you are talking about

Given the way the spikes are detected for the HH model, what you describe seems indeed very unreasonable, but I cannot say much without seeing it.
But it is true that apart from fixing the incorrect refractoriness, I don't think anyone worked on this in details... so there might be some lingering problems...

In the same way, adding refractoriness, though I indeed think it is questionable for the HH model, is expected to change the dynamics, so again, it is difficult to comment without seeing what is happening in your code.

[JonThom]

Please see the updated comment above (updated again with PyNEST script)

[Silmathoron]

Sorry, maybe my comment was unclear: please post code directly, not text, it will save us time.

[heplesser]

590 fixed #473, but this only means that hh_cond_exp_traub now handles the refractory time as configurable parameter and handles it correctly wrt changes in resolution. It is still a fixed refractory time.

The reference provided by @yan1982 does not provide the full model description. The documentation of the NEST implementation contains the following "Todo" remark: "Only the channel variables m,h,n are implemented. The original contains variables called y,s,r,q and \chi." The latter variables are not mentioned inthe Biol Cybern paper. So we would really need access to the original Traub & Miles (1991) definition of the model.

[JonThom]

Update:

• I used a spike generator to feed a single neuron several spikes at 0.1ms intervals, starting at t = 10 ms. The voltage and conductance traces show that the incoming spikes change the conductances but apparently not the voltage dynamics during the spike.

Without having examined the neuron's code, it seems to me that, as @astoeckel first noted, the neuron records an event for every timestep where it is above the threshold potential. With a resolution of 0.1 ms, (at least in a balanced network as the one posted above), we tend to get five events per spike. Presumably during a spike, the membrane is above threshold for approximately 0.5 ms, resulting in five events being registered (at 0.1 ms resolution). A temporary quick fix may be to set a 0.5 ms refractory period, resulting in each spike producing a single event. Since (I suppose) this corresponds to the duration of the part of the spike where the neuron is above threshold, it seems to me to be a reasonable constraint in any case.

[heplesser]

@JonThom Thank you for your analysis! I took another look at the code now, and it seems to me that the spike-detection method is problematic unless one has long refractory periods. The documentation for the model states that:

 (2) Spike Detection
 Spike detection is done by a combined threshold-and-local-maximum search: if
 there is a local maximum above a certain threshold of the membrane potential,
 it is considered a spike.

This is in principal a fine idea, but the implementation is problematic. For each time step, the update() method first stores the old membrane potential

    V_.U_old_ = S_.y_[ State_::V_M ];

then advances the state to the end of the time step and then checks for spikes (unless the neuron ise refractory) using

      if ( S_.y_[ State_::V_M ] >= P_.V_T + 30. && V_.U_old_ > S_.y_[ State_::V_M ] )

Upon a spike, the membrane potential is not reset. The first term ensures membrane potential is above threshold, the second term ensures that membrane potential has fallen during the time step, i.e., that we have passed a maximum. If the refractory time is short, it can happen we have several time steps with falling V_m with values all above the threshold, and then we get a spike for each one of them.

This is clearly wrong, and I assume the model was originally designed/used with refractory times so long that Vm would be below threshold before the end of the refractory time. Then, you would only get a single spike. At the very least here, we need a check that ensures that no new spike is fired before membrane potential has fallen below threshold. I am also no happy with the `P.V_T + 30expression for the threshold.V_Tis a reference voltage used in the channel dynamics; I think the model should have an independent paramterV_th, as threshold for "calling" a spike: the first drop in membrane potential afterV_th` has been crossed.

[jougs]

@ingablundell, @DimitriPlotnikov: Can you please have a look at this in the context of NESTML? The comment by @heplesser explains how things should be. Thanks!

[DimitriPlotnikov]

@jougs I believe in NESTML it also depends on how the model is written and we cannot do any magic to force the modeler to do it some way and forbid another.

[Silmathoron]

@DimitriPlotnikov No magic, but a check could be made to test whether the modeler included a reset condition on the mambrane potential, and if not, warn about potential problems for spike detection... And anyway, if we want a correct implementation for the HH model in NEST(ML), we need to correct the spike detection ;)

[DimitriPlotnikov]

@Silmathoron Yes, it is doable.

[heplesser]

@astoeckel @Silmathoron @JonThom @diesmann I have now looked more thoroughly at this model (see also #944) and the situation is somewhat complicated.

1. The model is called ..._traub, which is a notation which we in NEST reserve for models obtain from specific papers or books (it should really be ..._traub_miles or ..._tm_1991 to be consistent with other model names.
2. The model in Traub and Miles (1991) is significantly more complex than the model implemented in NEST at present. The original Traub-Miles model
   1. is a multicompartment model
   2. has two variables governing I_K activation/inactivation, NEST has only one
   3. has an I_Ca, governed by two activation/inactivation variables
   4. has an I_K[Ca], governed by one variable
   5. synapses are modelled as ds/dt = S_in(t) - s(t), i.e., as exponential decay, but the input variable S_in(t) is not a delta-function as in NEST, but a function that is set to 1 for a certain period of time after arrival of a spike; that period can be 1 ms, 3 ms, or 40 ms (for slow inhibition).
   6. emits output spikes according to the following rule (p 106): "The cell sends an output if its depolarized more than a threshold amount (20 mV depolarized relative to resting potential) and if no output has been sent for 3 ms".
      • The current implementation in NEST was created as part of NEST's participation in the Brette et al (2007) simulator comparison paper to support the hh_coba benchmark, see Appendix B of the paper. The paper does not discuss the simplifications made compared to the Traub-Miles book. Worse, it does not specify the rule for spike detection used.
      • Code for the benchmarks in the Brette-paper is available on ModelDB, where the spiking rule used is most easily found from the Brian code:

        P=NeuronGroup(4000,model=eqs,
         threshold=EmpiricalThreshold(threshold=-20*mV,refractory=3*ms),
         implicit=True,freeze=True,compile=False)

        This is quite close to the Traub-Miles model, except that the threshold is -20 mV instead of 20 mV above resting potential.

I think we should do two things:

1. Add to NEST a new model, e.g., hh_cond_exp, with adaptive stepsize integration and a proper spike detection mechanism along the lines sketched by @astoeckel previously. Users should use that model if they need HH-dynamics in NEST.
2. Rename the current hh_cond_exp_traub model to something like hh_cond_exp_2007_review with spike detection and refractoriness as in the original paper. This model would only be retained to support the Brette et al hh_coba benchmark case. We should postpone this change to NEST 3, since the renaming of hh_cond_exp_traub may break user code.

[Silmathoron]

Thanks for the nice summary @heplesser ! I agree with both options. Regarding the implementation of the spike detection mechanism, however, I'm not sure @astoeckel proposal is sufficient, though it is better than current one.

Problems could occur for high firing rates where the neuron might spike several times without going below V_th (unless we find a "good" value for it, but this would still be risky). It might be necessary to either:

• check for a local minimum before V_th has been crossed (tricky if there is noise...)
• implement a 2-step threshold (V_th_low, V_th_high), where crossing V_th_low out of refractory triggers a spike and crossing V_th_high during refractory period also does.

[Silmathoron]

An update on the issue of spike detection with the HH model, I have been playing with BRIAN lately and noticed that they have an example of HH where they use the sodium gate variable m to detect spike events and handle refractoriness. This intuitively sounds like an interesting way of doing it, so I'm adding it here for reference: https://brian2.readthedocs.io/en/2.0rc/examples/compartmental.hh_with_spikes.html

[heplesser]

When addressing this issue, also check other affected models, e.g., hh_cond_beta_gap_traub.

[heplesser]

See also #944.

[clinssen]

Some of the suggestions on this thread have already been merged, while others have not.

I suggest I create a PR to rename this model (hh_cond_exp_traub) to hh_cond_exp_2007_review as suggested by @heplesser, and to add documentation to the models (also hh_cond_beta_gap_traub) to document rather than change the behaviour, because these are not strictly speaking bugs but implementation details/quirks.

Document:

• Hard-coded threshold and note discrepancy with original paper
• Membrane potential is not reset during refractory period and synaptic integration continues during refractory period

NEST3 milestone would be good because the rename breaks the API.

[heplesser]

@clinssen Good points! We should discuss the underlying question in general: To which degree should we implement models from the literature "as is", and to which degree should we create well-designed models with solid numerics based on models from the literature.

I believe that we will do neuroscience a service by creating well designed and implemented models, but for reproducibility checks one might occasionally want to have the original in place as well. But maybe we should put those into a special place?

Concerning naming, if we change it I would go for hh_brette_2007 for some consistency with iaf_chs_2007 and iaf_chxk_2008.

[clinssen]

@heplesser: if we rename "hh_cond_exp_traub" to "hh_brette_2007", what should we rename "hh_cond_beta_gap_traub" to? It's a variant of "hh_cond_exp_traub" with beta function synaptic kernels and gap junctions added. "hh_brette_2007_beta_gap"? "hh_beta_gap_brette_2007"?

[heplesser]

@clinssen It is unfortunate that any neat scheme will collide with reality sooner rather than later. I would want to take this topic up in a developer meeting. One approach could be to try to find out who created and where it was first published, so that we could base our naming on that paper. Otherwise, I'd prefer hh_brette_2007_beta_gap, because the model is a derivative of hh_brette_2007.

A side note, @sarakonradi: The documentation of the ...beta_gap... variant is rather too close to the original model. I am quite certain that this model was not created in connection with the Brette et al (2007) paper, but much later when gap junctions were introduced.

[github-actions[bot]]

Issue automatically marked stale!