Different outputs of ode solvers inside for loop
1 view (last 30 days)
Show older comments
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?
3 Comments
Answers (1)
Torsten
on 18 Aug 2022
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,:)
%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
4 Comments
Torsten
on 18 Aug 2022
Edited: Torsten
on 18 Aug 2022
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.
See Also
Categories
Find more on Matrix Computations 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!