Disparity between connection probability equation and estimated connection probability

[oliverdibb]

Description of problem

Hi,

I’m using Brian to explore a model originally implemented in NEST, and I’m wondering if a subtle difference between connection probability interpretations might be contributing to some problems I’m having.

See the connection probability equations from the original paper:

Note - here both excitatory and inhibitory neurons are arranged together randomly (uniformly) in a grid (with a fixed distance between neurons along each axis), so each neuron has three grid coordinates
(x, y, z).

Parameter values determining connection probabilities for each group are shown in table 1, below:

The connection probability equation looks to be in exactly the form required for straightforward use with Brian. However, when taking the parameter values from the above table, estimates of connection probability taken from the generated synaptic connectivity matrix are many times higher than desired. This is using an estimate function that correctly estimates average connectivities of synapses generated with p=c for constant c.

Initially I thought little of it and just tuned the parameters manually to produce the desired connection probabilities, but I’d like to remove any unnecessary departures from the original model, in case they cause problems down the road.

I’ve come back to this issue a number of times and have read relevant documentation of both simulators fairly extensively. However, I haven’t seen anything to suggest a clear difference in connection probability interpretation between the two simulators.

Any help would be greatly appreciated!

Best wishes,
Oli

[mstimberg]

Hi @oliverdibb. Could you give a link to the original paper (or state some more details, in particular how many neurons are distributed over the grid with what distances between them)? And ideally some minimal Brian code (no need for the correct neuron model, etc.) that shows unexpected results?

At the first sight, I don’t see anything that should be interpreted differently here. I would have expected something like

syn.connect(p='c*exp(-sqrt((x_post-x_pre)**2 + (y_post-y_pre)**2 + (z_post - z_pre)**2)/lambda)')`

to reproduce the connection pattern.

[oliverdibb]

Hi @mstimberg,

Thank you for your quick reply!

The original paper is here: STDP Forms Associations between Memory Traces in Networks of Spiking Neurons - PMC
Relevant details are on pages 3 and 4.

The neurons are arranged in a 6 x 6 x 15 grid (432 excitatory, 108 inhibitory, 540 total), and the positions are unit free. Unfortunately, the connection probability equation you gave is exactly the equation I’m using in my model!

I’ll get working on some minimal code to reproduce the connections - thought I’d send the paper through in the meantime in case it would help.

Thanks again for the help!

[schmitts]

Hi!

A quick implementation of the connection in Python also gave me higher than quoted connection probabilities (see output at the bottom).

Cheers,

Sebastian

#!/usr/bin/env python3
"""
https://brian.discourse.group/t/disparity-between-connection-probability-equation-and-estimated-connection-probability/667

Sebastian Schmitt, 2022
"""

import itertools

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

indices = np.arange(15*6*6)
# shuffle to randomly distribute inh/exc neurons
np.random.shuffle(indices)

# a column of 15x6x6 neurons
x = indices % 6
y = (indices // 6) % 6
z = indices // (6*6)

"""
# to check that grid is correct
fig = plt.figure()
ax = fig.add_subplot(projection="3d")

ax.scatter(x, y, z)
plt.show()
"""

pos = np.column_stack([x, y, z])
N_exc = 432
pos_exc = pos[:N_exc]
pos_inh = pos[N_exc:]

parameters = {"E->I" : {"c" : 2e5,
                        "l" : 0.25,
                        "pos_pre" : pos_exc,
                        "pos_post" : pos_inh},
              "I->E" : {"c" : 4e5,
                        "l" : 0.25,
                        "pos_pre" : pos_inh,
                        "pos_post" : pos_exc},
              "I->I" : {"c" : 1e5,
                        "l" : 0.25,
                        "pos_pre" : pos_inh,
                        "pos_post" : pos_inh}
              }

for connection_type, params in parameters.items():

    n_connections = 0
    n_pairs = 0

    c = params["c"]
    l = params["l"]

    pos_pre = params["pos_pre"]
    pos_post = params["pos_post"]

    for (x_pre, y_pre, z_pre), (x_post, y_post, z_post) in itertools.product(pos_pre, pos_post):

        d = np.sqrt((x_post-x_pre)**2 + (y_post-y_pre)**2 + (z_post-z_pre)**2)

        p = c*np.exp(-d/l)

        if p > np.random.uniform():
            n_connections += 1

        n_pairs += 1

    print("Connection {}, Conn. prob.: {:.1%} ({}/{}) ".format(connection_type, n_connections / n_pairs, n_connections, n_pairs))

Connection E->I, Conn. prob.: 16.2% (7542/46656) 
Connection I->E, Conn. prob.: 18.1% (8461/46656) 
Connection I->I, Conn. prob.: 14.3% (1665/11664) 

[oliverdibb]

Hi all,

@schmitts - thanks for your reply! You beat me to it, and I appreciate you taking the time.

@mstimberg, here is the (hopefully) runnable code that should demonstrate the issue, and I’ve pasted the output below that.

def avg_conn_prob(syn):
    source = syn.source
    target = syn.target
    pres = syn.i

    nos, vals = np.histogram(pres, bins=np.arange(0, len(source), 1))
    prob_array = nos / len(target)
    avg_prob = np.mean(prob_array)

    return avg_prob
from brian2 import *

x_lim = 6
y_lim = 6
z_lim = 15
n_inh = 108
n_exc = 432
n_tot = n_exc + n_inh

poss_coordinates = tuple(
    [x, y, z]
    for x in range(0, x_lim)
    for y in range(0, y_lim)
    for z in range(0, z_lim)
)

import random
gen = random.Random()
neuron_coordinates = gen.sample(poss_coordinates, n_tot)

x_val = []
y_val = []
z_val = []
for coord_vect in neuron_coordinates:
    x_val.append(coord_vect[0])
    y_val.append(coord_vect[1])
    z_val.append(coord_vect[2])

exc_x = x_val[0:n_exc]
exc_y = y_val[0:n_exc]
exc_z = z_val[0:n_exc]
inh_x = x_val[n_exc:]
inh_y = y_val[n_exc:]
inh_z = z_val[n_exc:]

neuron_eqs = """
X: 1 (constant)
Y: 1 (constant)
Z: 1 (constant)
"""
exc = NeuronGroup(n_exc, neuron_eqs)
exc.X = exc_x
exc.Y = exc_y
exc.Z = exc_z
inh = NeuronGroup(n_inh, neuron_eqs)
inh.X = inh_x
inh.Y = inh_y
inh.Z = inh_z


param_dict = {'E_I': {'c': 2*10**5},
              'I_E': {'c': 4*10**5},
              'I_I': {'c': 1*10**5}}

general_conn_prob_eq = "{c}*exp(-(sqrt((X_pre - X_post)**2 + (Y_pre - Y_post)**2 + (Z_pre - Z_post)**2))/{l_scale})"
lda = 0.25

E_I = Synapses(exc, inh, '')
I_E = Synapses(inh, exc, '')
I_I = Synapses(inh, inh, '')

for syn, label, goal_prob in zip((E_I, I_E, I_I), ('E_I', 'I_E', 'I_I'), (0.04, 0.05, 0.04)):
    conn_eq = general_conn_prob_eq.format(c=param_dict[label]['c'], l_scale=lda)
    syn.connect(p=conn_eq)
    est_conn_prob = avg_conn_prob(syn)
    print(f'{label}: Observed connection prob. {est_conn_prob}, goal connection prob. {goal_prob}')

Output:

E_I: Observed connection prob. 0.1617470138351809, goal connection prob. 0.04
I_E: Observed connection prob. 0.18429820006922812, goal connection prob. 0.05
I_I: Observed connection prob. 0.14096573208722743, goal connection prob. 0.04

Thanks in advance!

[schmitts]

Hi @oliverdibb!

Now the case is pretty clear. Brian2 and the manual implementation agree. Either we have the same error in the interpretation of the formulas or there’s a typo or some other missing information in the paper.

Did you find the original source code online? If not, we could ask the authors for it.

Cheers,

Sebastian

[mstimberg]

I agree with @schmitts that it looks like Brian is interpreting things correctly and that there is either a bug in their code, on some information in the paper is missing or incorrect. E.g. if I divide the probability by 100 (i.e. let p(d) denote a probability in %), then the values are of a similar magnitude as those reported in the paper (around 5–6%). I don’t think we can resolve the mismatch without having access to their code.

[oliverdibb]

Hi both,

Thank you for your replies. I got in touch with the author and there is indeed a typo in the paper - \lambda = 0.15 was used, instead of the stated value \lambda = 0.25. Changing this corrects the erroneous mean probabilities.

I’ve also noticed that the refractory period is implemented differently to my version of the model. In theirs, the refractory period for each neuron is redrawn at each spike from a given gamma distribution, identical for every neuron in each group (for excitatory neurons, t \sim \Gamma(2, 5), for inhibitory neurons t \sim \Gamma(2, 1.5)). Initially, I implemented it as follows:

refractory = "(-log(1 - rand()) -log(1 - rand())) / rate" to draw from a gamma distribution with (shape, scale) = (2, 1/rate). This raises no errors, but testing the implementation raised some questions…

To test, I gave a NeuronGroup a variable with exactly the same string representation as refractory= so I could record it in a StateMonitor and confirm the average behaviour over the run. However, on running I get an error stating that rand() is a stateful function and hence can’t be used twice in a string. The same error was not raised when the same expression was used as the refractory parameter of the NeuronGroup. I assume though that the string implementation of the refractory period will still run into the same issue?

If this won’t work, I figured I could implement the refractory period as a variable and use a TimedArray refreshed every so often with a network operation. Does this sound along the right tracks, or am I missing something better?

Thanks again for your help.
Best wishes,
Oli

[mstimberg]

Glad you could figure this out!

@oliverdibb

To test, I gave a NeuronGroup a variable with exactly the same string representation as refractory= so I could record it in a StateMonitor and confirm the average behaviour over the run. However, on running I get an error stating that rand() is a stateful function and hence can’t be used twice in a string. The same error was not raised when the same expression was used as the refractory parameter of the NeuronGroup. I assume though that the string implementation of the refractory period will still run into the same issue?

If this won’t work, I figured I could implement the refractory period as a variable and use a TimedArray refreshed every so often with a network operation. Does this sound along the right tracks, or am I missing something better?

This is a tricky issue and currently not handled very well. The reason why we don’t allow two rand() calls in a single expression is that we process mathematical expressions with sympy, and it does not have a concept of stateful functions. It will therefore consider -log(rand()) - log(rand()) to be equivalent to -2*log(rand()). Now, this makes sense for expressions that are part of equations, but general statements (e.g. a reset) or expressions (like a threshold, or the refractory condition) are not treated “mathematically”, otherwise statements like u = u + d wouldn’t make any sense. But unfortunately the current code generation pipeline is a bit convoluted, and we end up converting “normal” expressions with sympy in places where we probably shouldn’t.

When you set refractory to your expression, this will actually not run into problems like the -2*log(rand()) “simplification” that I mentioned earlier, but it will not do what you expect it to do. The refractory expression is verified every time step, not only when a spike occurs (which is necessary for refractoriness conditions like v > v_threshold used e.g. in HH models, but is less obvious for expressions that evaluate to a time span). Long story short, the correct way to implement it would be to calculate it as part of the reset, i.e.:

• add a parameter like ref_time : second to the equations
• set refractory = 'ref_time'
• and add ref_time = (-log(rand() - log(rand()) ) / rate to your reset definition.

This will draw a new random time span for each spike and then use this as the refractory time.

Now, alas, this will run into the error about two rand() calls again… The workaround to deal with this is to “hide” that fact that you are using the same function twice. One way to do this is to define two new functions with different names, which both refer to the same function under the hood. This can be as easy as defining:

rand1 = DEFAULT_FUNCTIONS['rand']
rand2 = DEFAULT_FUNCTIONS['rand']

and then using ref_time = (-log(rand1() - log(rand2()) ) / rate in your reset.

To check whether this is doing what you think it is doing I’d do a more straightforward check: use a NeuronGroup that spikes every time step, except when it is refractory. Putting together what I said above this would be:

rand1 = DEFAULT_FUNCTIONS['rand']
rand2 = DEFAULT_FUNCTIONS['rand']
group = NeuronGroup(1, 'ref_time : second', threshold='True',
                    reset='ref_time = (-log(rand1() - log(rand2()) ) / rate',
                    refractory='ref_time')

You can record its spike with a SpikeMonitor, run it for a while and then make a histogram of the ISIs (to get the ISIs in ms, use e.g. np.diff(spike_mon.t/ms)).

Let us know whether this works as expected!