Possion-process Based Matlab Code for Detecting Deviated Activity From Neuronal (Spike) Data

function [popdev,devseq] = devFromPoisson(acti, sr, binw, stp, eventT )
% suppose you recorded neuronal activities during a long-term free-behavior
% sequence (e.g., 1 hour), and you also recorded the start time for several types of behaviors
% （e.g, licking, saccading, shaking and sleeping). you want to know whether
% the neurons you recorded are mostly correlated with a specific behavior
% type (e.g., licking)，you can refer
% https://doi.org/10.1016/j.cell.2024.02.032. 
% in this paper, Selmaan N. Chettih et al. (Cell 2024) find the behavior type corresponding to the recorded neurons based on a possion process method.
% this matlab code shows how to repeat the method used by Selmann.

% INPUT:
% acti: neurons x activity sequence
% sr: sampling rate of neuron activities, in seconds
% binw: bin width for counting spikes or computing Calcium signal, e.g.,
% 0.5 s
% stp: for sliding window, e.g., 0.1s
% eventT: a cell, each one corresponding to a behavior type. each cell contains rows for the offset time of behavior events, in seconds, relative to the
% begining of the sequence.
% OUTPUT:
% popdev: a cell, each one corresponding to a behavior type. the elements
% corresponding to population deviation values.
% devseq: neurons x deviation. deviation index for each neuron, unit:std (see the reference paper)

% WARNING: I didn't test this matlab code with my neuron data. if you want
% to use this code, please debug it first. But I think I followed the mathmatical thought in the paper
% correctly.

dur = size(acti,2)/sr;% recording time
times = 0:1/sr:dur;% time ticks
devseq = [];% for each neuron， deviation sequence
lamda = [];% lamada for each neuron

for n=1:size(acti,1)
    lamda(n) = sum(acti(n,:))/dur;
end

for t=0:stp:dur
    idxs = find(times>=t&times<t+binw);
    for n=1:size(acti,1)
        prob = poisspdf(sum(acti(n,idxs))/binw,lamda(n));
        devseq(n,t/stp+1) = -log(prob);
    end
end
% plot(times,mean(dev,1),'k-'); hold on; % deviation of pupulation activity
popdev = cell(length(eventT),1);
for e=1:length(eventT)
    wins = [eventT{e}-0.5 eventT{e}+0.5]; % window of each event for computing dev
    dev = [];
    for w=1:length(wins)
        idxs = find(times>=wins(w,1)&times<=wins(w,2));
        for n=1:size(acti,1)
            prob = poisspdf(sum(acti(n,idxs)),lamda(n));
            dev(n,w) = -log(prob);
        end
    end
    popdev{e} = mean(dev,1);
end

% figure;  % figure1E in the paper
% for e=1:length(eventT)
%     x = popdev{e};
%     cdfplot(x); hold on;
% end


end

Neurons fire spikes during free behavior sequence, how to know the neurons are more correlated with a certain type of behavior events. I recommend you to read the article by Selmaan N. Chettih et al. (Cell 2024) . I repeated their analysis method in this Matlab code. Hope this can help some students.