%
% networkpr.m
%
% solve the pinsky-rinzel 2 compartment Ca3 model
%   for a network of num cells
%
% the state is y = [Vs Vd h n s c q Ca]
%
%  usage   networkpr(T,Is,num,M)     e.g.,  networkpr(50,[1 -.5],2,[0 0;0 1])
%
%  where   T = duration of simulation
%         Is = somatic current injection
%         num = number of cells in network
%           M = synaptic input matrix
% 
%              cellnum   
% M matrix = [recieves input | ... ]
% 
%  figure produces Vs and Vd vs time
%

function [t y] = networkpr(T,Is,num,M)

y0 = [-4.6 -4.5 0.999 0.001 0.009 0.007 0.01 0.2 0 0]; %S and W

y0 = repmat(y0,1,num);

for i = 1:length(M)
    if nnz(M(i,:)) == 0
        wsyn(i) = 1;
    else
        wsyn(i) = 20/nnz(M(i,:)); %W(i) = synaptic connections for ith cell
        % we need 20 cells input to effectively excite a cell. TEMPORARY
    end
end
M = sparse(M);

[t,y] = ode23(@red,[0 T],y0,[],Is,num,M,wsyn);

% for i = 1:round(num/4):num
%     subplot(5,1,round(4*i/num)+1)
%     plot(t,y(:,i*10+1-10),t,y(:,i*10+2-10))
% end
%
% subplot(5,1,5)
% plot(t,y(:,num*10+1-10),t,y(:,num*10+2-10))

plot(t,y(:,1:10:end),t,y(:,2:10:end))
% keyboard
legend('V_s','V_d')
xlabel('time  (ms)','fontsize',16)
ylabel('Membrane Potential  (mV)','fontsize',16)
% figure
% plot(t,y(:,7),t,y(:,8)/450)
% legend('q','Ca')

return

%
% the reduced traub model of Pinsky & Rinzel
%

function dy = red(t,y,Is,num,M,wsyn)

gL = 0.1; gNa = 30; gKDR = 15; gCa = 10; gKAHP = 0.8; gKC = 15;

VNa = 120; VCa = 140; VK = -15; VL = 0; gAMPA=0.0045; gNMDA = 0*1.75/125; Smax = 125;

% Is = -0.5;

IAMPA = gAMPA*y(10:10:(num-1)*10+10).*(y(2:10:(num-1)*10+2)-60);
INMDA = gNMDA*y(9:10:(num-1)*10+9).*(y(2:10:(num-1)*10+2)-60)./(1+0.28*exp(-0.062*(y(2:10:(num-1)*10+2)-60)));
Isyn = wsyn'.*(INMDA + IAMPA);

if y(9:10:(num-1)*10+9) > 125 
    disp('Isyn > 125')
end

Id = 0;
gc = 2.1;
p = 0.5;
Cm = 3;

dy = zeros(10*num,1);

Ils = gL*(y(1:10:(num-1)*10+1)-VL);
minf = am(y(1:10:(num-1)*10+1))./(am(y(1:10:(num-1)*10+1))+bm(y(1:10:(num-1)*10+1)));
INa = gNa*minf.^2.*y(3:10:(num-1)*10+3).*(y(1:10:(num-1)*10+1)-VNa);
IKDR = gKDR*y(4:10:(num-1)*10+4).*(y(1:10:(num-1)*10+1)-VK);

dy(1:10:(num-1)*10+1) = (-Ils - INa - IKDR + (gc/p)*(y(2:10:(num-1)*10+2)-y(1:10:(num-1)*10+1)) + Is'/p)/Cm;

Ils = gL*(y(2:10:(num-1)*10+2)-VL);
ICa = gCa*y(5:10:(num-1)*10+5).^2.*(y(2:10:(num-1)*10+2)-VCa);

IKC = gKC*y(6:10:(num-1)*10+6).*min(y(8:10:(num-1)*10+8)/250,ones(num,1)).*(y(2:10:(num-1)*10+2)-VK);
IKAHP = gKAHP*y(7:10:(num-1)*10+7).*(y(2:10:(num-1)*10+2)-VK);

dy(2:10:(num-1)*10+2) = (-Ils-ICa-IKAHP-IKC-Isyn/(1-p)+(gc/(1-p))*(y(1:10:(num-1)*10+1)-y(2:10:(num-1)*10+2))+Id/(1-p))/Cm;

dy(3:10:(num-1)*10+3) = ah(y(1:10:(num-1)*10+1)).*(1-y(3:10:(num-1)*10+3)) - bh(y(1:10:(num-1)*10+1)).*y(3:10:(num-1)*10+3);

dy(4:10:(num-1)*10+4) = an(y(1:10:(num-1)*10+1)).*(1-y(4:10:(num-1)*10+4)) - bn(y(1:10:(num-1)*10+1)).*y(4:10:(num-1)*10+4);

dy(5:10:(num-1)*10+5) = as(y(2:10:(num-1)*10+2)).*(1-y(5:10:(num-1)*10+5)) - bs(y(2:10:(num-1)*10+2)).*y(5:10:(num-1)*10+5);

dy(6:10:(num-1)*10+6) = ac(y(2:10:(num-1)*10+2)).*(1-y(6:10:(num-1)*10+6)) - bc(y(2:10:(num-1)*10+2)).*y(6:10:(num-1)*10+6);

dy(7:10:(num-1)*10+7) = aq(y(8:10:(num-1)*10+8)).*(1-y(7:10:(num-1)*10+7)) - bq(y(8:10:(num-1)*10+8)).*y(7:10:(num-1)*10+7);

dy(8:10:(num-1)*10+8) = -0.13*ICa - 0.075*y(8:10:(num-1)*10+8);

for i = 1:num
    dy((i-1)*10+9) = sum(heaviside(M(i,:)'.*y(1:10:(num-1)*10+1)-10))-y((i-1)*10+9)/150;
    dy(i*10) = sum(heaviside(M(i,:)'.*y(1:10:(num-1)*10+1)-20))-y((i-1)*10+10)/2;
end
% dy(9:10:(num-1)*10+9) = sum(heaviside(y(1:10:(num-1)*10+1)-10))-y(9:10:(num-1)*10+9)/150 - heav(y(1:10:(num-1)*10+1)-10)';
%     % I_NMDA
% dy(10:10:(num-1)*10+10) = sum(heaviside(y(1:10:(num-1)*10+1)-20))-y(10:10:(num-1)*10+10)/2- heav(y(1:10:(num-1)*10+1)-20)';
%     %I_AMPA

return

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

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

function val = an(v)
val = 0.016*(35.1-v)./(exp((35.1-v)/5)-1);

function val = bn(v)
val = 0.25*exp(0.5-0.025*v);

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

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

function val = as(v)
val = 1.6./(1+exp(-0.072*(v-65)));

function val = bs(v)
val = 0.02*(v-51.1)./(exp((v-51.1)/5)-1);

function val = ac(v)
tf= v <= 50;
val = tf.*(exp((v-10)/11-(v-6.5)/27)/18.975);
val = val +(~tf.*((2*exp((6.5-v)/27))));

function val = bc(v)
tf = v <= 50;
val = tf.*(2*exp((6.5-v)/27) - exp((v-10)/11-(v-6.5)/27)/18.975);

function val = aq(v)
val = min((0.00002)*v,0.01*ones(length(v),1));

function val = bq(v)
val = 0.001;