# How to find the time between two spikes?

14 views (last 30 days)
Haya Ali on 23 Jun 2022
Answered: Karim on 23 Jun 2022
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

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