Векторизованный подход к вычислению PSTH (гистограмма времени перистимула) в MATLAB

Question

У меня есть вектор периодов спайка (потенциалы действия от нейрона) и вектор временных меток события стимула. Я хочу создать PSTH, чтобы увидеть, влияет ли стимул на скорость всплеска нейрона. Я могу сделать это, перейдя через каждое событие стимула (см. Простой пример ниже), но это очень медленно для длительных экспериментов, где происходит более 30 000 событий стимула и регистрируется много нейронов.

Как это можно сделать без петель?

Пример медленного пути:

% 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

Лучший способ это, вероятно, просто использовать существующие функции гистограммы. Они очень быстро и должны предоставить вам всю необходимую информацию. Конечно, это предполагает, что ящики не перекрываются. Учитывая ваш пример данных:

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 или 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 

Поскольку контейнеры включают в себя время между (T(i) + postStimTime) и (T(i+1) - preStimTime) мы должны взять любой другой бункер:

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 

Полученные результаты согласуются с тем, что мы Первоначально: в бункере 1 и 2 пики в 3-х местах.

Answer 2

Вот решение с одной петлей над всеми шипами и двух предположений:

• стимулы раз в фиксированные интервалы
• стимулы интервалы больше, чем интервал 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, в presicion 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