Dear Charl,

Thank you very much for your follow up.
I agree that using a linear system of ODEs would be preferable. I made a proxy version that uses a beta kernel (inspired from iaf_cond_beta_neuron.nestml). The rise and decay time constants, and a modulatory parameter to scale the conductance, may be tuned to approximate the polynomials derived by Zhang and colleagues. In the attachment, my result is provided from a simple grid search, comparing the proxy with two exponentials and the polynomial from Zhang and colleagues.    

It may be interesting indeed to perform a ReScience paper comparing the beta kernel version in NESTML with the version from Zhang and colleagues (for your information, their matlab code is available on github). However, in order to do so, another part of their implementation puzzles me. Their total model takes into account the polynomial component, but also multiplied by an inhibition efficiency that relies on the firing rate of the ChI population. In their code, from my understanding, they compute the inhibition at time i + k via a multiplication (efficiency due to ChI firing at time i X inhibition along time via the value of the polynomial at time k), but keep only the maximum value along i and k, rather than actually summing the effect of spikes at distinct times. I am not sure to understand the reason between this modelling choice. Do you think it may be interesting as well to compare this with a more traditional approach in NESTML, where the efficiency is computed online at the time t of the simulation, and the beta kernel allows to integrate all past spikes? Even if we integrate a new modulatory parameter for scaling the entire inhibition response and fit better the average behavior of the model from Zhang and colleagues．      

Please do not hesitate to let me know your thoughts about this.
Thank you,
Kind regards,
Jeanne Barthelemy   
From: Charl Linssen <nest-users@turingbirds.com>
Sent: Thursday, March 5, 2026 1:47 AM
To: users@nest-simulator.org <users@nest-simulator.org>
Subject: [NEST Users] Re: Heaviside function in a Kernel?
 
Dear Jeanne,

Much obliged for the link to the paper and clarification. If a postsynaptic response is well described by a polynomial, it makes total sense to use that as part of model fitting. From a simulation perspective, it may be more efficient though to describe the response in the form of a (preferably linear) system of ODEs. Therefore if you could use the alpha or beta functions for the postsynaptic response, that would indeed be an ideal solution.

That being said, please let me know if this could create an interesting opportunity for a ReScience paper comparing the responses of the alpha versus the polynomial; we could then look into doing a precise simulation of the polynomial more deeply.

My colleagues also kindly attended me to the possibility of using the so-called Linear Chain Trick in the context of delay differential equations, which could be an attractive alternative method for fitting and simulating postsynaptic responses.

Cheers, with kind regards,
Charl


On Tue, Mar 3, 2026, at 06:04, Jeanne Barthelemy wrote:
Dear Charl,

Thank you very much for this relevant answer.
I was trying to reproduce with NESTML the model from this paper: Zhang, Yan-Feng, et al. "An axonal brake on striatal dopamine output by cholinergic interneurons." Nature Neuroscience 28.4 (2025): 783-794.".  Though, I am starting to think that approximating their model with an alpha kernel, without any need for an heaviside function, may be a better option in my case.  

Thank you very much,
Kind regards,
Jeanne Barthelemy


From: Charl Linssen <nest-users@turingbirds.com>
Sent: Monday, March 2, 2026 11:47 PM
To: users@nest-simulator.org <users@nest-simulator.org>
Subject: [NEST Users] Re: Heaviside function in a Kernel?
 
Dear Jeanne,

Thank you for writing in and thank you for this interesting question. I'm sorry that the advanced kernel shape you need is not supported by the numerical backend. The central issue with this kernel is that it cannot be described by a linear and time-invariant system. As such it makes adding contributions from different spikes much harder, because you cannot simply add the effects of different spikes together (something that only makes sense mathematically for LTI systems).

Actually, there are two components that make this a tricky kernel to work with: the factor t^2 as well as the Heaviside function. May I ask how you obtained the shape and parameters of this kernel? Would it be possible to replace it with an LTI kernel like the “alpha” or “beta” function? (You can have a look at the neuron models in NESTML for a definition of both.)

As for the Heaviside function: I tried to think of an approach to implement this behaviour. One thing you could do is to create an additional buffer in the model that corresponds to the 200 ms time interval. When a spike comes in (let's say at time t_sp), it would then not only trigger the kernel response (ignoring the Heaviside term), but also add an entry to the buffer, so that 200 ms later, a negative contribution would be added to the state, cancelling out the value at that time to reset the value back to zero. The magnitude of this negative contribution would be of equal magnitude to the kernel at time t_sp + 200 ms, but with an opposite sign, so that it cancels out. I believe this would work well, even with multiple spikes arriving all the time.

If you can confirm that you are interested in making an implementation of this, I would be happy to look at it together.

Much obliged,
Charl


On Fri, Feb 27, 2026, at 10:20, Jeanne Barthelemy wrote:
Dear all, 

This is a question about best practices for NESTML. I would like to build in the equations block a Kernel, such as   for instance:

"kernel inh_time_end = ((-0.00002*t*t)/(pow(tau,2.0)) + 0.0087*t/tau + 0.006)/tau"

, but modulated by the heaviside function H depending on 200 - t (i.e., H(200.0 ms -t)). To sum up, I would like to get the following kernel:

 "kernel inh_time_end_modulated = H(200.0 ms -t)*((-0.00002*t*t)/(pow(tau,2.0)) + 0.0087*t/tau + 0.006)/tau"

There does not appear to be a direct implementation of the heaviside function H in nestml.  I tried to build my own function to convert booleans to 1.0 or 0.0 in my NestML model, upon the condition "200.0 ms > t", and use that instead of the heaviside, but the model is really slow to build on my laptop. I tried other versions with if conditions in the update block, or ternary operators, but they all led to bugs. So, I am unsure of the correct practice to implement an heaviside-modulated polynomial kernel. If anyone has a good and fast solution to propose I would be very happy to learn about it !

Thank you very much,
Kind regards,
Jeanne BARTHELEMY
_______________________________________________
NEST Users mailing list -- users@nest-simulator.org
To unsubscribe send an email to users-leave@nest-simulator.org

CAUTION: This email originated from outside the organization. Do not click links or open attachments unless you recognize the sender and know the content is safe.
 
_______________________________________________
NEST Users mailing list -- users@nest-simulator.org
To unsubscribe send an email to users-leave@nest-simulator.org


CAUTION: This email originated from outside the organization. Do not click links or open attachments unless you recognize the sender and know the content is safe.