Implement Brunel network with conductance-based synapses

Support
,
1

Description of problem

I am trying to implement a lIF network with conductance-based synapses, supposed to resemble a similar lIF network with current-based synapses as e.g. done in Meffin, et al. (2004), and Cavallari et al. (2014).
To be able to reproduce the Meffin network, I tried adapting the Brunel network (as implemented in the Brian2 examples) to conductance-based synapses. However, I can’t quite seem te grasp how to do this properly because of the connection between V and the conductance parameter.

I realise that there have been a lot of questions surrounding this topic already, but none of them seem to help me solve my issue in this specific model (in fact, the conversion from current-based to conductance-based seems easier in the models with more complex synapses as then there is a conductance/current parameter explicitly defined and updated.)
So I am hoping that someone in this community can help me.

Minimal code to reproduce problem

I added the following lines to the code from the example (with only parameters set B):

## Meffin to COBA method
Vbar = (theta+V_r)/2.
Ve_ = (0+60) * mV
Vi_ = (-80+60) * mV
be = J
bi = -g*J
ae = be / (Vbar-Ve_)
ai = bi / (Vbar-Vi_)

What you have aready tried

I tried several things (not being convinced either is the correct way to do this and none seem to produce correct activity ).

  1. Changing the effect of the synapses on_pre="v += ae*(Ve_-v_post)" and on_pre="v += ai*(Vi_-v_post)"
    and the Poisson input to P = PoissonGroup(C_ext, rates=nu_ext); SynP = Synapses(P, neurons, on_pre="v += ae*(Ve_ - v_post)"); SynP.connect()
    → No activity at all
  2. Changing the neuron and synapse equations to the following (but then there is something wrong with the units and as expected I get the DimensionMismatchError: Inconsistent units in differential equation defining variable 'v' :
neurons = NeuronGroup(N,
"""dv/dt = -v/tau - ainh*(v-Vi_) - aex*(v-Ve_): volt (unless refractory)
  ainh : 1
  aex : 1
  """,  threshold="v > theta", reset="v = V_r", refractory=tau_rp,  method="rk2")
exc_synapses = Synapses(excitatory_neurons, target=neurons, on_pre="aex += ae", delay=D)
inhib_synapses = Synapses(inhibitory_neurons, target=neurons, on_pre="ainh += ai", delay=D)
external_poisson_input = PoissonInput(target=neurons, target_var="aex", N=C_ext, rate=nu_ext, weight=ae)
  1. Although I cannot seem to find a reference in the paper by Meffin to a good time constant to use to solve the units problem I randomly tried dividing the whole expression for dv/dt by tau, but that just leads to plumetting voltage values and accompanying warnings
`WARNING    (string):19: RuntimeWarning: overflow encountered in multiply
 [py.warnings]
WARNING    (string):20: RuntimeWarning: invalid value encountered in subtract
 [py.warnings]
WARNING    (string):20: RuntimeWarning: invalid value encountered in multiply
 [py.warnings]
WARNING    (string):20: RuntimeWarning: overflow encountered in multiply
 [py.warnings]
WARNING    (string):19: RuntimeWarning: invalid value encountered in multiply
 [py.warnings]
WARNING    neurongroup's variable 'v' has NaN, very large values, or encountered an error in numerical integration. This is usually a sign that an unstable or invalid integration method was chosen. [brian2.groups.group.invalid_values]
2

Hi @ArankaS. Unfortunately, I do not have the time to reply to this in much detail right now, but hopefully I can give you some pointers (and/or someone else from the community can step in). One source of misunderstanding is that the Brunel network is not actually using current-based synapses (there is no synaptic current in the model), but so called “delta synapses”, where each spike directly makes the post-synaptic membrane potential jump up by a certain amount. For a current-based model, you’d have a synaptic current in the equations (dv/dt = ... + I_syn... : volt), which decays back to zero in the absence of spikes (dI_syn/dt = -I_syn / tau_syn : amp), and gets increased by every spike (in the synapse: on_pre="I_syn_post += ..."). For a conductance-based model, the synaptic current depends on the difference between the current membrane potential and the reversal potential of the synapse (I_syn = g_syn*(E_syn - v) : amp), the synaptic conductance decays back to 0 in the absence of spikes (dg_syn/dt = -g_syn/tau_syn : siemens), and each spike increases the synaptic conductance (in the synapse: on_pre="g_syn_post += ...").

Pages 17–19 of this presentation illustrate this a bit: brian-material/2022-TD-Brian-Sorbonne/neural_simulation_BIP.pdf at d1d74fc78c04c703a99730c1edaab9f677552a1e · brian-team/brian-material · GitHub

For an example of implementations of the three types of synapses, see the end of this Jupyter Notebook: Jupyter Notebook Viewer

Hope that gets you going!

4

Hi @mstimberg,
Thank you for the clarification, and the useful material! I had been reading up on synapse types (but had somehow missed the ‘delta-synapse’ concept) but the jupyter notebook gives a nice overview of the effects of different synapse types (and a nice sanbox to play around with)!

In the paper by Meffin et al. delta-synapses are being used, but they define both a current (I_\text{syn}^\text{current} = \sum_{k=e,i,x} c b_k \sum_n\sum_m \delta(t-t_n^m-D)), with reference to Brunel (2000) and a conductance type(I_\text{syn}^\text{conductance} = \sum_{k=e,i,x} c a_k (V-V_k) \sum_n\sum_m \delta(t-t_n^m-D)).

So I am assuming a correct interpretation of their conductance type synapse is as delta pulses with a variable pulse strength (dependent on the current membrane voltage).
In which case my point (1) of the things I had tried might be the correct solution.
I tested this on the example in the Jupiter notebook, and this gives me delta-synapses with a behaviour more like the conductance-based synapses (whereas the delta pulses in the notebook behave similar to the current-based synapses).

afbeelding
afbeelding997×251 8.58 KB

Based on this, I might feel safe assuming that my issue lies somewhere else than in the synapse implementation. Because I am still getting zero-activity when trying to apply the correspondence of Meffin, et al. to the Brunel implementation in Brian.

On the note of correspondences between models, I was wondering if there is a method of determining the correct update values for a current/conductance based network based on the delta pulses. In the notebook all three weights for the different synapses seem to be chosen to have corresponding behaviour. The formula of Meffin, et al. gives me approximately the update values for g_i, g_e based on I_e, I_i, is there some similar relation you can apply to the update value of v_post to obtain I_e, I_i?

5

Hi again. I had a quick look at the Meffin et al. paper, and I think your interpretation is correct. The terminology is a bit fuzzy in this area, what I called “delta synapses” is sometimes simply called current-based synapses, mostly to distinguish it from conductance-based models, and this corresponds to their definition. I have to say I haven’t seen something like their definition of “conductance-based” synapses before, but I agree that I’d implement it as a jump in v_post, dependent on the current membrane potential.

I am not 100% sure whether I understand your question correctly, but if you are asking about having parameters that give you the same behaviour for different type of synapse models, then the short answer is that there is no single “equivalent” value (for the notebook, I think I chose some values by hand, but not quite sure anymore). In some situations, it is useful to chose values that give you the same maximum value of the membrane potential after a spike. This will give equivalent behaviour with regard to some questions, e.g. how much do I have to scale up a synapse so that a single spike leads to a spike, or how many spikes have to arrive simultaneously to make a neuron spike that was at rest before? But, in particular with slow synapses, this does make less and less sense, since they sum up over much longer times than in a model with delta synapses – so just looking at the maximum does not tell the whole story.

6

Hi, thank you for the additional information!

In the mean time I think I found the issues with my code that caused the v_post-dependent jumps not to work.
On the one hand I made a stupid sign-error in calculating the values for a_. They should be calculated as follows:

be = J
bi = -g*J
ae = be / (Ve_ - Vbar)
ai = bi / (Vi_ - Vbar)

And secondly, my replacement for the Poisson-input was sending the same Poisson trains to each neuron, causing some weird perfectly synchroneous behaviour.
Making the function PoissonInput (which does produce independent input for each neuron, as I understand from the Brian2-documentation) dependent on v_post does not seem possible.
I have currently implemented this Poisson input as follows:

P = PoissonGroup((N_E+N_I)*C_ext, rates=nu_ext)
SynP = Synapses(P, neurons, on_pre="v += ae*(Ve_ - v_post)") 
SynP.connect(i='j*C_ext+k for k in range(0,C_ext)')

This leads to a huge increase in runtime though. So I was wondering if there is a more efficient way to implement this?

7

Hi, indeed having a big PoissonGroup is quite expensive, and it is doing a lot of unnecessary computation if all that ever happens is summing up all these spikes. I don’t quite remember, but it might actually be that in the Brunel paper they did not use P = PoissonGroup((N_E+N_I)*C_ext, rates=nu_ext) but rather P = PoissonGroup(N_E+N_I, rates=C_ext*nu_ext) which is very similar (and much faster) for low rates – it breaks down for high rates, though, since you cannot have more than one spike per time step.

All that said, I think you should be able to use PoissonInput with a small workaround. Instead of directly updating v_post, you’d include a variable like exc_spikes : integer in the neuron equations and use a combination of PoissonInput and a run_regularly operation:

# Update exc_spikes to number of spikes coming in this time step
poisson_inp = PoissonInput(neurons, "exc_spikes", C_ext, nu_ext,
                           when='before_synapses') 
# Update v based on the number of spikes
neurons.run_regularly("v += exc_spikes * ae * (Ve_ - v)", when="synapses")