How to find the time between two spikes?

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
  

Answer 1

[Karim]

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