
% Traub.Miles active soma Inhibitory cell
% modified FengXiao code

function traubi
close all

% reversal potentials

vL = -15;
vna = 100;
vk = -30;

% parameters

gL = .01;
gna = 3.32;
gk = 3.98;

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

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

% initial values 

v = zeros(size(t));
ina = v;
ik = 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)));

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)));
    
    % gating conductances
    
    gpna = gna*m^3*h;
    gpk = gk*n^4;
    
    % the 1e-3 in the next eqn converts nanoamps to microamps
    
    top =  tv + (dt/Cm)*(Is+gL*vL+gpna*vna+gpk*vk);
    
    bot = 1 + (dt/Cm)*(gL + gpna + gpk);
    
    v(j+1) = top/bot;    

    ina(j+1) = gpna*(v(j+1)-vna);
    ik(j+1) = gpk*(v(j+1)-vk);
    
    j = j + 1;
    
end

figure(1)
plot(t,v)
xlabel('t (ms)','fontsize',16)
ylabel('v (mV)','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);
