All-optical control#
In this tutorial we’ll see how to
configure and inject a microscope for two-photon imaging and stimulation
inject a calcium indicator for neurons visible to the microscope
set up two-photon optogenetics by
injecting an opsin with an extended action spectrum and
injecting a 2P light source targeting neurons imaged by the microscope
Preamble#
import brian2 as b2
from brian2 import np
import cleo
from cleo import opto, imaging, light
import matplotlib.pyplot as plt
# for reproducibility
rng = np.random.default_rng(92)
np.random.seed(92)
cleo.utilities.style_plots_for_docs()
Brian network setup#
All we need are LIF neurons that will spike in response to photostimulation for us to see the resulting fluorescence traces.
ng = b2.NeuronGroup(
200,
"""dv/dt = (-(v - E_L) + Rm*Iopto) / tau_m : volt
Iopto : amp""",
threshold="v > -50*mV",
reset="v=E_L",
namespace={
"tau_m": 20 * b2.ms,
"Rm": 500 * b2.Mohm,
"E_L": -70 * b2.mV,
},
)
ng.v = -70 * b2.mV
cleo.coords.assign_coords_rand_rect_prism(ng, [-0.2, 0.2], [-0.2, 0.2], [0, 0.4])
sim = cleo.CLSimulator(b2.Network(ng))
Microscope configuration#
By default a scope selects neurons based on focus depth and assigns them a signal-to-noise ratio (SNR) based on soma size.
(Larger cells’ SNR decays more slowly with distance from the focal plane.)
focus_depth and soma_radius are taken from scope parameters but can be overridden on injection.
We’ll use GCaMP6f as our indicator and leave the default parameters, though those can be changed in the function call gcamp6f().
scope = imaging.Scope(
focus_depth=100 * b2.um,
img_width=500 * b2.um,
sensor=imaging.gcamp6f(),
)
Heterogeneous expression is simulated by providing a rho_rel_generator function that takes the number of neurons as an argument and returns a vector of relative expression levels.
This is done on scope injection (as opposed to sensor injection afterward) because the scope needs this information to select ROIs with sufficiently high SNR.
expr_level_gen = lambda n: rng.lognormal(0, 0.2, n)
sim.inject(scope, ng, rho_rel_generator=expr_level_gen) # uses scope's parameters
# optional overrides (we'll use a small soma to intentionally get few targets)
sim.inject(scope, ng, focus_depth=200 * b2.um, soma_radius=5 * b2.um, rho_rel_generator=expr_level_gen)
CLSimulator(io_processor=None, devices={Scope(sim=..., name='Scope', save_history=True, sensor=GECI(sim=None, name='gcamp6f', save_history=True, model='\n dCa/dt = -gamma * (Ca - Ca_rest) / (1 + kappa_S + kappa_B) : mmolar (clock-driven)\n kappa_B = B_T * K_d / (Ca + K_d)**2 : 1\n\n CaB_active = Ca_rest + b : mmolar # add tiny bit to avoid /0\n db/dt = beta : mmolar (clock-driven)\n lam = 1/tau_off + 1/tau_on : 1/second\n kap = 1/tau_off : 1/second\n dbeta/dt = ( # should be M/s/s\n A * (lam - kap) * (Ca - Ca_rest) # M/s/s\n - (kap + lam) * beta # M/s/s\n - kap * lam * b # M/s/s\n ) : mmolar/second (clock-driven)\n \nexc_factor = 1 : 1\n\n dFF_baseline = 1 / (1 + (K_d / Ca_rest) ** n_H) : 1\n dFF = exc_factor * rho_rel * dFF_max * (\n 1 / (1 + (K_d / CaB_active) ** n_H)\n - dFF_baseline\n ) : 1\n rho_rel : 1\n ', on_pre='Ca += dCa_T / (1 + kappa_S + kappa_B)', sigma_noise=0.03748181818181818, dFF_1AP=0.09775500000000001, location='cytoplasm', cal_model=DynamicCalcium(on_pre='Ca += dCa_T / (1 + kappa_S + kappa_B)', model='\n dCa/dt = -gamma * (Ca - Ca_rest) / (1 + kappa_S + kappa_B) : mmolar (clock-driven)\n kappa_B = B_T * K_d / (Ca + K_d)**2 : 1', Ca_rest=50. * nmolar, gamma=292.3 * hertz, B_T=200. * umolar, kappa_S=110, dCa_T=7.6 * umolar), bind_act_model=DoubExpCalBindingActivation(model='\n CaB_active = Ca_rest + b : mmolar # add tiny bit to avoid /0\n db/dt = beta : mmolar (clock-driven)\n lam = 1/tau_off + 1/tau_on : 1/second\n kap = 1/tau_off : 1/second\n dbeta/dt = ( # should be M/s/s\n A * (lam - kap) * (Ca - Ca_rest) # M/s/s\n - (kap + lam) * beta # M/s/s\n - kap * lam * b # M/s/s\n ) : mmolar/second (clock-driven)\n ', A=7.61251 * khertz, tau_on=1.17164616 * second, tau_off=10.14020867 * msecond, Ca_rest=50. * nmolar), exc_model=NullExcitation(model='exc_factor = 1 : 1'), fluor_model='\n dFF_baseline = 1 / (1 + (K_d / Ca_rest) ** n_H) : 1\n dFF = exc_factor * rho_rel * dFF_max * (\n 1 / (1 + (K_d / CaB_active) ** n_H)\n - dFF_baseline\n ) : 1\n rho_rel : 1\n ', K_d=290. * nmolar, n_H=2.7, dFF_max=25.2), img_width=0.5 * mmetre, focus_depth=100. * umetre, location=array([0., 0., 0.]) * metre, direction=array([0., 0., 1.]), soma_radius=10. * umetre, snr_cutoff=1)})
Signal strength (ΔF/F for 1 spike) per ROI is thus defined as dFF_1AP * rho_rel.
While Cleo adjusts noise levels according to distance from the focal plane automatically, other adjustments can be made manually by calling target_neurons_in_plane() to reflect factors such as cell heterogeneity or indicator expression or increased noise with depth.
Neurons with SNR (defined as dFF_1AP / sigma_noise) below the scope’s snr_cutoff parameter are discarded.
i_targets, noise_focus_factor, focus_coords = scope.target_neurons_in_plane(ng, focus_depth=300 * b2.um, soma_radius=5 * b2.um)
# scale noise std σ randomly to simulate biological variability
std_noise = scope.sensor.sigma_noise * noise_focus_factor * rng.normal(1, .2, len(i_targets))
sim.inject(scope, ng, focus_depth=None, i_targets=i_targets, sigma_noise=std_noise, rho_rel_generator=expr_level_gen)
cleo.viz.plot(ng, colors=['#c500cc'], sim=sim)
(<Figure size 640x480 with 1 Axes>,
<Axes3D: xlabel='x (um)', ylabel='y (um)', zlabel='z (um)'>)
We can see targets off the visualized plane, resulting from our injections at 200 and 300 μm depths, resembling a multi-plane imaging experiment.
When focus_depth is set to None, corresponding to sculpted holographic imaging, the scope will select all neurons in the volume, or the user can specify a list of neurons to select via i_targets on injection.
After all targets are specified, the sensor protein must also be injected. Doing it this way allows for efficiency: the sensor protein is only simulated for the imaged neurons. Expression levels are as sampled previously on scope injection.
scope.inject_sensor_for_targets()
2p stimulation configuration#
Future work will address how to better quantify the relationship between 2P and the 1P action spectra, but for now we’ll pretend that Vf-Chrimson activation at 1060 nm is 1/100th that of peak 1P activation for the same power density.
opsin = opto.vfchrimson_4s()
opsin.spectrum.append((1060, .01))
sim.inject(opsin, ng)
CLSimulator(io_processor=None, devices={GECI(sim=..., name='gcamp6f', save_history=True, model='\n dCa/dt = -gamma * (Ca - Ca_rest) / (1 + kappa_S + kappa_B) : mmolar (clock-driven)\n kappa_B = B_T * K_d / (Ca + K_d)**2 : 1\n\n CaB_active = Ca_rest + b : mmolar # add tiny bit to avoid /0\n db/dt = beta : mmolar (clock-driven)\n lam = 1/tau_off + 1/tau_on : 1/second\n kap = 1/tau_off : 1/second\n dbeta/dt = ( # should be M/s/s\n A * (lam - kap) * (Ca - Ca_rest) # M/s/s\n - (kap + lam) * beta # M/s/s\n - kap * lam * b # M/s/s\n ) : mmolar/second (clock-driven)\n \nexc_factor = 1 : 1\n\n dFF_baseline = 1 / (1 + (K_d / Ca_rest) ** n_H) : 1\n dFF = exc_factor * rho_rel * dFF_max * (\n 1 / (1 + (K_d / CaB_active) ** n_H)\n - dFF_baseline\n ) : 1\n rho_rel : 1\n ', on_pre='Ca += dCa_T / (1 + kappa_S + kappa_B)', sigma_noise=0.03748181818181818, dFF_1AP=0.09775500000000001, location='cytoplasm', cal_model=DynamicCalcium(on_pre='Ca += dCa_T / (1 + kappa_S + kappa_B)', model='\n dCa/dt = -gamma * (Ca - Ca_rest) / (1 + kappa_S + kappa_B) : mmolar (clock-driven)\n kappa_B = B_T * K_d / (Ca + K_d)**2 : 1', Ca_rest=50. * nmolar, gamma=292.3 * hertz, B_T=200. * umolar, kappa_S=110, dCa_T=7.6 * umolar), bind_act_model=DoubExpCalBindingActivation(model='\n CaB_active = Ca_rest + b : mmolar # add tiny bit to avoid /0\n db/dt = beta : mmolar (clock-driven)\n lam = 1/tau_off + 1/tau_on : 1/second\n kap = 1/tau_off : 1/second\n dbeta/dt = ( # should be M/s/s\n A * (lam - kap) * (Ca - Ca_rest) # M/s/s\n - (kap + lam) * beta # M/s/s\n - kap * lam * b # M/s/s\n ) : mmolar/second (clock-driven)\n ', A=7.61251 * khertz, tau_on=1.17164616 * second, tau_off=10.14020867 * msecond, Ca_rest=50. * nmolar), exc_model=NullExcitation(model='exc_factor = 1 : 1'), fluor_model='\n dFF_baseline = 1 / (1 + (K_d / Ca_rest) ** n_H) : 1\n dFF = exc_factor * rho_rel * dFF_max * (\n 1 / (1 + (K_d / CaB_active) ** n_H)\n - dFF_baseline\n ) : 1\n rho_rel : 1\n ', K_d=290. * nmolar, n_H=2.7, dFF_max=25.2), Scope(sim=..., name='Scope', save_history=True, sensor=GECI(sim=..., name='gcamp6f', save_history=True, model='\n dCa/dt = -gamma * (Ca - Ca_rest) / (1 + kappa_S + kappa_B) : mmolar (clock-driven)\n kappa_B = B_T * K_d / (Ca + K_d)**2 : 1\n\n CaB_active = Ca_rest + b : mmolar # add tiny bit to avoid /0\n db/dt = beta : mmolar (clock-driven)\n lam = 1/tau_off + 1/tau_on : 1/second\n kap = 1/tau_off : 1/second\n dbeta/dt = ( # should be M/s/s\n A * (lam - kap) * (Ca - Ca_rest) # M/s/s\n - (kap + lam) * beta # M/s/s\n - kap * lam * b # M/s/s\n ) : mmolar/second (clock-driven)\n \nexc_factor = 1 : 1\n\n dFF_baseline = 1 / (1 + (K_d / Ca_rest) ** n_H) : 1\n dFF = exc_factor * rho_rel * dFF_max * (\n 1 / (1 + (K_d / CaB_active) ** n_H)\n - dFF_baseline\n ) : 1\n rho_rel : 1\n ', on_pre='Ca += dCa_T / (1 + kappa_S + kappa_B)', sigma_noise=0.03748181818181818, dFF_1AP=0.09775500000000001, location='cytoplasm', cal_model=DynamicCalcium(on_pre='Ca += dCa_T / (1 + kappa_S + kappa_B)', model='\n dCa/dt = -gamma * (Ca - Ca_rest) / (1 + kappa_S + kappa_B) : mmolar (clock-driven)\n kappa_B = B_T * K_d / (Ca + K_d)**2 : 1', Ca_rest=50. * nmolar, gamma=292.3 * hertz, B_T=200. * umolar, kappa_S=110, dCa_T=7.6 * umolar), bind_act_model=DoubExpCalBindingActivation(model='\n CaB_active = Ca_rest + b : mmolar # add tiny bit to avoid /0\n db/dt = beta : mmolar (clock-driven)\n lam = 1/tau_off + 1/tau_on : 1/second\n kap = 1/tau_off : 1/second\n dbeta/dt = ( # should be M/s/s\n A * (lam - kap) * (Ca - Ca_rest) # M/s/s\n - (kap + lam) * beta # M/s/s\n - kap * lam * b # M/s/s\n ) : mmolar/second (clock-driven)\n ', A=7.61251 * khertz, tau_on=1.17164616 * second, tau_off=10.14020867 * msecond, Ca_rest=50. * nmolar), exc_model=NullExcitation(model='exc_factor = 1 : 1'), fluor_model='\n dFF_baseline = 1 / (1 + (K_d / Ca_rest) ** n_H) : 1\n dFF = exc_factor * rho_rel * dFF_max * (\n 1 / (1 + (K_d / CaB_active) ** n_H)\n - dFF_baseline\n ) : 1\n rho_rel : 1\n ', K_d=290. * nmolar, n_H=2.7, dFF_max=25.2), img_width=0.5 * mmetre, focus_depth=100. * umetre, location=array([0., 0., 0.]) * metre, direction=array([0., 0., 1.]), soma_radius=10. * umetre, snr_cutoff=1), 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), (1060, 0.01)], 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')})
Now, to target these neurons with 2P laser power, we use tp_light_from_scope() to create a Light object with a 2P laser GaussianEllipsoid light profile centered on each of the scope’s targets:
laser = light.tp_light_from_scope(scope, wavelength=1060*b2.nmeter, name='laser')
sim.inject(laser, ng)
cleo.viz.plot(ng, colors=['#c500cc'], sim=sim)
(<Figure size 640x480 with 1 Axes>,
<Axes3D: xlabel='x (um)', ylabel='y (um)', zlabel='z (um)'>)
Simulating an all-optical experiment#
We’ll stimulate each neuron with a different number of laser pulses and see the resulting calcium traces.
from cleo.ioproc import LatencyIOProcessor
# for seeing ground-truth spikes
i_all_targets = scope.i_targets_for_neuron_group(ng)
smon = b2.SpikeMonitor(ng, record=i_all_targets)
sim.network.add(smon)
amplitude_mW = 2.5
pulse_width_ms = 2
interpulse_ms = 20
num_pulses = np.arange(1, scope.n + 1)
train_ends_ms = interpulse_ms * num_pulses
class IOProc(LatencyIOProcessor):
def process(self, state_dict, time_ms):
t_in_cycle = time_ms % interpulse_ms
on = (0 <= t_in_cycle) & (t_in_cycle < pulse_width_ms) & (time_ms < (interpulse_ms * num_pulses))
return {"laser": on * amplitude_mW}, time_ms
sim.set_io_processor(IOProc(1))
sim.reset()
tmax = 1000
sim.run(tmax * b2.ms)
INFO No numerical integration method specified for group 'neurongroup', using method 'exact' (took 0.18s). [brian2.stateupdaters.base.method_choice]
INFO No numerical integration method specified for group 'syn_VfChrimson_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_gcamp6f_neurongroup', using method 'euler' (took 0.02s, trying other methods took 0.05s). [brian2.stateupdaters.base.method_choice]
Now we’ll plot the result:
import matplotlib.pyplot as plt
sep = .5
n2plot = len(i_all_targets)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6), sharex=True)
spt = np.empty((0,))
spi = np.empty((0,))
for i_target, i_ng in enumerate(i_all_targets):
single_spt = smon.t[smon.i == i_ng]
spt = np.concatenate((spt, single_spt))
spi = np.concatenate((spi, np.full_like(single_spt, i_target)))
ax1.scatter(spt / b2.ms, spi, marker='|', label='spikes')
t_pulse_ms = np.arange(0, num_pulses.max() * interpulse_ms, interpulse_ms)
ax1.scatter(t_pulse_ms, np.full_like(t_pulse_ms, -1), marker='|', color='r', label='pulses')
ax1.set(ylabel='neuron index', title='ground-truth spikes', xlim=(0, tmax))
ax1.legend()
ax2.plot(scope.t_ms[:], np.array(scope.dFF)[:, :n2plot] + sep*np.arange(n2plot), lw=.2, c='w')
y_scale = 2
y0 = 0
x0 = -20
ax2.set(xlabel='time (ms)', title='ΔF/F₀ per ROI', xlim=(x0, tmax))
ax2.annotate('', xy=(x0, y0), xytext=(x0, y0 + y_scale), arrowprops=dict(arrowstyle='-'))
ax2.annotate(f'{y_scale * 100}% ΔF/F', xy=(x0, y0), va='bottom', ha='right', rotation=90)
ax2.yaxis.set_visible(False)
ax2.spines['left'].set_visible(False)
Results are as expected:
We see one spike for every pulse, though for longer trains we might start to see missed spikes due to the opsin activation plateauing.
If neurons were closer together, we’d expect to see some off-target effects; if there are any here, they are apparently subthreshold.
Signal strength is proportional to the number of spikes, as expected, but varies somewhat due to heterogeneous expression.
We also see that noise levels vary from ROI to ROI due to varying distances from the focal plane.