# Zero division problem

**URL:** <https://brian.discourse.group/t/zero-division-problem/893>\
**Category:** Support\
**Tags:** synapses, equations\
**Created:** [30 January 2023 15:02 UTC](https://brian.discourse.group/t/zero-division-problem/893 "2023-01-30T15:02:41Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Arina\_Belova](https://avatars.discourse-cdn.com/v4/letter/a/4bbf92/32.png) [@Arina\_Belova](https://brian.discourse.group/u/Arina_Belova)\
**Post date:** [30 January 2023 15:02 UTC](https://brian.discourse.group/t/zero-division-problem/893/1 "2023-01-30T15:02:41Z")

</div>

# 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=15_second, 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 = 1_second,  
tau\_2 = 3_second,  
proba=0.2  
)

ei\_synapses = get\_synapses(  
“ei\_synapses”,  
exc\_neurons,  
inh\_neurons,  
J=-0.045,  
tau\_1 = 1_second, # ???  
tau\_2 = 2_second,  
proba=0.5  
)

ie\_synapses = get\_synapses(  
“ie\_synapses”,  
inh\_neurons,  
exc\_neurons,  
J=0.014,  
tau\_1 = 1_second, # ???  
tau\_2 = 3_second,  
proba=0.5  
)

ii\_synapses = get\_synapses(  
“ii\_synapses”,  
inh\_neurons,  
inh\_neurons,  
J=-0.057,  
tau\_1 = 1_second, # ???  
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.

# Expected output (if relevant)

# Actual output (if relevant)

# Full traceback of error (if relevant)

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [31 January 2023 17:40 UTC](https://brian.discourse.group/t/zero-division-problem/893/2 "2023-01-31T17:40:13Z")

</div>

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](https://brian.discourse.group/t/zerodivisionerror-float-division/519/2)

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [31 January 2023 23:30 UTC](https://brian.discourse.group/t/zero-division-problem/893/3 "2023-01-31T23:30:28Z")

</div>

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…

```python
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

```python
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`

```python
synapses.J = J
synapses.tau_1 = tau_1
synapses.tau_2 = tau_2

```

Does it cause the error?

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [1 February 2023 17:08 UTC](https://brian.discourse.group/t/zero-division-problem/893/4 "2023-02-01T17:08:52Z")

</div>

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:

```auto
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!