classdef Neuron < handle
    properties
        input_data
        params
        u0
        index
        n_out
        idx_out
        enabled
    end
    methods
        function self = Neuron(varargin)
            self.params = parse_inputs(varargin{:});
            self.index = indices();
            self.u0 = initial_conditions(self.index);
            self.enabled = true(size(self.u0));
            [self.idx_out, self.n_out] = output_indices();
        end
        function [du, varargout] = rhs(self, t, u, NO_k, R)
            t = t(:).';
            p = self.params;
            idx = self.index;
            
    %% State variables
            
            % Elshin neuron model
            v_sa = u(idx.v_sa, :);      % membrane potential of soma/axon, mV
            v_d = u(idx.v_d, :);        % membrane potential of dendrite, mV
            K_sa = u(idx.K_sa, :);      % K+ concentration of soma/axon, mM
            Na_sa = u(idx.Na_sa, :);    % Na+ concentration of soma/axon, mM
            K_d = u(idx.K_d, :);        % K+ concentration of dendrite, mM
            Na_d = u(idx.Na_d, :);      % Na+ concentration of dendrite, mM
            K_e = u(idx.K_e, :);        % K+ concentration of ECS, mM
            Na_e = u(idx.Na_e, :);      % Na+ concentration of ECS, mM
            Buff_e = u(idx.Buff_e, :);  % Buffer concentration for K+ buffering in ECS, mM
 
            O2 = u(idx.O2, :);          % Tissue oxygen concentration, mM
            CBV = u(idx.CBV, :);        % BOLD CBV
            HbR = u(idx.HbR, :);        % BOLD deoxyhemoglobin concentration
            
            % Gating variables
            m1 = u(idx.m1, :);          % Activation gating variable, soma/axon NaP channel (Na+)
            m2 = u(idx.m2, :);          % Activation gating variable, soma/axon KDR channel (K+)
            m3 = u(idx.m3, :);          % Activation gating variable, soma/axon KA channel (K+)
            m4 = u(idx.m4, :);          % Activation gating variable, dendrite NaP channel (Na+)
            m5 = u(idx.m5, :);          % Activation gating variable, dendrite NMDA channel (Na+)
            m6 = u(idx.m6, :);          % Activation gating variable, dendrite KDR channel (K+)
            m7 = u(idx.m7, :);          % Activation gating variable, dendrite KA channel (K+)
            m8 = u(idx.m8, :);          % Activation gating variable, soma/axon NaT channel (Na+)
            
            h1 = u(idx.h1, :);          % Inactivation gating variable, soma/axon NaP channel (Na+)
            h2 = u(idx.h2, :);          % Inactivation gating variable, soma/axon KA channel (K+)
            h3 = u(idx.h3, :);          % Inactivation gating variable, dendrite NaP channel (Na+)
            h4 = u(idx.h4, :);          % Inactivation gating variable, dendrite NMDA channel (Na+)
            h5 = u(idx.h5, :);          % Inactivation gating variable, dendrite KA channel (K+)
            h6 = u(idx.h6, :);          % Inactivation gating variable, soma/axon NaT channel (Na+)
            
            % NO pathway 
            Ca_n = u(idx.Ca_n, :);                  
            nNOS_act_n = u(idx.nNOS_act_n, :);
            NO_n = u(idx.NO_n, :);
            
    %% Algebraic Expressions
            
            %% Elshin neuron model
            %% Nernst potential for Na,K ions in soma and dendrite (Cl constant)
            E_Na_sa     = p.ph * log(Na_e ./ Na_sa);
            E_K_sa      = p.ph * log(K_e ./ K_sa);
            E_Na_d      = p.ph * log(Na_e ./ Na_d);
            E_K_d       = p.ph * log(K_e ./ K_d);
            
            %% Ion fluxes
            % Leak fluxes of Na,K,Cl in soma and dendrite using HH
            J_Naleak_sa = p.gNaleak_sa * (v_sa - E_Na_sa);
            J_Kleak_sa  = p.gKleak_sa * (v_sa - E_K_sa);
            J_Naleak_d  = p.gNaleak_d * (v_d - E_Na_d);
            J_Kleak_d   = p.gKleak_d * (v_d - E_K_d);

            % Na flux through NaP channel in soma using GHK
            m1alpha     = 1 ./ (6 * (1 + exp(-((0.143 * v_sa) + 5.67))));
            m1beta      = exp(-((0.143 * v_sa) + 5.67)) ./ (6 * (1 + exp(-((0.143 * v_sa) + 5.67))));
            h1alpha     = 5.12e-8 * exp(-((0.056 * v_sa) + 2.94));
            h1beta      = 1.6e-6 ./ (1 + exp(-(((0.2 * v_sa)) + 8)));
            J_NaP_sa    = (m1.^2 .* h1 .* p.gNaP_GHk * p.Farad .* v_sa .* (Na_sa - (exp(-v_sa / p.ph) .* Na_e))) ./ (p.ph * (1 - exp(-v_sa / p.ph)));

            % Na flux through NaT channel in soma using GHK
            m8alpha     = 0.32 * ((-v_sa - 51.9) ./ (exp(-(0.25 * v_sa + 12.975)) - 1));
            m8beta      = 0.28 * ((v_sa + 24.89) ./ (exp(0.2 * v_sa + 4.978) - 1));
            h6alpha     = 0.128 * exp(-(0.056 * v_sa + 2.94));
            h6beta      = 4 ./ (1 + exp(-(0.2 * v_sa + 6)));
            J_NaT_sa    = (m8.^3 .* h6 .* p.gNaT_GHk * p.Farad .* v_sa .* (Na_sa - (exp(-v_sa / p.ph) .* Na_e))) ./ (p.ph * (1 - exp(-v_sa / p.ph)));
            
            % K flux through KDR channel in soma using GHK
            m2alpha     = 0.016 * ((v_sa + 34.9) ./ (1 - exp(-((0.2 * v_sa) + 6.98))));
            m2beta      = 0.25 * exp(-((0.025 * v_sa) + 1.25));
            J_KDR_sa    =(m2.^2 .* p.gKDR_GHk * p.Farad .* v_sa .* (K_sa - (exp(-v_sa / p.ph) .* K_e))) ./ (p.ph * (1 - exp(-v_sa / p.ph)));

            % K flux through KA channel in soma using GHKinput_current
            m3alpha     = 0.02 * ((v_sa + 56.9) ./ (1 - exp(-((0.1 * v_sa) + 5.69))));
            m3beta      = 0.0175 * ((v_sa + 29.9) ./ (exp(((0.1 * v_sa) + 2.99)) - 1));
            h2alpha     = 0.016 * exp(-((0.056 * v_sa) + 4.61));
            h2beta      = 0.5 ./ (1 + exp(-((0.2 * v_sa) + 11.98)));
            J_KA_sa     = (m3.^2 .* h2 .* p.gKA_GHk * p.Farad .* v_sa .* (K_sa - (exp(-v_sa / p.ph) .* K_e))) ./ (p.ph * (1 - exp(-v_sa / p.ph)));

            % Na flux through NaP channel in dendrite using GHK
            m4alpha     = 1 ./ (6 * (1 + exp(-((0.143 * v_d) + 5.67))));
            m4beta      = exp(-((0.143 * v_d) + 5.67)) ./ (6 * (1 + exp(-((0.143 * v_d) + 5.67))));
            h3alpha     = 5.12e-8 * exp(-((0.056 * v_d) + 2.94));
            h3beta      = 1.6e-6 ./ (1 + exp(-(((0.2 * v_d)) + 8)));
            J_NaP_d     = (m4.^2 .* h3 .* p.gNaP_GHk * p.Farad .* v_d .* (Na_d - (exp(-v_d / p.ph) .* Na_e))) ./ (p.ph * (1 - exp(-v_d / p.ph)));

            % Na/K flux through NMDA channel in dendrite using GHK
            m5alpha     = 0.5 ./ (1 + exp((13.5 - K_e) / 1.42));
            m5beta      = 0.5 - m5alpha;
            h4alpha     = 1 ./ (2000 * (1 + exp((K_e - 6.75) / 0.71)));
            h4beta      = 5e-4 - h4alpha;
            J_NMDA_K_d    = ( (m5 .* h4 .* p.gNMDA_GHk * p.Farad .* v_d .* (K_d - (exp(-v_d / p.ph) .* K_e))) ./ (p.ph * (1 - exp(-v_d / p.ph))) ) ./ (1 + 0.33 * p.Mg * exp(-(0.07 * v_d + 0.7)));
            J_NMDA_Na_d    = ( (m5 .* h4 .* p.gNMDA_GHk * p.Farad .* v_d .* (Na_d - (exp(-v_d / p.ph) .* Na_e))) ./ (p.ph * (1 - exp(-v_d / p.ph))) ) ./ (1 + 0.33 * p.Mg * exp(-(0.07 * v_d + 0.7)));

            
            % K flux through KDR channel in dendrite using GHK
            m6alpha     = 0.016 * ((v_d + 34.9) ./ (1 - exp(-((0.2 * v_d) + 6.98))));
            m6beta      = 0.25 * exp(-((0.025 * v_d) + 1.25));
            J_KDR_d     = (m6.^2 .* p.gKDR_GHk * p.Farad .* v_d .* (K_d - (exp(-v_d / p.ph) .* K_e))) ./ (p.ph * (1 - exp(-v_d / p.ph)));

            % K flux through KA channel in dendrite using GHK
            m7alpha     = 0.02 * ((v_d + 56.9) ./ (1 - exp(-((0.1 * v_d) + 5.69))));
            m7beta      = 0.0175 * ((v_d + 29.9) ./ (exp(((0.1 * v_d) + 2.99)) - 1));
            h5alpha     = 0.016 * exp(-((0.056 * v_d) + 4.61));
            h5beta      = 0.5 ./ (1 + exp(-((0.2 * v_d) + 11.98)));
            J_KA_d      = (m7.^2 .* h5 .* p.gKA_GHk * p.Farad .* v_d .* (K_d - (exp(-v_d / p.ph) .* K_e))) ./ (p.ph * (1 - exp(-v_d / p.ph)));
            
            % flux through the NaK-ATPase pump
            J_pump1_sa      = (1 + (p.K_init_e ./ K_e)).^(-2) .* (1 + (p.Na_init_sa ./ Na_sa)) .^ (-3);
            J_pump1init_sa  = 0.0312; %(1 + (p.K_init_e / p.K_init_e)).^(-2) .* (1 + (p.Na_init_sa / p.Na_init_sa)).^(-3);
            J_pump1_d       = (1 + (p.K_init_e ./ K_e)).^(-2) .* (1 + (p.Na_init_d ./ Na_d)).^(-3);
            J_pump1init_d   = 0.0312; %(1 + (p.K_init_e / p.K_init_e)).^(-2) .* (1 + (p.Na_init_d / p.Na_init_d)).^(-3);
            
            % Determine whether there is limited oxygen: O2switch=0 ATP is plentiful, O2switch=1 ATP is limited (oxygen-limited regime)
            O2_p            = p.O2_0 * (1 - p.O2switch) + O2 * p.O2switch;
            J_pump2         = 2 * (1 + p.O2_0 ./ (((1 - p.alpha_O2) * O2_p) + p.alpha_O2 * p.O2_0)).^(-1);     

            J_pump_sa   = p.Imax * J_pump1_sa .* J_pump2;
            J_pump_d    = p.Imax * J_pump1_d .* J_pump2;

            J_Napump_sa = 3 * J_pump_sa;
            J_Kpump_sa  = -2 * J_pump_sa;
            J_Napump_d  = 3 * J_pump_d;
            J_Kpump_d   = -2 * J_pump_d;  
            
            %% Total ion fluxes
            %Total ion fluxes in soma
            J_Na_tot_sa = J_NaP_sa + J_Naleak_sa + J_Napump_sa + J_NaT_sa;
            J_K_tot_sa  = J_KDR_sa + J_KA_sa + J_Kleak_sa + J_Kpump_sa;
            J_leak_tot_sa = p.gleak_sa * (v_sa - p.E_Cl_sa); % leak 
            
            % Total ion fluxes in dendrite
            J_Na_tot_d  = J_NaP_d + J_Naleak_d + J_Napump_d + J_NMDA_Na_d;
            J_K_tot_d   = J_KDR_d + J_KA_d + J_Kleak_d + J_Kpump_d + J_NMDA_K_d;
            J_leak_tot_d  = p.gleak_d * (v_d - p.E_Cl_d); % leak
                        
            % Total ion fluxes in soma and dendrite
            J_tot_sa    = J_Na_tot_sa + J_K_tot_sa + J_leak_tot_sa;
            J_tot_d     = J_Na_tot_d + J_K_tot_d + J_leak_tot_d;

            %% Tissue oxygen 
            J_pump2_0       = 0.0952;   % 2 * (1 + p.O2_0 ./ (((1 - p.alpha_O2) * 0) + p.alpha_O2 * p.O2_0)).^(-1);
            J_pump2_O2_0    = 1;        % 2 * (1 + p.O2_0 ./ (((1 - p.alpha_O2) * p.O2_0) + p.alpha_O2 * p.O2_0)).^(-1);
            P_02            = (J_pump2 - J_pump2_0 ) ./ ( J_pump2_O2_0 - J_pump2_0);            
            CBF             = p.CBF_init * (R.^4 / p.R_init^4);
            J_O2_vascular   = CBF .* ((p.O2_b - O2) ./ (p.O2_b - p.O2_0));
            J_O2_background = p.CBF_init * P_02 * (1 - p.gamma_O2);
            J_O2_pump       = p.CBF_init * P_02 * p.gamma_O2 .* ((J_pump1_sa + J_pump1_d) ./ (J_pump1init_sa + J_pump1init_d));
            
            %% BOLD
            f_out           = CBV.^(1/p.d) + p.tau_TAT * (1/(p.tau_MTT + p.tau_TAT) .* ( CBF/p.CBF_init  - CBV.^(1/p.d) ));
            CMRO2           = J_O2_background + J_O2_pump;
            CMRO2_init      = p.CBF_init * P_02;
            OEF             = CMRO2 * p.E_0 ./ CBF;
            %BOLD            = p.V_0 * ( p.a_1 * (1 - HbR) - p.a_2 * (1 - CBV) );
            
            %% NO pathway
            
            % Glutamate input: vesicle released when the extracellular K+ is over 5.5 mM (Ke_switch)
            Glu = self.shared(t, u);
            
            % NO pathway
            w_NR2A = Glu ./ (p.K_mA + Glu); %[-] 
            w_NR2B = Glu ./ (p.K_mB + Glu); %[-]
          
            I_Ca = (-4 * p.v_n * p.G_M * p.P_Ca_P_M * (p.Ca_ex / p.M)) / (1 + exp(-80 * (p.v_n + 0.02))) * (exp(2 * p.v_n * p.F / (p.R_gas * p.T))) / (1 - exp(2 * p.v_n * p.F / (p.R_gas * p.T))); %[fA]
            I_Ca_tot = I_Ca .* (p.n_NR2A * w_NR2A + p.n_NR2B * w_NR2B); %[fA]
            
            CaM = Ca_n / p.m_c; %[uM]
            tau_nk = p.x_nk ^ 2 ./  (2 * p.D_cNO); %[s]
            
            p_NO_n = p.NOswitch * ( nNOS_act_n * p.V_max_NO_n * p.O2_n / (p.K_mO2_n + p.O2_n) * p.LArg_n / (p.K_mArg_n + p.LArg_n) ); %[uM/s]
            c_NO_n = p.k_O2_n * NO_n.^2 * p.O2_n; %[uM/s]
            d_NO_n = (NO_k - NO_n) ./ tau_nk; %[uM/s]
            
            % ECS electrodiffusion TEST            
            Delta_K_e = 3.493 - K_e;
            average_K_e = (3.493 + K_e)/2;
            Delta_Na_e = 150 - Na_e;
            average_Na_e = (150 + Na_e)/2;
            Delta_v_e = -(8.315 * 300) / 96.485 * ( 1 * 3.6e-9 * Delta_K_e + 1 * 2.4e-9 * Delta_Na_e ) / ( 1 * 3.6e-9 * average_K_e + 1 * 2.4e-9 * average_Na_e );
            flu_diff_Fick_K = (3.6e-9 / (1.24e-4)^2) * (Delta_K_e);
            flu_diff_Fick_Na = (2.4e-9 / (1.24e-4)^2) * (Delta_Na_e);
            flu_diff_el_K 	= (3.6e-9 / (1.24e-4)^2) * ( ((1 * 96.485)/(8.315 * 300)) * ( average_K_e * Delta_v_e ) );
            flu_diff_el_Na = (2.4e-9 / (1.24e-4)^2) * (((1 * 96.485)/(8.315 * 300)) * ( average_Na_e * Delta_v_e ) );
            
            flu_diff_K 	= flu_diff_Fick_K + flu_diff_el_K;
            flu_diff_Na = flu_diff_Fick_Na + flu_diff_el_Na;
            
            
      %% Conservation equations                                                         
            
            %% Elshin neuron model
            
            % change in membrane potential
            du(idx.v_sa, :)     = 1/p.Cm * ( -J_tot_sa + 1 / (2 * p.Ra * p.dhod.^2) * (v_d - v_sa) + self.input_current(t) );
            du(idx.v_d, :)      = 1/p.Cm * (-J_tot_d + 1 / (2 * p.Ra * p.dhod.^2) * (v_sa - v_d));

            %change in concentration of Na,K in the soma
            du(idx.Na_sa, :)    = -p.As / (p.Farad * p.Vs) * J_Na_tot_sa + p.D_Na * (p.Vd + p.Vs) ./ (2 * p.dhod.^2 * p.Vs) * (Na_d - Na_sa);
            du(idx.K_sa, :)     = -p.As / (p.Farad * p.Vs) * J_K_tot_sa + p.D_K * (p.Vd + p.Vs) ./ (2 * p.dhod.^2 * p.Vs) * (K_d - K_sa);

            %change in concentration of Na,K in the dendrite
            du(idx.Na_d, :)     = -p.Ad / (p.Farad * p.Vd) * J_Na_tot_d + p.D_Na * (p.Vs + p.Vd) ./ (2 * p.dhod.^2 * p.Vd) * (Na_sa - Na_d);
            du(idx.K_d, :)      = -p.Ad / (p.Farad * p.Vd) * J_K_tot_d + p.D_K * (p.Vs + p.Vd) ./ (2 * p.dhod.^2 * p.Vd) * (K_sa - K_d);

            % change in buffer for K+ in the extracellular space
            du(idx.Buff_e, :)   = p.Mu * K_e .* (p.B0 - Buff_e) ./ (1 + exp(-((K_e - 5.5) ./ 1.09))) - (p.Mu * Buff_e);

            % change in concentration of Na,K in the extracellular space 
            du(idx.Na_e, :)     = 1/(p.Farad * p.fe) * (((p.As * J_Na_tot_sa) / p.Vs) + ((p.Ad * J_Na_tot_d) / p.Vd)) + flu_diff_Na;
            du(idx.K_e, :)      = 1/(p.Farad * p.fe) * (((p.As * J_K_tot_sa) / p.Vs)  + ((p.Ad * J_K_tot_d) / p.Vd)) - du(idx.Buff_e) + flu_diff_K;
            
            % change in tissue oxygen
            du(idx.O2, :)       = J_O2_vascular - J_O2_background - J_O2_pump;
            
            % change in BOLD response
            du(idx.CBV, :)       = 1/(p.tau_MTT + p.tau_TAT) .* ( CBF/p.CBF_init  - CBV.^(1/p.d) ); 
            du(idx.HbR, :)       = 1/p.tau_MTT * ( CMRO2./CMRO2_init - HbR./CBV .* f_out );

            % Change in activation gating variables m
            du(idx.m1, :)       = 1000 * ((m1alpha .* (1 - m1)) - (m1beta .* m1));
            du(idx.m2, :)       = 1000 * ((m2alpha .* (1 - m2)) - (m2beta .* m2));
            du(idx.m3, :)       = 1000 * ((m3alpha .* (1 - m3)) - (m3beta .* m3));
            du(idx.m4, :)       = 1000 * ((m4alpha .* (1 - m4)) - (m4beta .* m4));
            du(idx.m5, :)       = 1000 * ((m5alpha .* (1 - m5)) - (m5beta .* m5));
            du(idx.m6, :)       = 1000 * ((m6alpha .* (1 - m6)) - (m6beta .* m6));
            du(idx.m7, :)       = 1000 * ((m7alpha .* (1 - m7)) - (m7beta .* m7));            
            du(idx.m8, :)       = 1000 * ((m8alpha .* (1 - m8)) - (m8beta .* m8));

            % Change in inactivation gating variables h
            du(idx.h1, :)       = 1000 * ((h1alpha .* (1 - h1)) - (h1beta .* h1));
            du(idx.h2, :)       = 1000 * ((h2alpha .* (1 - h2)) - (h2beta .* h2));
            du(idx.h3, :)       = 1000 * ((h3alpha .* (1 - h3)) - (h3beta .* h3));
            du(idx.h4, :)       = 1000 * ((h4alpha .* (1 - h4)) - (h4beta .* h4));
            du(idx.h5, :)       = 1000 * ((h5alpha .* (1 - h5)) - (h5beta .* h5));
            du(idx.h6, :)       = 1000 * ((h6alpha .* (1 - h6)) - (h6beta .* h6));
            
            % NO pathway            
            du(idx.Ca_n, :) = (I_Ca_tot / (2 * p.F * p.V_spine) - (p.k_ex * (Ca_n - p.Ca_rest))) / (1 + p.lambda_buf); %[uM/s]
            du(idx.nNOS_act_n, :) = p.V_maxNOS * CaM ./ (p.K_actNOS + CaM) - p.mu2_n * nNOS_act_n; %[uM/s]
            du(idx.NO_n, :) = p_NO_n - c_NO_n + d_NO_n; %[uM/s]
            
            
            du = bsxfun(@times, self.enabled, du);
            if nargout == 2
            	Uout = zeros(self.n_out, size(u, 2));
                Uout(self.idx_out.Glu, :) = Glu; %self.input_Glu(t);
            	Uout(self.idx_out.w_NR2A, :) = w_NR2A;
            	Uout(self.idx_out.w_NR2B, :) = w_NR2B;
            	Uout(self.idx_out.I_Ca, :) = I_Ca;
            	Uout(self.idx_out.I_Ca_tot, :) = I_Ca_tot;
            	Uout(self.idx_out.CaM, :) = CaM;
            	Uout(self.idx_out.tau_nk, :) = tau_nk;
            	Uout(self.idx_out.p_NO_n, :) = p_NO_n;
            	Uout(self.idx_out.c_NO_n, :) = c_NO_n;
            	Uout(self.idx_out.d_NO_n, :) = d_NO_n;
                Uout(self.idx_out.current, :) = self.input_current(t);
                Uout(self.idx_out.CBF, :) = CBF;
                Uout(self.idx_out.OEF, :) = OEF;
                Uout(self.idx_out.CMRO2, :) = CMRO2;
                Uout(self.idx_out.CMRO2_init, :) = CMRO2_init;
            	Uout(self.idx_out.J_KDR_sa, :) = J_KDR_sa;
            	Uout(self.idx_out.J_KA_sa, :) = J_KA_sa;
                Uout(self.idx_out.J_Kleak_sa, :) = J_Kleak_sa;
                Uout(self.idx_out.J_Kpump_sa, :) = J_Kpump_sa;
                Uout(self.idx_out.J_KDR_d, :) = J_KDR_d;
                Uout(self.idx_out.J_KA_d, :) = J_KA_d;
                Uout(self.idx_out.J_Kleak_d, :) = J_Kleak_d;
                Uout(self.idx_out.J_Kpump_d, :) = J_Kpump_d;
                Uout(self.idx_out.J_NMDA_K_d, :) = J_NMDA_K_d;
                Uout(self.idx_out.flu_diff_K, :) = flu_diff_K;
                Uout(self.idx_out.flu_diff_Fick_K, :) = flu_diff_Fick_K;
                Uout(self.idx_out.flu_diff_el_K, :) = flu_diff_el_K;
            	varargout = {Uout};
            end
        end
        
       function [Glu, J_K_NEtoSC, J_Na_NEtoSC, NO_n, O2] = shared(self, t, u)    %shared variables
            t = t(:).';
            p = self.params;
            idx = self.index;
            NO_n = u(idx.NO_n, :);
            Buff_e = u(idx.Buff_e, :);  % Buffer concentration for K+ buffering in ECS, mM
            v_sa = u(idx.v_sa, :);      % membrane potential of soma/axon, mV
            v_d = u(idx.v_d, :);        % membrane potential of dendrite, mV
            K_sa = u(idx.K_sa, :);      % K+ concentration of soma/axon, mM
            Na_sa = u(idx.Na_sa, :);    % Na+ concentration of soma/axon, mM
            K_d = u(idx.K_d, :);        % K+ concentration of dendrite, mM
            Na_d = u(idx.Na_d, :);      % Na+ concentration of dendrite, mM
            K_e = u(idx.K_e, :);        % K+ concentration of ECS, mM
            Na_e = u(idx.Na_e, :);      % Na+ concentration of ECS, mM
            m1 = u(idx.m1, :);          % Activation gating variable, soma/axon NaP channel (Na+)
            m2 = u(idx.m2, :);          % Activation gating variable, soma/axon KDR channel (K+)
            m3 = u(idx.m3, :);          % Activation gating variable, soma/axon KA channel (K+)
            m4 = u(idx.m4, :);          % Activation gating variable, dendrite NaP channel (Na+)
            m5 = u(idx.m5, :);          % Activation gating variable, dendrite NMDA channel (Na+)
            m6 = u(idx.m6, :);          % Activation gating variable, dendrite KDR channel (K+)
            m7 = u(idx.m7, :);          % Activation gating variable, dendrite KA channel (K+)
            m8 = u(idx.m8, :);          % Activation gating variable, soma/axon NaT channel (Na+)
            h1 = u(idx.h1, :);          % Inactivation gating variable, soma/axon NaP channel (Na+)
            h2 = u(idx.h2, :);          % Inactivation gating variable, soma/axon KA channel (K+)
            h3 = u(idx.h3, :);          % Inactivation gating variable, dendrite NaP channel (Na+)
            h4 = u(idx.h4, :);          % Inactivation gating variable, dendrite NMDA channel (Na+)
            h5 = u(idx.h5, :);          % Inactivation gating variable, dendrite KA channel (K+)
            h6 = u(idx.h6, :);          % Inactivation gating variable, soma/axon NaT channel (Na+)
            O2 = u(idx.O2, :);
            
            %% Glutamate function dependent on neuron membrane potential
            Glu = p.GluSwitch * 0.5 * p.Glu_max * ( 1 + tanh( (K_e - p.Ke_switch) / p.Glu_slope) );  % based on extracellular K+ (Kager2000) - smooth
            
            %% Fluxes for the ODEs - (unfortunately) must be reproduced in this function to share with Astrocyte.m
            K = K_e;             
            E_Na_sa     = p.ph * log(Na_e ./ Na_sa);
            E_K_sa      = p.ph * log(K_e ./ K_sa);
            E_Na_d      = p.ph * log(Na_e ./ Na_d);
            E_K_d       = p.ph * log(K_e ./ K_d);
            
            J_Kleak_sa  = p.gKleak_sa * (v_sa - E_K_sa);
            J_KDR_sa    =(m2.^2 .* p.gKDR_GHk * p.Farad .* v_sa .* (K_sa - (exp(-v_sa / p.ph) .* K))) ./ (p.ph * (1 - exp(-v_sa / p.ph)));
            J_KA_sa     = (m3.^2 .* h2 .* p.gKA_GHk * p.Farad .* v_sa .* (K_sa - (exp(-v_sa / p.ph) .* K))) ./ (p.ph * (1 - exp(-v_sa / p.ph)));
            J_Kleak_d   = p.gKleak_d * (v_d - E_K_d);
            J_NMDA_K_d  = ( (m5 .* h4 .* p.gNMDA_GHk * p.Farad .* v_d .* (K_d - (exp(-v_d / p.ph) .* K))) ./ (p.ph * (1 - exp(-v_d / p.ph))) ) ./ (1 + 0.33 * p.Mg * exp(-(0.07 * v_d + 0.7)));
            J_KDR_d     = (m6.^2 .* p.gKDR_GHk * p.Farad .* v_d .* (K_d - (exp(-v_d / p.ph) .* K))) ./ (p.ph * (1 - exp(-v_d / p.ph)));
            J_KA_d      = (m7.^2 .* h5 .* p.gKA_GHk * p.Farad .* v_d .* (K_d - (exp(-v_d / p.ph) .* K))) ./ (p.ph * (1 - exp(-v_d / p.ph)));
   
            J_NaP_sa    = (m1.^2 .* h1 .* p.gNaP_GHk * p.Farad .* v_sa .* (Na_sa - (exp(-v_sa / p.ph) .* Na_e))) ./ (p.ph * (1 - exp(-v_sa / p.ph)));
            J_Naleak_sa = p.gNaleak_sa * (v_sa - E_Na_sa);
            J_NaT_sa    = (m8.^3 .* h6 .* p.gNaT_GHk * p.Farad .* v_sa .* (Na_sa - (exp(-v_sa / p.ph) .* Na_e))) ./ (p.ph * (1 - exp(-v_sa / p.ph)));
            J_NaP_d     = (m4.^2 .* h3 .* p.gNaP_GHk * p.Farad .* v_d .* (Na_d - (exp(-v_d / p.ph) .* Na_e))) ./ (p.ph * (1 - exp(-v_d / p.ph)));
            J_Naleak_d  = p.gNaleak_d * (v_d - E_Na_d);
            J_NMDA_Na_d    = ( (m5 .* h4 .* p.gNMDA_GHk * p.Farad .* v_d .* (Na_d - (exp(-v_d / p.ph) .* Na_e))) ./ (p.ph * (1 - exp(-v_d / p.ph))) ) ./ (1 + 0.33 * p.Mg * exp(-(0.07 * v_d + 0.7)));
            
            % Oxygen dependent ATP pump
            O2_p        = p.O2_0 * (1 - p.O2switch) + O2 * p.O2switch;
            J_pump2     = 2 * (1 + p.O2_0 ./ (((1 - p.alpha_O2) * O2_p) + p.alpha_O2 * p.O2_0)).^(-1); 
            J_pump1_sa  = (1 + (p.K_init_e ./ K)).^(-2) .* (1 + (p.Na_init_sa ./ Na_sa)) .^ (-3);
            J_pump1_d   = (1 + (p.K_init_e ./ K)).^(-2) .* (1 + (p.Na_init_d ./ Na_d)).^(-3);
            J_pump_sa   = p.Imax * J_pump1_sa .* J_pump2;
            J_pump_d    = p.Imax * J_pump1_d .* J_pump2;
            J_Kpump_sa  = -2 * J_pump_sa;
            J_Kpump_d   = -2 * J_pump_d;  
            J_Napump_sa = 3 * J_pump_sa;
            J_Napump_d = 3 * J_pump_d;

            J_K_tot_sa  = J_KDR_sa + J_KA_sa + J_Kleak_sa + J_Kpump_sa;
            J_K_tot_d   = J_KDR_d + J_KA_d + J_Kleak_d + J_Kpump_d + J_NMDA_K_d;
            J_Na_tot_sa = J_NaP_sa + J_Naleak_sa + J_Napump_sa + J_NaT_sa;
            J_Na_tot_d  = J_NaP_d + J_Naleak_d + J_Napump_d + J_NMDA_Na_d;
            
            Delta_K_e = 3.493 - K_e;
            average_K_e = (3.493 + K_e)/2;
            Delta_Na_e = 150 - Na_e;
            average_Na_e = (150 + Na_e)/2;
            Delta_v_e = -(8.315 * 300) / 96.485 * ( 1 * 3.6e-9 * Delta_K_e + 1 * 2.4e-9 * Delta_Na_e ) / ( 1 * 3.6e-9 * average_K_e + 1 * 2.4e-9 * average_Na_e );
            flu_diff_Fick_K = (3.6e-9 / (1.24e-4)^2) * (Delta_K_e);
            flu_diff_el_K 	= (3.6e-9 / (1.24e-4)^2) * ( ((1 * 96.485)/(8.315 * 300)) * ( average_K_e * Delta_v_e ) );
            flu_diff_K 	= flu_diff_Fick_K + flu_diff_el_K;
            
            dKedt = 1/(p.Farad * p.fe) * (((p.As * J_K_tot_sa) / p.Vs)  + ((p.Ad * J_K_tot_d) / p.Vd)) - (p.Mu * K_e .* (p.B0 - Buff_e) ./ (1 + exp(-((K_e - 5.5) ./ 1.09))) - (p.Mu * Buff_e)) + flu_diff_K;
            J_K_NEtoSC = p.SC_coup * dKedt;
            
            dNaedt = 1/(p.Farad * p.fe) * (((p.As * J_Na_tot_sa) / p.Vs) + ((p.Ad * J_Na_tot_d) / p.Vd));
%             J_Na_NEtoSC = p.SC_coup * dNaedt;
            J_Na_NEtoSC = -p.SC_coup * dKedt;
       end
       
              %% Current input to neuron
       function current = input_current(self, t) 
            p = self.params;
            
            if p.CurrentType == 1
                
            	current = p.Istrength * rectpuls(t - (p.t_0 + p.lengthpulse/2), p.lengthpulse);  

            elseif p.CurrentType == 2
                
                current = p.Istrength * rectpuls(t - (p.t_0 + p.lengthpulse/2), p.lengthpulse) ...
                        + p.Istrength * rectpuls(t - (p.t_0 + p.lengthpulse + 8 + 1/2), 1);
           
            elseif p.CurrentType == 3 || p.CurrentType == 4
                
            % Load experimental data and use as a scaled current input         
                dt = p.dt; XLIM2 = p.XLIM2;
                T = 0:dt:XLIM2;
                neural_data = self.input_data;
                index_t = round(t/dt + 1);  % Find index corresponding to time t
                current = p.Istrength * neural_data(index_t);                 
            end
       end        
        
        function f = input_ECS(self, t)
            % Input of K+ into the ECS, if you want to use set t0 and tend
            % here and turn off other inputs
            p = self.params;
            f = zeros(size(t));
            lengtht = p.ECS_length;
                    f(p.t0_ECS <= t & t < p.t0_ECS + lengtht ) = p.ECS_input;
        end
        
        function names = varnames(self)
            names = [fieldnames(self.index); fieldnames(self.idx_out)];
        end
    end
end    
        
function idx = indices()    %for state variables
    idx.CBV = 1;
    idx.HbR = 2;   
    idx.v_sa = 3;
    idx.v_d = 4;
    idx.K_sa = 5;
    idx.Na_sa = 6;
    idx.K_d = 7;
    idx.Na_d = 8;
    idx.K_e = 9;
    idx.Na_e = 10;
    idx.Buff_e = 11;
    idx.O2 = 12;
    idx.m1 = 13;
    idx.m2 = 14;
    idx.m3 = 15;
    idx.m4 = 16;
    idx.m5 = 17;
    idx.m6= 18;
    idx.m7 = 19;
    idx.m8 = 20;
    idx.h1 = 21;
    idx.h2 = 22;
    idx.h3 = 23;
    idx.h4 = 24;
    idx.h5 = 25;
    idx.h6 = 26;
    idx.Ca_n = 27;                  
    idx.nNOS_act_n = 28;
    idx.NO_n = 29;
end        

function [idx, n] = output_indices()    %for variables in nargout loop
    idx.Glu = 3; 
    idx.w_NR2A = 4; 
    idx.w_NR2B = 5; 
    idx.I_Ca = 6; 
    idx.I_Ca_tot = 7; 
    idx.CaM = 10; 
    idx.tau_nk = 11; 
    idx.p_NO_n = 12;    
    idx.c_NO_n = 13; 
    idx.d_NO_n = 14;       
    idx.current = 15;
    idx.CBF = 16;
    idx.OEF = 17;
    idx.CMRO2 = 19;
    idx.CMRO2_init = 20;
    idx.J_KDR_sa = 21;
    idx.J_KA_sa = 22;
    idx.J_Kleak_sa = 23;
    idx.J_Kpump_sa = 24;
    idx.J_KDR_d = 25;
    idx.J_KA_d = 26;
    idx.J_Kleak_d = 27;
    idx.J_Kpump_d = 28;
    idx.J_NMDA_K_d = 18;
    idx.flu_diff_K = 19;
    idx.flu_diff_Fick_K = 20;
    idx.flu_diff_el_K = 21;
    n = numel(fieldnames(idx));     
end
      
function params = parse_inputs(varargin)
    parser = inputParser();
    
    parser.addParameter('CurrentType', 1);  % Type of current input. 1: normal, 2: two stimulations, 3: obtained from data, 4: data + LC pain pathway
    
    
    parser.addParameter('GluSwitch', 1); 
    parser.addParameter('KSwitch', 1); 
    parser.addParameter('NOswitch', 1); 
    parser.addParameter('O2switch', 1); 
    
    % Glutamate parameters
    parser.addParameter('Glu_max', 1846);   % [uM] (one vesicle, Santucci2008)
    parser.addParameter('Glu_slope', 0.1);  % [mM] Slope of sigmoidal (M.E.)
    parser.addParameter('Ke_switch', 5.5);   % [mM] Threshold past which glutamate vesicle is released (M.E.)
    
    % Elshin model constants
    parser.addParameter('Istrength', 0.022);    % Current input strength
    parser.addParameter('SC_coup', 11.5);        % [-] Coupling scaling factor, values can range between 1.06 to 14.95 according to estimation from experimental data of Ventura and Harris, Maximum dilation at 9.5
       
    % BOLD constants
    parser.addParameter('tau_MTT', 3);              % [s] Transit time
    parser.addParameter('tau_TAT', 20);        
    parser.addParameter('d', 0.4); 
    parser.addParameter('a_1', 3.4);              % Buxton
    parser.addParameter('a_2', 1); 
    parser.addParameter('V_0', 0.03); 
    parser.addParameter('E_0', 0.4); 
    
    % global constants
    parser.addParameter('F', 9.65e4); %C mol^-1; Faraday's constant
    parser.addParameter('R_gas', 8.315); %J mol^-1 K^-1; Gas constant
    parser.addParameter('T', 300); % K; Temperature
    
    % input 
    parser.addParameter('startpulse', 200);    
    parser.addParameter('lengthpulse', 200);
    parser.addParameter('lengtht1', 10);
        
    % Elshin neuron model
    parser.addParameter('E_Cl_sa', -70);        % Nernst potential for Cl- in soma      [mV]
    parser.addParameter('E_Cl_d', -70);         % Nernst potential for Cl- in dendrite   [mV]
    parser.addParameter('Ra', 1.83e5);          % Input resistance of dendritic tree [ohms]
    parser.addParameter('dhod', 4.5e-2);        % Half length of dendrite [cm]
    parser.addParameter('As', 1.586e-5);        % Surface area of soma [cm2]
    parser.addParameter('Ad', 2.6732e-4);       % Surface area of dendrite [cm2]
    parser.addParameter('Vs', 2.16e-9);         % Volume of soma [cm3]
    parser.addParameter('Vd', 5.614e-9);        % Volume of dendrite [cm3]
    parser.addParameter('fe', 0.15);            % ECS to neuron ratio
    parser.addParameter('Cm', 7.5e-7);          % Membrance capacitance [s/ohms cm2]    *******
    parser.addParameter('Farad', 96.485);           % Faraday constant [C / mmol]
    parser.addParameter('ph', 26.6995);         % RT/F
    parser.addParameter('Mu', 8e-4);            % [m/s]
    parser.addParameter('Mu2', 8e-4);            % [m/s]
    parser.addParameter('B0', 500);             % Effective total buffer concentration [mM]
    
    parser.addParameter('gNaP_GHk', 2e-6);
    parser.addParameter('gKDR_GHk', 10e-5);
    parser.addParameter('gKA_GHk', 1e-5);
    parser.addParameter('gNMDA_GHk', 1e-5);   
    parser.addParameter('gNaT_GHk', 10e-5); 
    
%    parser.addParameter('gNaleak_sa', 9.5999e-6);   
%    parser.addParameter('gKleak_sa', 3.4564e-5);
%    parser.addParameter('gleak_sa', 10*9.5999e-6);
%    parser.addParameter('gNaleak_d', 1.0187e-5);   
%    parser.addParameter('gKleak_d', 3.4564e-5); 
%   parser.addParameter('gleak_d', 10*1.0187e-5);  
%   parser.addParameter('Imax', 0.013);  
     
%     parser.addParameter('gNaleak_sa', 4.1333e-5);   
%     parser.addParameter('gKleak_sa', 1.4623e-4);
%     parser.addParameter('gleak_sa', 10*4.1333e-5);
%     parser.addParameter('gNaleak_d', 4.1915e-5);   
%     parser.addParameter('gKleak_d', 1.4621e-4); 
%     parser.addParameter('gleak_d', 10*4.1915e-5);  
%     parser.addParameter('Imax', 0.013*4);  
    
     parser.addParameter('gNaleak_sa', 6.2378e-5);   
     parser.addParameter('gKleak_sa', 2.1989e-4);
     parser.addParameter('gleak_sa', 10*6.2378e-5);
     parser.addParameter('gNaleak_d', 6.2961e-5);   
     parser.addParameter('gKleak_d', 2.1987e-4); 
     parser.addParameter('gleak_d', 10*6.2961e-5);  
     parser.addParameter('Imax', 0.013*6);  

    parser.addParameter('O2_0', 0.02);          % [mM]
    parser.addParameter('alpha_O2',  0.05);         % percentage of ATP production independent of O2
    parser.addParameter('D_Na', 1.33e-5);       % [cm2 / s]
    parser.addParameter('D_K', 1.96e-5);        % [cm2 / s]
    
    parser.addParameter('K_init_e', 2.9); 
    parser.addParameter('Na_init_sa', 10); 
    parser.addParameter('Na_init_d', 10); 
    
    parser.addParameter('R_init', 1.9341e-5); 
    parser.addParameter('CBF_init', 3.2e-2); 
    parser.addParameter('O2_b', 4e-2);   
    parser.addParameter('gamma_O2', 0.10);  
    parser.addParameter('Mg', 1.2);  

    % NO pathway 
    parser.addParameter('m_c', 4);              % [-] Number of Ca2+ bound per calmodulin (approximated as parameter, originally an algebraic variable that changed from 3.999 to 4)
    
    parser.addParameter('K_mA', 650);           % [uM] - fit to Santucci2008
    parser.addParameter('K_mB', 2800);          % [uM] - fit to Santucci2008

    parser.addParameter('v_n', -0.04);          % [V] ; the neuronal membrane potential , assumed to be approx constant in this model
    parser.addParameter('G_M', 46000);          % [fS]! was 46 pS! ; the conductance of the NMDA channel to Ca2+ compaired  
    parser.addParameter('P_Ca_P_M', 3.6);       % [-] ; the relative conductance of the NMDA channel to Ca2+ compared to monovalent ions
    parser.addParameter('Ca_ex', 2e3);          % [microM] ; the external calcium concentration (in Comerford+David2008: 1.5 mM!)
    parser.addParameter('M', 1.3e5);            % [microM] ; the concentration of monovalent ions in the neuron

    parser.addParameter('n_NR2A', 0.63);        % [-] ; average number of NR2A NMDA receptors per synapse (Santucci2008)
    parser.addParameter('n_NR2B', 11);          % [-] ; average number of NR2B NMDA receptors per synapse (Santucci2008)
    parser.addParameter('V_max_NO_n', 4.22);    % [s^-1] ; maximum catalytic rate of NO production (Chen2006) - obtained from fig 6 & equ 17 & 18
    parser.addParameter('O2_n', 200);           % [uM] ; tissue O2 concentration in the neuron (M.E.)
    parser.addParameter('K_mO2_n', 243);        % [uM] ; Chen2006
    parser.addParameter('LArg_n', 100);         % [uM] ; 
    parser.addParameter('K_mArg_n', 1.5);       % [uM] ; 
    
    parser.addParameter('k_O2_n', 9.6e-6);      % [uM^-2 s^-1] ; % (Kavdia2002)
    parser.addParameter('x_nk', 25);            % [um] ;  (M.E.)
    parser.addParameter('D_cNO', 3300);         % [um^2 s^-1] ; Diffusion coefficient NO (Malinski1993)
    
    parser.addParameter('V_spine', 8e-8);       % [nL] ; volume of the neuronal dendritic spine Santucci2008
    parser.addParameter('k_ex', 1600);          % [s^-1] ; decay rate constant of internal calcium concentration Santucci2008
    parser.addParameter('Ca_rest', 0.1);        % [uM] ; resting calcium concentration (in Comerford+David2008: 2.830 mM; in Santucci2008P: 0.1 \muM)
    parser.addParameter('lambda_buf', 20);      % [-] ; buffer capacity Santucci2008

    parser.addParameter('V_maxNOS', 25e-3);     % [] ; M.E.
    parser.addParameter('K_actNOS', 9.27e-2);   % [uM] ; 
    parser.addParameter('mu2_n', 0.0167);       % [s^-1] ; rate constant at which the nNOS is deactivated Comerford2008

    parser.addParameter('dt', 0.001);
    parser.addParameter('XLIM2', 150);
    parser.addParameter('input_data', linspace(0, 1200, 1000)); 
    
    parser.parse(varargin{:})
    params = parser.Results;
    params.t_0 = params.startpulse;
    params.t_1 = params.t_0 + params.lengtht1;
    params.t_2 = params.t_0 + params.lengthpulse;
    params.t_3 = params.t_1 + params.lengthpulse;
end

function u0 = initial_conditions(idx)
    u0 = zeros(length(fieldnames(idx)), 1);

    u0(idx.CBV) = 1.3167;       % [-]
    u0(idx.HbR) = 0.6665;       % [-]
    u0(idx.v_sa) = -70.0337;    % [mV]
    u0(idx.v_d) = -70.0195;     % [mV]
    u0(idx.K_sa) = 134.1858;    % [mM]
    u0(idx.Na_sa) = 9.2691;     % [mM]
    u0(idx.K_d) = 134.4198;     % [mM]
    u0(idx.Na_d) = 9.3203;      % [mM]
    u0(idx.K_e) = 3.493;        % [mM]
    u0(idx.Na_e) = 150;         % [mM]
    u0(idx.Buff_e) = 165.9812;  % [mM]
    u0(idx.O2) = 0.0281;        % [mM]    
    u0(idx.m1) = 0.01281;       % [-]
    u0(idx.m2) = 0.001209;      % [-]
    u0(idx.m3) = 0.1190;        % [-]
    u0(idx.m4) = 0.01284;       % [-]
    u0(idx.m5) = 0.000869;      % [-]
    u0(idx.m6) = 0.001213;      % [-]
    u0(idx.m7) = 0.1191;        % [-]
    u0(idx.m8) = 0.004962;      % [-]
    u0(idx.h1) = 0.9718;        % [-]
    u0(idx.h2) = 0.1214;        % [-]
    u0(idx.h3) = 0.9718;        % [-]
    u0(idx.h4) = 0.9899;        % [-]
    u0(idx.h5) = 0.1210;        % [-]
    u0(idx.h6) = 0.9961;        % [-]
    u0(idx.Ca_n) = 0.1;         % [uM]
    u0(idx.nNOS_act_n) = 0.318; % [uM]
    u0(idx.NO_n) = 0.1671;      % [uM]
end

