guide for using NEMOS for fitting a GLM to firing rate during different behaviors

[WeissShahaf]

behavior responsive cells pose an interesting challange. a normal GLM cannot handle this, maybe if we turn the behavior epochs to binomial continuous variables (on /off for every FR bin).

can you publish a guide on how to use NEMOS in this case?

[BalzaniEdoardo]

Hello and thanks for your message!

I may need a bit more context to understand what type of model you would want to setup, but I can try to guess.

We have ways to model responses to individual events. Modeling a temporal predictor that model a response to the event "start epoch" and "end epoch" can be done by creating a one-hot vector that is 1 when the event happens and 0 otherwise. Then, this vector should be convolved with a basis to create the model design matrix. This would be equivalent of using the spike times as predictors, as you first count the spikes, and then convolve the counts, see this tutorial for example).

One can also create the boolean predictor that you are describing and convolve that instead. Unfortunately we don't have a tutorial showing how to do that, but it is fairly simple to do with pyanapple. An example like that would be a good addition, and it is probably worth an how-to-guide. For now, I just post below.

import numpy as np
import pynapple as nap
import nemos as nmo
import matplotlib.pyplot as plt

# some start/end of stimuli
stimulus_a = nap.IntervalSet([0, 5, 10], [3, 8, 13])
stimulus_b = nap.IntervalSet([1, 7, 11], [2, 8, 12])

# simulate some spike times for a neuron
spk = nap.TsGroup({1: nap.Ts(np.sort(np.random.uniform(0, 15, 20)))})

# count spikes
counts = spk.count(0.01)
n_samples, n_neurons = counts.shape

# initialize predictors
predictors = nap.TsdFrame(t=counts.t, d=np.zeros((n_samples, 2)), columns = ['stimulus A', 'stimulus B'])

# Check whether there is a flash within a given bin of spikes
idx_a = stimulus_a.in_interval(counts)
idx_b = stimulus_b.in_interval(counts)

# Replace everything that is not nan with 1 in the corresponding column
predictors.d[~np.isnan(idx_a), 0] = 1
predictors.d[~np.isnan(idx_b), 1] = 1

plt.subplot(211)
plt.title("stimulus")
plt.plot(predictors)
plt.subplot(212)
plt.title("counts")
plt.plot(counts)
plt.tight_layout()
plt.show()


# define a basis and convolve
basis = nmo.basis.RaisedCosineLogConv(4, window_size=80)

# design
X = basis.compute_features(predictors)

Let me know if I interpreted correctly, otherwise could you provide an synthetic data example of a stimulus and explain a bit more on how you would go about constructing the glm design?

Thanks again,
Edoardo

[WeissShahaf]

thanks for the reply.
I have a dozen or more annotated behaviors from BORIS per session.

let's say i only take two: escape and pursuit.
then the code in your example will look like this:

import numpy as np
import pynapple as nap
import nemos as nmo
import matplotlib.pyplot as plt
def main():
"""
Main function to run the GLM analysis on experimental neural data using the nemos package.
The script performs model selection using L1 regularization to identify the minimum
set of behavioral predictors for a neuron's firing rate.
"""
# 1. Load and preprocess experimental data
# This section is adapted from your original script to load your specific
# dataset and prepare it for the GLM analysis.

print('Step 1: Loading and Preprocessing Data...')

# --- Data Loading ---
# Note: Please ensure these file paths are correct for your system.
try:
    # These paths point to your data files. The script will fail if they are not found.
    firing_rates = np.load(f"//gpfs.corp.brain.mpg.de/stem/data/project_hierarchy/Figures/20250509/afm16924/concat_PSTHs/firing_rates_240522.npy")
    behaviour_df = pd.read_pickle(f"//gpfs.corp.brain.mpg.de/stem/data/project_hierarchy/Figures/20250509/afm16924/concat_PSTHs/behaviors_240522.pkl")
    print('Data loaded successfully.')
    spk_times=np.load((r'D:\GitHub\NewGit\Project_hierarchy_stem\Analysis\neuron_tuning\GLM\example_spikes_344.npz'))
    spk_times=spk_times['spk_times']
    
except FileNotFoundError as e:
    print(f"Error: Data file not found at the specified path.")
    print(f"Details: {e}")
    print("Please ensure the file paths are correct and accessible.")
    return # Exit the script if data cannot be loaded

# --- Data Preprocessing ---
print('Preprocessing Data...')
vframerate = 50
dt = 1.0 / vframerate
i_neuron = 195  # The specific neuron to analyze

# Select the firing rate for the chosen neuron. Ensure i_neuron is a valid index.
if i_neuron >= firing_rates.shape[0]:
    print(f"Warning: i_neuron index {i_neuron} is out of bounds. Using neuron 0 instead.")
    i_neuron = 0
firing_rate = firing_rates[i_neuron]
num_bins = len(firing_rate)
time_ax = np.linspace(0, num_bins * dt, num_bins)

# Convert firing rate (Hz) to spike counts for the Poisson model.
# The nemos GLM works with spike counts, which are modeled as a Poisson process.
spike_train = (firing_rate * dt).astype(np.int32)
spk = nap.TsGroup({1: nap.Ts(np.sort(spk_times))})

# Create binary predictors for each unique behavior from the 'behaviours' column
unique_behaviours = behaviour_df['behaviours'].unique()
binary_behaviors = {}
print(f"Creating binary predictors for behaviors: {unique_behaviours}")

predictors_list=[]

for beh in unique_behaviours:
    binary_signal = np.zeros(len(time_ax), dtype=int)
    predictor_dic={}
    start_times=[]
    end_times=[]
    
    # Filter the dataframe for the current behavior
    subdf = behaviour_df[behaviour_df["behaviours"] == beh].sort_values("frames_s").reset_index(drop=True)
    
    i = 0
    while i < len(subdf):
        row = subdf.iloc[i]
        event_type = row["start_stop"]
        event_time = row["frames_s"]
        
        start_time = event_time
        
        if event_type == "POINT":
            # For point events, mark a very short duration (one time bin)
            end_time = event_time + dt
            i += 1
        elif event_type == "START":
            # For interval events, find the corresponding STOP
            if i + 1 < len(subdf) and subdf.iloc[i + 1]["start_stop"] == "STOP":
                end_time = subdf.iloc[i + 1]["frames_s"]
                i += 2 # Move past the START and STOP pair
            else:
                # If no STOP is found, treat as a single time-bin event
                end_time = start_time + dt
                i += 1
        else: # Unpaired STOP or other event type
            i += 1
            continue
        start_times.append(start_time)
        end_times.append(end_time)

        # Set the binary signal to 1 for the duration of the event
        idx = np.where((time_ax >= start_time) & (time_ax < end_time))[0]
        
        binary_signal[idx] = 1
    
    binary_behaviors[beh] = binary_signal
    predictor_dic[beh]={'binary_signal':binary_signal,
                        'start_times': start_times,
                        'end_times': end_times}
    
    predictors_list.append(predictor_dic)

# --- Build the final design matrix (X) and response vector (y) ---
# Create a DataFrame to hold all variables, then extract for nemos.
predictor_df = pd.DataFrame(binary_behaviors)

# Design matrix X: contains all predictors (behaviors)
X = predictor_df.values

# count spikes
counts = spk.count(0.01)
n_samples, n_neurons = counts.shape

# initialize predictors
#predictors = nap.TsdFrame(t=counts.t, d=np.zeros((n_samples, 2)), columns = unique_behaviours) # all behaviors(14)
predictors = nap.TsdFrame(t=counts.t, d=np.zeros((n_samples, 2)), columns = ['escape','pursuit']) # sample subset

# count spikes
counts = spk.count(0.01)
n_samples, n_neurons = counts.shape

# initialize predictors
escape=predictors_list[10]#escape
pursuit=predictors_list[2]#pursuit   

# some start/end of stimuli
start_times= escape['escape']['start_times']
end_times = escape['escape']['end_times']
stimulus_a = nap.IntervalSet([start_times], [end_times])

start_times= pursuit['pursuit']['start_times']
end_times = pursuit['pursuit']['end_times']
stimulus_b = nap.IntervalSet([start_times], [end_times])

# Check whether there is a flash within a given bin of spikes
idx_a = stimulus_a.in_interval(counts)
idx_b = stimulus_b.in_interval(counts)

# Replace everything that is not nan with 1 in the corresponding column
predictors.d[~np.isnan(idx_a), 0] = 1
predictors.d[~np.isnan(idx_b), 1] = 1

plt.subplot(211)
plt.title("stimulus")
plt.plot(predictors)
plt.subplot(212)
plt.title("counts")
plt.plot(counts)
plt.tight_layout()
plt.show()


# define a basis and convolve
basis = nmo.basis.RaisedCosineLogConv(4, window_size=80)

# design
X = basis.compute_features(predictors)

if name == 'main':
main()

## but how do i go from there?

in the tutorial the basis and model are refereing to the firing. and here the basis is refereing to the behavior