Clear Filters
Clear Filters

How to find the time between two spikes?

19 views (last 30 days)
Please help me to to find the time between two spikes i.e. inter-spike interval T in matlab. Below is the code for spike generation.
close all
clear
clc
global Vr
tic
%flagJ = 1;
% Amplitude of step function [100e-6 A.cm-2]
J0 = 10e-6;
% Simulation time [40e-3 s]
tMax = 40e-3;
% Time stimulus applied [5e-3 s]
tStart = 5e-3;
% Number of grid points [8001]
num = 8001;
Jext = zeros(num,1); % external current density [A.cm^-2]
t = linspace(0,tMax,num);
dt = t(2) - t(1);
num1 = find(t > tStart, 1 ); % index for stimulus ON
Jext(num1:num) = J0 % external stimulus current
% FIXED PARAMETERS ====================================================
sf = 1e3; % conversion of V to mV
VR = -65e-3; % resting voltage [V]
Vr = VR*1e3; % resting voltage [mV]
VNa = 50e-3; % reversal voltage for Na+ [V]
VK = -77e-3; % reversal voltage for K+ [V]
Cm = 1e-6; % membrane capacitance/area [F.cm^-2]
gKmax = 36e-3; % K+ conductance [S.cm^-2]
gNamax = 120e-3; % Na+ conductance [S.cm.-2)]
gLmax = 0.3e-3; % max leakage conductance [S.cm-2]
T = 20; % temperature [20 deg C]
% SETUP ===============================================================
JNa = zeros(num,1); % Na+ current density (A.cm^-2)
JK = zeros(num,1); % K+ current density (A.cm^-2)
JL = zeros(num,1); % leakage current density (A.cm^-2)
Jm = zeros(num,1); % membrane current (A.cm^-2)
V = zeros(num,1); % membrane potential (V)
gNa = zeros(num,1); % Na+ conductance
gK = zeros(num,1); % K+ conductance
gL = ones(num,1); % gL conductance
n = zeros(num,1); % K+ gate parameter
m = zeros(num,1); % Na+ gate parameter
h = zeros(num,1); % Na+ gate parameter
% Initial Values -----------------------------------------------------
V(1) = VR; % initial value for membrane potential
[An, Bn, Am, Bm, Ah, Bh] = AlphaBeta(V(1), T);
n(1) = An / (An + Bn);
m(1) = Am / (Am + Bm);
h(1) = Ah / (Ah + Bh);
gK(1) = gKmax * n(1)^4;
gNa(1) = gNamax * m(1)^3 * h(1);
gL = gLmax .* gL;
JK(1) = gK(1) * (V(1) - VK);
JNa(1) = gNa(1) * (V(1) - VNa);
Vadj = (Jext(1) - JK(1) - JNa(1))/gL(1);
JL(1) = gL(1) * (V(1) - VR + Vadj);
Jm(1) = JNa(1) + JK(1) + JL(1);
V(1) = VR + (dt/Cm) * (-JK(1) - JNa(1) - JL(1) + Jext(1));
% Solve ODE -----------------------------------------------------------
for cc = 1 : num-1
[An, Bn, Am, Bm, Ah, Bh] = AlphaBeta(V(cc), T);
An = sf * An; Am = sf * Am; Ah = sf * Ah;
Bn = sf * Bn; Bm = sf * Bm; Bh = sf * Bh;
n(cc+1) = n(cc) + dt * (An *(1-n(cc)) - Bn * n(cc));
m(cc+1) = m(cc) + dt * (Am *(1-m(cc)) - Bm * m(cc));
h(cc+1) = h(cc) + dt * (Ah *(1-h(cc)) - Bh * h(cc));
gK(cc+1) = n(cc+1)^4 * gKmax;
gNa(cc+1) = m(cc+1)^3 * h(cc+1) * gNamax;
JK(cc+1) = gK(cc+1) * (V(cc) - VK);
JNa(cc+1) = gNa(cc+1) * (V(cc) - VNa);
JL(cc+1) = gL(cc+1) * (V(cc) - VR + Vadj);
Jm(cc+1) = JNa(cc+1) + JK(cc+1) + JL(cc+1);
V(cc+1) = V(cc) + (dt/Cm) * (-JK(cc+1) - JNa(cc+1) - JL(cc+1) + Jext(cc+1));
end
Vdot = (Jext - Jm)./Cm;
% GRAPHICS ============================================================
% Width & height of Figure Window
w = 0.3; d = 0.30;
% Position of Figure window (x,y)
x1 = 0.02; y1 = 0.45;
x2 = 0.35; y2 = 0.05;
x3 = 0.67;
% Membrane Potential
figure(2)
set(gcf,'units','normalized');
set(gcf,'position',[x2 y1 w d]);
set(gcf,'color','w');
LW = 2;
x = t.*sf; y = V.*sf;
plot(x,y,'b','linewidth',LW); % membrane voltage
xlabel('time t [ ms ]'); ylabel('membrane voltage V_m [ mV ]');
set(gca,'fontsize',12)
grid on
box on
hold on
%if flagJ == 1
%tm = sprintf('J_0 = %2.0f \\muA.cm^{-2} ',J0.*1e6);
%title(tm)
%end
toc
% FUNCTIONS =========================================================
function [An, Bn, Am, Bm, Ah, Bh] = AlphaBeta(V, T)
global Vr
V = V*1000;
dV = (V - Vr);
phi = 3^((T-6.3)/10);
An = phi * (eps + 0.10 - 0.01 .* dV) ./ (eps + exp(1 - 0.1 .* dV) - 1);
Am = phi * (eps + 2.5 - 0.1 .* dV) ./ (eps + exp(2.5 - 0.1 .* dV) - 1);
Ah = phi * 0.07 .* exp(-dV ./ 20);
Bn = phi * 0.125 .* exp(-dV ./ 80);
Bm = phi * 4 .* exp(-dV/18);
Bh = phi * 1 ./ (exp(3.0 - 0.1 .* dV) + 1);
end

Accepted Answer

Karim
Karim on 23 Jun 2022
hello, you can use the "findpeaks" routine do achieve this, see below for the code
best regards
% first run the main code
MainCode();
Elapsed time is 0.394297 seconds.
% in the main script, the voltage is saved in variable "y"
% use findpeaks to find the voltage peaks and the indexes
[Voltage_pks,locs] = findpeaks(y);
% in the main script, the time vector is saved in the variable "x"
% used the above indices to find the coressponding timestamps
Time_pks = x(locs)
Time_pks = 1×7
6.5900 11.4350 16.2700 21.1050 25.9400 30.7700 35.6050
% now find the difference between the peaks
Time_diff = diff(Time_pks)
Time_diff = 1×6
4.8450 4.8350 4.8350 4.8350 4.8300 4.8350

More Answers (0)

Categories

Find more on Mathematics in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!