Different outputs of ode solvers inside for loop

I attached my matlab code to this question.

In brief, I am trying to generate hodgkin and huxley neuron model by solving odes, the problem is that I am getting different values for LastSpikeInterval(inside the matrix variable) with FOR loop (the name of the file is: Plot_traces_new_new.m) at different ranges (for example: k=0:0.2:18 and k=14:0.2:18) which makes the outputs of the model inaccurate. Does someone know how to fix this issue?

[Torsten]

Could you explain in plain words what LastSpikeInterval is ? The length of the time interval between t=340 and the last peak in the array "voltage_model" that exceeds -20 ?

[Mounir Arab]

True LastPeakinterval is the interval between the timing of last peak and 340 ms

[Mounir Arab]

I think I am obtaining different solutions (voltage_model) for the same values of k if i changed the ranges of k

Answer 1

[Torsten]

This code gives reproducible results; I don't know exactly why your code didn't work.

matrix=zeros(90,2);

%Model Trace;

t2_stim=340;

time=0:0.02:500;

K = 14:0.2:18;
%K = 0:0.2:18;

for i=1:numel(K)
    k = K(i);
    
    
    [voltage_model]= Plot_Full(time,k);
    %plot(time,voltage_model)
    
    
    
    
    [pks locs]=findpeaks(voltage_model,'MinPeakHeight',-20);
    SpikeTiming=locs*0.02;
    ISI=diff(SpikeTiming);
    
    if isempty(pks)==0
        LastSpikeInterval=t2_stim-SpikeTiming(end);
    else
        LastSpikeInterval=0;
        
    end
    matrix(i,1)=k;
    matrix(i,2)=LastSpikeInterval;
    
end
matrix(1:end,:)

ans = 90×2

14.0000 1.7000 14.2000 6.2400 14.4000 1.0000 14.6000 17.6000 14.8000 86.4800 15.0000 4.8000 15.2000 4.0600 15.4000 72.2600 15.6000 24.9600 15.8000 77.0000

%matrix(71:end,:)

function [v_i_model] = Plot_Full(time,k)

%function [v_i_model] = Plot_Full(time,v_i)

thetan=-13;%
sigman=-14;%

thetae=-60;
sigmae=5;

thetaw=-43;
sigmaw=8;
thetaz=-67;
sigmaz=7.3;

thetay=-45;
sigmay=5;

% Initial Conditons

V1=-78;
v_i=V1;

alphaH_i = 0.128*exp(-(v_i+50)/18);
betaH_i= 4/(1+exp(-(v_i+27)/5));
h_i = alphaH_i/(alphaH_i+betaH_i);

n_i = 1./(1+exp((v_i-thetan)/sigman));
Caid_i=0.103;
Cai_i=0.103;
ed_i = 1./(1+exp((v_i-thetae)/sigmae));

wd_i =1-1/(1+exp((v_i-thetaw)/sigmaw));
zd_i = 1/(1+exp((v_i-thetaz)/sigmaz));

w_i =1-1/(1+exp((v_i-thetaw)/sigmaw));
z_i = 1/(1+exp((v_i-thetaz)/sigmaz));

yd_i = 1/(1+exp((v_i-thetay)/sigmay));
y_i=yd_i;
vd_i=v_i;

init=[v_i;h_i;n_i;y_i;vd_i;Caid_i;wd_i;zd_i;ed_i;yd_i;Cai_i;w_i;z_i;];


%[~,output]=ode113('DiffEquations_Full',time,init); %113
[~, output] = ode113(@(time, init) DiffEquations_Full( init,time, k), time, init);


v_i_model = output(:,1);
end


function [ output ] = DiffEquations_Full( init,time, k)

%function [output] = DiffEquations_Full(time,init)


t1_stim=40;
t2_stim=340;

iapp= 305;

Cm=15;
Cd=15;
gc=50;

gNa=350;
VNa=50;
thetam=-32; %
sigmam=-6.8; %
tauh=1;

gK=1800;
VK=-90;
taunbar=10;
thetan=-13;%
sigman=-14;%

gCaL=5;
gCaLd=1;
Ca_ex=2.5;
RTF=26.7;
thetas=-20;
sigmas=-0.05;

gSK=14; %14 %14.6 10.8
gSKd=k;
ks=0.5;
f=0.1;
eps=0.0015;
kCa=0.3;

gAd=0;
thetaa=-20;
sigmaa=-10;
thetae=-60;
sigmae=5;
taue=20;

gD=0;
gDd=20;
thetaw=-43;
sigmaw=8;
thetaz=-67;
sigmaz=7.3;
tauw=1;
tauz=3500;

gM=0;
gMd=0;
thetay=-45;
sigmay=5;
tauy=30;

gL=0.1;
gLd=0.1;
VL=-70;


V=init(1);
h=init(2);
n=init(3);

yy=init(4);
Vd=init(5);
Caid=init(6);
wd=init(7);
zd=init(8);
ed=init(9);
yd=init(10);
Cai=init(11);
w=init(12);
z=init(13);

%% Soma

%Sodium Current (Na+)
minf = 1./(1+exp((V-thetam)/sigmam));
alphah = 0.128*exp(-(V+50)/18);
betah = 4/(1+exp(-(V+27)/5));
hinf = alphah/(alphah+betah);
iNa = gNa*(minf^3)*h*(V-VNa);


%Potassium Current (K+)
ninf = 1./(1+exp((V-thetan)/sigman));
taun = taunbar./cosh((V-thetan)/(2*sigman));
iK = gK*(n^4)*(V-VK);


% High-threshold L-type Ca 2+ Current (Ca-L)
sinf = 1./(1+exp((V-thetas)/sigmas));

iCaL = gCaL*(sinf^2)*V*(Ca_ex./(1-exp((2*V)./RTF)));


% Ca 2+ -dependent K+ current (SK)
kinf = (Cai^2)/(Cai^2+(ks^2));

iSK= gSK*kinf*(V-VK);


% D-type K+ Current (D)
winf=1-1/(1+exp((V-thetaw)./sigmaw));
zinf=1/(1+exp((V-thetaz)/sigmaz));
iD=gD*w*z*(V-VK);

%M-type K+ current (M)
yinf = 1./(1+exp(-(V-thetay)./sigmay));
iM=gM*yy*(V-VK);


%Leak Current
iL=gL*(V-VL);


%% Dendrites

% High-threshold L-type Ca 2+ Current (Ca-L)

sinfd = 1./(1+exp((Vd-thetas)/sigmas));
iCaLd = gCaLd*(sinfd^2)*Vd*(Ca_ex./(1-exp((2*Vd)./RTF)));


% Ca 2+ -dependent K+ current (SK)
kinfd = (Caid^2)/(Caid^2+(ks^2));

iSKd= gSKd*kinfd*(Vd-VK);

% A-type K+ Current (A)
ainfd = 1./(1+exp((Vd-thetaa)/sigmaa));
einfd = 1./(1+exp((Vd-thetae)/sigmae));

iAd=gAd*ainfd*ed*(Vd-VK);

% D-type K+ Current (D)
winfd=1-1/(1+exp((Vd-thetaw)/sigmaw));
zinfd=1/(1+exp((Vd-thetaz)/sigmaz));
iDd=gDd*wd*zd*(Vd-VK);

%M-type K+ current (M)
yinfd = 1./(1+exp(-(Vd-thetay)./sigmay));
iMd=gMd*yd*(Vd-VK);


%Leak Current
iLd=gLd*(Vd-VL);
% Applied Current


if time>t1_stim && time<t2_stim
    istim=iapp;
else
    istim=0;    
end

output(1,1)=(-iNa-iK-iCaL-iSK-iM-iD-iL+gc*(Vd-V)+istim)/Cm;
output(2,1)=(hinf-h)/tauh;
output(3,1)=(ninf-n)/taun;

output(4,1)=(yinf-yy)/tauy;
output(5,1)=(-iCaLd-iSKd-iMd-iDd-iAd-iLd+gc*(V-Vd))/Cd;
output(6,1)=-f*(eps*(iCaLd)+ kCa*(Caid-0.1));
output(7,1)=(winfd-wd)/tauw;
output(8,1)=(zinfd-zd)/tauz;
output(9,1)=(einfd-ed)/taue;
output(10,1)=(yinfd-yd)/tauy;
output(11,1)=-f*(eps*(iCaL)+ kCa*(Cai-0.1));
output(12,1)=(winf-w)/tauw;
output(13,1)=(zinf-z)/tauz;

end

[Torsten]

@Mounir Arab comment moved here:

Would you get the same value of LastSpikeInterval fork between 14 and 18 if you changed k to k=0:0.2:18?

[Torsten]

Why don't you test it ? I did ...

[Mounir Arab]

I tried it a million times, I am still getting different values.

[Torsten]

Running the code above and changing

K = 14:0.2:18;
%K = 0:0.2:18;

to

%K = 14:0.2:18;
K = 0:0.2:18;

and

matrix(1:end,:)
%matrix(71:end,:)

to

%matrix(1:end,:)
matrix(71:end,:)

gives at least the same values for k from 14 to 15.8 for both cases. Here, you already encountered discrepancies in your two Excel sheets. I cannot see the following values.