[beginner question]step through brian2 to understand inner workings of neuronal models

[akatav]

Hi, Beginner here. I don't have academic background in computational neuro theory or the various neuronal/synaptic models.I tried to step through the HH example HH.py in neurodynex3 which calls brian2.Network.run() function.I think this part is code generated using the math equations object given.How to step through the code that actually loops through the inner working of the HH model ? i think this code is in C++.can you please point me to this code ?

[mstimberg]

Hi @akatav, and sorry for the late reply. For this type of question, please rather go to the discussion forum at https://brian.discourse.group, we try to keep the issue tracker about bugs and feature discussions.

But to reply to your question: the code that gets generated by Brian is not necessarily very easy to read/follow, so I am not sure that stepping through it is the best approach. If you are interested in how the equations are numerically integrated, you can apply the integration method to the equations, and it will give you the step-by-step update of the variables. For the example of the HH model, this would look like this:

from brian2 import Equations, exponential_euler
eqs = Equations("""
I_e = input_current(t,i) : amp
membrane_Im = I_e + gNa*m**3*h*(ENa-vm) + 
    gl*(El-vm) + gK*n**4*(EK-vm) : amp
alphah = .07*exp(-.05*vm/mV)/ms    : Hz
alpham = .1*(25*mV-vm)/(exp(2.5-.1*vm/mV)-1)/mV/ms : Hz
alphan = .01*(10*mV-vm)/(exp(1-.1*vm/mV)-1)/mV/ms : Hz
betah = 1./(1+exp(3.-.1*vm/mV))/ms : Hz
betam = 4*exp(-.0556*vm/mV)/ms : Hz
betan = .125*exp(-.0125*vm/mV)/ms : Hz
dh/dt = alphah*(1-h)-betah*h : 1
dm/dt = alpham*(1-m)-betam*m : 1
dn/dt = alphan*(1-n)-betan*n : 1
dvm/dt = membrane_Im/C : volt
""")
print(exponential_euler(eqs))

which will give you

_BA_h = 0.07*exp(-0.05*vm/mV)/(ms*(-1.0/(ms + 20.0855369231877*ms*exp(-0.1*vm/mV)) - 0.07*exp(-0.05*vm/mV)/ms))
_h = -_BA_h + (_BA_h + h)*exp(dt*(-1.0/(ms + 20.0855369231877*ms*exp(-0.1*vm/mV)) - 0.07*exp(-0.05*vm/mV)/ms))
_BA_m = (2.5*mV/(-mV*ms + 12.1824939607035*mV*ms*exp(-0.1*vm/mV)) - 0.1*vm/(-mV*ms + 12.1824939607035*mV*ms*exp(-0.1*vm/mV)))/(-2.5*mV/(-mV*ms + 12.1824939607035*mV*ms*exp(-0.1*vm/mV)) + 0.1*vm/(-mV*ms + 12.1824939607035*mV*ms*exp(-0.1*vm/mV)) - 4*exp(-0.0556*vm/mV)/ms)
_m = -_BA_m + (_BA_m + m)*exp(dt*(-2.5*mV/(-mV*ms + 12.1824939607035*mV*ms*exp(-0.1*vm/mV)) + 0.1*vm/(-mV*ms + 12.1824939607035*mV*ms*exp(-0.1*vm/mV)) - 4*exp(-0.0556*vm/mV)/ms))
_BA_n = (0.1*mV/(e*mV*ms/e**(0.1*vm/mV) - mV*ms) - 0.01*vm/(e*mV*ms/e**(0.1*vm/mV) - mV*ms))/(-0.1*mV/(e*mV*ms/e**(0.1*vm/mV) - mV*ms) + 0.01*vm/(e*mV*ms/e**(0.1*vm/mV) - mV*ms) - 0.125/(e**(0.0125*vm/mV)*ms))
_n = -_BA_n + e**(dt*(-0.1*mV/(e*mV*ms/e**(0.1*vm/mV) - mV*ms) + 0.01*vm/(e*mV*ms/e**(0.1*vm/mV) - mV*ms) - 0.125/(e**(0.0125*vm/mV)*ms)))*(_BA_n + n)
_BA_vm = (EK*gK*n**4/C + ENa*gNa*h*m**3/C + El*gl/C + input_current(t, i)/C)/(-gK*n**4/C - gNa*h*m**3/C - gl/C)
_vm = -_BA_vm + (_BA_vm + vm)*exp(dt*(-gK*n**4/C - gNa*h*m**3/C - gl/C))
h = _h
m = _m
n = _n
vm = _vm

This will then be translated into the actual programming language code, which could be in Python (using numpy), Cython, or C++ (when using C++ standalone mode). If you actually want to step through the final generated code, then I think you have two options:

1. Force python code generation by setting the Python/numpy code generation target:

   from brian2 import prefs
   prefs.codegen.target = 'numpy'

   you can then use normal Python debugging tools (e.g. those of your editor).

2. You use the C++ standalone mode by setting

   from brian2 import set_device
   set_device('cpp_standalone', directory='debug', debug=True)

   which will create a C++ project in the "debug" folder. You can then (outside of Python) run the generated main binary with a C debugger, e.g. gdb.

I'll close the issue here, but feel free to add further comments if this does not help you. Or, as mentioned earlier, join us on the discussion forum.