在 MATLAB 中计算 PSTH（周刺激时间直方图）的矢量化方法

Question

我有一个尖峰时间向量（来自神经元的动作电位）和一个刺激事件时间戳向量。我想创建一个 PSTH 来查看刺激是否影响神经元的尖峰率。我可以通过遍历每个刺激事件来做到这一点（参见下面的简单示例），但是对于有超过 30,000 个刺激事件并且正在记录许多神经元的长期实验来说，这非常慢。

如果没有 for 循环，如何做到这一点？

慢速方式示例：

% set variables
spikeTimes = [0.9 1.1 1.2 2.5 2.8 3.1];
stimTimes = [1 2 3 4 5];        
preStimTime = 0.2;
postStimTime = 0.3;
for iStim = 1:length(stimTimes)
    % find spikes within time window
    inds = find((spikeTimes > (stimTimes(iStim) - preStimTime)) & (spikeTimes < (stimTimes(iStim) + postStimTime)));
    % align spike times to stimulus onset
    stimONtimes = spikeTimes(inds) - stimTimes(iStim);
    % store times in array for plotting
    PSTH_array(iStim,1:length(stimONtimes)) = stimONtimes;
end

Answer 1

这是一个解决方案，在所有尖峰和两个假设上使用一个循环：

刺激时间间隔固定
刺激间隔大于 PSTH 间隔

假设刺激时间是固定的时间间隔：

delta_times = mean(diff(stimTimes));
assert(max(abs(diff(stimTimes)-delta_times))<1e-3);

现在我们将尖峰时间与第一个刺激之前的 preStimTime 对齐：

spikeTimes0 = spikeTimes - stimTimes(1) + preStimTime;

现在我们希望使用第二个假设为每个尖峰计算其计数的刺激：

assert((postStimTime-preStimTime)<dekta_times);
stimuli_index = floor(spikeTimes0 / delta_times);

相对于该刺激进行计算：

spike_time_from_stimuli = spikeTimes0 - stimuli_index*delta_times;

现在让我们构建 PSTH，精度为 0.01（与所有其他时间的单位相同）：

dt = 0.01;
times_around_stimuli = preStimTime:dt:postStimTime;
n_time_bins = length(times_around_stimuli);
n_stimuli = length(stimTimes);
PSTH = zeros(n_stimuli, n_time_bins)
for i=1:length(spikeTimes)
     time_index = ceil(spike_time_from_stimuli(i) / dt);
     % Ignore time-bins far from the event
     if time_index > n_time_bins
           continue;
     end
     PSTH(stimuli_index(i),time_index) = PSTH(stimuli_index(i),time_index) + 1;
end

Answer 2

最好的方法可能是只使用现有的直方图函数。它们非常很快，应该可以为您提供所需的所有信息。当然，这是假设 bin 不重叠。鉴于您的示例数据：

spikeTimes = [0.9 1.1 1.2 2.5 2.8 3.1];
stimTimes = [1 2 3 4 5];        
preStimTime = 0.2;
postStimTime = 0.3;

您可以像这样构造垃圾箱：

bins = sort([stimTimes - preStimTime, stimTimes + postStimTime])

或

bins = [stimTimes - preStimTime; stimTimes + postStimTime];
bins = bins(:).'

bins =
   0.80000   1.30000   1.80000   2.30000   2.80000   3.30000   3.80000   4.30000   4.80000   5.30000

然后您可以使用 histcounts, discretize, or histc，具体取决于您想要的结果和您拥有的 MATLAB 版本。我将使用 histc（因为我没有那么多花哨的东西）但是所有三个函数的输入都是相同的。 histcounts（edges，对我们没用）有一个额外的输出，discretize（实际计数）少一个输出。

[N, IDX] = histc(spikeTimes, bins)

N =    
   3   0   0   1   2   0   0   0   0   0

IDX =    
   1   1   1   4   5   5

由于 bins 包括 (T(i) + postStimTime) 和 (T(i+1) - preStimTime) 之间的时间，我们需要每隔一个 bin：

N = N(1:2:end)

N =
   3   0   2   0   0

同样，我们只对发生在奇数时隙的峰值感兴趣，我们需要调整索引以匹配新的 IDX:

v = mod(IDX, 2)

v =
   1   1   1   0   1   1

IDX = ((IDX+1)/2).*v

IDX =
   1   1   1   0   3   3

结果与我们最初得到的一致：bin 1 中有 3 个尖峰，bin 3 中有 2 个尖峰。

在 MATLAB 中计算 PSTH（周刺激时间直方图）的矢量化方法

Vectorized approach to compute PSTH (peristimulus time histogram) in MATLAB

matlab

vectorization

neuroscience