Video visualization#

In this tutorial we’ll see how to inject a video visualizer into a simulation.

Preamble:

from brian2 import *
import matplotlib.pyplot as plt

from cleo import *

utilities.style_plots_for_docs()

# numpy faster than cython for lightweight example
prefs.codegen.target = 'numpy'
# for reproducibility
np.random.seed(1866)

c_exc = 'xkcd:tomato'
c_inh = 'xkcd:cerulean blue'

Set up the simulation#

Network#

We’ll use excitatory and inhibitory populations of exponential integrate-and-fire neurons.

n_e = 400
n_i = n_e // 4
def eif(n, name):
    ng = NeuronGroup(
        n,
        """
        dv/dt = (-(v - E_L) + Delta_T*exp((v-theta)/Delta_T) + Rm*I) / tau_m : volt
        I : amp
        """,
        threshold="v>30*mV",
        reset="v=-55*mV",
        namespace={
            "tau_m": 20 * ms,
            "Rm": 500 * Mohm,
            "theta": -50 * mV,
            "Delta_T": 2 * mV,
            "E_L": -70*mV,
        },
        name=name,
    )
    ng.v = -70 * mV
    return ng

exc = eif(n_e, "exc")
inh = eif(n_i, "inh")
W = 250
p_S = 0.3
S_ei = Synapses(exc, inh, on_pre="v_post+=W*mV/n_e")
S_ei.connect(p=p_S)
S_ie = Synapses(inh, exc, on_pre="v_post-=W*mV/n_i")
S_ie.connect(p=p_S)
S_ee = Synapses(exc, exc, on_pre="v_post+=W*mV/n_e")
S_ee.connect(condition='abs(i-j)<=20')

mon_e = SpikeMonitor(exc)
mon_i = SpikeMonitor(inh)

net = Network(exc, inh, S_ei, S_ie, S_ee, mon_e, mon_i)

Coordinates and optogenetics#

Here we configure the coordinates and optogenetic stimulation. For more details, see the “Optogenetic stimulation” tutorial. Note that we save the arguments used in the plotting function for reuse later on when generating the video.

from cleo.coords import assign_coords_uniform_cylinder
from cleo.viz import plot

r = 1
assign_coords_uniform_cylinder(
    exc, xyz_start=(0, 0, 0.3), xyz_end=(0, 0, 0.4), radius=r
)
assign_coords_uniform_cylinder(
    inh, xyz_start=(0, 0, 0.3), xyz_end=(0, 0, 0.4), radius=r
)

from cleo.opto import chr2_4s, fiber473nm, Light

opsin = chr2_4s()
fibers = Light(
    coords=[[0, 0, 0], [700, 0, 0]] * umeter,
    light_model=fiber473nm(),
    max_Irr0_mW_per_mm2=30,
    save_history=True,
)

plotargs = {
    "colors": [c_exc, c_inh],
    "zlim": (0, 600),
    "scatterargs": {"s": 20},  # to adjust neuron marker size
    "axis_scale_unit": umeter,
}

plot(
    exc,
    inh,
    **plotargs,
    devices=[fibers],
)

(<Figure size 640x480 with 1 Axes>,
 <Axes3D: xlabel='x (um)', ylabel='y (um)', zlabel='z (um)'>)
../_images/video_visualization_5_1.png

Simulator, optogenetics injection#

Here we create the simulator and inject the OptogeneticIntervention.

sim = CLSimulator(net)
sim.inject(opsin, exc, Iopto_var_name='I')
sim.inject(fibers, exc)

CLSimulator(io_processor=None, devices={Light(brian_objects=set(), sim=..., name='Light', value=array([0., 0.]), save_history=True, light_model=OpticFiber(R0=100. * umetre, NAfib=0.37, wavelength=0.473 * umetre, K=125. * metre ** -1, S=7370. * metre ** -1, ntis=1.36), coords=array([[0. , 0. , 0. ],
       [0.7, 0. , 0. ]]) * mmetre, direction=array([0., 0., 1.]), max_Irr0_mW_per_mm2=30, max_Irr0_mW_per_mm2_viz=None, default_value=array([0., 0.])), FourStateOpsin(brian_objects={NeuronGroup(clock=Clock(dt=100. * usecond, name='defaultclock'), when=start, order=0, name='light_agg_ChR2_exc'), Synapses(clock=Clock(dt=100. * usecond, name='defaultclock'), when=start, order=0, name='opto_syn_ChR2_exc')}, sim=..., name='ChR2', action_spectrum=[(400, 0.34), (422, 0.65), (460, 0.96), (470, 1), (473, 1), (500, 0.57), (520, 0.22), (540, 0.06), (560, 0.01)], required_vars=[('Iopto', amp), ('v', volt)], g0=114. * nsiemens, gamma=0.00742, phim=2.33e+23 * (second ** -1) / (meter ** 2), k1=4.15 * khertz, k2=0.868 * khertz, p=0.833, Gf0=37.3 * hertz, kf=58.1 * hertz, Gb0=16.1 * hertz, kb=63. * hertz, q=1.94, Gd1=105. * hertz, Gd2=13.8 * hertz, Gr0=0.33 * hertz, E=0. * volt, v0=43. * mvolt, v1=17.1 * mvolt, model="\n        dC1/dt = Gd1*O1 + Gr0*C2 - Ga1*C1 : 1 (clock-driven)\n        dO1/dt = Ga1*C1 + Gb*O2 - (Gd1+Gf)*O1 : 1 (clock-driven)\n        dO2/dt = Ga2*C2 + Gf*O1 - (Gd2+Gb)*O2 : 1 (clock-driven)\n        C2 = 1 - C1 - O1 - O2 : 1\n\n        Theta = int(phi_pre > 0*phi_pre) : 1\n        Hp = Theta * phi_pre**p/(phi_pre**p + phim**p) : 1\n        Ga1 = k1*Hp : hertz\n        Ga2 = k2*Hp : hertz\n        Hq = Theta * phi_pre**q/(phi_pre**q + phim**q) : 1\n        Gf = kf*Hq + Gf0 : hertz\n        Gb = kb*Hq + Gb0 : hertz\n\n        fphi = O1 + gamma*O2 : 1\n        # TODO: get this voltage dependence right \n        # v1/v0 when v-E == 0 via l'Hopital's rule\n        # fv = (1 - exp(-(V_VAR_NAME_post-E)/v0)) / -2 : 1\n        fv = f_unless_x0(\n            (1 - exp(-(V_VAR_NAME_post-E)/v0)) / ((V_VAR_NAME_post-E)/v1),\n            V_VAR_NAME_post - E,\n            v1/v0\n        ) : 1\n\n        IOPTO_VAR_NAME_post = -g0*fphi*fv*(V_VAR_NAME_post-E)*rho_rel : ampere (summed)\n        rho_rel : 1", extra_namespace={'f_unless_x0': <brian2.core.functions.Function object at 0x7fd1dc9ab4f0>})})

Processor#

And we set up open-loop optogenetic stimulation:

from cleo.ioproc import LatencyIOProcessor

fibers.update([5, 0])
class OpenLoopOpto(LatencyIOProcessor):
    def process(self, state_dict, time_ms):
        # random walk stimulation
        fiber_intensities = fibers.value + np.random.randn(2)*.5
        fiber_intensities[fiber_intensities < 0] = 0
        if time_ms > 50:
            fiber_intensities = [0, 5]
        return ({"Light": fiber_intensities}, time_ms)

sim.set_io_processor(OpenLoopOpto(sample_period_ms=1))

CLSimulator(io_processor=<__main__.OpenLoopOpto object at 0x7fd1dc4be470>, devices={Light(brian_objects=set(), sim=..., name='Light', value=array([5, 0]), save_history=True, light_model=OpticFiber(R0=100. * umetre, NAfib=0.37, wavelength=0.473 * umetre, K=125. * metre ** -1, S=7370. * metre ** -1, ntis=1.36), coords=array([[0. , 0. , 0. ],
       [0.7, 0. , 0. ]]) * mmetre, direction=array([0., 0., 1.]), max_Irr0_mW_per_mm2=30, max_Irr0_mW_per_mm2_viz=None, default_value=array([0., 0.])), FourStateOpsin(brian_objects={NeuronGroup(clock=Clock(dt=100. * usecond, name='defaultclock'), when=start, order=0, name='light_agg_ChR2_exc'), Synapses(clock=Clock(dt=100. * usecond, name='defaultclock'), when=start, order=0, name='opto_syn_ChR2_exc')}, sim=..., name='ChR2', action_spectrum=[(400, 0.34), (422, 0.65), (460, 0.96), (470, 1), (473, 1), (500, 0.57), (520, 0.22), (540, 0.06), (560, 0.01)], required_vars=[('Iopto', amp), ('v', volt)], g0=114. * nsiemens, gamma=0.00742, phim=2.33e+23 * (second ** -1) / (meter ** 2), k1=4.15 * khertz, k2=0.868 * khertz, p=0.833, Gf0=37.3 * hertz, kf=58.1 * hertz, Gb0=16.1 * hertz, kb=63. * hertz, q=1.94, Gd1=105. * hertz, Gd2=13.8 * hertz, Gr0=0.33 * hertz, E=0. * volt, v0=43. * mvolt, v1=17.1 * mvolt, model="\n        dC1/dt = Gd1*O1 + Gr0*C2 - Ga1*C1 : 1 (clock-driven)\n        dO1/dt = Ga1*C1 + Gb*O2 - (Gd1+Gf)*O1 : 1 (clock-driven)\n        dO2/dt = Ga2*C2 + Gf*O1 - (Gd2+Gb)*O2 : 1 (clock-driven)\n        C2 = 1 - C1 - O1 - O2 : 1\n\n        Theta = int(phi_pre > 0*phi_pre) : 1\n        Hp = Theta * phi_pre**p/(phi_pre**p + phim**p) : 1\n        Ga1 = k1*Hp : hertz\n        Ga2 = k2*Hp : hertz\n        Hq = Theta * phi_pre**q/(phi_pre**q + phim**q) : 1\n        Gf = kf*Hq + Gf0 : hertz\n        Gb = kb*Hq + Gb0 : hertz\n\n        fphi = O1 + gamma*O2 : 1\n        # TODO: get this voltage dependence right \n        # v1/v0 when v-E == 0 via l'Hopital's rule\n        # fv = (1 - exp(-(V_VAR_NAME_post-E)/v0)) / -2 : 1\n        fv = f_unless_x0(\n            (1 - exp(-(V_VAR_NAME_post-E)/v0)) / ((V_VAR_NAME_post-E)/v1),\n            V_VAR_NAME_post - E,\n            v1/v0\n        ) : 1\n\n        IOPTO_VAR_NAME_post = -g0*fphi*fv*(V_VAR_NAME_post-E)*rho_rel : ampere (summed)\n        rho_rel : 1", extra_namespace={'f_unless_x0': <brian2.core.functions.Function object at 0x7fd1dc9ab4f0>})})

Inject VideoVisualizer#

A VideoVisualizer is an InterfaceDevice like recorders and stimulators and needs to be injected in order to properly interact with the Brian network. Keep in mind the following:

  • It must be injected after all other devices for the devices='all' argument to work as expected.

  • Similarly to recording and stimulation, you must specify the target neuron groups (to display, in this case) on injection

  • The dt argument makes a huge difference on the amount of time it takes to generate the video. You may want to keep this high while experimenting and only lower it when you are ready to generate a high-quality video since the process is so slow.

from cleo.viz import VideoVisualizer

vv = VideoVisualizer(dt=1 * ms, devices="all")
sim.inject(vv, exc, inh)

CLSimulator(io_processor=<__main__.OpenLoopOpto object at 0x7fd1dc4be470>, devices={VideoVisualizer(brian_objects={<SpikeMonitor, recording from 'spikemonitor_2'>, <SpikeMonitor, recording from 'spikemonitor_3'>}, sim=..., name='VideoVisualizer', devices=[Light(brian_objects=set(), sim=..., name='Light', value=array([5, 0]), save_history=True, light_model=OpticFiber(R0=100. * umetre, NAfib=0.37, wavelength=0.473 * umetre, K=125. * metre ** -1, S=7370. * metre ** -1, ntis=1.36), coords=array([[0. , 0. , 0. ],
       [0.7, 0. , 0. ]]) * mmetre, direction=array([0., 0., 1.]), max_Irr0_mW_per_mm2=30, max_Irr0_mW_per_mm2_viz=None, default_value=array([0., 0.]))], dt=1. * msecond, fig=None, ax=None, neuron_groups=[NeuronGroup(clock=Clock(dt=100. * usecond, name='defaultclock'), when=start, order=0, name='exc'), NeuronGroup(clock=Clock(dt=100. * usecond, name='defaultclock'), when=start, order=0, name='inh')], _spike_mons=[<SpikeMonitor, recording from 'spikemonitor_2'>, <SpikeMonitor, recording from 'spikemonitor_3'>], _num_old_spikes=[0, 0], _value_per_device_per_frame=[], _i_spikes_per_ng_per_frame=[]), Light(brian_objects=set(), sim=..., name='Light', value=array([5, 0]), save_history=True, light_model=OpticFiber(R0=100. * umetre, NAfib=0.37, wavelength=0.473 * umetre, K=125. * metre ** -1, S=7370. * metre ** -1, ntis=1.36), coords=array([[0. , 0. , 0. ],
       [0.7, 0. , 0. ]]) * mmetre, direction=array([0., 0., 1.]), max_Irr0_mW_per_mm2=30, max_Irr0_mW_per_mm2_viz=None, default_value=array([0., 0.])), FourStateOpsin(brian_objects={NeuronGroup(clock=Clock(dt=100. * usecond, name='defaultclock'), when=start, order=0, name='light_agg_ChR2_exc'), Synapses(clock=Clock(dt=100. * usecond, name='defaultclock'), when=start, order=0, name='opto_syn_ChR2_exc')}, sim=..., name='ChR2', action_spectrum=[(400, 0.34), (422, 0.65), (460, 0.96), (470, 1), (473, 1), (500, 0.57), (520, 0.22), (540, 0.06), (560, 0.01)], required_vars=[('Iopto', amp), ('v', volt)], g0=114. * nsiemens, gamma=0.00742, phim=2.33e+23 * (second ** -1) / (meter ** 2), k1=4.15 * khertz, k2=0.868 * khertz, p=0.833, Gf0=37.3 * hertz, kf=58.1 * hertz, Gb0=16.1 * hertz, kb=63. * hertz, q=1.94, Gd1=105. * hertz, Gd2=13.8 * hertz, Gr0=0.33 * hertz, E=0. * volt, v0=43. * mvolt, v1=17.1 * mvolt, model="\n        dC1/dt = Gd1*O1 + Gr0*C2 - Ga1*C1 : 1 (clock-driven)\n        dO1/dt = Ga1*C1 + Gb*O2 - (Gd1+Gf)*O1 : 1 (clock-driven)\n        dO2/dt = Ga2*C2 + Gf*O1 - (Gd2+Gb)*O2 : 1 (clock-driven)\n        C2 = 1 - C1 - O1 - O2 : 1\n\n        Theta = int(phi_pre > 0*phi_pre) : 1\n        Hp = Theta * phi_pre**p/(phi_pre**p + phim**p) : 1\n        Ga1 = k1*Hp : hertz\n        Ga2 = k2*Hp : hertz\n        Hq = Theta * phi_pre**q/(phi_pre**q + phim**q) : 1\n        Gf = kf*Hq + Gf0 : hertz\n        Gb = kb*Hq + Gb0 : hertz\n\n        fphi = O1 + gamma*O2 : 1\n        # TODO: get this voltage dependence right \n        # v1/v0 when v-E == 0 via l'Hopital's rule\n        # fv = (1 - exp(-(V_VAR_NAME_post-E)/v0)) / -2 : 1\n        fv = f_unless_x0(\n            (1 - exp(-(V_VAR_NAME_post-E)/v0)) / ((V_VAR_NAME_post-E)/v1),\n            V_VAR_NAME_post - E,\n            v1/v0\n        ) : 1\n\n        IOPTO_VAR_NAME_post = -g0*fphi*fv*(V_VAR_NAME_post-E)*rho_rel : ampere (summed)\n        rho_rel : 1", extra_namespace={'f_unless_x0': <brian2.core.functions.Function object at 0x7fd1dc9ab4f0>})})

Run simulation and visualize#

Here we display a quick plot before generating the video:

T = 100
sim.run(T * ms)

WARNING    'T' is an internal variable of group 'light_prop_ChR2_exc', but also exists in the run namespace with the value 100. The internal variable will be used. [brian2.groups.group.Group.resolve.resolution_conflict]
INFO       No numerical integration method specified for group 'exc', using method 'euler' (took 0.01s, trying other methods took 0.04s). [brian2.stateupdaters.base.method_choice]
INFO       No numerical integration method specified for group 'inh', using method 'euler' (took 0.01s, trying other methods took 0.02s). [brian2.stateupdaters.base.method_choice]

INFO       No numerical integration method specified for group 'opto_syn_ChR2_exc', using method 'euler' (took 0.02s, trying other methods took 0.08s). [brian2.stateupdaters.base.method_choice]

fig, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True)
stim_vals = np.array(fibers.values)
sptexc = mon_e.spike_trains()
ax1.eventplot([t/ms for t in sptexc.values()], lineoffsets=list(sptexc.keys()), color=c_exc)
ax1.set(ylabel="neuron index", title="exc spiking")
sptinh = mon_i.spike_trains()
ax2.eventplot([t/ms for t in sptinh.values()], lineoffsets=list(sptinh.keys()), color=c_inh)
ax2.set(ylabel="neuron index", title="inh spiking")
ax3.plot(fibers.t_ms, stim_vals[:, 0], c="#72a5f2", lw=2, label='fiber 1')
ax3.plot(fibers.t_ms, stim_vals[:, 1], c="#72c5f2", lw=2, label='fiber 2')
ax3.legend()
ax3.set(ylabel=r"$Irr_0$ (mm/mW$^2$)", title="optogenetic stimulus", xlabel="time (ms)");
../_images/video_visualization_14_0.png 
fibers.t_ms[0:5]

[0.0, 0.0, 1.0, 2.0, 3.0]

The VideoVisualizer stores the data it needs during the simulation, but hasn’t yet produced any visual output. We first use the generate_Animation(), plugging in the arguments we used for the original plot.

Also, we set the max_Irr0_mW_per_mm2_viz attribute of the optogenetic intervention. This effectively scales how bright the light appears in the visualization. That is, a high maximum irradiance makes the stimulus values small in comparison and produces a faint light, while a low ceiling makes the values relatively large and produces a bright light in the resulting video.

fibers.max_Irr0_mW_per_mm2_viz = np.max(stim_vals)
ani = vv.generate_Animation(plotargs, slowdown_factor=10)
../_images/video_visualization_17_0.png

The generate_Animation() function returns a matplotlib FuncAnimation object, which you can then use however you want. You will probably want to save a video.

Note that at this point the video still hasn’t been rendered; that happens when you try and save or visualize the animation. This step takes a while if your temporal resolution is high, so we suggest you do this only after your experiment is finalized and after you’ve experimented with low framerate videos to finalize video parameters.

Here we embed the video using HTML so you can see the output:

from matplotlib import rc
rc('animation', html='jshtml')

ani