Changing stimuli programmatically in C++ standalone

Following up the discussion with @DavidKing2020 in User-provided function in run_regularly, here’s a little example showing that you can include rather complex programming logic in C++ standalone mode, with only basic Brian functions (i.e. without writing C++ code yourself). In the below example (full code here), I use a TimedArray to store 10 stimuli for 5 neurons, each 100ms long (each of the stimuli here is just a constant rate for each of the neurons, but it could also be time-varying).

# Each "input" is one combination of rates across the 5 neurons
n_stimuli = 10
rates = TimedArray(np.random.rand(n_stimuli, 5)*100*Hz, dt=100*ms)
stimulus_length = 100*ms

The group that transforms the rates into spikes is a NeuronGroup, which gives us more flexibility than a PoissonGroup. It does not simply replay the stimulus over time, but instead uses a start_t variable (set to t at the beginning of each block) and a stimulus_idx to shift the time so that it reads the stimulus from the correct block. It also multiplies the rate with a factor, so that we can scale it up if needed.

max_repetitions = 5
input_spikes = NeuronGroup(5, '''
                   rate = stimulus_strength * rates(t - start_t + stimulus_length*stimulus_idx, i) : second**-1
                   start_t : second  # Start time of the current trial
                   stimulus_idx : integer  # Index of the stimulus to show
                   stimulus_strength : 1  # Factor to scale up stimulus
                   repetitions : integer  # Number of times the stim has been presented''',
                           threshold='rand()<rate*dt')
input_spikes.stimulus_strength = 1

The target neurons are very simple lif neurons with an exponential synaptic current (connected 1–to–1 with the input neurons, but this could be different):

tau = 10*ms; tau_syn = 5*ms
target_neurons = NeuronGroup(5, '''dv/dt = (-v + I)/tau : 1
                                   dI/dt = -I/tau_syn : 1  # synaptic current''',
                             threshold='v > 1', reset='v = 0')

connections = Synapses(input_spikes, target_neurons, on_pre='I += 1')
connections.connect(j='i')

Now, we’ll have to be a bit creative. Our criterion for a “good trial” is that we have a certain number of spikes across all neurons for this block. To make this number easily accessible, we create a dummy group that stores this value and increase the spike counter for each spike via a synapse. We also reset the counter after each block:

# A dummy "synapse" and a dummy group. Used to keep track of the total number of spikes in
# the target group
spike_counter = NeuronGroup(1, 'spike_counter : integer')
counter_synapse = Synapses(target_neurons, spike_counter, on_pre='spike_counter += 1')
counter_synapse.connect()
# Reset the spike counter every 100ms
counter_synapse.run_regularly('spike_counter = 0', when='end', dt=stimulus_length)

Finally, we can put everything together via a synapse that connects this spike counter to the input neurons, and changes the stimulus_idx and stimulus_strength according to our rule: here, we want to switch to the next stimulus if at least 5 spikes were recorded or if we have exceeded the maximum number of trials. If this is not the case, we keep the same stimulus, and multiply its strength by 1.5. Since we still don’t have proper support for some kind of if operator in Brian code, we’ll have to work around by multiplying with 0s or 1s:

 Another "synapse", that uses the information in the spike counter to update the stimulus
# and or 
stimulus_control = Synapses(spike_counter, input_spikes)
stimulus_control.connect()
stimulus_control.run_regularly('''
switch_stimulus = spike_counter_pre > 5 or repetitions_post >= max_repetitions

repetitions_post = int(not switch_stimulus)*(repetitions_post + 1)
stimulus_idx_post += int(switch_stimulus)
stimulus_strength_post = 1*int(switch_stimulus) + 1.5*stimulus_strength_post*int(not switch_stimulus)

start_t = t
stimulus_strength_post = int(not stimulus_idx_post >= n_stimuli)*stimulus_strength_post
''', dt=stimulus_length, when='before_end')

The model is now complete and should do what we want. Unfortunately, we cannot dynamically change the runtime of a simulation in C++ standalone mode, so we have to run for the maximum possible time:

run(max_repetitions*n_stimuli*stimulus_length)

This should not be that bad in practice, since we switch off the stimulus by setting its strength to 0 if we go beyond the last stimulus. Simulating the unnecessary time should therefore be relatively quick.

Here’s a little post that shows that everything seems to be working as expected (note that stimulus_idx continues to increase in the final period, but this does not affect anything).

3 Likes

I made a the following changes to the above example:

  input_spikes = NeuronGroup(5, '''
                     rate: second**-1
                     start_t : second  # Start time of the current trial
                     stimulus_idx : integer  # Index of the stimulus to show
                     stimulus_strength : 1  # Factor to scale up stimulus
                     repetitions : integer  # Number of times the stim has been presented''',
                             threshold='rand()<rate*dt')
  input_spikes.stimulus_strength = 1
  input_spikes.rate='stimulus_strength * rates(t - start_t + stimulus_length*stimulus_idx, i)'

The result is somehow different from the original output though I expected it to be the same

image

Hi @DavidKing2020 . When you set

input_spikes.rate='stimulus_strength * rates(t - start_t + stimulus_length*stimulus_idx, i)'

you perform an initialization at the beginning of the simulation (i.e. using t=0ms, etc.). The value of input_spikes.rate will not change over time. You can compare this to a simulation where you set
the neuron’s membrane potential to a random value at the start of the simulation, e.g.:

neuron.v = 'v_reset + rand()*(v_threshold - v_reset)

You wouldn’t expect that this expression gets evaluated repeatedly for each time step during the simulation. The same mechanism is at play in your example, the expression you give for the rate gets only evaluated once.

Thanks Marcel. I’m trying to add the following neuron and synapse to control the amount of idle time of input spikes.

    silence_spikes = SpikeGeneratorGroup(1, [0], [stimulus_length - 10*ms], period=stimulus_length)
    silence_input_syn = Synapses(silence_spikes, input_spikes, on_pre='rate_post = 0*Hz')
    silence_input_syn.connect()
    silence_input_syn.delay = 0*ms

The error is:

SyntaxError: Illegal line 'rate_post = 0*Hz' in abstract code. Cannot write to subexpression 'rate_post'.

The addition works fine with the attribute:

    rate: second**-1

but fails with

 rate = stimulus_strength * rates(t - start_t + stimulus_length*stimulus_idx, i) : second**-1

Ok, now I get what you are trying to do. Some explanation first: if you write rate = ..., you define what we call a “subexpression”. This basically means that wherever you write rate, it will be substituted by its right-hand-side. So in my above example, you could also directly write threshold='rand()<stimulus_strength * rates(t - start_t + stimulus_length*stimulus_idx, i)*dt' instead of defining rate. A subexpression is not something you can give a value. The semantics of this would be unclear: if you define rate as a function of time but then set it to 0Hz, should it stay at that value forever? I think in your case, you rather want to silence the last 10ms of every stimulus, but it should become active again for the next stimulus. There are various ways to do this:

  1. You could directly include the “silent time” into the array that you wrap in the TimedArray.
  2. You can multiply the rate by 0 during the silent time, by using something like
rate = int(t - start_t < stimulus_length - 10*ms) * stimulus_strength * rates(...) : second**-1

This will multiply the definition of the rate by 1, except during the last 10ms where it will be multiplied by 0.
3. In some cases, something closer to the approach you tried could also be useful. You could add a new variable and refer to it in the rate definition:

is_quiet : boolean
rate = int(not is_quiet) * stimulus_strength * rates : second**-1

You could then target the is_quiet variable by a synapse and switch it between True and False when appropriate.

1 Like

Thanks. That’s solved the problem.