Sunday, November 25, 2012

Finding AHPs



AHP = average y-distance between the blue and the green dots.

Since apparently in my last action potentials post it was not clear what I was looking for, I decided to post the above graph.

I found voltage thresholds (blue dots) and action potentials' troughs (green squares) to find the average difference between them (after-hyperpolarization potential). I did not find a good automated way of finding voltage threshold, so I just hard-coded the minimal value for dv/dt: once the trace crossed it, the voltage at which it was crossed was designated as voltage threshold. If you look at the graph, you can see that it worked rather well, but I would still like to find it in a more elegant way (I tried finding local third-derivative maxima, but had some problems with the function).

The troughs were obtained through findpeaks function (Signal Processing Toolbox package), such that:
[p, i_p] = findpeaks (f);
saves the values of peaks as p and their indexes as i_p. Since in my case I was looking for troughs, I inverted my trace (and then inverted the peaks back):
[troughs, i_t] = findpeaks (trace);
troughs = –troughs;
If you want to plot the troughs, you can just write:
plot (i_t, troughs, 'og'); %green circles
plot (trace, 'r'); %red trace
dv/dt was obtained through function diff, such that:
dn_dx = diff (f, n); %produces n-th derivative of f
In order to analyze the data, I first had to separate out individual action potentials:

function [output] = findAPs (input_trace)

    barrage = input_trace; %select APs to analyze from the trace

    %[peaks, i_p] = findpeaks(barrage); %if you want to analyze peaks

    [~, i_t] = findpeaks (-barrage); %obtain x-values for troughs

    APs = cell(length(i_t), 1); %allocate empty cell space

    %adding first and last pts of barrage to minima:
    minima = [1; i_t; length(barrage)];

    %separating APs:
    for i = 1:length(minima)-1
        APs{i} = barrage(minima(i):minima(i+1));
    end

    output = APs; %return a cell of AP traces
 end
Then I found dv/dt for each AP and selected for the first value above the accepted threshold level to find the threshold:

function [output] = find_v_th (APs)
    v_th = struct([]);

    for i = 1:length(APs)
        dv = diff (APs{i}, 1); %first derivative

        %find first instance of dv/dt>0.001:
        v_th_i = find (dv>0.001, 1, 'first')-1; %offset by one b/c of diff

        if (v_th_i) %if v-threshold was found, enter voltage and x vals
            v_th(i).v = APs{i}(v_th_i);
            v_th(i).x = v_th_i;
        else
            v_th(i).v = [];
            v_th(i).x = [];
        end
    end

    output = v_th;
end
Plotting dv/dt was the same as plotting troughs. To find AHP, one can create an array of trough values and threshold values and subtract one from another, or just average each array and subtract the averages. To plot the AHPs and dv/dt's took a little more effort, since I had to put the APs and dv/dt's with AHP values back together in one barrage, making sure to preserve the x-values. I recommend doing all this in object-oriented manner.

No comments: