This code calculate m n h values in For loop, and m n h values are probability values, they must be between 0 and 1, So m n h values must be between 0 and 1 when the for loop is over. But I cant When I run the code, I can not get these values.Can you

% === simulation time=== 
simulationtime= 1000; % Milisaniye cinsinden 
deltaT=0.01;
t=0:deltaT:simulationtime;

% === specify external current I === 
changeTimes = [0]; %milliseconds% 
CurrentLevel = [50]; %Change this to see the effect of the different currents on the voltage. (Suggested values: 3, 20, 50, 1000)

%Set the externally applied current over time.
%Here, the first 500 time-outs are now at 50, the next 1500 hours.
%zero current (zero neuron rest potential).
%and the remaining timeouts are constant current.
I(1:500) = CurrentLevel; 
I(501:2000) = 0; 
I(2001:numel(t)) = CurrentLevel;
%Interpret the line above and subtract it for the following line constant current.
%and the effect of the voltage time

%===fixed parameters===%
g_K=36;    
g_Na=120; 
g_L=0.3;
E_K =-12;   %-12
E_Na=115;   %115
E_L=10.6;   %10,6
C=1;

%===initial situations===%
V=0; % Başlangıç hattı gerilimiz
alfa_n = 0.01 * ( (V+55) / (1-exp(-(V+55)/10)) ); 
beta_n = 0.125*exp(-(V+65)/80);

alfa_m = 0.1*( (V+40) / (1-exp(-(V+40)/10)) ); 
beta_m = 4*exp(-(V+65)/20);

alfa_h = 0.07*exp(-(V+65)/20); 
beta_h = 1/(exp(-(V+35)/10)+1); 

n(1) = alfa_n/(alfa_n+beta_n); 
m(1) = alfa_m/(alfa_m+beta_m); 
h(1) = alfa_h/(alfa_h+beta_h); 

for i=1:numel(t)-1 %Coefficients, currents and derivative calculations in each step.

      %---calculate coefficients---%
      %The equations in this section are the same as above, only calculated in each step.
      alfa_n(i) = 0.01 * ( (10-V(i)) / (exp((10-V(i))/10)-1) );
      beta_n(i) = 0.125*exp(-V(i)/80);
      alfa_m(i) = 0.1*( (25-V(i)) / (exp((25-V(i))/10)-1) );
      beta_m(i) = 4*exp(-V(i)/18);
      alfa_h(i) = 0.07*exp(-V(i)/20);
      beta_h(i) = 1/(exp((30-V(i))/10)+1);

      %---calculate currents---%
      I_Na = (m(i)^3) * g_Na * h(i) * (V(i)-E_Na); 
      I_K = (n(i)^4) * g_K * (V(i)-E_K); 
      I_L = g_L *(V(i)-E_L); 
      I_iyon = 0 - I_K - I_Na - I_L; 

      %---Euler calculates derivative using first order approximation.---%
      V(i+1) = V(i) + deltaT*I_iyon/C;
      gm=sqrt(((2*alfa_m(i)*beta_m(i))/(60*128*(alfa_m(i)+beta_m(i)))))*randn;
      gn=sqrt(((2*alfa_n(i)*beta_n(i))/(18*128*(alfa_n(i)+beta_n(i)))))*randn;
      gh=sqrt(((2*alfa_h(i)*beta_h(i))/(60*128*(alfa_h(i)+beta_h(i)))))*randn;
      n(i+1) = n(i) + deltaT*(alfa_n(i) *(1-n(i)) - beta_n(i) * n(i))+ gm; 
      m(i+1) = m(i) + deltaT*(alfa_m(i) *(1-m(i)) - beta_m(i) * m(i))+ gn; 
      h(i+1) = h(i) + deltaT*(alfa_h(i) *(1-h(i)) - beta_h(i) * h(i))+ gh; 

end
E_Na_mu = 115;
E_Na_std = 5;
E_Na = normrnd(E_Na_mu, E_Na_std);

V = V-70; %Set resting potential to -70mv
figure
%===plot Voltage===%
plot(t,V,'LineWidth',3)
%hold on
legend({'Voltaj'})
ylabel('Voltaj (mv)')
xlabel('Zaman (ms)')
title('')

%===plot Conductance===%
%figure
%p1 = plot(t,g_K*n.^4,'LineWidth',2);
%hold on
%p2 = plot(t,g_Na*(m.^3).*h,'r','LineWidth',2);
%legend([p1, p2], 'Conductivity for Potassium', 'Conductivity for Sodium')
%ylabel('Conductivity')
%xlabel('Time (ms)')
%title('Conductivity for Potassium and Sodium Ions')

Answer 1

[Torsten]

Use

disp(n(i+1))
disp(m(i+1))
disp(h(i+1))

at the end of the for-loop.

Best wishes

Torsten.

[ayten kurnaz]

I did what you said.This time I got this error 'Index exceeds matrix dimensions.' Error in This Line: disp(n(i+1)) = n(i) + deltaT*(alfa_n(i) *(1-n(i)) - beta_n(i) * n(i))+ gm;

[Torsten]

What if you initialize

V(1)=0

instead of

V=0

?

Best wishes

Torsten.

[ayten kurnaz]

No,ıt doesnt work.By the way thank you for your answers.

[Torsten]

Try

V(1)=0; % Baslangiç hatti gerilimiz
alfa_n(1) = 0.01 * ( (V(1)+55) / (1-exp(-(V(1)+55)/10)) ); 
beta_n(1) = 0.125*exp(-(V(1)+65)/80);

alfa_m(1) = 0.1*( (V(1)+40) / (1-exp(-(V(1)+40)/10)) ); 
beta_m(1) = 4*exp(-(V(1)+65)/20);

alfa_h(1) = 0.07*exp(-(V(1)+65)/20); 
beta_h(1) = 1/(exp(-(V(1)+35)/10)+1); 

n(1) = alfa_n(1)/(alfa_n(1)+beta_n(1)); 
m(1) = alfa_m(1)/(alfa_m(1)+beta_m(1)); 
h(1) = alfa_h(1)/(alfa_h(1)+beta_h(1));

instead of

V=0; % Başlangıç hattı gerilimiz
alfa_n = 0.01 * ( (V+55) / (1-exp(-(V+55)/10)) ); 
beta_n = 0.125*exp(-(V+65)/80);

alfa_m = 0.1*( (V+40) / (1-exp(-(V+40)/10)) ); 
beta_m = 4*exp(-(V+65)/20);

alfa_h = 0.07*exp(-(V+65)/20); 
beta_h = 1/(exp(-(V+35)/10)+1); 

n(1) = alfa_n/(alfa_n+beta_n); 
m(1) = alfa_m/(alfa_m+beta_m); 
h(1) = alfa_h/(alfa_h+beta_h);

[ayten kurnaz]

ıt gives the same error.

[Torsten]

You included the "disp" commands as

for i=1:numel(t)-1 
      %Coefficients, currents and derivative calculations in each step.

      %---calculate coefficients---%
      %The equations in this section are the same as above, only calculated in each step.
      alfa_n(i) = 0.01 * ( (10-V(i)) / (exp((10-V(i))/10)-1) );
      beta_n(i) = 0.125*exp(-V(i)/80);
      alfa_m(i) = 0.1*( (25-V(i)) / (exp((25-V(i))/10)-1) );
      beta_m(i) = 4*exp(-V(i)/18);
      alfa_h(i) = 0.07*exp(-V(i)/20);
      beta_h(i) = 1/(exp((30-V(i))/10)+1);
      %---calculate currents---%
      I_Na = (m(i)^3) * g_Na * h(i) * (V(i)-E_Na); 
      I_K = (n(i)^4) * g_K * (V(i)-E_K); 
      I_L = g_L *(V(i)-E_L); 
      I_iyon = 0 - I_K - I_Na - I_L; 
      %---Euler calculates derivative using first order approximation.---%
      V(i+1) = V(i) + deltaT*I_iyon/C;
      gm=sqrt(((2*alfa_m(i)*beta_m(i))/(60*128*(alfa_m(i)+beta_m(i)))))*randn;
      gn=sqrt(((2*alfa_n(i)*beta_n(i))/(18*128*(alfa_n(i)+beta_n(i)))))*randn;
      gh=sqrt(((2*alfa_h(i)*beta_h(i))/(60*128*(alfa_h(i)+beta_h(i)))))*randn;
      n(i+1) = n(i) + deltaT*(alfa_n(i) *(1-n(i)) - beta_n(i) * n(i))+ gm; 
      m(i+1) = m(i) + deltaT*(alfa_m(i) *(1-m(i)) - beta_m(i) * m(i))+ gn; 
      h(i+1) = h(i) + deltaT*(alfa_h(i) *(1-h(i)) - beta_h(i) * h(i))+ gh; 
      disp(n(i+1))
      disp(m(i+1))
      disp(h(i+1))

end

didn't you ?

Best wishes

Torsten.

[ayten kurnaz]

It worked but m,n and h are not between 0 and 1.

[Torsten]

"randn" can be quite large and negative ...

You won't be able to guarantee that n,m and h remain between 0 and 1 if you use a random variable in the updating.

Best wishes

Torsten.

[ayten kurnaz]

Thank you. Best wishes Ayten