Overview#

Introduction#

Who is this package for?#

Cleo (Closed-Loop, Electrophysiology, and Optophysiology experiment simulation testbed) is a Python package developed to bridge theory and experiment for mesoscale neuroscience. We envision two primary uses cases:

  1. For prototyping closed-loop control of neural activity in silico. Animal experiments are costly to set up and debug, especially with the added complexity of real-time intervention—our aim is to enable researchers, given a decent spiking model of the system of interest, to assess whether the type of control they desire is feasible and/or what configuration(s) would be most conducive to their goals.

  2. The complexity of experimental interfaces means it’s not always clear what a model would look like in a real experiment. Cleo can help anyone interested in observing or manipulating a model while taking into account the constraints present in real experiments. Because Cleo is built around the Brian simulator, we especially hope this is helpful for existing Brian users who for whatever reason would like a convenient way to inject recorders (e.g., electrodes or 2P imaging) or stimulators (e.g., optogenetics) into the core network simulation.

What is closed-loop control?

In short, determining the inputs to deliver to a system from its outputs. In neuroscience terms, making the stimulation parameters a function of the data recorded in real time.

Structure and design#

Cleo wraps a spiking network simulator and allows for the injection of stimulators and/or recorders. The models used to emulate these devices are often non-trivial to implement or use in a flexible manner, so Cleo aims to make device injection and configuration as painless as possible, requiring minimal modification to the original network.

Cleo also orchestrates communication between the simulator and a user-configured IOProcessor object, modeling how experiment hardware takes samples, processes signals, and controls stimulation devices in real time.

For an explanation of why we choose to prioritize spiking network models and how we chose Brian as the underlying simulator, see Design rationale.

Why closed-loop control in neuroscience?

Fast, real-time, closed-loop control of neural activity enables intervention in processes that are too fast or unpredictable to control manually or with pre-defined stimulation, such as sensory information processing, motor planning, and oscillatory activity. Closed-loop control in a reactive sense enables the experimenter to respond to discrete events of interest, such as the arrival of a traveling wave or sharp wave ripple, whereas feedback control deals with driving the system towards a desired point or along a desired state trajectory. The latter has the effect of rejecting noise and disturbances, reducing variability across time and across trials, allowing the researcher to perform inference with less data and on a finer scale. Additionally, closed-loop control can compensate for model mismatch, allowing it to reach more complex targets where open-loop control based on imperfect models is bound to fail.

Installation#

Make sure you have Python ≥ 3.7, then use pip: pip install cleosim.

Caution

The name on PyPI is cleosim since cleo was already taken, but in code it is still used as import cleo. The other Cleo appears to actually be a fairly well developed package, so I’m sorry if you need to use it along with this Cleo in the same environment. In that case, there are workarounds.

Or, if you’re a developer, install poetry and run poetry install from the repository root.

Usage#

Brian network model#

The starting point for using Cleo is a Brian spiking neural network model of the system of interest. For those new to Brian, the docs are a great resource. If you have a model built with another simulator or modeling language, you may be able to import it to Brian via NeuroML.

Perhaps the biggest change you may have to make to an existing model to make it compatible with Cleo’s optogenetics and electrode recording is to give the neurons of interest coordinates in space. See the tutorials or the cleo.coords module for more info.

You’ll need your model in a Brian Network object before you move on. E.g.,:

import brian2 as b2
ng = b2.NeuronGroup( # a simple population of 100 LIF neurons
    500,
    """dv/dt = (-v - 70*mV + (500*Mohm)*Iopto + 2*xi*sqrt(tau_m)*mvolt) / tau_m : volt
    Iopto : amp""",
    threshold='v>-50*mV',
    reset='v=-70*mV',
    namespace={'tau_m': 20*b2.ms},
)
ng.v = -70*b2.mV
net = b2.Network(ng)

Your neurons need x, y, and z coordinates for Cleo’s spatially defined recording and stimulation models to work:

import cleo
cleo.coords.assign_coords_rand_rect_prism(
    ng, xlim=(-0.2, 0.2), ylim=(-0.2, 0.2), zlim=(0.2, 1), unit=b2.mm
)

CLSimulator#

Once you have a network model, you can construct a CLSimulator object:

sim = cleo.CLSimulator(net)

The simulator object wraps the Brian network and coordinates device injection, processing input and output, and running the simulation.

Recording#

Recording devices take measurements of the Brian network. To use a Recorder, you must inject it into the simulator via cleo.CLSimulator.inject(). The recorder will only record from the neuron groups specified on injection, allowing for such scenarios as singling out a cell type to record from. Some extremely simple implementations (which do little more than wrap Brian monitors) are available in the cleo.recorders module. See the Electrode recording and All-optical control tutorials for more detail on how to record from a simulation more realistically, but here’s a quick example of how to record multi-unit spiking activity with an electrode:

# configure and inject a 32-channel shank
coords = cleo.ephys.linear_shank_coords(
    array_length=1*b2.mm, channel_count=32, start_location=(0, 0, 0.2)*b2.mm
)
mua = cleo.ephys.MultiUnitSpiking(
    r_perfect_detection=20 * b2.um,
    r_half_detection=40 * b2.um,
)
probe = cleo.ephys.Probe(coords, signals=[mua])
sim.inject(probe, ng)

---------------------------------------------------------------------------
DimensionMismatchError                    Traceback (most recent call last)
Cell In[4], line 10
      5 mua = cleo.ephys.MultiUnitSpiking(
      6     r_perfect_detection=20 * b2.um,
      7     r_half_detection=40 * b2.um,
      8 )
      9 probe = cleo.ephys.Probe(coords, signals=[mua])
---> 10 sim.inject(probe, ng)

File ~/checkouts/readthedocs.org/user_builds/cleosim/envs/v0.14.2/lib/python3.12/site-packages/cleo/base.py:383, in CLSimulator.inject(self, device, *neuron_groups, **kwparams)
    381         device.sim = self
    382         device.init_for_simulator(self)
--> 383     device.connect_to_neuron_group(ng, **kwparams)
    384 for brian_object in device.brian_objects:
    385     if brian_object not in self.network.objects:

File ~/checkouts/readthedocs.org/user_builds/cleosim/envs/v0.14.2/lib/python3.12/site-packages/cleo/ephys/probes.py:158, in Probe.connect_to_neuron_group(self, neuron_group, **kwparams)
    146 """Configure probe to record from given neuron group
    147 
    148 Will call :meth:`Signal.connect_to_neuron_group` for each signal
   (...)
    155     Passed in to signals' connect functions, needed for some signals
    156 """
    157 for signal in self.signals:
--> 158     signal.connect_to_neuron_group(neuron_group, **kwparams)
    159     self.brian_objects.update(signal.brian_objects)

File ~/checkouts/readthedocs.org/user_builds/cleosim/envs/v0.14.2/lib/python3.12/site-packages/cleo/ephys/spiking.py:209, in MultiUnitSpiking.connect_to_neuron_group(self, neuron_group, **kwparams)
    201 def connect_to_neuron_group(self, neuron_group: NeuronGroup, **kwparams) -> None:
    202     """Configure signal to record from specified neuron group
    203 
    204     Parameters
   (...)
    207         group to record from
    208     """
--> 209     neuron_channel_dtct_probs = super(
    210         MultiUnitSpiking, self
    211     ).connect_to_neuron_group(neuron_group, **kwparams)
    212     if self._dtct_prob_array is None:
    213         self._dtct_prob_array = neuron_channel_dtct_probs

File ~/checkouts/readthedocs.org/user_builds/cleosim/envs/v0.14.2/lib/python3.12/site-packages/cleo/ephys/spiking.py:100, in Spiking.connect_to_neuron_group(self, neuron_group, **kwparams)
     98 distances = np.sqrt(dist2) * meter  # since units stripped
     99 # probs is n_neurons by n_channels
--> 100 probs = self._detection_prob_for_distance(distances)
    101 # cut off to get indices of neurons to monitor
    102 # [0] since nonzero returns tuple of array per axis
    103 i_ng_to_keep = np.nonzero(np.all(probs > self.cutoff_probability, axis=1))[0]

File ~/checkouts/readthedocs.org/user_builds/cleosim/envs/v0.14.2/lib/python3.12/site-packages/cleo/ephys/spiking.py:149, in Spiking._detection_prob_for_distance(self, r)
    147 h = b - a
    148 with np.errstate(divide="ignore"):
--> 149     decaying_p = h / (r - c)
    150 decaying_p[decaying_p == np.inf] = 1  # to fix NaNs caused by /0
    151 p = 1 * (r <= a) + decaying_p * (r > a)

File ~/checkouts/readthedocs.org/user_builds/cleosim/envs/v0.14.2/lib/python3.12/site-packages/brian2/units/fundamentalunits.py:1193, in Quantity.__array_ufunc__(self, uf, method, *inputs, **kwargs)
   1190 elif uf.__name__ in UFUNCS_MATCHING_DIMENSIONS + UFUNCS_COMPARISONS:
   1191     # Ok if dimension of arguments match (for reductions, they always do)
   1192     if method == "__call__":
-> 1193         fail_for_dimension_mismatch(
   1194             inputs[0],
   1195             inputs[1],
   1196             error_message=(
   1197                 "Cannot calculate {val1} %s {val2}, the units do not match"
   1198             )
   1199             % uf.__name__,
   1200             val1=inputs[0],
   1201             val2=inputs[1],
   1202         )
   1203     if uf.__name__ in UFUNCS_COMPARISONS:
   1204         return uf_method(*[np.asarray(i) for i in inputs], **kwargs)

File ~/checkouts/readthedocs.org/user_builds/cleosim/envs/v0.14.2/lib/python3.12/site-packages/brian2/units/fundamentalunits.py:266, in fail_for_dimension_mismatch(obj1, obj2, error_message, **error_quantities)
    264         raise DimensionMismatchError(error_message, dim1)
    265     else:
--> 266         raise DimensionMismatchError(error_message, dim1, dim2)
    267 else:
    268     return dim1, dim2

DimensionMismatchError: Cannot calculate [[253.59191569 234.93521353 ... 834.37801816 865.74091715]
 [268.86480928 236.73975294 ... 700.22711039 732.46869875]
 ...
 [622.11971536 591.11944268 ... 404.99082745 434.55930455]
 [775.16460481 743.06894225 ... 210.53690407 240.89840169]] mm^2 subtract 0. m, the units do not match (units are m^2 and m).

Stimulation#

Stimulator devices manipulate the Brian network, and are likewise inject()ed into the simulator, into specified neuron groups. Optogenetics (1P and 2P) is the main stimulation modality currently implemented by Cleo. This requires injection of both a light source and an opsin—see the Optogenetic stimulation and All-optical control tutorials for more detail.

fiber = cleo.light.Light(
    coords=(0, 0, 0.5)*b2.mm,
    light_model=cleo.light.fiber473nm(),
    wavelength=473*b2.nmeter
)
chr2 = cleo.opto.chr2_4s()
sim.inject(fiber, ng).inject(chr2, ng)

Note

Markov opsin kinetics models require target neurons to have membrane potentials in realistic ranges and an Iopto term defined in amperes. If you need to interface with a model without these features, you may want to use the simplified ProportionalCurrentOpsin. You can find more details, including a comparison between the two model types, in the optogenetics tutorial.

Note

Recorders and stimulators need unique names which serve as keys to access/set values in the IOProcessor’s input/output dictionaries. The name defaults to the class name, but you can specify it on construction.

I/O Processor#

Just as in a real experiment where the experiment hardware must be connected to signal processing equipment and/or computers for recording and control, the CLSimulator must be connected to an IOProcessor. If you are only recording, you may want to use the RecordOnlyProcessor. Otherwise you will want to implement the LatencyIOProcessor, which not only takes samples at the specified rate, but processes the data and delivers input to the network after a user-defined delay, emulating the latency inherent in real experiments. This is done by creating a subclass and defining the process() function:

class MyProcessor(cleo.ioproc.LatencyIOProcessor):
    def process(self, state_dict, sample_time_ms):
        # state_dict contains a {'recorder_name': value} dict of network.
        i_spikes, t_ms_spikes, y_spikes = state_dict['Probe']['MultiUnitSpiking']
        # on-off control
        irr0_mW_per_mm2 = 5 if len(i_spikes) < 10 else 0
        # output is a {'stimulator_name': value} dict and output time
        return {'Light': irr0_mW_per_mm2}, sample_time_ms + 3  # (3 ms delay)

sim.set_io_processor(MyProcessor(sample_period_ms=1))

The On-off control, PI control, and LQR optimal control using ldsctrlest tutorials give examples of closed-loop control ranging from simple to complex.

Visualization#

cleo.viz.plot() allows you to easily visualize your experimental configuration:

cleo.viz.plot(ng, colors=['#c500cc'], sim=sim, zlim=(200, 1000))

Cleo also features some video visualization capabilities.

Running experiments#

Use cleo.CLSimulator.run() function with the desired duration. This wrap’s Brian’s brian2.core.network.Network.run() function:

sim.run(50 * b2.ms)  # kwargs are passed to Brian's run function

Use reset() to restore the default state (right after initialization/injection) for the network and all devices. This could be useful for running a simulation multiple times under different conditions.

To facilitate access to data after the simulation, devices offer a cleo.InterfaceDevice.save_history option on construction, by default True. If true, that object will store relevant variables as attributes. For example:

import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.scatter(mua.t_ms, mua.i, marker='.', c='white', s=2)
ax1.set(ylabel='channel index', title='spikes')
ax2.step(fiber.t_ms, fiber.values, c='#72b5f2')
ax2.set(xlabel='time (ms)', ylabel='irradiance (mW/mm²)', title='photostimulation')

Design rationale#

Why not prototype with more abstract models?#

Cleo aims to be practical, and as such provides models at the level of abstraction corresponding to the variables the experimenter has available to manipulate. This means models of spatially defined, spiking neural networks.

Of course, neuroscience is studied at many spatial and temporal scales. While other projects may be better suited for larger segments of the brain and/or longer timescales (such as HNN or BMTK’s PopNet or FilterNet), this project caters to finer-grained models because they can directly simulate the effects of alternate experimental configurations. For example, how would the model change when swapping one opsin for another, using multiple opsins simultaneously, or with heterogeneous expression? How does recording or stimulating one cell type vs. another affect the experiment? Would using a more sophisticated control algorithm be worth the extra compute time, and thus later stimulus delivery, compared to a simpler controller?

Questions like these could be answered using an abstract dynamical system model of a neural circuit, but they would require the extra step of mapping the afore-mentioned details to a suitable abstraction—e.g., estimating a transfer function to model optogenetic stimulation for a given opsin and light configuration. Thus, we haven’t emphasized these sorts of models so far in our development of Cleo, though they should be possible to implement in Brian if you are interested. For example, one could develop a Poisson linear dynamical system (PLDS), record spiking output, and configure stimulation to act directly on the system’s latent state.

And just as experiment prototyping could be done on a more abstract level, it could also be done on an even more realistic level, which we did not deem necessary. That brings us to the next point…

Why Brian?#

Brian is a relatively new spiking neural network simulator written in Python. Here are some of its advantages:

  • Flexibility: allowing (and requiring!) the user to define models mathematically rather than selecting from a pre-defined library of cell types and features. This enables us to define arbitrary models for recorders and stimulators and easily interface with the simulation

  • Ease of use: it’s all just Python

  • Speed (at least when compiled to C++ or GPU—see Brian2GENN, Brian2CUDA)

NEST is a popular alternative to Brian also strong in point neuron simulations. However, it appears to be less flexible, and thus harder to extend. NEURON is another popular alternative to Brian. Its main advantage is its first-class support of detailed, morphological, multi-compartment neurons. In fact, strong alternatives to Brian for this project were BioNet (docs, paper) and NetPyNE (docs, paper), which already offer a high-level interface to NEURON with extracellular potential recording. Optogenetics could be incorporated with pre-existing .hoc code, though the light model would need to be implemented. From brief examination of the source code of BioNet, it appears that closed-loop stimulation would not be too difficult to add. It is unclear for NetPyNE.

In the end, we chose Brian since our priority was to model circuit/population-level dynamics over molecular/intra-neuron dynamics. Also, Brian does have support for multi-compartment neurons, albeit less fully featured, if that is needed.

PyNN would have been an ideal interface to build around, as it supports multiple simulator backends. The difficulty is, implementing objects not natively supported by SNN simulators (e.g., opsins, calcium indicators, and light source) has required bespoke, idiosyncratic code applicable only to one simulator. To do so in a native, efficient way, as we have attempted to do with Brian, would require significant work. A collaborative effort extending a multi-simulator framework such as PyNN for this purpose may be of value if there is enough community interest in expanding the open-source SNN experiment simulation toolbox.

Future development#

Here are some features which are missing but could be useful to add:

  • Electrode microstimulation

  • A more accurate LFP signal (only usable for morphological neurons) based on the volume conductor forward model as in LFPy or Vertex

  • The Mazzoni-Lindén LFP approximation for LIF point-neuron networks

  • Voltage indicators

  • An expanded calcium indicator library—currently the parameter set for the full NAOMi model is only available for GCaMP6f. The phenomenological S2F model should be easy to implement and fit to data, has parameters for several of the latest and greatest GECIs (jGCaMP7 and jGCaMP8 varieties), and should be cheaper to simulate as well.