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

14 visualizaciones (últimos 30 días)
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)
  7 comentarios

Iniciar sesión para comentar.

Respuesta aceptada

Star Strider
Star Strider el 2 de Feb. de 2021
Editada: Star Strider el 2 de Feb. de 2021
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.
.
  4 comentarios

Iniciar sesión para comentar.

Más respuestas (1)

Mathieu NOE
Mathieu NOE el 2 de Feb. de 2021
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
  3 comentarios
Roksana Khalid
Roksana Khalid el 2 de Feb. de 2021
It gave these values for width and half-width:
Xwidth = 59.9300
APhalfwidth = 29.9650
Mathieu NOE
Mathieu NOE el 3 de Feb. de 2021
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')

Iniciar sesión para comentar.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by