Implementing condition v value in STN (AdEx)

Science, Projects, and Showcases Projects
1

I am building an STN neuron model using Brian2 to simulate burst activity. The burst mechanism involves resetting the membrane potential (v) after a spike using the rule:

v_reset+max⁡(u−15,20) (when u<0)

I tried implementing this condition in the reset argument of NeuronGroup within the create_neurons function. However, Brian2 does not seem to support max() directly. What is the proper way to apply this reset condition in Brian2? I followed the instructions in the image below and applied.

image
image750×110 19.1 KB

class STN(NeuronModel):
    def __init__(self, N, params):
        super().__init__(N, params)
        self.neurons = None

    def create_neurons(self):
        eqs = AdEx.eqs
        reset = '''
        v = vr + ((u < 0) * ((u - 15*pA) * (u - 15*pA > 20*pA) + 20*pA * (u - 15*pA <= 20*pA)));
        u += d
        '''
        self.neurons = NeuronGroup(
            self.N, eqs, threshold='v > th', reset=reset, method='euler'
        )

        self.neurons.g_L = self.params['g_L']['value'] * eval(self.params['g_L']['unit'])
        self.neurons.E_L = self.params['E_L']['value'] * eval(self.params['E_L']['unit'])
        self.neurons.Delta_T = self.params['Delta_T']['value'] * eval(self.params['Delta_T']['unit'])
        self.neurons.vt = self.params['vt']['value'] * eval(self.params['vt']['unit'])
        self.neurons.vr = self.params['vr']['value'] * eval(self.params['vr']['unit'])
        self.neurons.tau_w = self.params['tau_w']['value'] * eval(self.params['tau_w']['unit'])
        self.neurons.th = self.params['th']['value'] * eval(self.params['th']['unit'])
        self.neurons.a = self.params['a']['value'] * eval(self.params['a']['unit'])
        self.neurons.d = self.params['d']['value'] * eval(self.params['d']['unit'])
        self.neurons.C = self.params['C']['value'] * eval(self.params['C']['unit'])

        return self.neurons```
2

Hi @keun, the absence of min/max functions comes up regularly, I think we should fix this ! Brian does support the more powerful clip function, though, which can replace both min and max. In your case, you should be able to use:

reset = 'vr + clip(w - 15*pA, 20*pA, inf*pA)'

which clips values below 20pA to 20pA. I’m a bit confused about the “with w < 0” part, though. What should happen if w ≥ 0? Also, if w < 0, then max(w - 15, 20) will always be 20, no?

3

Hi @mstimberg , thank you for your reply.

If w<0, I want to select the larger of (w-15, 20), but it’s hard to apply that condition with clip because of the unit issue. I tried with the code below, but it doesn’t work the way I want because CLIP is not comparing the specified values, but specifying the position of the min, max values. Is there any way to modify it differently?

4

Hi @keun.

This is the part that I did not understand: if w < 0, then w-15 will be always lower than -15, so 20 will always be the larger value.

I am sorry, but I do not understand this either. What do you mean by “CLIP is not comparing the specified values, but specifying the position of the min, max values”? The clip functions clips a value between a minimum and a maximum value. If you replace the minimum or maximum value by -inf or inf respectively, it will be exactly equivalent to min or max. Regarding units, all three elements should have the same units, so e.g. use clip(w - 15*pA, 20*pA, inf*pA) if w has units of Ampère or use clip(w -15, 20, inf) if w is unitless.

5

Thank you for your kind response. According to the paper, I initially understood that w could also take values less than zero depending on the conditions, but it seems I misunderstood. Since v is reset based on max(w-15, 20) regardless of the range of w, I plan to use clip, as you explained.

However, I have one question: The STN model uses the AdEx model, and the units are in pA (amperes). Would there be any issues in calculations when performing operations with values that do not have explicitly multiplied units?

Once again, I really appreciate your help. It was very useful!

image
image740×114 19.1 KB

6

Hi @keun

From the code you posted it already looks as if you did take care of the units correctly? Everything has to be consistent, so if your model equations describe e.g. w as having units of amp, then you’ll have to write things as clip(w - 15*pA, 20*pA, inf*pA) so that everything matches. But you can also make your variables unitless, and use conventions as in many papers (voltages are in mV, currents are in pA, etc.). The only dimension that you cannot remove is the unit of time, which means that e.g. time constants would still have to be quantities in seconds (or ms, etc.).

7

Hi @mstimberg,

Thank you for your thoughtful response. Your answer was extremely helpful and has allowed me to resolve my issue.

1 Like