Integrate specific ODEs with `integrate_odes()`

[clinssen]

Feature request: specify which ODEs to integrate via parameters to the integrate_odes() function. Note that there might be interdependencies between the ODEs that need to be respected/checked.

[Silmathoron]

I don't understand this... does it imply some ODEs would not be integrated? Could you illustrate your point?

[jougs]

The basic idea was that we want to support the case where you have multiple sets of coupled ODEs that are mutually independent.

If you for example had two coupled equations U and V and another equation W, which would be independent of the first two, you could specify integrate_odes(U,V) and integrate(W) at different points in the update block or in different conditional blocks. This is currently not possible. If either U or V was missing from the first call to integrate_odes() that would result in an error.

[Silmathoron]

Mhh... this sounds pretty risky and error prone... do you have a specific use case in mind? I this related to the synapses?

[jougs]

I think this was just a "could-be-done" type of idea rather than being driven by a specific use-case. And we're not currently working on it.

[clinssen]

The following potential use case for this feature just came by on the Brian Simulator mailing list: a situation can occur where an integrate-and-fire neuron is driven so strongly by an injected current, that its membrane potential reaches float infinity during its refractory state. Ideally, integration of the membrane potential dynamics should be paused in the refractory state, but other dynamics might need to continue to be integrated. This would be supported in NESTML if we could write integrate_odes("V_m") seperately from, for example, integrate_odes("I_in, I_ex, I_adap").

[Silmathoron]

I'm not sure... I'm pretty convinced I would want my models to crash if the current I'm injecting leads the membrane potential over float infinity within refractory time...

[clinssen]

The issue of how to handle refractoriness came up also in https://github.com/nest/nest-simulator/issues/1097 The issue here is that within an integration timestep, the membrane potential can move away from its clamped value, because numerical integration of V_m continues during this interval. This can be an issue particularly when there is a strong external current stimulating the neurons, and we are using conductance-based synapses so the instantaneous value of V_m determines the synaptic current.

This issue holds for NESTML as well, as refractoriness is typically implemented as

update:
    integrate_odes()
    if r != 0: # neuron is absolute refractory
        r =  r - 1
        V_m = V_reset    # clamp potential **at the end** of each timestep
    elif V_m >= V_th:
        r = RefractoryCounts
        V_m = V_reset    # clamp potential when going into refractory state
        emit_spike()
    end
end

Integrating only specifc ODEs is a potential solution, as this code could be rewritten as:

update:
    if r != 0: # neuron is absolute refractory
        r =  r - 1
        integrate_odes(g_in, g_ex)     # do not integrate V_m: no need to clamp V_m any more
    elif V_m >= V_th:
        r = RefractoryCounts
        V_m = V_reset    # clamp potential when going into refractory state
        emit_spike()
    else:
        integrate_odes()   # g_in, g_ex and V_m
    end
end

Note that an alternative could be to re-formulate the ODE for V_m from:

V_m' = ( -I_leak - I_syn_exc - I_syn_inh + I_e + I_stim ) / C_m

to:

V_m' = r != 0 ? 0. : ( -I_leak - I_syn_exc - I_syn_inh + I_e + I_stim ) / C_m

However, this is not currently supported. (Is that a bug?)

[jougs]

I actually find your last form of writing the equation quite interesting. I would, however, go for a more pythonic instead of the C-style syntax you proposed:

V_m' = 0.0 if r != 0 else (-I_leak - I_syn_exc - I_syn_inh + I_e + I_stim) / C_m

I'm not aware of the existence of this in NESTML, but am not opposed to adding it. @Silmathoron: what do you think?

[Silmathoron]

This actually sounds quite reasonable since we're dealing mostly with hybrid systems... and it would ensure that clamping is performed correctly inside whatever numerical integrator is used. I would also go for the pythonic syntax.