>> clear all;
GmaxNa=120; %mS/cm^2,max conductance
EkNa=50;%mV, Nerst potential
GmaxK=36; %mS/cm^2,max conductance
EkK=-80;%mV, Nerst potential
GvCl=0.3;%mS/cm^2,constant conductance
EkCl=-56;%mV, Nerst potential
C=1;%uF/cm^2, capacitance of membrane
%----------initial-------------%
V=-67.3;% initial membrane potential
%m = 0.08;%@-48mV
%n = 0.32;%@-48mV
%h = 0.59;%@-48mV
%m = 0.065;%@-57.6mV
%n = 0.32;%@-57.6mV
%h = 0.59;%@-57.6mV
Er=-67.3;%mV, rest potential
Vr=V-Er;% potential relative to Er
alpha_n = 0.01*(-Vr+10)/(exp((-Vr+10)/10)-1);
beta_n = 0.125*exp(Vr/80);
alpha_m = 0.1*(-Vr+25)/(exp((-Vr+25)/10)-1);
beta_m = 4*exp(-Vr/18);
alpha_h = 0.07*exp(-Vr/20);
beta_h = 1/(exp((-Vr+30)/10)+1);
m=alpha_m/(alpha_m+beta_m);
n=alpha_n/(alpha_n+beta_n);
h=alpha_h/(alpha_h+beta_h);
GvNa = GmaxNa*m^3*h;
GvK = GmaxK*n^4;
I_Na=GvNa*(V-EkNa);%uA
I_K=GvK*(V-EkK);%uA
I_Cl=GvCl*(V-EkCl);%uA
%Er=(GvNa*EkNa+GvK*EkK+GvCl*EkCl)/(GvNa+GvK+GvCl);
dt=0.001;%msec, time step
%----------initial end--------------%
num_idx=0;
%-----------------------------------------begin for-----------------------------------------%
for time=0:dt:4
Vr=V-Er;
alpha_n = 0.01*(-Vr+10)/(exp((-Vr+10)/10)-1);
beta_n = 0.125*exp(-Vr/80);
alpha_m = 0.1*(-Vr+25)/(exp((-Vr+25)/10)-1);
beta_m = 4*exp(-Vr/18);
alpha_h = 0.07*exp(-Vr/20);
beta_h = 1/(exp((-Vr+30)/10)+1);
GvNa = GmaxNa*m^3*h;
GvK = GmaxK*n^4;
I_Na=GvNa*(V-EkNa);%uA
I_K=GvK*(V-EkK);%uA
I_Cl=GvCl*(V-EkCl);
%-----------------------------------changes after dt---------------------------------%
dm = (alpha_m*(1-m)-beta_m*m)*dt;
dn = (alpha_n*(1-n)-beta_n*n)*dt;
dh = (alpha_h*(1-h)-beta_h*h)*dt;
m = m+dm;
n = n+dn;
h = h+dh;
%-------------------injected current---------------------%
if time<0.1
Icm=100;%uA/cm^2
else
Icm=0;
end
%-------------------injected current---------------------%
dV=(Icm-(I_Na+I_K+I_Cl))*dt/C;%mV
V=V+dV;
%-----------------------------------changes after dt---------------------------------%
num_idx=num_idx+1;
disp_V(num_idx)=V;
disp_time(num_idx)=time;
end
%-----------------------------------------end for--------------------------------------------%
plot(disp_time,disp_V);