% Pinsky_Rinzel.m

% This is for the dynamics of the two-compartment Pinsky-Rinzel model.

function [t,y] = call_RHS();
   clear;
   tmax = 1000;
   tspan = [0,tmax];

%  initial conditions
   Vs0=-70;
   Vd0=-70;
   h0=1;
   n0=0.;
   s0=0;
   c0=0;
   q0=0;
   Ca0=0.2;
   
%  initial condition vector
   y0 = [Vs0 Vd0 h0 n0 s0 c0 q0 Ca0];

%  Current pulse timing
   ton=100;
   toff=900;

%  Is=2.1, gc=8, p=0.3 give some bursts in Vs

%  Applied current parameter 
   Is=2.1; % default is 2.1 

%  Coupling conductance 
   gc=10; % book uses 0.15; default is 10

%  Fraction of area in soma
   p=0.5;  % default is 0.5
   
%  Call ODE solver, pass some parameters
   [t,y] = ode45(@(t,y)RHS(t,y,ton,toff,Is,gc,p),tspan,y0);

%  change defaults for plotting
   set(0,'DefaultAxesFontSize',16)

%  set up variable vectors
   Vs=y(:,1);
   Vd=y(:,2);
   h=y(:,3);
   n=y(:,4);
   s=y(:,5);
   c=y(:,6);
   q=y(:,7);
   Ca=y(:,8);
   
%  applied current pulse for plotting
   Nsizevec=size(n);
   Nsize=Nsizevec(1);
   tme=linspace(0,tmax,Nsize);
   for j=1:Nsize
       if tme(j) > ton & tme(j) < toff
           Iap(j)=Is;
       else
           Iap(j)=0;
       end
   end

   figure(1);
%  plot Iap time course
   ha=plot(t,Iap,'black');
   set(ha,'LineWidth',2);
   axis([0 tmax 0 0.2]);
   box off;
   xlabel('t (ms)');
   ylabel('I_s');
   
   figure(2);
%  plot soma V time course
   ha=plot(t,Vs,'black');
   set(ha,'LineWidth',2);
   axis([0 tmax -70 30]);
   box off;
   xlabel('t (ms)');
   ylabel('V_s (mV)');
   
   figure(3);
%  plot dendrite V time course
   ha=plot(t,Vd,'black');
   set(ha,'LineWidth',2);
   axis([0 tmax -70 30]);
   box off;
   xlabel('t (ms)');
   ylabel('V_d (mV)');

   figure(4);
   % plot AHP activation variable activation variable
   ha=plot(t,q,'black','Linewidth',2);
   axis([0 tmax 0 0.5]);
   box off;
   xlabel('t (ms)');
   ylabel('q');
   
   figure(5);
   % plot Ca concentration 
   ha=plot(t,Ca,'black','Linewidth',2);
   axis([0 tmax 0 300]);
   box off;
   xlabel('t (ms)');
   ylabel('Ca');
   
   figure(6);
   % plot another activation variable 
   ha=plot(t,c,'black','Linewidth',2);
   axis([0 tmax 0 1]);
   box off;
   xlabel('t (ms)');
   ylabel('c');
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Define the RHS function for the ODEs
   function dydt = RHS(t,y,ton,toff,Is,gc,p);
 
%     set up the variable vectors
      Vs=y(1);
      Vd=y(2);
      h=y(3);
      n=y(4);
      s=y(5);
      c=y(6);
      q=y(7);
      Ca=y(8);

%     Shape parameters
%     Soma
      V1=Vs+46.9;
      V2=Vs+19.9;
      V3=Vs+24.9;
      V4=Vs+40;
      V5=Vs+20;

%     Dendrite
      V6=Vd+8.9;
      V7=Vd+53.5;
      V8=Vd+50;

%     Rate coefficients
%     Soma
      alpha_m = -0.32*V1/(exp(-V1/4)-1);
      beta_m = 0.28*V2/(exp(V2/5)-1);
      alpha_h = 0.128*exp(-(43+Vs)/18);
      beta_h = 4/(exp(-V5/5)+1);
      alpha_n = -0.016*V3/(exp(-V3/5)-1);
      beta_n = 0.25*exp(-V4/40);

%     Dendrite
      alpha_s = 1.6/(exp(-0.072*(Vd-5))+1);
      beta_s = 0.02*V6/(exp(V6/5)-1);
      if Vd <= -10
         alpha_c = 0.0527*(exp(V8/11)-exp(V7/27));
         beta_c = 2*exp(-V7/27)-alpha_c;
      else
         alpha_c = 2*exp(-V7/27);
         beta_c = 0;
      end
      alpha_q = min(0.00002*Ca,0.01);
      beta_q = 0.001;

%  infinity functions 
      minf = alpha_m / (alpha_m + beta_m);
      xinf = min(Ca/250,1);

%  parameters
      Cm = 3; Isyn=0; 
      gNabar = 30; gLbar = 0.1; gDRbar = 15; gCabar = 10;
      gAHPbar = 0.8; gKCabar = 15;
      VNa = 60; VK = -75; VCa = 80; VL=-70; % paper uses VL=0

%  calculate conductances
      gNa = gNabar*(minf^2)*h;
      gDR = gDRbar*n;
      gCa = gCabar*s^2;
      gKCa = gKCabar*c*xinf;
      gAHP = gAHPbar*q;

%  calculate currents
%     Soma
      ILs=gLbar*(Vs-VL);
      INa=gNa*(Vs-VNa);
      IDR=gDR*(Vs-VK);
 
%     Dendrite
      ILd=gLbar*(Vd-VL);
      ICa=gCa*(Vd-VCa);
      IAHP=gAHP*(Vd-VK);
      IKCa=gKCa*(Vd-VK);
      
%     Coupling        
      Ic=gc*(Vd-Vs);

%     Current applied to soma
      if t > ton & t < toff
         Iin = Is;
      else
         Iin = 0;
      end

%     Differential equations
      Vsdot = (-ILs-INa-IDR+Ic/p+Iin/p)/Cm; 
      Vddot = (-ILd-ICa-IAHP-IKCa-Ic/(1-p)+Isyn/(1-p))/Cm; 
      hdot = alpha_h*(1-h)-beta_h*h;
      ndot = alpha_n*(1-n)-beta_n*n;
      sdot = alpha_s*(1-s)-beta_s*s;
      cdot = alpha_c*(1-c)-beta_c*c;
      qdot = alpha_q*(1-q)-beta_q*q;
      Cadot = -0.13*ICa-0.075*Ca;
      dydt=[Vsdot;Vddot;hdot;ndot;sdot;cdot;qdot;Cadot];

   end
end 
