% Traub.Miles active soma
% modified FengXiao code

function traub1
close all

% reversal potentials

vL = -15;
vna = 115;
vna = 100;
vk = -15;
vk = -30;
vca = 140;

% parameters

gL = .01;
gna = 3.32;
gca = 6.64;
gk = 3.98;
gkca = 0.1;  % uS

Cm = 3*3320*1e-5;	%nF

c = 5200; %mmol/megacoulomb
A = 3320; %um^2
d = 5e-4; %um
D = c/(A*d);	% mM/ms/muA

xib = 0.1;  % ms^-1

j = 1; dt = 0.025;  T=60;
t = 0:dt:T;

% initial values 

v = zeros(size(t));
xi = v;
ina = v;
ik = v;
ica = v;
ikca = v;

v(1) = vL;

m = am(v(1))/(am(v(1))+bm(v(1)));
h = ah(v(1))/(ah(v(1))+bh(v(1)));
n = an(v(1))/(an(v(1))+bn(v(1)));
y = ay(v(1))/(ay(v(1))+by(v(1)));
s = as(v(1))/(as(v(1))+bs(v(1)));
r = ar(v(1))/(ar(v(1))+br(v(1)));

q = aq(v(1),xi(1))/(aq(v(1),xi(1))+bq(v(1),xi(1)));

while j*dt<=T
   
    Is = 0;
    if j*dt > 5 & j*dt < 25
        Is = 0.5; 		%nA
    end
    
    % gating vars
    
    tv = v(j);
    tmp = am(tv);
    m = (dt*tmp + m)/(1+dt*(tmp+bm(tv)));
    
    tmp = ah(tv);
    h = (dt*tmp+h)/(1+dt*(tmp+bh(tv)));
    
    tmp = an(tv);
    n = (dt*tmp+n)/(1+dt*(tmp+bn(tv)));
    
    tmp = ay(tv);
    y = (dt*tmp+y)/(1+dt*(tmp+by(tv)));
    
    tmp = as(tv);
    s = (dt*tmp+s)/(1+dt*(tmp+bs(tv)));
    
    tmp = ar(tv);
    r = (dt*tmp+r)/(1+dt*(tmp+br(xi(j))));
    
    tmp = aq(tv,xi(j));
    q = (dt*tmp+q)/(1+dt*(tmp+bq(tv,xi(j))));
    
    % gating conductances
    
    gpna = gna*m^3*h;
    gpk = gk*n^4*y;
    gpca = gca*s^5*r;
    gpkca = gkca*q;
    
    % the 1e-3 in the next eqn converts nanoamps to microamps
    
    xi(j+1) = (xi(j) - dt*D*(1e-3)*gpca*(tv-vca))/(1+dt*xib);
        
    top =  tv + (dt/Cm)*(Is+gL*vL+gpna*vna+gpk*vk+gpca*vca+gpkca*vk);
    
    bot = 1 + (dt/Cm)*(gL + gpna + gpk + gpca + gpkca);
    
    v(j+1) = top/bot;    

    ina(j+1) = gpna*(v(j+1)-vna);
    ik(j+1) = gpk*(v(j+1)-vk);
    ica(j+1) = gpca*(v(j+1)-vca);
    ikca(j+1) = gpkca*(v(j+1)-vk);
    
    j = j + 1;
    
end

figure(1)
plot(t,v)
xlabel('t (ms)','fontsize',16)
ylabel('v (mV)','fontsize',16)

figure(2)
plot(t,xi)
xlabel('t (ms)','fontsize',16)
ylabel('Ca (mM)','fontsize',16)

figure(3)
plot(t,ina,t,ik,t,ica)
legend('I_{Na}','I_K','I_{Ca}')
xlabel('t (ms)','fontsize',16)
ylabel('Ionic Currents  (nA)','fontsize',16)

function val =  am(v)
val = 0.32*(13-v)./(exp((13-v)./4)-1);

function val = bm(v)
val = 0.28*(v-40)./(exp((v-40)./5)-1);

function val = ah(v)
val = 0.128*exp((17-v)./18);

function val = bh(v)
val = 4./(exp((40-v)./5)+1);

function val = an(v)
val = 0.032*(15-v)./(exp((15-v)./5)-1);

function val = bn(v)
val = 0.5*exp((10-v)./40);

function val = ay(v)
val = 0.028*exp((15-v)./15) + 2/(exp((85-v)./10)+1);

function val = by(v)
val = 0.4./(exp((40-v)./10)+1);

function val = as(v)
val = 0.04*(60-v)./(exp((60-v)./10)-1);

function val = bs(v)
val = 0.005*(v-45)./(exp((v-45)./10)-1);

function val = ar(v)
val = 0.005;

function val = br(xi)
val = 0.025*(200-xi)./(exp((200-xi)./20)-1);

function val = aq(v,xi)
val = exp(v/27)*0.005*(200-xi)./(exp((200-xi)./20)-1);

function val = bq(v,xi)
val = 0.002;