Skip to content

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

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

Pynapple

Data structures and preparation

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

Download Jupyter notebook: 01_current_injection_code.ipynb

Gallery generated by mkdocs-gallery