Bridging model and experiment in systems neuroscience with Cleo: the Closed-Loop, Electrophysiology, and Optophysiology simulation testbed

Systems neuroscience has experienced an explosion of new tools for reading and writing neural activity, enabling exciting new experiments such as all-optical or closed-loop control that effect powerful causal interventions. At the same time, improved computational models are capable of reproducing behavior and neural activity with increasing fidelity. Unfortunately, these advances have drastically increased the complexity of integrating different lines of research, resulting in the missed opportunities and untapped potential of suboptimal experiments. Experiment simulation can help bridge this gap, allowing model and experiment to better inform each other by providing a low-cost testbed for experiment design, model validation, and methods engineering. Specifically, this can be achieved by incorporating the simulation of the experimental interface into our models, but no existing tool integrates optogenetics, two-photon calcium imaging, electrode recording, and flexible closed-loop processing with neural population simulations. To address this need, we have developed Cleo: the Closed-Loop, Electrophysiology, and Optophysiology experiment simulation testbed. Cleo is a Python package enabling injection of recording and stimulation devices as well as closed-loop control with realistic latency into a Brian spiking neural network model. It is the only publicly available tool currently supporting two-photon and multi-opsin/wavelength optogenetics. To facilitate adoption and extension by the community, Cleo is open-source, modular, tested, and documented, and can export results to various data formats. Here we describe the design and features of Cleo, validate output of individual components and integrated experiments, and demonstrate its utility for advancing optogenetic techniques in prospective experiments using previously published systems neuroscience models.


Light superposition
Here we describe how we model dynamics for light of two different wavelengths, given photon flux ϕ λ1 and ϕ λ2 .Bansal et al. [1] use a weighted sum of activation functions However, these are sublinear functions, so adding them results in exaggerated activation.In the extreme case, imagine two light sources that are just 1 nm apart in wavelength, with ε = 1.Thus: which contradicts what we expect for sublinear G: A more accurate approach instead assumes the following where φ λ2 is the standard-wavelength equivalent ("effective flux") of ϕ λ2 .
Demonstrating with G a1 of the four-state model [2], we solve for φ in terms of ϕ and ε: We then compute our activation functions as Unfortunately, this doesn't yield a simple constant conversion factor.However, if we approximate G as linear, we can use a weighted sum of fluxes: Plotting G for multiple opsins shows this linear approximation does yield a lower activation curve than the Bansal et al. approach and is qualitatively close to the true φ derived above.See Extended Data Fig. ?? for a comparison of the three methods.

Action spectrum normalization
Some action spectra are measured with equal power density/pulse width across wavelengths, while others are reported with equal photon flux.We store ours as equal power density spectra since they seem to be more common and allow for the opsin model to use both accurate power density and photon flux values.We use ε to represent sensitivity relative to the peak-sensitivity wavelength, ε ϕ and ε P representing the equal photon flux and power density versions, respectively.
For a given wavelength λ, We let G λ represent the response at wavelength λ, while G(ϕ) represents the response at the peak wavelength for the same flux ϕ.We will assume G is a linear function, as above, using C to represent a constant: Then we make either photon flux or irradiance (power density) constant: λ Thus, we can convert from ε ϕ (λ) to ε P (λ) by dividing by λ and normalizing.

GECI convolution simulation as an ODE
Song et al. convolve the intracellular calcium trace Ca 2+ with a double exponential kernel to capture variable rise and decay times in the fluorescence signal.To simplify simulation (to not have to keep a buffer of past calcium values), we can represent this convolution as an ODE.Let c(t) and b(t) be the free and bound calcium concentrations and h(t) be the kernel function where u(t) is the unit step function, included to ensure the kernel is causal.
We can represent the convolution as multiplication in the Laplace domain:
Now we get a common denominator and rearrange: Now we use the s 2 and s terms to convert to a second-order ODE, using the fact that the Laplace transform of b ′′ (t) and b ′ (t) are s 2 B(s) − sb(0 − ) − b(0 − ) and sB(s) − b(0 − ), respectively.Also, we assume that b(0) = b ′ (0) = 0 to avoid undefined δ(t) and δ ′ (t) terms after the inverse Laplace transform; this appears to have only a minor effect.
Rearranging to a first-order ODE system by introducing Special thanks to DinosaurEgg on Math Stack Exchange for helping solve this problem.

Hippocampal epilepsy model validation
Theta band power was computed using SciPy's spectrogram function with a Tukey window of width 3.906 sec, α = 0.25, and overlap width of 3.809 sec.

Prospective experiment 3
The reference signal was generated by delivering a 1 nA square wave input from 100 to 300 ms to entorhinal cortex using the original model's I ext term, without noise added.Training data was generated by running the system for 13 seconds with alternating on and off periods of length T ∼ N (200, 50) µs.During "on" intervals, I ext ∼ |N (0, Irr max /3) |.Irr max was 75 mW/mm 2 , which is described as an upper safety limit for 473 nm light delivery to the brain [3].Gaussian process noise was generated with mean µ = 0.167 nA and using the exponentiated quadratic kernel with σ = 0.083 nA, l = 30 ms and was added to input current I ext in each control scenario.The parameters of light delivery were altered from the 473 nm optic fiber defaults to allow for greater propagation-K and S were both divided by 10.The training data was fit using the ldsCtrlEst library's SSID and EM fitting methods with latent dimensionality n x = 4 [4].ldsCtrlEst's Gaussian linear quadratic regulator (LQR) controller was used with a gain computed from Q = C T C, R = 0.001 state and input penalties, with a simulated 3 ms of latency.Model-predictive control (MPC) was implemented with a control and prediction horizon of 6 and 36 control periods (each of which was 3 ms long) respectively.The standard, quadratic cost function utilized constant error weights Q = C T C, R = 10 −6 , and was optimized using OSQP via the JuMP optimization interface [5,6].MPC was simulated with 6 ms latency.