How do I set the threshold value on the action potential on the upward and downward slope to define width?

Hi,

I am doing action potential analysis and at the moment I am defining my threshold value manually using ginput for analysing width, however I would like to set identical y values for analysis of width of the curve. Right now, I zoom into the figure to try and set the x and y points to get the width of the AP.

My code for producing the figure and setting ginput values manually for the defined variables:

% Roksana's script to plot .smr AP data 
% 20210127

% For troubleshooting analysis - cleans the workspace 
clear
close all

%Import data
load('C:\xyz.mat')

%Define channel to plot
data = xyz;
vm = data.values;

%scale x-axis by interval
t = linspace(0, length(vm)*data.interval,length(vm));

plot(t,vm)
xlabel('Time (s)')
ylabel('Voltage (mV)')
title('Membrane potential')

%find max Vm value between 2 points selected on the figure
[x, y] = ginput(2);

%Report values towards AP width
xfirst = x(1,1)
xsecond = x(2,1)

%Round values for calculations of parameters
xpt1 = round(x(1));
xpt2 = round(x(2));
ypt1 = round(y(1));
ypt2 = round(y(2));

%Find a point between the first and second points on x-axis
isin = (t >= xpt1 & t <= xpt2);

%Identify the peak between the two points on the x-axis
%peak = max(vm(isin))

peak = max(vm)

%Display trough (mV); from baseline
APtrough = min(vm(isin))

%%Display AP amplitude (mV)
APamp = (peak + (ypt1*(-1)))

%%Display half-amplitude(mV)
APhalfamp = ((peak-APtrough)*0.5)

%%%Area under the curve and width values from threshold value - require zoom%%%
%pause script until button press for zooming in
zoom on;
pause ();
zoom off;

%find max Vm value between 2 points selected on the figure
[a, b] = ginput(2);

%Report values towards AP width
thresholdt1 = a(1,1)
thresholdt2 = a(2,1)
thresholdValue1 = b(1)
thresholdValue2 = b(2)

%%Display width(ms)
width = (thresholdt2 - thresholdt1)*1000  %%%%needs to be defined by threshold

%%Display half-width(ms)
APhalfwidth = (width/2)

[darova]

YOu want to detect those points automatically. Is this correct?

Can you please show some more pictures?

[Roksana Khalid]

Yes, those points.

I tried using a grid for more accuracy but it takes so long and it's not perfect. For example, the last run produced these values:

thresholdt1 = 0.0156

thresholdt2 = 0.0179

[Roksana Khalid]

Ultimately, being able to get the threshold value at t1 and it's equivalent at t2, I think I would also be able to find area under the curve.

[KALYAN ACHARJYA]

Has the issue been fixed?

[Roksana Khalid]

not yet

[Mathieu NOE]

hello Roksana

can you share your data file pls ?

tx

[Roksana Khalid]

Here is the .mat file I was using.

Answer 1

[Star Strider]

I am not certain how robust this is, since we only have one record. However it could likely be adapted easily to other records if they are similar.

D = load('xyz.mat');
F = fieldnames(D);
C = struct2cell(D);

Ch1s = C{1}.values;                                         % Signal Vector
L1 = numel(Ch1s);                                           % Length
Ch1t = linspace(0, L1, L1).'*C{1}.interval;                 % Time Vector

[CPL,m,b] = ischange(Ch1s, 'linear', 'Threshold',10);       % Requires R2017b Or Later
cp1 = find((m>1), 1, 'first');
cp2 = find((m<-1), 1, 'last');

figure
plot(Ch1t, Ch1s)
hold on
% plot(Ch1t, m)
plot(Ch1t(cp1), Ch1s(cp1), '^r')
plot(Ch1t(cp2), Ch1s(cp2), 'vr')
hold off
grid
xlabel('Time (s)')
ylabel('Voltage (mV)')
title('Membrane potential')
pw = Ch1t(cp2) - Ch1t(cp1);                                 % Pulse Width
text(Ch1t(cp1), Ch1s(cp1), '\rightarrow',  'HorizontalAlignment','right', 'VerticalAlignment','middle', 'FontWeight','bold')
text(Ch1t(cp2), Ch1s(cp2), sprintf('\\leftarrow Width = %.4f', pw),  'HorizontalAlignment','left', 'VerticalAlignment','middle', 'FontWeight','bold')

producing this plot:

EDIT — (2 Feb 2021 at 19:28)

My code using the newly-posted ‘AP2.mat’:

My code is otherwise unchanged, except for importing the new data file, and increasing the precision of the sprintf numeric field.

.

[Roksana Khalid]

So if I zoom into the AP, the firing threshold can actually be seen at the kink. It's an abrupt change in the slope and this point on the y-axis is what I want to set on the rising and falling sides of the curve. It's approximately -60.8mV.

[Star Strider]

I doubt that any function could be that accurate in locating that particular point, especially being robust enough to do so on any number of depolarisations.

If you want to define the firing threshold as being -60 mV, it is then straightforward to find the times on both sides of the action potential that correspond to that value:

D = load('AP2.mat'); 
F = fieldnames(D); 
C = struct2cell(D); 

Ch1s = C{1}.values;                                         % Signal Vector 
L1 = numel(Ch1s);                                           % Length 
Ch1t = linspace(0, L1, L1).'*C{1}.interval;                 % Time Vector

Ch1sf = sgolayfilt(Ch1s,  3, 49);                           % Essentially A Lowpass Filter

[max,idx] = max(Ch1s);                                      % Peak & Index
thrshld = -60;                                              % Defined Threshold
tc(1) = interp1(Ch1sf(1:idx), Ch1t(1:idx), thrshld);        % Start Deflection Time
tc(2) = interp1(Ch1sf(idx:end), Ch1t(idx:end), thrshld);    % Stop Deflection Time
pw = tc(2)-tc(1);                                           % Pulse Width

figure
plot(Ch1t, Ch1s)
hold on
plot(tc(1), -60, '+r')
plot(tc(2), -60, '+r')
hold off
grid
xlabel('Time (s)')
ylabel('Voltage (mV)')
title('Membrane potential')
text(tc(1), -60, '\rightarrow',  'HorizontalAlignment','right', 'VerticalAlignment','middle', 'FontWeight','bold')
text(tc(2), -60, sprintf('\\leftarrow Width = %.5f', pw),  'HorizontalAlignment','left', 'VerticalAlignment','middle', 'FontWeight','bold')

With that, the pulse width is about 0.00226 seconds.

.

[Roksana Khalid]

This is super ^^

Thanks so much!

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.

Answer 2

[Mathieu NOE]

this is my first attempt to make it fully automatic

based on this excellent submission (attached function) : Piecewise linear interpolation

so I try to fit a 5 segment of linear segments to your data, and the rest comes automatically

other methods could also work (e.g. based on derivatives , but this can be sensitive to noise)

% Roksana's script to plot .smr AP data 
% 20210127
% For troubleshooting analysis - cleans the workspace 
clear
close all

%Import data
load('xyz.mat')

%Define channel to plot
% data = xyz;
% vm = data.values;

% S2C1_17082020_WM_apomorphine_20uM_Ch1 = 
% 
%   struct with fields:
% 
%        title: 'voltage'
%      comment: 'No comment'
%     interval: 4.0000e-05
%        scale: 0.0153
%       offset: 0
%        units: 'mV'
%        start: 0
%       length: 1122
%       values: [1122×1 double]
      
vm = S2C1_17082020_WM_apomorphine_20uM_Ch1.values;
samples = length(vm);
dt = S2C1_17082020_WM_apomorphine_20uM_Ch1.interval;
t_final = samples*dt;

%scale x-axis by interval
% t = linspace(0, length(vm)*data.interval,length(vm));
t = linspace(0, t_final,samples);

figure(1);
plot(t,vm)
xlabel('Time (s)')
ylabel('Voltage (mV)')
title('Membrane potential')

% five linear segments approximation 
    % first select 3 of the 5 points : start , peak (max) , end
    [vm_max,ind] =max(vm);
    x_vm_max = t(ind);

    % Init of fminsearch 
    xdata = t(:);ydata = vm(:);
    a = (max(xdata)-min(xdata))/4;    b = a*2;
    global xdata ydata x_vm_max vm_max
    x0 = [a,b];
    x = fminsearch(@obj_fun, x0);
    a_sol = x(1);    b_sol = x(2);
    XI = [min(xdata),a_sol,x_vm_max,b_sol,max(xdata)]; % vector of 1-D look-up table "x" points
    YI = lsq_lut_piecewise( xdata, ydata, XI ); % obtain vector of 1-D look-up table "y" points
    % plot fit
    figure(2);
    plot(xdata,ydata,'.',XI,YI,'+-')
    legend('experimental data (x,y(x))','LUT points (XI,YI)')
    title('Piecewise 1-D look-up table least square estimation')

% selected threshold = max value of the y values corresponding to x = a_sol and x =  b_sol
vm1 = interp1(XI,YI,a_sol);
vm2 = interp1(XI,YI,b_sol);
threshold = ceil(max(vm1,vm2)); % let's round it to the upper integer value

% finally, selected y values above "threshold" are 
ind = find(vm>threshold);
vm_select = vm(ind);
t_select = t(ind);

figure(3);
plot(t,vm,t_select,vm_select,'+-')
legend('experimental data','selected data')
xlabel('Time (s)')
ylabel('Voltage (mV)')
title('Membrane potential')

%%Display AP amplitude (mV)
APamp = (vm_max - (vm1+vm2)/2)  % compute difference between peak and "baseline" at same x location =  average of vm1 and vm2
%%Display half-amplitude(mV)
% APhalfamp = ((peak-APtrough)*0.5) % not sure what to do here

% %%%Area under the curve and width values from threshold value - require zoom%%%
Area = trapz(t_select,vm_select-min(vm_select))   % units ? mV*s ; NB : to compute area, y must be offset so that y min value is zero


%%Display width(ms)
% this is the x distance between a_sol and b_sol
Xwidth = 1000*(-a_sol+b_sol)  % units :  ms

% width = (thresholdt2 - thresholdt1)*1000  %%%%needs to be defined by threshold
% %%Display half-width(ms)
% APhalfwidth = (width/2)
APhalfwidth = (Xwidth/2)


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function err = obj_fun(x)
global xdata ydata x_vm_max vm_max
    XI = [min(xdata),x(1),x_vm_max,x(2),max(xdata)];  % vector of 1-D look-up table "x" points
    YI = lsq_lut_piecewise( xdata, ydata, XI ); % obtain vector of 1-D look-up table "y" points
    Yint = interp1(XI,YI,xdata);
    err = sum((Yint-ydata).^2);
end

[Roksana Khalid]

This works well on the xyz.mat action potential, thanks!

The threshold and width make sense from generated figures, but I think the threshold should be -61.14mV not -62mV, based on the phase plot associated with the same curve.

Also, I tried running it on an additional dataset but the LUT points are messed up. :/

[Roksana Khalid]

It gave these values for width and half-width:

Xwidth = 59.9300

APhalfwidth = 29.9650

[Mathieu NOE]

hello

I tested again my code with the LUT with your new data and got "normal" behaviour not like your

??

anyway, I ended up doing very simple math , no fancy functions needed at all anymore

I found that the best results IMO are obtained simply by using a rule on a fixed threshold on the first derivative of the data ( no smoothing anymore , neither LUT , etc ...) , maybe this fixed value (around 4e3 mV/s) has some physical meaning for you;

I tested with the two supplied datasets and I'm pretty happy with the results !

you only need the firstsecondderivatives sub function

code below :

%Import data
% load('xyz.mat')
load('AP2.mat')

     
vm = S2C1_17082020_WM_apomorphine_20uM_Ch1.values;
samples = length(vm);
dt = S2C1_17082020_WM_apomorphine_20uM_Ch1.interval;
t_final = samples*dt;

%scale x-axis by interval
t = linspace(0, t_final,samples);

% compute first derivative on data
[dvm, ddvm] = firstsecondderivatives(t,vm);

figure(1);
subplot(2,1,1),plot(t,vm,'b');
ylabel('V (mV)');
subplot(2,1,2),plot(t,dvm,'b');
ylabel('dV/dt (mV/s)');
xlabel('Time (s)');
title('Membrane potential')

% 1st knee point has slope approx 1/N th of max  slope  (positive threshold, positive slope)
% 2nd knee point has slope approx 1/N th of max  slope  (negative threshold, positive slope)
% a_pos1 = 150;
% pos_slope_max = max(dvm(dvm>=0));
% pos_slope_trigger = pos_slope_max/a_pos1;   % 4e3
pos_slope_trigger =  4e3;

% a_pos2 = 25;
% neg_slope_max = min(dvm(dvm<=0));
% neg_slope_trigger = neg_slope_max/a_pos2;   % -4e3
neg_slope_trigger =  -4e3;

indpos = find(dvm>=pos_slope_trigger)
t0_pos1 = t(indpos(1));

indneg = find(dvm<=neg_slope_trigger)
t0_pos2 = t(indneg(end));

vm_knee1 = interp1(t,vm,t0_pos1,'linear');
vm_knee2 = interp1(t,vm,t0_pos2,'linear');

figure(2);
plot(t,vm,'b',t0_pos1,vm_knee1,'+r',t0_pos2,vm_knee2,'+g')
ylabel('V (mV)')
xlabel('Time (s)')
title('Membrane potential')