Multi-channel, bidirectional optogenetics#

Cleo supports the simultaneous use of multiple Light devices, multiple channels per device, and multiple opsins per neuron group. Here’ we’ll see how to use all these features.

# boilerplate
%load_ext autoreload
%autoreload 2
from brian2 import *
import matplotlib.pyplot as plt

import cleo
from cleo import *

cleo.utilities.style_plots_for_docs()

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

# colors
c = {
    'main': '#C500CC',
    '473nm': '#72b5f2',
    '590nm': (1, .875, 0),
}

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

Network setup#

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

n = 500
ng = NeuronGroup(
    n,
    """
    dv/dt = (
        -(v - E_L)
        + Delta_T*exp((v-theta)/Delta_T)
        + Rm*(I_exc + I_inh + I_bg)
    ) / tau_m + sigma*sqrt(2/tau_m)*xi: volt
    I_exc : amp
    I_inh : amp
    """,
    threshold="v > -10*mV",
    reset="v=E_L",
    namespace={
        "tau_m": 20 * ms,
        "Rm": 500 * Mohm,
        "theta": -50 * mV,
        "Delta_T": 2 * mV,
        "E_L": -70*mV,
        "sigma": 5 * mV,
        "I_bg": 60 * pamp,
    },
)
ng.v = np.random.uniform(-70, -50, n) * mV

W = 1 * mV
p_S = 0.3
n_neighbors = 40
S_ee = Synapses(ng, model='w: 1', on_pre="v_post+=W*w/sqrt(N)")
S_ee.connect(condition='abs(i-j)<=n_neighbors and i!=j')
S_ee.w = np.exp(np.random.randn(int(S_ee.N - 0))) * np.random.choice([-1, 1], int(S_ee.N - 0))

spike_mon = SpikeMonitor(ng)
# for visualizing currents
Imon = StateMonitor(ng, ('I_exc', 'I_inh'), record=range(50))
net = Network(ng, S_ee, spike_mon, Imon)
sim = cleo.CLSimulator(net)

from cleo.coords import assign_coords_grid_rect_prism
xmax_mm = 2
assign_coords_grid_rect_prism(ng, xlim=(0, xmax_mm), ylim=(-.15, .15), zlim=(0, .3), shape=(20, 5, 5))
cleo.viz.plot(ng, colors=[c['main']])

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

Injecting a multi-channel Light#

A Light device can have multiple channels; all the user needs is to specify the coordinates (and optionally direction) of each light source (channel). A LightModel (e.g., that of an optical fiber) defines how light propagates from each source.

Here we inject 590 nm light for activating Vf-Chrimson. Lacking a more rigorous quantification, we assume absorption and scattering coefficients of 590 nm light in the brain are roughly 0.8 times that of 470 nm light (see Jacques 2013 for some justification).

from cleo.light import Light, OpticFiber
from cleo.opto import vfchrimson_4s

n_fibers = 4
coords = np.zeros((n_fibers, 3))
end_space = 1 / (2 * n_fibers)
coords[:, 0] = np.linspace(end_space, 1 - end_space, n_fibers) * xmax_mm
coords[:, 1] = -0.025
amber_fibers = Light(
    name="amber fibers",
    coords=coords * mm,
    light_model=OpticFiber(K=0.125 * 0.8 / mm, S=7.37 * 0.8 / mm),
    wavelength=590 * nmeter,
)
sim.inject(amber_fibers, ng)
vfc = vfchrimson_4s()
sim.inject(vfc, ng, Iopto_var_name="I_exc")
cleo.viz.plot(ng, colors=[c["main"]], sim=sim)

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

Bidirectional control via a second opsin#

Here we will demonstrate increasing and decreasing activity in the same experiment by injecting an inhibitory opsin with a minimally overlapping activation spectrum. Alternatively, we could achieve bidirectional control with excitatory opsins on excitatory and inhibitory neurons.

We will use the anion channel GtACR2 which is maximally activated by 470 nm light.

Let’s first visualize how the action spectra for the two opsins overlap:

gtacr2 = cleo.opto.gtacr2_4s()
cleo.light.plot_spectra(vfc, gtacr2)

(<Figure size 640x480 with 1 Axes>,
 <Axes: title={'center': 'Action/excitation spectra'}, xlabel='λ (nm)', ylabel='ε'>)
../_images/multi_opto_8_1.png

We see that GtACR2 will be totally unaffected by amber (590 nm) light, while Vf-Chrimson will be slightly activated by blue (470 nm) light.

coords[:, 1] = 0.025
blue_fibers = Light(
    name="blue fibers",
    coords=coords * mm,
    light_model=cleo.light.fiber473nm(),
)
sim.inject(blue_fibers, ng)
sim.inject(gtacr2, ng, Iopto_var_name="I_inh")
cleo.viz.plot(ng, colors=[c["main"]], sim=sim)

WARNING    /home/kyle/Dropbox (GaTech)/projects/cleo/cleo/light/light_dependence.py:107: UserWarning: λ = 590.0 nm is outside the range of the action spectrum data for GtACR2. Assuming ε = 0.
  warnings.warn(
 [py.warnings]

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

Open-loop stimulation#

We will now design a stimulus pattern to demonstrate bidirectional control segregated by channel.

from cleo.ioproc import LatencyIOProcessor

class OpenLoopOpto(LatencyIOProcessor):
    def __init__(self):
        super().__init__(sample_period_ms=1)

    # since this is open-loop, we don't use state_dict
    def process(self, state_dict, time_ms):
        amplitude_mW_mm2 = 1
        time_offsets = np.array([0, -20, -40, -60])
        t = time_ms + time_offsets
        amber = ((t >= 20) & (t < 40)) * amplitude_mW_mm2
        blue = ((t >= 60) & (t < 63)) * amplitude_mW_mm2
        
        # return output dict and time
        return ({"amber fibers": amber, "blue fibers": blue}, time_ms)

sim.set_io_processor(OpenLoopOpto())

CLSimulator(io_processor=<__main__.OpenLoopOpto object at 0x7fe758a361d0>, devices={BansalFourStateOpsin(sim=..., name='GtACR2', save_history=True, on_pre='', spectrum=[(400, 0.4), (410, 0.49), (420, 0.56), (430, 0.65), (440, 0.82), (450, 0.88), (460, 0.88), (470, 1.0), (480, 0.91), (490, 0.67), (500, 0.41), (510, 0.21), (520, 0.12), (530, 0.06), (540, 0.02), (550, 0.0), (560, 0.0)], required_vars=[('Iopto', amp), ('v', volt)], Gd1=17. * hertz, Gd2=10. * hertz, Gr0=0.58 * hertz, g0=44. * nsiemens, phim=2.e+23 * (second ** -1) / (meter ** 2), k1=40. * khertz, k2=20. * khertz, Gf0=1. * hertz, Gb0=3. * hertz, kf=1. * hertz, kb=5. * hertz, gamma=0.05, p=1, q=0.1, E=-69.5 * 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\n        IOPTO_VAR_NAME_post = -g0*fphi*(V_VAR_NAME_post-E)*rho_rel : ampere (summed)\n        rho_rel : 1'), Light(sim=..., name='amber fibers', save_history=True, value=array([0., 0., 0., 0.]), light_model=OpticFiber(R0=100. * umetre, NAfib=0.37, K=100. * metre ** -1, S=5896. * metre ** -1, ntis=1.36), wavelength=0.59 * umetre, direction=array([0., 0., 1.]), max_Irr0_mW_per_mm2=None, max_Irr0_mW_per_mm2_viz=None), BansalFourStateOpsin(sim=..., name='VfChrimson', save_history=True, on_pre='', spectrum=[(470.0, 0.4123404255319149), (490.0, 0.593265306122449), (510.0, 0.7935294117647058), (530.0, 0.8066037735849055), (550.0, 0.8912727272727272), (570.0, 1.0), (590.0, 0.9661016949152542), (610.0, 0.7475409836065574), (630.0, 0.4342857142857143)], required_vars=[('Iopto', amp), ('v', volt)], Gd1=0.37 * khertz, Gd2=175. * hertz, Gr0=0.667 * mhertz, g0=17.5 * nsiemens, phim=1.5e+22 * (second ** -1) / (meter ** 2), k1=3. * khertz, k2=200. * hertz, Gf0=20. * hertz, Gb0=3.2 * hertz, kf=10. * hertz, kb=10. * hertz, gamma=0.05, p=1, q=1, E=0. * volt, 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\n        IOPTO_VAR_NAME_post = -g0*fphi*(V_VAR_NAME_post-E)*rho_rel : ampere (summed)\n        rho_rel : 1'), Light(sim=..., name='blue fibers', save_history=True, value=array([0., 0., 0., 0.]), light_model=OpticFiber(R0=100. * umetre, NAfib=0.37, K=125. * metre ** -1, S=7370. * metre ** -1, ntis=1.36), wavelength=0.473 * umetre, direction=array([0., 0., 1.]), max_Irr0_mW_per_mm2=None, max_Irr0_mW_per_mm2_viz=None)})

Run simulation and plot results#

sim.reset()
sim.run(200*ms)

INFO       No numerical integration method specified for group 'neurongroup', using method 'euler' (took 0.15s, trying other methods took 0.00s). [brian2.stateupdaters.base.method_choice]
INFO       No numerical integration method specified for group 'syn_GtACR2_neurongroup', using method 'euler' (took 0.02s, trying other methods took 0.06s). [brian2.stateupdaters.base.method_choice]
INFO       No numerical integration method specified for group 'syn_VfChrimson_neurongroup', using method 'euler' (took 0.01s, trying other methods took 0.01s). [brian2.stateupdaters.base.method_choice]

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)

ax1.plot(spike_mon.t / ms, spike_mon.i[:], ",")
ax1.set(ylabel="neuron index", title="spiking")

ax2.step(
    blue_fibers.t_ms, blue_fibers.values + np.arange(n_fibers) * 1.3 + 0.1, c=c["473nm"]
)
ax2.step(amber_fibers.t_ms, amber_fibers.values + np.arange(n_fibers) * 1.3, c=c["590nm"])
ax2.set(
    yticks=[],
    ylabel="intensity per channel",
    title="photostimulation",
    xlabel="time (ms)",
)

[[],
 Text(0, 0.5, 'intensity per channel'),
 Text(0.5, 1.0, 'photostimulation'),
 Text(0.5, 0, 'time (ms)')]
../_images/multi_opto_15_1.png

As you can see, Vf-Chrimson has fast dynamics, enabling high-frequency control. GtACR2, on the other hand, continues acting long after the light is removed due to its slow deactivation kinetics. We can confirm this by plotting the current from each of the opsins in the first segment of neurons:

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.plot(Imon.t / ms, Imon.I_exc.T / namp, lw=0.2, c=c["590nm"])
ax1.set(title="Vf-Chrimson current", ylabel="$I_{exc}$ (nA)")
ax2.plot(Imon.t / ms, Imon.I_inh.T / namp, lw=0.2, c=c["473nm"])
ax2.set(title="GtACR2 current", ylabel="$I_{inh}$ (nA)")

[Text(0.5, 1.0, 'GtACR2 current'), Text(0, 0.5, '$I_{inh}$ (nA)')]
../_images/multi_opto_17_1.png

We can also confirm that blue light has a non-negligible effect on Vf-Chrimson.

Conclusion#

As a recap, in this tutorial we’ve seen how to:

  • configure a multi-channel Light device,

  • use more than one of them simultaneously, and

  • inject multiple opsins of overlapping action spectra into the same network.