% Derrick Roos
%
% P-Cell simulator
%
% vSoma = p_cell(dt,Tcutoff,Tfin,I0,plot_option)
%
% plot_option = 0 plots ion and gating channels
%               1 plots alpha and tau of gating functions in addition,
%               2 no plots come up
%
% ex: p_cell(.01,50,80,1,1)

function vSoma = p_cell(dt,Tcutoff,Tfin,I0,plot_option)

V = -100.5:1:50;
if plot_option ==1
% plot gating function
M_tau = 1./(am(V)+bm(V));
M_inf = am(V).*M_tau;
subplot(2,2,1)
plotyy(V,M_inf,V,M_tau)
title('m channel')
ylabel('M_inf')
H_tau = 1./(ah(V)+bh(V));
H_inf = ah(V).*H_tau;
subplot(2,2,2)
plotyy(V,H_inf,V,H_tau)
title('h channel')
ylabel('H_inf')
subplot(2,2,3)
N_tau = 1./(an(V)+bn(V));
N_inf = an(V).*N_tau;
plotyy(V,N_inf,V,N_tau)
title('n channel')
ylabel('N_inf')
subplot(2,2,4)
Q_tau = 1./(aq(V)+bq(V));
Q_inf = aq(V).*Q_tau;
plotyy(V,Q_inf,V,Q_tau)
title('q channel')
ylabel('Q_inf')
end

% nuimber of compartments
comps = 4;

% for P-Cell
Gm = 0.0032; %soma
Cm = 0.032; % soma
Gm(2:comps) = 0.0096;
Cm(2:comps) = 0.288;
Gcore = 0.04;
GNa = 1; % mewS
vNa = 40;
vK = -70;
GK = .5;
vCa = 150;
GCa = 0;
GKCA = 0.0017;
Gleak = Gm;
vLeak = -50;
APinflux = 4;
APdecay = 30;
v1 = -50;


N = ceil(Tfin/dt);

t = zeros(N,1);   % allocate space for long vectors
v = zeros(N,comps);
n = t;
m = t;
h = t;
q = t;
CaAP = t;
ICaAP = t;
IK = t;
INa = t;
Istim = zeros(comps,1);

n(1) = an(v1)/(an(v1)+bn(v1));
m(1) = am(v1)/(am(v1)+bm(v1));
h(1) = ah(v1)/(ah(v1)+bh(v1));
q(1) = aq(v1)/(aq(v1)+bq(v1));
CaAP(1)=0;
v(1,:) = v1;

IK(1) = GK*n(1)^4*(vK-v1);
INa(1) = GNa*m(1)^3*h(1)*(vNa-v1);
ICaAP(1) = 0;

for j=2:N,
    t(j) = (j-1)*dt;

    % pick compartment for injected current
    Istim(1) = I0*(t(j)>2)*(t(j)<Tcutoff);

    % Update Soma
    n(j) = ( n(j-1) + dt*an(v(j-1,1)) )/(1 + dt*(an(v(j-1,1))+bn(v(j-1,1))) );
    m(j) = ( m(j-1) + dt*am(v(j-1,1)) )/(1 + dt*(am(v(j-1,1))+bm(v(j-1,1))) );
    h(j) = ( h(j-1) + dt*ah(v(j-1,1)) )/(1 + dt*(ah(v(j-1,1))+bh(v(j-1,1))) );
    q(j) = ( q(j-1) + dt*aq(v(j-1,1)) )/(1 + dt*(aq(v(j-1,1))+bq(v(j-1,1))) );
    CaAP(j) = ( dt*(vCa - v(j-1,1))*APinflux*q(j)^5 + CaAP(j-1,1) )/( 1+dt*APdecay );

    top = Cm(1)*v(j-1,1)+dt*(Gleak(1)*vLeak + GNa*m(j)^3*h(j)*vNa + GK*n(j)^4*vK + GCa*q(j)^5*vCa +GKCA*vK*CaAP(j)+ Istim(1));

    bot = Cm(1) + dt*(Gleak(1) + GNa*m(j)^3*h(j) + GK*n(j)^4 + GCa*q(j)^5 + GKCA*CaAP(j));

    % Add compartmental current
    top = top + dt*Gcore*v(j-1,2);
    bot = bot + dt*Gcore;
    v(j,1) = top/bot;

    for k = 2:comps

        top = Cm(k)*v(j-1,k)+dt*(Gleak(k)*vLeak + Istim(k));
        bot = Cm(k) + dt*Gleak(k);

        % check for right neighbor
        if k < comps
            top = top + dt*Gcore*v(j-1,k+1);
            bot = bot + dt*Gcore;
        end
        % left neighbor compartment
        top = top + dt*Gcore*v(j-1,k-1);
        bot = bot + dt*Gcore;

        % update voltage
        v(j,k) = top/bot;

    end

    IK(j) = GK*n(j)^4*(vK-v(j,1));

    INa(j) = GNa*m(j)^3*h(j)*(vNa - v(j,1));

    ICaAP(j) = GKCA*CaAP(j)*(vK-v(j,1));

end

vSoma = v(:,1);

if plot_option ~= 2
figure;
subplot(2,3,1)
plot(t,v(:,1),t,v(:,comps))
legend('Soma','end dendrite')
ylabel('v (mV)','fontsize',16)

subplot(2,3,2)
plot(t,n)
ylabel('n','fontsize',16)
subplot(2,3,3)
plot(t,m)
xlabel('t (ms)','fontsize',16)
ylabel('m','fontsize',16)
subplot(2,3,4)
plot(t,h)
ylabel('h','fontsize',16)
xlabel('t (ms)','fontsize',16)
subplot(2,3,5)
plot(t,q)
ylabel('q','fontsize',16)
xlabel('t (ms)','fontsize',16)
subplot(2,3,6)
plot(t,CaAP)
ylabel('CaAP','fontsize',16)
xlabel('t (ms)','fontsize',16)
figure;
plot(t,IK,t,INa,t,ICaAP)
xlabel('t (ms)','fontsize',16)
ylabel('I (\muA/cm^2)','fontsize',16)
legend('I_K','I_{Na}','I_{CaAP}')
end

% functions for P-Cellfunction
function val = am(v)
A = 0.2; B = -40; C = 1;
val = A*(v-B)./(1-exp((B-v)/C));

function val = bm(v)
A= 0.06 ; B = -49 ; C = 20;
val = A*(B-v)./(1-exp((v-B)/C));

function val = ah(v)
A= 0.08 ; B = -40; C = 1;
val = A*(B-v)./(1-exp((v-B)/C));

function val = bh(v)
A= 0.4 ; B = -36; C = 2;
val = A./(1+exp((B-v)/C));

function val = an(v)
A= 0.02 ; B = -15; C = 0.8;
val = A*(v-B)./(1-exp((B-v)/C));

function val = bn(v)
A= 0.04 ; B = -40; C = 0.4;
val = A*(B-v)./(1-exp((v-B)/C));

function val = aq(v)
% A=  0.08; B = -10; C = 11;
A = 0.08; B = -25; C = 1;
val = A*(v-B)./(1-exp((B-v)/C));

function val = bq(v)
% A= 0.001 ; B = -10; C = 0.5;
A = 0.005; B = -20; C = 20;
val = A*(B-v)./(1-exp((v-B)/C));