Note

Click here to download the full example code

# Fit injected current

## Learning objectives

- Learn how to explore spiking data and do basic analyses using pynapple
- Learn how to structure data for nemos
- Learn how to fit a basic Generalized Linear Model using nemos
- Learn how to retrieve the parameters and predictions from a fit GLM for intrepetation.

```
# Import everything
import jax
import os
import matplotlib.pyplot as plt
import nemos as nmo
import nemos.glm
import numpy as np
import pynapple as nap
import workshop_utils
# configure plots some
plt.style.use(workshop_utils.STYLE_FILE)
# Set the default precision to float64, which is generally a good idea for
# optimization purposes.
jax.config.update("jax_enable_x64", True)
```

## Data Streaming

- Stream the data. Format is Neurodata Without Borders (NWB) standard

```
path = workshop_utils.data.download_data("allen_478498617.nwb", "https://osf.io/vf2nj/download",
'../data')
```

## Pynapple

### Data structures and preparation

- Open the NWB file with pynapple

```
data = nap.load_file(path)
print(data)
```

`units`

: dictionary of neurons, holding each neuron's spike timestamps.`epochs`

: start and end times of different intervals, defining the experimental structure, specifying when each stimulation protocol began and ended.`stimulus`

: injected current, in Amperes, sampled at 20k Hz.`response`

: the neuron's intracellular voltage, sampled at 20k Hz.

```
trial_interval_set = data["epochs"]
# convert current from Ampere to pico-amperes, to match the above visualization
# and move the values to a more reasonable range.
current = data["stimulus"] * 1e12
spikes = data["units"]
```

`trial_interval_set`

: dictionnary of start and end times of different intervals, defining the experimental structure, specifying when each stimulation protocol began and ended.

```
trial_interval_set.keys()
```

`Noise 1`

: epochs of random noise

```
noise_interval = trial_interval_set["Noise 1"]
noise_interval
```

- Let's focus on the first epoch.

```
noise_interval = noise_interval.loc[[0]]
noise_interval
```

`current`

: Tsd (TimeSeriesData) : time index + data

```
current
```

`restrict`

: restricts a time series object to a set of time intervals delimited by an IntervalSet object

```
current = current.restrict(noise_interval)
current
```

`TsGroup`

: a custom dictionnary holding multiple`Ts`

(timeseries) objects with potentially different time indices.

```
spikes
```

We can index into the `TsGroup`

to see the timestamps for this neuron's
spikes:

```
spikes[0]
```

Let's restrict to the same epoch `noise_interval`

:

```
spikes = spikes.restrict(noise_interval)
print(spikes)
spikes[0]
```

Let's visualize the data from this trial:

```
fig, ax = plt.subplots(1, 1, figsize=(8, 2))
ax.plot(current, "grey")
ax.plot(spikes.to_tsd([-5]), "|", color="k", ms = 10)
ax.set_ylabel("Current (pA)")
ax.set_xlabel("Time (s)")
```

### Basic analyses

The Generalized Linear Model gives a predicted firing rate. First we can use pynapple to visualize this firing rate for a single trial.

`count`

: count the number of events within`bin_size`

```
# bin size in seconds
bin_size = 0.001
count = spikes.count(bin_size)
count
```

Let's convert the spike counts to firing rate :

`smooth`

: convolve with a Gaussian kernel

```
# the inputs to this function are the standard deviation of the gaussian and
# the full width of the window, given in bins. So std=50 corresponds to a
# standard deviation of 50*.001=.05 seconds
firing_rate = count.smooth(std=50, size=1000)
# convert from spikes per bin to spikes per second (Hz)
firing_rate = firing_rate / bin_size
```

Note that firing_rate is a `TsdFrame`

!

```
print(type(firing_rate))
```

```
# we're hiding the details of the plotting function for the purposes of this
# tutorial, but you can find it in the associated github repo if you're
# interested:
# https://github.com/flatironinstitute/nemos-workshop-feb-2024/blob/binder/src/workshop_utils/plotting.py
workshop_utils.plotting.current_injection_plot(current, spikes, firing_rate)
```

What is the relationship between the current and the spiking activity?
`compute_1d_tuning_curves`

: compute the firing rate as a function of a 1-dimensional feature.

```
tuning_curve = nap.compute_1d_tuning_curves(spikes, current, nb_bins=15)
tuning_curve
```

Let's plot the tuning curve of the neuron.

```
workshop_utils.plotting.tuning_curve_plot(tuning_curve)
```

## Nemos

### Preparing data

- Get data from pynapple to nemos-ready format:
- predictors and spikes must have same number of time points

```
binned_current = current.bin_average(bin_size)
print(f"current shape: {binned_current.shape}")
# rate is in Hz, convert to KHz
print(f"current sampling rate: {binned_current.rate/1000.:.02f} KHz")
print(f"\ncount shape: {count.shape}")
print(f"count sampling rate: {count.rate/1000:.02f} KHz")
```

- predictors must be 2d, spikes 1d

```
# add singleton dimensions for axis 1.
predictor = np.expand_dims(binned_current, 1)
# grab the spike counts for our one neuron
count = count[:, 0]
# check that the dimensionality matches nemos expectation
print(f"predictor shape: {predictor.shape}")
print(f"count shape: {count.shape}")
```

- predictors and spikes must be jax arrays

```
predictor = jax.numpy.asarray(predictor.values)
count = jax.numpy.asarray(count.values)
```

### Fitting the model

- GLM objects need regularizers and observation models

```
model = workshop_utils.model.GLM(regularizer=nmo.regularizer.UnRegularized(solver_name="LBFGS"))
```

- call fit and retrieve parameters

```
model.fit(predictor, count)
```

```
print(f"firing_rate(t) = exp({model.coef_} * current(t) + {model.intercept_})")
```

```
print(f"coef_ shape: {model.coef_.shape}")
print(f"intercept_ shape: {model.intercept_.shape}")
```

- generate and examine model predictions.

```
predicted_fr = model.predict(predictor)
# convert units from spikes/bin to spikes/sec
predicted_fr = predicted_fr / bin_size
# let's reintroduce the time axis by defining a TsdFrame.
# we must convert the firing rate to a numpy array (from jax.numpy) to make it
# pynapple compatible
predicted_fr = nap.TsdFrame(t=binned_current.t, d=np.asarray(predicted_fr))
# and let's smooth the firing rate the same way that we smoothed the smoothed
# spike train
smooth_predicted_fr = predicted_fr.smooth(50, 1000)
# and plot!
workshop_utils.plotting.current_injection_plot(current, spikes, firing_rate,
# plot the predicted firing rate that has
# been smoothed the same way as the
# smoothed spike train
predicted_firing_rate=smooth_predicted_fr)
```

- what do we see?

```
# compare observed mean firing rate with the model predicted one
print(f"Observed mean firing rate: {np.mean(count) / bin_size} Hz")
print(f"Predicted mean firing rate: {np.mean(predicted_fr)} Hz")
```

- examine tuning curve -- what do we see?

```
tuning_curve_model = nap.compute_1d_tuning_curves_continuous(predicted_fr, current, 15)
fig = workshop_utils.plotting.tuning_curve_plot(tuning_curve)
fig.axes[0].plot(tuning_curve_model, color="tomato", label="glm")
fig.axes[0].legend()
```

### Finishing up

- Finally, let's look at spiking and scoring/metrics

```
spikes = jax.random.poisson(jax.random.PRNGKey(0), predicted_fr.values)
```

```
log_likelihood = model.score(predictor, count, score_type="log-likelihood")
print(f"log-likelihood: {log_likelihood}")
```

```
model.score(predictor, count, score_type='pseudo-r2-Cohen')
```

## Further Exercises

- what else can we do?

### Other stimulation protocols

### Train and test sets

### Model extensions

## Citation

The data used in this tutorial is from the Allen Brain Map, with the following citation:

Contributors: Agata Budzillo, Bosiljka Tasic, Brian R. Lee, Fahimeh Baftizadeh, Gabe Murphy, Hongkui Zeng, Jim Berg, Nathan Gouwens, Rachel Dalley, Staci A. Sorensen, Tim Jarsky, Uygar Sümbül Zizhen Yao

Dataset: Allen Institute for Brain Science (2020). Allen Cell Types Database -- Mouse Patch-seq [dataset]. Available from brain-map.org/explore/classes/multimodal-characterization.

Primary publication: Gouwens, N.W., Sorensen, S.A., et al. (2020). Integrated morphoelectric and transcriptomic classification of cortical GABAergic cells. Cell, 183(4), 935-953.E19. https://doi.org/10.1016/j.cell.2020.09.057

Patch-seq protocol: Lee, B. R., Budzillo, A., et al. (2021). Scaled, high fidelity electrophysiological, morphological, and transcriptomic cell characterization. eLife, 2021;10:e65482. https://doi.org/10.7554/eLife.65482

Mouse VISp L2/3 glutamatergic neurons: Berg, J., Sorensen, S. A., Miller, J., Ting, J., et al. (2021) Human neocortical expansion involves glutamatergic neuron diversification. Nature, 598(7879):151-158. doi: 10.1038/s41586-021-03813-8

**Total running time of the script:** ( 0 minutes 0.000 seconds)

Download Python source code: 01_current_injection_code.py