Neo export#

Cleo allows you to export data from the whole CLSimulator or from individual InterfaceDevices objects to Neo objects. This facilitates data analysis and visualization with packages that take Neo as an input format, as well as export to various open and proprietary ephys data file formats, which could be useful for integration with existing pipelines developed for experimental data.

Setup#

import brian2 as b2
from brian2 import np
import cleo
from cleo import opto, ephys, light
import neo

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

cleo.utilities.style_plots_for_docs()

We’ll set up a basic simulation to see how export works for different devices and the whole simulation.

ng = b2.NeuronGroup(
    100,
    """dv/dt = ((-70*mV - v) + (500*Mohm)*Iopto) / (20*ms) : volt
    Iopto : amp""",
    threshold="v > -50*mV",
    reset="v = -70*mV",
)
ng.v = -70 * b2.mV
syn = b2.Synapses(ng, on_pre="v += 3*mV")
cleo.coords.assign_coords_rand_rect_prism(ng, (-.2, .2), (-.2, .2), (0, .4))

light = light.Light(
    coords=[[-0.1, 0, 0], [0.1, 0, 0]] * b2.mm, light_model=light.fiber473nm()
)
probe = ephys.Probe(
    coords=[[0, 0, 50], [0, 0, 150], [0, 0, 250]] * b2.um,
    signals=[
        ephys.TKLFPSignal(),
        ephys.MultiUnitSpiking(
            r_perfect_detection=50 * b2.um, r_half_detection=100 * b2.um
        ),
        ephys.SortedSpiking(
            r_perfect_detection=50 * b2.um, r_half_detection=100 * b2.um
        ),
    ],
)
sim = cleo.CLSimulator(b2.Network(ng))

class IOProc(cleo.ioproc.LatencyIOProcessor):
    def process(self, state_dict, t_samp_ms):
        return {'Light': .5 * np.random.rand(light.n)}, t_samp_ms
sim.set_io_processor(IOProc(1))
sim.inject(light, ng).inject(probe, ng, tklfp_type="exc")
sim.inject(opto.chr2_4s(), ng)

cleo.viz.plot(ng, colors=['orange'], sim=sim)

(<Figure size 640x480 with 1 Axes>,
 <Axes3D: xlabel='x (mm)', ylabel='y (mm)', zlabel='z (mm)'>)
../_images/999bb9a6d2cc878a24ae29805fa4ac2892be0f811328c42cb72c90249af0c849.png 
sim.run(50 * b2.ms)

INFO       No numerical integration method specified for group 'neurongroup', using method 'exact' (took 0.23s). [brian2.stateupdaters.base.method_choice]
INFO       No numerical integration method specified for group 'syn_ChR2_neurongroup', using method 'euler' (took 0.02s, trying other methods took 0.06s). [brian2.stateupdaters.base.method_choice]

Exporting individual devices#

Depending on whether time values are regularly spaced, data is exported to either AnalogSignal or IrregularlySampledSignal objects, with appropriate annotations and per-channel annotations.

light.to_neo()

AnalogSignal with 2 channels of length 50; units mW/mm**2; datatype float64 
name: 'Light'
description: 'Exported from Cleo Light device'
annotations: {'export_datetime': datetime.datetime(2023, 9, 6, 9, 59, 33, 181618)}
sampling rate: 1.0 1/ms
time: 0.0 ms to 50.0 ms

light.to_neo().array_annotations

{'x': array([-0.1,  0.1]) * mm,
 'y': array([0., 0.]) * mm,
 'z': array([0., 0.]) * mm,
 'direction_x': array([0., 0.]),
 'direction_y': array([0., 0.]),
 'direction_z': array([1., 1.]),
 'i_channel': array([0, 1])}

Different signals from the same device are grouped together:

probe_neo = probe.to_neo()
probe_neo

Group with 2 groups, 1 analogsignals
name: 'Probe'
description: 'Exported from Cleo Probe device'

probe_neo.children

(AnalogSignal with 3 channels of length 50; units uV; datatype float64 
 name: 'Probe.TKLFPSignal'
 description: 'Exported from Cleo TKLFPSignal object'
 annotations: {'export_datetime': datetime.datetime(2023, 9, 6, 9, 59, 33, 239616)}
 sampling rate: 1.0 1/ms
 time: 0.0 ms to 50.0 ms,
 Group with 3 spiketrains
 name: 'Probe.MultiUnitSpiking'
 description: 'Exported from Cleo MultiUnitSpiking object'
 annotations: {'export_datetime': datetime.datetime(2023, 9, 6, 9, 59, 33, 241341)},
 Group with 55 spiketrains
 name: 'Probe.SortedSpiking'
 description: 'Exported from Cleo SortedSpiking object'
 annotations: {'export_datetime': datetime.datetime(2023, 9, 6, 9, 59, 33, 265177)})

Whole-simulation export#

Exporting the whole CLSimulator object generates a Block object with a single Segment object containing all the data from the simulation.

sim.to_neo()

Block with 1 segments, 1 groups
description: 'Exported from Cleo simulation'
annotations: {'export_datetime': datetime.datetime(2023, 9, 6, 9, 59, 33, 320150)}
# segments (N=1)
0: Segment with 2 analogsignals, 58 spiketrains
   # analogsignals (N=2)
   0: AnalogSignal with 2 channels of length 50; units mW/mm**2; datatype float64 
      name: 'Light'
      description: 'Exported from Cleo Light device'
      annotations: {'export_datetime': datetime.datetime(2023, 9, 6, 9, 59, 33, 320674)}
      sampling rate: 1.0 1/ms
      time: 0.0 ms to 50.0 ms
   1: AnalogSignal with 3 channels of length 50; units uV; datatype float64 
      name: 'Probe.TKLFPSignal'
      description: 'Exported from Cleo TKLFPSignal object'
      annotations: {'export_datetime': datetime.datetime(2023, 9, 6, 9, 59, 33, 321221)}
      sampling rate: 1.0 1/ms
      time: 0.0 ms to 50.0 ms

Trial structure#

Using the convention that a Segment represents a single “trial,” there are a couple ways we can generate trial-structured Neo data:

  1. Resetting the simulation and exporting the data from each trial individually, merging them into a single Block object.

  2. Structuring multiple trials into a single call to run() and segmenting the data on export. This is not implemented, but could be in the future if there is sufficient demand.

Let’s try the first of these options:

block1 = sim.to_neo()
block1.segments[0].name = 'trial 1'
sim.reset()
sim.run(50 * b2.ms)
block2 = sim.to_neo()
block2.segments[0].name = 'trial 2'

block = neo.Block()
block.segments.extend(block1.segments)
block.segments.extend(block2.segments)
block.groups.extend(block1.groups)
block.groups.extend(block2.groups)
block

Block with 2 segments, 2 groups
# segments (N=2)
0: Segment with 2 analogsignals, 58 spiketrains
   name: 'trial 1'
   # analogsignals (N=2)
   0: AnalogSignal with 2 channels of length 50; units mW/mm**2; datatype float64 
      name: 'Light'
      description: 'Exported from Cleo Light device'
      annotations: {'export_datetime': datetime.datetime(2023, 9, 6, 9, 59, 33, 375717)}
      sampling rate: 1.0 1/ms
      time: 0.0 ms to 50.0 ms
   1: AnalogSignal with 3 channels of length 50; units uV; datatype float64 
      name: 'Probe.TKLFPSignal'
      description: 'Exported from Cleo TKLFPSignal object'
      annotations: {'export_datetime': datetime.datetime(2023, 9, 6, 9, 59, 33, 376227)}
      sampling rate: 1.0 1/ms
      time: 0.0 ms to 50.0 ms
1: Segment with 2 analogsignals, 63 spiketrains
   name: 'trial 2'
   # analogsignals (N=2)
   0: AnalogSignal with 2 channels of length 50; units mW/mm**2; datatype float64 
      name: 'Light'
      description: 'Exported from Cleo Light device'
      annotations: {'export_datetime': datetime.datetime(2023, 9, 6, 9, 59, 34, 138750)}
      sampling rate: 1.0 1/ms
      time: 0.0 ms to 50.0 ms
   1: AnalogSignal with 3 channels of length 50; units uV; datatype float64 
      name: 'Probe.TKLFPSignal'
      description: 'Exported from Cleo TKLFPSignal object'
      annotations: {'export_datetime': datetime.datetime(2023, 9, 6, 9, 59, 34, 139149)}
      sampling rate: 1.0 1/ms
      time: 0.0 ms to 50.0 ms

The more attractive approach suggested by Neo’s API of block1.merge(block2) does not work, but could potentially be fixed in the future.

Example analysis use case#

Let’s try elephant, following their statistics tutorial:

import matplotlib.pyplot as plt
from elephant.statistics import instantaneous_rate, time_histogram, mean_firing_rate
import quantities as pq

fig, ax = plt.subplots()
st1 = block.segments[0].spiketrains[0]
histogram_rate = time_histogram([st1], 5*pq.ms, output='rate')
inst_rate = instantaneous_rate(st1, sampling_period=1*pq.ms)

# plotting the original spiketrain
ax.plot(
    st1,
    [0] * len(st1),
    "r",
    marker=2,
    ms=25,
    markeredgewidth=2,
    lw=0,
    label="poisson spike times",
)

# mean firing rate
ax.hlines(
    mean_firing_rate(st1),
    xmin=st1.t_start,
    xmax=st1.t_stop,
    linestyle="--",
    label="mean firing rate",
)

# time histogram
ax.bar(
    histogram_rate.times,
    histogram_rate.magnitude.flatten(),
    width=histogram_rate.sampling_period,
    align="edge",
    alpha=0.3,
    label="time histogram (rate)",
)

# instantaneous rate
ax.plot(
    inst_rate.times.rescale(pq.ms),
    inst_rate.rescale(histogram_rate.dimensionality).magnitude.flatten(),
    label="instantaneous rate",
)

# axis labels and legend
ax.set(
    xlabel=f"time [{st1.times.dimensionality.latex}]",
    ylabel=f"firing rate [{histogram_rate.dimensionality.latex}]",
    xlim=(st1.t_start, st1.t_stop),
)
ax.legend()

<matplotlib.legend.Legend at 0x7f527c923e20>
../_images/3ff0b20d39c6a16c111b08a76671ed26740a6f840f2c93e7f80d484984ae2e25.png