Zero division problem

Support
,
1

Description of problem

Trying to implement uniform network structure as in "Slow dynamics and high variability in balanced cortical networks with clustered connections " by Litwin-Kumar and Doiron we got zero division error in implementing the differential equations for synapses.
The full code for reproducibility of the problem is below:

Minimal code to reproduce problem

def get_population(name, N, tau_mem, mu):
“”“Get population of neurons
name – name of population
N – number of neurons
tau_mem – time constant of neuron’s membrane
mu – ???
“””

v_thresh = 1 # units? 
v_reset = 0 

neurons = NeuronGroup(
    N,
    """
    tau_mem : second
    mu : 1
    v_reset : 1 
    v_thresh : 1
    I_syn_ee_synapses : 1 / second
    I_syn_ei_synapses : 1 / second
    I_syn_ie_synapses : 1 / second
    I_syn_ii_synapses : 1 / second
    dv/dt = (mu - v)/tau_mem + (I_syn_ee_synapses +
                          I_syn_ei_synapses +
                          I_syn_ie_synapses +
                          I_syn_ii_synapses) : 1 (unless refractory)
    """,
    threshold="v > v_thresh",
    reset="v = v_reset",
    refractory=5*second,
    method="exact",
    name=name,
)

neurons.mu = mu
neurons.tau_mem = tau_mem # different for exc and inh neurons 15 and 10 ms 
neurons.v_thresh = v_thresh
neurons.v_reset = v_reset

return neurons

def get_synapses(name, source, target, J, tau_1, tau_2, proba):
“”"Construct connections and retrieve synapses

name -- name of synapses
source -- source of connections
target -- target of connections
J -- weight of the synapse
tau_1 -- time constant of source group
tau_2 -- time constant of target group
proba -- probability of connection between groups
"""

# may need to multiply J by 15mV to dimensionise it
# for now dimensionless as I do not have a clue how units 
# work in Brian


synapses_eqs = """
J : 1
tau_1 : second
tau_2 : second

dF/dt = -F/tau_2 + x/(tau_1 * tau_2) : 1 / second (clock-driven)
dx/dt = -x/tau_1 : 1 (clock-driven)

I_syn_{}_post = F : 1 / second (summed)
""".format(name)

synapses_action = """x += J"""

synapses = Synapses(
    source,
    target,
    model=synapses_eqs,
    on_pre=synapses_action,
    method="exact",
    name=name,
)
synapses.connect(p = proba)

# N_syn = len(synapses)

# different values depending on types of pre/post synaptic neurons
# TODO: scale J appropriately to its neuron group size!!!!!
synapses.J = J

return synapses

start_scope()

configure neuron populations

exc_neurons = get_population(“exc_neurons”, N=4000, tau_mem=15second, mu=np.random.uniform(1.1, 1.2))
inh_neurons = get_population(“inh_neurons”, N=1000, tau_mem=10second, mu=np.random.uniform(1.0, 1.05))

configure synapses

ee_synapses = get_synapses(
“ee_synapses”,
exc_neurons,
exc_neurons,
J=0.024,
tau_1 = 1second,
tau_2 = 3second,
proba=0.2
)

ei_synapses = get_synapses(
“ei_synapses”,
exc_neurons,
inh_neurons,
J=-0.045,
tau_1 = 1second, # ???
tau_2 = 2second,
proba=0.5
)

ie_synapses = get_synapses(
“ie_synapses”,
inh_neurons,
exc_neurons,
J=0.014,
tau_1 = 1second, # ???
tau_2 = 3second,
proba=0.5
)

ii_synapses = get_synapses(
“ii_synapses”,
inh_neurons,
inh_neurons,
J=-0.057,
tau_1 = 1second, # ???
tau_2 = 2second,
proba=0.5
)

run for burnin time to settle network activity

defaultclock.dt = 0.1 * second
burnin = 500 * second
run(burnin)

record from now on

spike_monitor_exc = SpikeMonitor(exc_neurons)
spike_monitor_inh = SpikeMonitor(inh_neurons)
state_monitor_ee = StateMonitor(ee_synapses, [“x”], record=True)

duration = 4200 * second
run(duration, report=“text”)

What you have aready tried

We tried to switch from cyton execution to numpy one, but there is still a problem with division by 0 somewhere. The equations are implemented the best way possible in our opinion.

Expected output (if relevant)

Actual output (if relevant)

Full traceback of error (if relevant)

2

Hi @Arina_Belova. Unfortunately, I did not yet have time to look into this issue in more detail, but sometimes the analytical solutions (generated by using method='exact') can contain terms that create a division by zero for some parameter values. A typical example would be when you have two time constants (e.g. tau_1 and tau_2 in your example) with the same values – or of course when one of the time constants is zero. I don’t see anything obvious like this in your equations, though.

As a workaround, you can change the integration method to e.g. method='rk2', this will probably make the error go away.

For debugging, the technique I mention in this comment could be helpful: ZeroDivisionError: float division - #2 by mstimberg

3

Hi @Arina_Belova , it seems I spotted some problem, but I’m not sure, so try to look at this more carefully. When you call get_synapses function, variables tau_1 and tau_2 exist only in the local context. At the moment when brian2 makes actually synapses (just before the simulation), these variables are gone…

def get_synapses(name, source, target, J, tau_1, tau_2, proba):
“”"Construct connections and retrieve synapses

name -- name of synapses
source -- source of connections
target -- target of connections
J -- weight of the synapse
tau_1 -- time constant of source group
tau_2 -- time constant of target group
proba -- probability of connection between groups
"""
# may need to multiply J by 15mV to dimensionise it
# for now dimensionless as I do not have a clue how units 
# work in Brian


synapses_eqs = """
J : 1
tau_1 : second
tau_2 : second

dF/dt = -F/tau_2 + x/(tau_1 * tau_2) : 1 / second (clock-driven)
dx/dt = -x/tau_1 : 1 (clock-driven)

I_syn_{}_post = F : 1 / second (summed)
""".format(name)

synapses_action = """x += J"""

synapses = Synapses(
    source,
    target,
    model=synapses_eqs,
    on_pre=synapses_action,
    method="exact",
    name=name,
)
synapses.connect(p = proba)

# N_syn = len(synapses)

# different values depending on types of pre/post synaptic neurons
# TODO: scale J appropriately to its neuron group size!!!!!
synapses.J = J

return synapses

Try to bake them into equations

synapses_eqs = """
J : 1
tau_1={} : second
tau_2 ={}: second

dF/dt = -F/tau_2 + x/(tau_1 * tau_2) : 1 / second (clock-driven)
dx/dt = -x/tau_1 : 1 (clock-driven)

I_syn_{}_post = F : 1 / second (summed)
""".format(tau_1,tau_2,name)

OR
explicitly set them in the same way as you set J

synapses.J     = J
synapses.tau_1 = tau_1
synapses.tau_2 = tau_2

Does it cause the error?

1 Like
4

Good catch, @rth :clap: , this indeed looks like the issue. Even if tau_1 and tau_2 were “visible” at the point of the run() call, the tau_1 and tau_2 parameters defined in the equation would take preference (you’d get a warning about this, though). Parameters have the value 0 if not initialized manually, so this would indeed explain the division by zero error.

As an alternative to “baking them into the equations” (you could directly replace tau_1 and tau_2 by their respective values, actually), or initializing them via synapses.tau_1 = ..., you could also use the Synapses namespace argument:

synapses_eqs = """
J : 1

dF/dt = -F/tau_2 + x/(tau_1 * tau_2) : 1 / second (clock-driven)
dx/dt = -x/tau_1 : 1 (clock-driven)

I_syn_{}_post = F : 1 / second (summed)
""".format(name)

synapses_action = """x += J"""

synapses = Synapses(
    source,
    target,
    model=synapses_eqs,
    on_pre=synapses_action,
    method="exact",
    name=name,
    namespace={'tau_1': tau_1, 'tau_2': tau_2}
)

(the same would work for J of course).
A final comment: if you keep them in the equations as parameters (and assign them values via synapses.tau_1 = ... , etc.), consider adding the (constant) flag, and – if you use a single value for all synapses – the (shared) flag, i.e. tau_1 : second (constant, shared) in the equations. This doesn’t change anything fundamentally, but can enable Brian to do some optimizations and also wastes less memory.

Let us know if this works for you!

1 Like