# 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=10
second, 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 = 3
second,
proba=0.2
)

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

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

ii_synapses = get_synapses(
“ii_synapses”,
inh_neurons,
inh_neurons,
J=-0.057,
tau_1 = 1second, # ???
tau_2 = 2
second,
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.

# Full traceback of error (if relevant)

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

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

Good catch, @rth , 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