Fit V1 cell

Learning objectives

Learn how to combine GLM with other modeling approach.
Review previous tutorials.

import jax
import math
import os
import matplotlib.pyplot as plt
import nemos as nmo
import numpy as np
import pynapple as nap
import requests
import tqdm
import workshop_utils

# required for second order methods (BFGS, Newton-CG)
jax.config.update("jax_enable_x64", True)
# configure plots some
plt.style.use(workshop_utils.STYLE_FILE)

Data Streaming

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

Pynapple

data = nap.load_file(path)

print(data)

epochs = data["epochs"]
spikes = data["units"]
stimulus = data["whitenoise"]

stimulus is white noise shown at 40 Hz
white noise is a good stimulus for mapping basic stimulus properties of V1 simple cells

fig, ax = plt.subplots(1, 1, figsize=(12,4))
ax.imshow(stimulus[0], cmap='Greys_r')
stimulus.shape

print(spikes)
spikes = spikes[[34]]

goal is to predict the neuron's response to this white noise stimuli
several ways we could do this, what do you think?

Spike-triggered average

compute spike-triggered average to visualize receptive field.

sta = nap.compute_event_trigger_average(spikes, stimulus, binsize=0.025,
                                        windowsize=(-0.15, 0.0))

sta

sta[1, 0]

visualize spike-triggered average and decide on our spatial filter.

fig, axes = plt.subplots(1, len(sta), figsize=(3*len(sta),3))
for i, t in enumerate(sta.t):
    axes[i].imshow(sta[i,0], vmin = np.min(sta), vmax = np.max(sta),
                   cmap='Greys_r')
    axes[i].set_title(str(t)+" s")

receptive_field = np.mean(sta.get(-0.125, -0.05), axis=0)[0]

fig, ax = plt.subplots(1, 1, figsize=(4,4))
ax.imshow(receptive_field, cmap='Greys_r')

use the spike-triggered average to preprocess our visual input.

# element-wise multiplication and sum
print((receptive_field * stimulus[0]).sum())
# dot product of flattened versions
print(np.dot(receptive_field.flatten(), stimulus[0].flatten()))

filtered_stimulus = np.einsum('t h w, h w -> t', stimulus, receptive_field)
# add the extra dimension for feature
filtered_stimulus = np.expand_dims(filtered_stimulus, 1)

fig, ax = plt.subplots(1, 1, figsize=(12,4))
ax.plot(filtered_stimulus)

Preparing data for nemos

get counts and filtered_stimulus into proper shape for nemos

# grab spikes from when we were showing our stimulus, and bin at 1 msec
# resolution
bin_size = .001
counts = spikes.restrict(filtered_stimulus.time_support).count(bin_size)
print(counts.rate)
print(filtered_stimulus.rate)

print(counts[:5])
print(filtered_stimulus[:5])

filtered_stimulus = workshop_utils.data.fill_forward(counts, filtered_stimulus)
filtered_stimulus

Set up the basis and prepare the temporal predictor for the GLM.

basis = nmo.basis.RaisedCosineBasisLog(8)
window_size = 100
time, basis_kernels = basis.evaluate_on_grid(window_size)
time *= bin_size * window_size
convolved_input = nmo.utils.convolve_1d_trials(basis_kernels, filtered_stimulus)
# convolved_input has shape (n_time_pts, n_features, n_basis_funcs), and
# n_features is the singleton dimension from filtered_stimulus, so let's
# squeeze it out:
convolved_input = np.squeeze(convolved_input)
# and, as also described in the head direction tutorial, when doing this we
# need to remove the first window_size time points from the neuron counts and
# the last time point from the convolved input:
counts = counts[window_size:]
convolved_input = convolved_input[:-1]
# and grab the counts for our single neuron
counts = counts[:, 0]

Fitting the GLM

Fit the GLM

model = workshop_utils.model.GLM(regularizer=nmo.regularizer.UnRegularized(solver_name="LBFGS"))
model.fit(convolved_input, counts)

Examine the resulting temporal filter

temp_weights = np.einsum('b, t b -> t', model.coef_, basis_kernels)
plt.plot(time, temp_weights)

Fit V1 cell

Learning objectives

Data Streaming

Pynapple

Spike-triggered average

Preparing data for nemos

Fitting the GLM

Further exercises