Different outputs of ode solvers inside for loop

1 次查看(过去 30 天)
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 个评论
Mounir Arab
Mounir Arab 2022-8-18
True LastPeakinterval is the interval between the timing of last peak and 340 ms
Mounir Arab
Mounir Arab 2022-8-18
I think I am obtaining different solutions (voltage_model) for the same values of k if i changed the ranges of k

请先登录,再进行评论。

回答(1 个)

Torsten
Torsten 2022-8-18
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
  4 个评论
Mounir Arab
Mounir Arab 2022-8-18
I tried it a million times, I am still getting different values.
Torsten
Torsten 2022-8-18
编辑:Torsten 2022-8-18
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.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Matrix Computations 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by