To run the code below:
SHIFT+ENTER
on your keyboard or press the play button
() in the toolbar above.Feel free to create new cells using the plus button
(), or pressing SHIFT+ENTER
while this cell
is selected.
This is an interactive verison of the idealized model of the smooth pursuit reflex. This version does not explain the model itself, but shows how Brian's "runtime mode" can be used to interact with a running simulation. In this mode, the generated code based on the model descriptions is seemlessly integrated with the Python environment and can execute arbitrary Python code at any point during the simulation via a specially annotated function, called a "network operation".
For a non-interactive version of this example which generates the article's figure see this notebook.
This notebook is based on matplotlib and ipympl, which enables quick updates of the plot in real-time. For a version based on plotly (as the other, non-interactive examples), see this notebook.
# Needs ipywidgets and ipympl
%matplotlib widget
import ipywidgets as widgets
import threading
from brian2 import *
plt.ioff()
The model itself (mostly identical to the non-interactive example, except that some of the constants are included as parameters in the equation and can therefore change during the simulation):
alpha = (1/(50*ms))**2 # characteristic relaxation time is 50 ms
beta = 1/(50*ms) # friction parameter
eqs_eye = '''
dx/dt = velocity : 1
dvelocity/dt = alpha*(x0-x)-beta*velocity : 1/second
dx0/dt = -x0/tau_muscle : 1
dx_object/dt = (noise - x_object)/tau_object: 1
dnoise/dt = -noise/tau_object + tau_object**-0.5*xi : 1
tau_object : second
tau_muscle : second
'''
eye = NeuronGroup(1, model=eqs_eye, method='euler')
taum = 20*ms
motoneurons = NeuronGroup(2, model= 'dv/dt = -v/taum : 1', threshold = 'v>1',
reset = 'v=0', refractory = 5*ms, method='exact')
motosynapses = Synapses(motoneurons, eye, model = 'w : 1', on_pre = 'x0+=w')
motosynapses.connect() # connects all motoneurons to the eye
motosynapses.w = [-0.5,0.5]
N = 20
width = 2./N # width of receptive field
gain = 4.
eqs_retina = '''
I = gain*exp(-((x_object-x_eye-x_neuron)/width)**2) : 1
x_neuron : 1 (constant)
x_object : 1 (linked) # position of the object
x_eye : 1 (linked) # position of the eye
dv/dt = (I-(1+gs)*v)/taum : 1
gs : 1 # total synaptic conductance
'''
retina = NeuronGroup(N, model = eqs_retina, threshold = 'v>1', reset = 'v=0', method='exact')
retina.v = 'rand()'
retina.x_eye = linked_var(eye, 'x')
retina.x_object = linked_var(eye, 'x_object')
retina.x_neuron = '-1.0 + 2.0*i/(N-1)'
sensorimotor_synapses = Synapses(retina, motoneurons, model = 'w : 1 (constant)', on_pre = 'v+=w')
sensorimotor_synapses.connect(j = 'int(x_neuron_pre > 0)')
sensorimotor_synapses.w = '20*abs(x_neuron_pre)/N_pre'
M = StateMonitor(eye, ('x', 'x0', 'x_object'), record = True)
S_retina = SpikeMonitor(retina)
S_motoneurons = SpikeMonitor(motoneurons)
We create an empty plot that will be updated during the run:
# Plot preparation
fig, (ax_spikes, ax_position) = plt.subplots(2, 1, gridspec_kw={'height_ratios': (2, 1)}, sharex=True)
h_retina = ax_spikes.plot([], [], '|k', markeredgecolor='k', label='retina')[0]
h_left = ax_spikes.plot([], [], '|', color='C0', markeredgecolor='C0', label='left motoneuron')[0]
h_right = ax_spikes.plot([], [], '|', color='C1', markeredgecolor='C1', label='right motoneuron')[0]
ax_spikes.set(yticks=[], ylabel='neuron index', xticks=[], xlim=(0, 10), ylim=(0, 22))
ax_spikes.spines['bottom'].set_visible(False)
ax_position.axhline(0, color='gray')
h_eye = ax_position.plot([], [], 'k', label='eye')[0]
h_object = ax_position.plot([], [], color='C2', label='object')[0]
ax_position.set(yticks=[-1, 1], yticklabels=['left', 'right'], xlabel='time (s)',
xticks=np.arange(11, 2), xticklabels=np.arange(11, 2)-10,
xlim=(0, 10), ylim=(-1, 1))
ax_position.legend(loc='upper right', bbox_to_anchor=(1.0, 2.0));
We now create interactive widgets that the user can use to start/stop the simulation, as well as for setting certain simulation parameters.
time_label = widgets.Label(value='Time: 0 s')
start_stop_button = widgets.Button(tooltip='Start simulation', icon='play')
tau_obj_slider = widgets.FloatSlider(orientation='horizontal', description='tau_object',
value=500, min=100, max=1000)
tau_muscle_slider = widgets.FloatSlider(orientation='horizontal', description='tau_muscle',
value=20, min=5, max=100)
weight_slider = widgets.FloatSlider(orientation='horizontal', description='w_muscle',
value=0.5, min=0, max=2)
sliders = widgets.VBox([widgets.HBox([time_label, start_stop_button]),
tau_obj_slider, tau_muscle_slider, weight_slider])
layout = widgets.HBox([fig.canvas, sliders])
We interact with the running simulation via a "network operation", a Python function that will be regularly called by Brian during the simulation run (here, every 100ms of biological time). This function can access arbitrary attributes of the model to get or set their values. We use this here to 1) update the plot with the data from the last second and 2) set parameters of the model to the values requested by the user.
should_stop = False
@network_operation(dt=100*ms)
def plot_output(t):
cutoff = (t - 10*second)
# Plot the data of the last 10 seconds
indices = S_retina.t > cutoff
h_retina.set_data((S_retina.t[indices] - cutoff)/second, S_retina.i[indices])
motoneuron_trains = S_motoneurons.spike_trains()
to_plot = motoneuron_trains[0][motoneuron_trains[0] > cutoff]
h_left.set_data((to_plot - cutoff)/second, np.ones(len(to_plot))*N)
to_plot = motoneuron_trains[1][motoneuron_trains[1] > cutoff]
h_right.set_data((to_plot - cutoff)/second, np.ones(len(to_plot))*(N+1))
indices = M.t > cutoff
h_eye.set_data((M.t[indices] - cutoff)/second, M.x[0][indices])
h_object.set_data((M.t[indices] - cutoff)/second, M.x_object[0][indices])
fig.canvas.draw_idle()
time_label.value = 'Time: {:.1f}s'.format(float(t[:]))
# Set the simulation parameters according to user settings
eye.tau_object = tau_obj_slider.value*ms
eye.tau_muscle = tau_muscle_slider.value*ms
motosynapses.w = [-weight_slider.value, weight_slider.value]
if should_stop:
net.stop()
We store the model and the "network operation" in a Network
object, and store its current state to allow for repeated execution.
net = Network(collect())
net.store()
We now define two helper functions used to start/stop simulations. The actual simulation will be run in a background thread so that the user interface stays reactive while the simulation is running:
def do_run(runtime):
net.restore()
net.run(runtime)
running = False
def button_pressed(b):
global running
global should_stop
if running:
should_stop = True
running = False
start_stop_button.tooltip = 'Start simulation'
start_stop_button.icon = 'play'
else:
should_stop = False
running = True
time_label.value = 'starting...'
start_stop_button.tooltip = 'Stop simulation'
start_stop_button.icon = 'stop'
thread = threading.Thread(target=do_run, args=(100*second, ))
thread.start()
start_stop_button.on_click(button_pressed)
We are now ready to display the plot and user interface, which can then be used to start the simulation and interact with the simulation parameters:
display(layout)