ZeroDivisionError when changing synaptic weights of large network with HH neurons

Description of problem

I aim to simulate a large network of 20.000 Hodgkin-Huxley neurons. This works, but when I change the synaptic weights to the desired value, I almost always get a ZeroDivionError: float division. This does not appear to occur in my code, but in something that brian does to process the code. I cannot seem to fix it. Can anyone help me?

Minimal code to reproduce problem

area = 300*umetre**2
Cm = (1*ufarad*cm**-2) * area
gl = (0.3*msiemens*cm**-2) * area

El = -65*mV
EK = -77*mV
ENa = 55*mV
g_na = (40*msiemens*cm**-2) * area
g_kd = (35*msiemens*cm**-2) * area
tau_ampa = 5*ms
tau_nmda = 150*ms
Eampa = 0.0*mV
Enmda = 0.0*mV

Mg2 = 1.

N_noise = 1 # number of inputs
N_rate = 5*Hz  # rate of each input
w_noise = 0.05*nS  #weight of the noise (PART OF PROBLEM)


N = 20000               #number of neurons
simtime = 1*second      #simulation time


eqs = Equations('''
dv/dt = (-gl*(v-El)-
         g_na*(m*m*m)*h*(v-ENa)-
         g_kd*(n*n*n*n)*(v-EK)-I_syn+I)/Cm  : volt
dm/dt = (xm-m)/tau_m : 1
dn/dt = (xn-n)/tau_n : 1
dh/dt = (xh-h)/tau_h : 1
I_syn = g_ampa*(v-Eampa) + g_nmda*(v-Enmda)/(1+Mg2*exp(-0.062*v/mV)/3.57) : amp
I :amp
dg_nmda/dt = -g_nmda*(1./tau_nmda) : siemens
dg_ampa/dt = -g_ampa*(1./tau_ampa) : siemens
alpha_m = 0.182*(v/mV+35)/(1-exp(-(v/mV+35)/9))/ms: Hz
beta_m = -0.124*(v/mV+35)/(1-exp((v/mV+35)/9))/ms: Hz
alpha_h = 0.25*exp(-(v/mV+90)/12)/ms : Hz
beta_h = 0.25*exp((v/mV+62)/6)/exp((v/mV+90)/12)/ms : Hz
alpha_n = 0.01*(v/mV-25)/(1-exp(-(v/mV-25)/9))/ms: Hz
beta_n = -0.002*(v/mV-25)/(1-exp((v/mV-25)/9))/ms: Hz
xm = alpha_m/(alpha_m+beta_m) :1
xn = alpha_n/(alpha_n+beta_n) :1
xh = alpha_h/(alpha_h+beta_h) :1
tau_m = 1/(alpha_m+beta_m) :second
tau_n = 1/(alpha_n+beta_n) :second
tau_h = 1/(alpha_h+beta_h) :second
''')



P = NeuronGroup(N, model=eqs, threshold='v>-20*mV', refractory=3*ms,
                method='exponential_euler')

Syn = Synapses(P, P, on_pre={'ampa': 'g_ampa+=0.02*nS',
                             'nmda': 'g_nmda+=0.02*nS'}) # make connections  (HERE THE WEIGHTS ARE ALSO THE PROBLEM)

Syn.connect(p=0.1)

Syn.ampa.delay = '(2.+randn())*ms'
Syn.nmda.delay = '(2.+randn())*ms'

P.v = -65*mV

Synnos = PoissonInput(P, 'g_ampa', N_noise, N_rate, w_noise, order=1) #synaptic noise 

run(simtime, report='text')

What you have aready tried

I’ve tried to reduce the number of neurons and number of connections. This only changes the values of the synaptic weights for which I get this error.

Expected output (if relevant)

Actual output (if relevant)

Full traceback of error (if relevant)

Traceback (most recent call last):
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "/home/Nina/Downloads/pycharm-2020.2.3/plugins/python/helpers/pydev/_pydev_bundle/pydev_umd.py", line 197, in runfile
    pydev_imports.execfile(filename, global_vars, local_vars)  # execute the script
  File "/home/Nina/Downloads/pycharm-2020.2.3/plugins/python/helpers/pydev/_pydev_imps/_pydev_execfile.py", line 18, in execfile
    exec(compile(contents+"\n", file, 'exec'), glob, loc)
  File "/home/Nina/PycharmProjects/pythonProject/main.py", line 89, in <module>
    run(simtime, report='text')
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/units/fundamentalunits.py", line 2392, in new_f
    result = f(*args, **kwds)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/core/magic.py", line 374, in run
    namespace=namespace, profile=profile, level=2+level)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/core/magic.py", line 232, in run
    namespace=namespace, profile=profile, level=level+1)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/core/base.py", line 280, in device_override_decorated_function
    return func(*args, **kwds)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/units/fundamentalunits.py", line 2392, in new_f
    result = f(*args, **kwds)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/core/network.py", line 1080, in run
    obj.run()
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/core/base.py", line 183, in run
    codeobj()
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/codegen/codeobject.py", line 102, in __call__
    return self.run()
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/codegen/runtime/cython_rt/cython_rt.py", line 153, in run
    return self.compiled_code.main(self.namespace)
  File "_cython_magic_aa3bcd64e641e8c37f7a0a182088943a.pyx", line 532, in _cython_magic_aa3bcd64e641e8c37f7a0a182088943a.main
ZeroDivisionError: float division

Hi. A division by zero error typically indicates a numerical problem during the integration of the differential equations. In general, it could be an inappropriate numerical integration method or a time step that is too big, but this does not seem to be the case here. I had a closer look at what happens when something starts to go wrong (the membrane potential of a single neuron goes to negative infinity), and in my run the issue came about when the membrane potential was very close to -35mV. The problem with the equations is that the equations for \alpha_m and \beta_m are not defined for a value of v exactly at -35mV. Hitting this value (and not, say, -35.1mV or -34.9mV) is not likely, but with many neurons getting a lot of input (and therefore spiking a lot) it will happen at some point. There is potentially the same issue for \alpha_n and \beta_n for the value 25mV. Since this \frac{\dots}{1-\exp(\dots)} expression is quite common in rate equations, we offer a special function exprel that can deal with this issue. It is defined as \mathrm{exprel}(x) = \frac{\exp(x) -1}{x}, except for being 1 for x=0.
You can use this to redefine the rate equations that use the problematic expression (also see the COBAHH example). E.g. for \alpha_m, your equation is equivalent to

alpha_m = 0.182*9/exprel((-(v/mV+35)/9))/ms : Hz

I was able to successfully run your network with this change. Since there are only excitatory synapses the network is firing like crazy, but I assume that your actual simulation will use inhibitory synapses as well (or weaker excitatory synapses).

Oh, one more thing I noticed:

Syn.ampa.delay = '(2.+randn())*ms'
Syn.nmda.delay = '(2.+randn())*ms'

This will lead to a number of negative delays, since the random number returned by randn() could well be below -2. I’m not 100% sure what Brian will do with these delays, but I’d rather avoid them. You can use the clip function to make sure they do not go below zero:

Syn.ampa.delay = 'clip(2.+randn(), 0, inf)*ms'
Syn.nmda.delay = 'clip(2.+randn(), 0, inf)*ms'

Thank you so much for the quick and accurate help!

I am actually trying to simulate signals from actual cultures of neuronal networks that only consist of excitatory upper layer cortical neurons (derived from hiPSCs). So there will be no inhibitory synapses. The in vitro networks actually don’t show crazy firing, but show really regular and slow network bursting behavior. So I’m curious to find out if I will eventually see this in the model when maybe adding plasticity or spike-frequency adaptation.

1 Like

I have a follow-up question, hope that is okay.

I’ve implemented the exprel function as you suggested. It seemed to work but I think that was just a coincidence. The voltage still blows up when it is close to -35 mV. However, the exprel function did not give an error. Now I am trying to find out where it goes wrong so I tried to record alpha_m with the StateMonitor and then I get the error:

''Error encountered with object named "statemonitor_1

  •       traceam = StateMonitor(P, 'alpha_m', record=True)
    

An error occurred when preparing an object. NotImplementedError: The “exprel” function needs C99 compiler support’’

I couldn’t figure out what this meant for the functioning of exprel. Does this mean it does not work? Or that I just cannot record it? Is there something I can do to fix this?

I hope the problem is clear. Thanks in advance!

Of course!

I am a bit confused by the error. The exprel function indeed needs C99 compiler support, but you should probably have a compatible compiler. From you earlier error message I assume that you are using Linux, is this still the case?
But the most confusing thing is that you get this error only when you record alpha_m – this would mean that when you are not recording it, your equations are not making use of alpha_m which is very odd. Could you give the equations you are using?

Of course! It is indeed very odd. I am testing with just one neuron which gets some external current with which I get the zero division error. The weird thing is, the gating variable m seems to go to -14, while it is perfectly between 0 and 1 before that. And the negative m then causes the voltage v to go to minus infinity. But the only reason I can think of that causes m to suddenly become negative, is when v is very negative, but that is only the case after m becomes negative. There should be no other possibility for m to become negative right? The minimal code to reproduce the error is as follows:

from brian2 import *
from brian2tools import *
import matplotlib.pyplot as plt
import numpy as np

#parameters 
area = 300*umetre**2  # doesn't matter because you multiply and divide by it
Cm = (1*ufarad*cm**-2) * area
gl = (0.3*msiemens*cm**-2) * area

El = -65*mV
EK = -77*mV
ENa = 55*mV
g_na = (40*msiemens*cm**-2) * area
g_kd = (35*msiemens*cm**-2) * area

simtime = 1*second      #simulation time

#The model
eqs = Equations('''
dv/dt = (-gl*(v-El)-g_na*(m*m*m)*h*(v-ENa)-g_kd*(n*n*n*n)*(v-EK)+I)/Cm  : volt
dm/dt = (xm-m)/tau_m : 1
dn/dt = (xn-n)/tau_n : 1
dh/dt = (xh-h)/tau_h : 1
I : amp
alpha_m = 0.182*9/exprel((-(v/mV+35)/9))/ms : Hz
beta_m = 0.124*9/exprel((v/mV+35)/9)/ms :Hz
alpha_h = 0.25*exp(-(v/mV+90)/12)/ms : Hz
beta_h = 0.25*exp((v/mV+62)/6)/exp((v/mV+90)/12)/ms : Hz
alpha_n = 0.01*9/exprel(-(v/mV-25)/9)/ms :Hz
beta_n = 0.002 *9/exprel((v/mV-25)/9)/ms :Hz
xm = alpha_m/(alpha_m+beta_m) :1
xn = alpha_n/(alpha_n+beta_n) :1
xh = alpha_h/(alpha_h+beta_h) :1
tau_m = 1/(alpha_m+beta_m) :second
tau_n = 1/(alpha_n+beta_n) :second
tau_h = 1/(alpha_h+beta_h) :second
''')

P = NeuronGroup(1, model=eqs, threshold='v>0*mV', refractory=1*ms, method='exponential_euler')

P.v = -65*mV
P.I[0] = 4e-12*amp

#Record variables
trace = StateMonitor(P, 'v', record=True)
#traceam = StateMonitor(P, 'alpha_m', record=True)
tracem = StateMonitor(P, 'm', record=True)
traceh = StateMonitor(P, 'h', record=True)
tracen = StateMonitor(P, 'n', record=True)

spikes = SpikeMonitor(P)

run(simtime, report='text')

I am indeed using Linux. The exprel function seems to work as expected in the command window by the way. Just not when I want to record alpha_m.
I hope this suffices.

Thanks for providing the code. Could it be that you are using an outdated version of Brian? We might have fixed a related bug at some point. The current version is 2.4.2 and if I run your code with it and plot v, I get the following which looks perfectly reasonable:
v_trace

PS: Not really important, but note that you can record several variables with a single monitor, e.g.

traces = StateMonitor(P, ['v', 'm', 'h', 'n'], record=True)
# Then use traces.v, traces.m, etc.

Oh I see that I am indeed using brian2 version 2.3. Would have never thought of that, thank you for the quick response again!

And thanks for the tip about the statemonitor, silly of me to not have done that.

1 Like

Unfortunately I again have a follow up question :upside_down_face:.

Updating brian worked on my windows laptop and the exact script like mentioned above (with 1 neuron) runs perfectly and I got the same output as you. However, I wish to work on a desktop with linux so I also updated brian2 there but there it doesn’t solve the problem. I am doing the exact thing I mentioned above which is exactly the same as on my windows laptop but I get the error:

ERROR      Brian 2 encountered an unexpected error. If you think this is a bug in Brian 2, please report this issue either to the discourse forum at <http://brian.discourse.group/>, or to the issue tracker at <https://github.com/brian-team/brian2/issues>. Please include this file with debug information in your report: /tmp/brian_debug_4i4_o4tl.log  Additionally, you can also include a copy of the script that was run, available at: /tmp/brian_script_qvhah2zk.py Thanks! [brian2]
Traceback (most recent call last):
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/core/network.py", line 897, in before_run
    obj.before_run(run_namespace)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/groups/group.py", line 1143, in before_run
    self.create_code_objects(run_namespace)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/groups/group.py", line 1136, in create_code_objects
    code_object = self.create_default_code_object(run_namespace)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/groups/group.py", line 1129, in create_default_code_object
    codeobj_class=self.codeobj_class
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/codegen/codeobject.py", line 450, in create_runner_codeobj
    compiler_kwds=compiler_kwds
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/devices/device.py", line 298, in code_object
    dtype=prefs['core.default_float_dtype'])
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/codegen/generators/base.py", line 263, in translate
    vector_statements)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/codegen/generators/cython_generator.py", line 208, in translate_statement_sequence
    kwds = self.determine_keywords()
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/codegen/generators/cython_generator.py", line 364, in determine_keywords
    user_func = self._add_user_function(varname, var, added)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/codegen/generators/cython_generator.py", line 220, in _add_user_function
    func_code = impl.get_code(self.owner)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/core/functions.py", line 280, in get_code
    self.availability_check()
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/codegen/cpp_prefs.py", line 279, in __call__
    'support'.format(self.name))
NotImplementedError: The "exprel" function needs C99 compiler support
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "/home/Nina/Downloads/pycharm-2020.2.3/plugins/python/helpers/pydev/_pydev_bundle/pydev_umd.py", line 197, in runfile
    pydev_imports.execfile(filename, global_vars, local_vars)  # execute the script
  File "/home/Nina/Downloads/pycharm-2020.2.3/plugins/python/helpers/pydev/_pydev_imps/_pydev_execfile.py", line 18, in execfile
    exec(compile(contents+"\n", file, 'exec'), glob, loc)
  File "/home/Nina/PycharmProjects/pythonProject/HH_single.py", line 141, in <module>
    run(simtime, report='text')
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/units/fundamentalunits.py", line 2434, in new_f
    result = f(*args, **kwds)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/core/magic.py", line 374, in run
    namespace=namespace, profile=profile, level=2+level)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/core/magic.py", line 232, in run
    namespace=namespace, profile=profile, level=level+1)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/core/base.py", line 278, in device_override_decorated_function
    return func(*args, **kwds)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/units/fundamentalunits.py", line 2434, in new_f
    result = f(*args, **kwds)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/core/network.py", line 1008, in run
    self.before_run(namespace)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/core/base.py", line 278, in device_override_decorated_function
    return func(*args, **kwds)
  File "/home/Nina/.conda/envs/brianenv/lib/python3.6/site-packages/brian2/core/network.py", line 899, in before_run
    raise BrianObjectException("An error occurred when preparing an object.", obj) from ex
brian2.core.base.BrianObjectException: Error encountered with object named "neurongroup_stateupdater".
Object was created here (most recent call only, full details in debug log):
  File "/home/Nina/PycharmProjects/pythonProject/HH_single.py", line 84, in <module>
    P = NeuronGroup(1, model=eqs, threshold='v>0*mV', refractory=2*ms, method='exponential_euler')
An error occurred when preparing an object. (See above for original error message and traceback.)

Exprel works perfectly fine if I use it in the command window. C99 compiler support should be fine. My laptop and the pc have all packages and versions the same. The difference is of course the system:
Laptop:
Python 3.6.12 |Anaconda, Inc.| (default, Sep 9 2020, 00:29:25) [MSC v.1916 64 bit (AMD64)]

desktop:
Python 3.6.12 |Anaconda, Inc.| (default, Sep 8 2020, 23:10:56) [GCC 7.3.0] on linux

And brian2genn is installed on the linux desktop like we discussed in the other topic, but I’m guessing that shouldn’t influence it.

Apologies for asking so much of you.

Is this also with Brian 2.4.2? There was a bug in earlier versions with respect to the detection of C99 support.

yes the python interpreter has brian2 2.4.2

Hmm, I might have an idea what is going wrong. I realized that our test for C99 support is not necessarily using the exact same compiler that will be used for actually compiling things, e.g. if you have one version of the compiler installed system-wide and one in your conda environment. Could you run the following two commands in a Linux terminal and post what they report (with the activated conda environment, so in your case after conda activate brianenv) ?

cc --version
$CC --version
(brianenv) [Nina@utwks62681 bin]$ cc --version
cc (GCC) 4.8.5 20150623 (Red Hat 4.8.5-39)
Copyright (C) 2015 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.


(brianenv) [Nina@utwks62681 bin]$ $CC --version
x86_64-conda_cos6-linux-gnu-cc (crosstool-NG 1.23.0.449-a04d0) 7.3.0
Copyright (C) 2017 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Ok, I think I understand now. We are using cc for the check of the compiler compatibilities, and on your machine that links to an old version of GCC which we identify as incompatible (we might be too strict here, actually). The version that will be used for actual compilation is the one given by $CC which is newer and has all the capabilities we need. This is something we need to fix in Brian itself, but this will take a moment.

Until then I see two options:

  1. you make cc on your machine link to a more recent gcc version. The disadvange of this is that this can potentially mess with other things on your system, so I would not recommend it.
  2. the other option is to add the following code to your script (somewhere near the beginning). It will make Brian skip the check to determine whether your compiler has C99 support or not and just pretend it has:
# Skip C99 support check
import brian2.codegen.cpp_prefs
brian2.codegen.cpp_prefs._compiler_supports_c99 = True

Okay! thank you so much once again.