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).
