Access postsynaptic variables in synapse equations

Support
1

Description of problem

I’m playing with voltage- and spike-time-dependent plasticity, aka Clopath et al. model. While LTD has a spike-train term with delta functions and can be moved to on_pre resetting, LTP has all continuous variables and should be computed in a clock-driven manner. Both LTP and LTD use linear rectificatopn for voltage and low-pass filtered voltage.

\frac{dw_{ij}}{dt} = -A^{-}S_j(t)\lfloor \bar{u}^-_i - \theta_{-}\rfloor_{+}+ A^+\bar{x}_j \lfloor \bar{u}^+_i - \theta_{-}\rfloor_{+} \lfloor v_i - \theta_{+}\rfloor_{+}
\begin{array}{rlrl} du^-_i/dt &=\left(v-u^-_i \right)/\tau^- & \tau^-&=10 \text{ms} \\ du^+_i/dt &=\left(v-u^+_i \right)/\tau^+ & \tau^+ &=7 \text{ms} \\ d\bar{x}_i/dt &=\left(S_i(t)-\bar{x}_i \right)/\tau^x & \tau^x &=15 \text{ms} \\ \end{array}

where w_{ij} is a synaptic conductance from j^{th} presynaptic neuron to i^{th} postsynaptic. v is a voltage and S(t) = \sum \delta(t-t') is a spike train.

To make code a little faster, I moved computations of low pass filtered voltages and ReLU functions to the neuron side. So, part of my neuron equations looks like that

nrn_equ="""
.... neuron dynamics ....
dvlp/dt = (v-vlp)/tlpf_p  : 1 # low-pass filter of the voltage for LTD
dvlm/dt = (v-vlm)/tlpf_m  : 1 # low-pass filter of the voltage for LTP

vup      = int(v>Tltp)*(v-Tltp)     : 1
vlpup    = int(vlp>Tltd)*(vlp-Tltd) : 1
vlmup    = int(vlm>Tltd)*(vlm-Tltd) : 1
"""

Synaptic model can be implemented like this

syn = Synapses(lgninput, neurons, model="""
dwsyn/dt = Altp*x*vup_post*vlpup_post/tausynscal: 1 (clock-driven) 
dx/dt = -x/tlpf_s    : 1 (clock-driven) # slow presynaptic FR 
dy/dt = -y/tfrest    : 1 (clock-driven) # low-pass filter for FR estimate
""",
on_pre='''
    x    += 1*ms/tlpf_s
    wsyn -= Altd*(y*second/tfrest/frRef)*vlmup_post
 .... synaptic dynamics .....
''',
on_post='''
    y    += 1*ms/tfrest''')

This code works and works well.

============================

However, to debug my model, I want to separate LTP and LTD components and record them. So I modified synaptic equations like this:

syn = Synapses(lgninput, neurons, model="""
dwsyn/dt = LTP/tausynscal: 1 (clock-driven) 
dx/dt = -x/tlpf_s               : 1 (clock-driven) # slow presynaptic FR 
dy/dt = -y/tfrest               : 1 (clock-driven) # low-pass filter for FR estimate
LTP=Altp*x*vup_post*vlpup_post  : 1
LTD                             : 1
""",
on_pre='''
    x    += 1*ms/tlpf_s
    LTD   = Altd*(y*second/tfrest/frRef)*vlmup_post
    wsyn -= LTD
 .... synaptic dynamics .....''',
on_post='''
    y    += 1*ms/tfrest''')

This code returns an error, which I don’t understand

ERROR      Brian 2 encountered an unexpected error. If you think this is a bug in Brian 2, please report this issue either to the discourse forum at <http://brian.discourse.group/>, or to the issue tracker at <https://github.com/brian-team/brian2/issues>. Please include this file with debug information in your report: /tmp/brian_debug_8n4izity.log  Additionally, you can also include a copy of the script that was run, available at: /tmp/brian_script_35ik7dpv.py Thanks! [brian2]
Traceback (most recent call last):
  File "PLS-cortex.py", line 448, in <module>
    onenrnpres_srec = StateMonitor(syn,"wsyn x y LTP LTD".split(),record=synids,dt=mth["/rec/cont/dt"]*ms)
  File "/home/rth/.local/lib/python3.8/site-packages/brian2/monitors/statemonitor.py", line 248, in __init__
    self.variables.add_reference(f'_source_{varname}',
  File "/home/rth/.local/lib/python3.8/site-packages/brian2/core/variables.py", line 1829, in add_reference
    self.add_referred_subexpression(name, group, var, index)
  File "/home/rth/.local/lib/python3.8/site-packages/brian2/core/variables.py", line 1768, in add_referred_subexpression
    self.add_referred_subexpression(new_name,
  File "/home/rth/.local/lib/python3.8/site-packages/brian2/core/variables.py", line 1773, in add_referred_subexpression
    self.add_reference(new_name,
  File "/home/rth/.local/lib/python3.8/site-packages/brian2/core/variables.py", line 1821, in add_reference
    raise TypeError(f"Cannot link variable '{name}' to '{varname}' in "
TypeError: Cannot link variable '___source_LTP_synapses_vup_post_synapses__vup_post_neurongroup_v' to '_vup_post_neurongroup_v' in group 'synapses' -- need to precalculate direct indices but index _postsynaptic_idx can change

It may be related to standalone OpenMP mode, I used to load all cores of my processor.

set_device('cpp_standalone')
prefs.devices.cpp_standalone.openmp_threads = npth
2

Hi!

I guess the culprit is setting only some synapses to be recorded in conjunction with standalone mode:

record=synids

Can you try to record all synapses (only as a check)?

Cheers,

Sebastian

3

Hi Sebastian, thank you for looking into it!
It doesn’t seem that the problem in recording, but in caches somewhere.
I add this line of code:

os.system("rm -fR tmp __pycache__ /home/rth/.cython/brian_extensions/_cython_magic_* "+standalone_dir)

and the problem is gone.
However, it’s a bit odd, so maybe something else is going on here.

4

Hi Sebastian,
After 1/2 hour investigation I should admit that you were right, and I overlooked a bug in my script :beetle:.

It is related to synaptic recording.

1 Like
5

Hi everyone,

huh, I don’t quite see how the cache plays a role here… The problem is a very technical one, and one that we could certainly handle better, but it is entangled deeply with the code generation machinery. Here’s a (minimal ?) example that generates the same error:

G = NeuronGroup(2, '''v : 1
                      x = 2*v : 1''')
S = Synapses(G, G, '''w : 1
                      y = w + x_post : 1''')
S.connect()
mon = StateMonitor(S, 'y', record=[0])

This fails with:

TypeError: Cannot link variable '___source_y_synapses_x_post_synapses__x_post_neurongroup_v' to '_x_post_neurongroup_v' in group 'synapses' -- need to precalculate direct indices but index _postsynaptic_idx can change

The problem is due to the use of the subexpression x from the postsynaptic cell as part of a subexpression in the synapse. Brian can handle these things, and its approach is using indices. E.g. here, we need the value of v[postsynaptic_idcs[record_idx]] when we loop over record_idx to record all the values we ask for. To make our code machinery not too complex, we try to limit this a bit, and pre-calculate the direct indices that we need. However, this errors out since it is worried about postysnaptic_idcs being “dynamic” – the StateMonitor wants to calculate these indices when it is initialized, but in particular in standalone mode this might not be possible since the postsynaptic indices have not been determined yet. Again, just to be clear, this is not handled well here and there are many situations (including yours, I think) where we could actually make this work. It is just not supported by the way things are done currently.

All that said, there are at least two workarounds (I’ll demonstrate them with my simple example, but hopefully they are easy to transfer).

  1. Do not use subexpressions, but record things individually and reconstruct them afterwards:
# ...group definitions
S_mon = StateMonitor(S, 'w', record=[0])
G_mon = StateMonitor(S, 'x', record=S.j[0])
# Result is  `S_mon.w + G_mon.x`
  1. Use (constant over dt) for the subexpressions
G = NeuronGroup(2, '''v : 1
                      x = 2*v : 1 (constant over dt)''')
S = Synapses(G, G, '''w : 1
                      y = w + x_post : 1 (constant over dt)''')

Using constant over dt means that these are no longer subexpressions, but e.g. x = 2*v is replaced by x: 1 together with a run_regularly operation that sets x = 2*v at the beginning of the time step. When such a variable is referred in a differential equation, this is not exactly equivalent (e.g. rk4 will only use one value instead of several slightly different values), but in practice this shouldn’t be an issue if the variable changes slowly with respect to the simulation time step.

Actually, using record=True is not possible in standalone mode. The reason is similar to what I described above: StateMonitor wants to determine the indices of the synapses to record from, but in standalone mode the synapses do not exist yet.

1 Like
6

@mstimberg, yes, it was my fault. I tried to debug and overlooked some flags in my script.

I see how it may work. Do you mean

G_mon = StateMonitor(G, 'x', record=S.j[0])
7

Sorry, yes. Alternatively (and more reasonable if you would otherwise record the same values several times), you could use

G_mon = StateMonitor(G, 'x', record=True)

and use the indices when you reconstruct the values, i.e. something like

y_values = S_mon.w + G_mon.x[S.j[:], :]
1 Like