%
% prsolve.m
%
% solve the pinsky-rinzel 2 compartment Ca3 model
%
% the state is y = [Vs Vd h n s c q Ca]
%
%  usage   prsolve(T,Is)     e.g.,  prsolve(50,1)
%
%  where   T = duration of simulation
%         Is = somatic current injection
%
%  figure produces Vs and Vd vs time
%

function prsolve(T,Is)

y0 = [-4.6 -4.5 0.999 0.001 0.009 0.007 0.01 0.2];

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

plot(t,y(:,1),t,y(:,2))
legend('V_s','V_d')
xlabel('time  (ms)','fontsize',16)
ylabel('Membrane Potential  (mV)','fontsize',16)

return

%
% the reduced traub model of Pinsky & Rinzel
%

function dy = red(t,y,Is)

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

VNa = 120; VCa = 140; VK = -15; VL = 0;

% Is = -0.5;
Isyn = 0;
Id = 0;
gc = 2.1;
p = 0.5;
Cm = 3;

dy = zeros(8,1);

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

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

Ils = gL*(y(2)-VL);
ICa = gCa*y(5).^2.*(y(2)-VCa);
IKC = gKC*y(6)*min(y(8)/250,1)*(y(2)-VK);
IKAHP = gKAHP*y(7).*(y(2)-VK);

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

dy(3) = ah(y(1))*(1-y(3)) - bh(y(1))*y(3);

dy(4) = an(y(1))*(1-y(4)) - bn(y(1))*y(4);

dy(5) = as(y(2))*(1-y(5)) - bs(y(2))*y(5);

dy(6) = ac(y(2))*(1-y(6)) - bc(y(2))*y(6);

dy(7) = aq(y(8))*(1-y(7)) - bq(y(8))*y(7);

dy(8) = -0.13*ICa - 0.075*y(8);

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)
if v <= 50
   val = exp((v-10)/11-(v-6.5)/27)/18.975;
else
   val = 2*exp((6.5-v)/27);
end

function val = bc(v)
if v <= 50
   val = 2*exp((6.5-v)/27) - exp((v-10)/11-(v-6.5)/27)/18.975;
else
   val = 0;
end

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

function val = bq(v)
val = 0.001;