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')

[Feng Du]

% 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;

Hi Mathieu,

Thanks for your last week's answer, I understand the math works.

But for the AP analyses, I still feel confused with your codes, can you explain how the part of the line came to be? And why is the value of a_pos1 defined as 150, and a_pos2 = 25？

Have a nice day!

Best regaeds,

Feng

[Mathieu NOE]

hello again @Feng Du

it took me a few minutes to remember what I did here ... it coming back to surface :)

% 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)

ok so the simple approach here was to find the knee points by first measuring the max positive and max negative slope of the signal then define by which factors (one positive and one negative) I have to divide these slope maxima values to find the right slope values at the knee points

Hope this clarifies these lines (it the same stuff explained a bit differently)

I admit that for sake of clarity , I forgot to mention that the "N" factor I was talking in the comments (above) will be named differently after in y code , which could cause some trouble to follow through :

(N=> a_pos1) a_pos1 = 150; % for the 1st knee point slope computation

and

(N=> a_pos2) a_pos2 = 25 ; % for the 2nd knee point slope computation

so for the first knee , we have a max positive slope of

pos_slope_max = 5.7554e+05 and to get the right threshold value I had to divide this by 150 (hence a_pos1 = 150) , as coded in : pos_slope_trigger = pos_slope_max/a_pos1;

same stuff for the second knee :

neg_slope_max = -1.0698e+05 , divided by a_pos2 = 25 , as in the code : neg_slope_trigger = neg_slope_max/a_pos2;

after this I noticed that pos_slope_trigger and neg_slope_trigger would have very similar values that I rounded to +4e3 and -4e3 respectively.

at the end the code can be simplified with those hard coded values and the unnecessary lines are commented just to remember how we get there.

% 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;

Of course there are other ways I could have figured out to find those treshold values

[Feng Du]

Hi Mathieu,

Thank you for your patient explanation, but I am still confused about the origin of 150 and 25. Please forgive my insufficient mathematical knowledge.

Can you give an example with my data?

Fs = 20000, The stimulation time at 1s. The recording time is 2s.

It is not an action potential, but the response of the GPe neuron after whisker stimulation is recorded.

I tried to find its onset point.

Best regards,

Feng

[Mathieu NOE]

hello again

I plotted your data , which points do you want to identify ?

[Feng Du]

Hi Mathieu,

I want to determine the point indicated by the second arrow. The first arrow is at 1 second, which represents the time of the stimulus. I want to know how long after the stimulus the response appears, so I call it the onset point.

Best,

Feng

[Mathieu NOE]

hello Feng

so what I did here is simply add a bit of gaussian smoothing on the derivative of your data

then I saw that the slope at the requested point location was around 200 units / s and then you got your point .

%Import data
load('My_signal.mat')
     
vm = M;
samples = length(vm);
Fs = 20000;
dt =1/Fs;
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);

% add some smoothing on the derivative
dvm = smoothdata(dvm,'gaussian',round(samples/100));

figure(1);
subplot(2,1,1),plot(t,vm,'b');
ylabel('y');
subplot(2,1,2),plot(t,dvm,'b');
ylabel('dy/dt');
xlabel('Time (s)');
% 1st knee point has slope approx 200
pos_slope_trigger =  200;
indpos = find(dvm>=pos_slope_trigger)
t0_pos1 = t(indpos(1));
vm_knee1 = interp1(t,vm,t0_pos1,'linear');
figure(2);
plot(t,vm,'b',t0_pos1,vm_knee1,'.r','markersize',25)
ylabel('y');
xlabel('Time (s)')

[Mathieu NOE]

if you prefer , you can also define the threshold as being a small fraction of the max value of the derivative of signal

so only one line to change :

% pos_slope_trigger =  200;
pos_slope_trigger =  max(dvm)/5;

[Star Strider]

Perhaps this —

imshow(imread('figure.png'))

LD = load('My_signal');
s = LD.M;
Fs = 2E+4;
L = numel(s);
t = linspace(0, L-1, L)/Fs;

Lvlm = islocalmax(s, 'MinProminence',5);            % Peaks
lmt = t(Lvlm);

Lvic = ischange(s, 'Threshold',10);                 % Changes
ict = t(Lvic);

tr = ict((ict > 1) & (ict < lmt(1)));               % Time Threshold

figure
plot(t,s)
hold on
plot(tr(1), s(t == tr(1)), '|r')
hold off
grid
text(tr(1), s(t ==tr(1)), sprintf('\\leftarrow t = %.4f, \\Delta(t) = %.4f',tr(1), 1-tr(1)))
xl1 = xline(1, '--r', 'Stimulus');
xl2 = xline(tr(1), '--r', 'Response');
xl1.LabelHorizontalAlignment = 'center';
xl2.LabelHorizontalAlignment = 'center';

.

[Feng Du]

Thank you guys, you helped me a lot!!

Best regards,

Feng

[Star Strider]

My pleasure!

Answer 3

[vignesh]

data = scale(data)

data.Time = (data.Time - data.Time(1))/1000;

data.X = 1.5*data.X;

data.X = data.X - mean(data.X);

data.Y = data.Y - mean(data.Y);

end