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?
回答(1 个)
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,:)
%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 个评论
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 Center 和 File Exchange 中查找有关 Matrix Computations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!