% adult_HVCX.m

% This MATLAB file contains the program for an adult HVC_X neuron,
% as described in "Neuronal Intrinsic Plasticity Changes During
% Development of a Learned Behavior", Matthew T. Ross, Diana Flores,
% Richard Bertram, Frank Johnson, and Richard L. Hyson, eNeuro,
% 4(5):e0297-17.2017, 2017.

% Provides a single pulse of applied current.
% Calls the ode file hvcmodel_ode.m 


close all;
clear all;
clc;

% parameters
% Applied Current
Iap_ini = 0;
del_Iap = 200;
ntraces = 2;
   ton  = 1500;
   toff = 1800;
data_sim = zeros(10000,3,3);

% Nerst potentials in mV
VK  =   -90;
VNa =   50;
VH  =   -43;  

% Parameters that change through development
VL  =   -63;
gL   = 5;
gNap =  1;
gCaL = 1; 
gNa  = 2300;
gK   = 120;
gA   = 5;
gM   = 15.4;
gCaT = 3.8;
gSK  = 2.1;
gH   = 2.25; 
kr = 0.17; 

C   = 260;
WtauH = 1;  
sigmaM   = -4;
thetaM   = -38;    
tauNbar  = 10;
tauZ     = 75;
thetaZ   = -39;
thetaN   = -30;
sigmaN   = -6;

% Parameters fixed through development
% Time constants in mS
tauHpbar = 1000;
tauE     = 20;
tauR0    = 200;
tauR1    = 87.5;

%Function constants in mV        
thetaE   = -60;
thetaAT  = -59;
thetaRs  = -87.7;
thetaMp  = -40;
thetaHp  = -48;
thetaA   = -20;
thetaS   = -13;
thetaBT  = .4;
thetaRT  = -67;
thetaRRT = 68;
thetaRf  = -87.7;

sigmaAT  = -6;
sigmaMp  = -6;
sigmaHp  = 6;
sigmaA   = -10;
sigmaE   = 5;
sigmaZ   = 5;
sigmaS   = -8.6;
sigmaBT  = -0.1;
sigmaRT  = 2;
sigmaRRT = 2.2;
sigmaRf = 6.4;
sigmaRs = 6.4;

Ca_ex = 2.5;
RTF   = 26.7;
ks    = 0.4;
f     = 0.01;
eps   = 0.0015;
kCa   = 0.3;
bCa   = 0.05;
prf   = 100;
  
% define initial conditions
    V_0  = -52.3;
	n_0  = 0.008;
	h_0  = 0.85;
    hp_0 = 0.67;
    e_0  = 0.18;
    z_0  = 0.06;
    rT_0 = 0.0;
    Cai_0 = 0.05;
    rf_0  = 0.004;
    rs_0  = 0.004;
	
init= [V_0 n_0 h_0 hp_0 e_0 z_0 rT_0 Cai_0 rf_0 rs_0];

% call solver routine 
t0 = 0;   %initial time [msec]
tf = 3e3; %final time [msec]
ts = 0.1; %time step [msec]


for i=1:ntraces
Iap = Iap_ini + del_Iap*(i-1);
parameters = [Iap VL VK VNa VH gL gNap gCaL gNa gK gA gM gCaT gSK gH kr C tauNbar tauHpbar tauE tauZ tauR0 tauR1 thetaM thetaN thetaMp thetaHp thetaA thetaE thetaZ thetaS thetaAT thetaBT thetaRT thetaRRT thetaRf thetaRs sigmaM sigmaN sigmaMp sigmaHp sigmaA sigmaE sigmaZ sigmaS sigmaAT sigmaBT sigmaRT sigmaRRT sigmaRf sigmaRs Ca_ex RTF ks f eps kCa bCa prf ton toff WtauH];

options = odeset('InitialStep',0.01,'MaxStep',0.01);
[t,y]= ode23(@(t,y) hvcmodel_ode(t, y, parameters), [t0 tf] ,init);

V   = y(:,1);
n   = y(:,2);
h   = y(:,3);
hp  = y(:,4);
e   = y(:,5);
z   = y(:,6);
rT  = y(:,7);
Cai = y(:,8);
rf  = y(:,9);
rs  = y(:,10);


ini_steady = [V(end) n(end) h(end) hp(end) e(end) z(end) rT(end) Cai(end) rf(end) rs(end)];
 
plot(t,V)
axis([1300 2300 -100 40 ])
hold on
end


