TorusNeuronMod.m

classdef TorusNeuronMod < handle % Modify class state by reference
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%                            NEURON PROPERTIES                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	properties(GetAccess = public, SetAccess = public)
	% % Neuron Biophysical Parameters % %
        % Parameters from Khosravi-Hashemi, 2012 and 2011 
            % Includes leak, delayed K-rectifier, and Na
            % Experimental verification from Chacron and Fortune, 2010 %%
            Cm = 1.0;   % microFarads/cm^2
            gNa = 30;   % microS
            gK = 10;    % microS
            gL = 0.18;  % microS
            ENa = 60;   % mV 
            EK = -85;   % mV
            EL = -65;   % mV (leak)
        % Describes Ih current from HCN channels %
            % Based on Migliore, 2012; Poolos, 2002; Neymotin, 2013 %
            % Note: we ignore the t & v independent Ilk current as of now
            f_h = 1;    % scaling factor to vary Ih strength (Neymotin goes between 0 and 2)
            Eh = -30;   % mV
            g_h = 7;    % maximal conductance density (This value from Migliore 2012 Table 1)
        % Calcium Channel (T-type) Parameters %%
            f_ca = 1; 
            ECa = 120;
            gT_Ca = 0.32; % microS
            tau_hCa = 30;
        % Bias Current Level %
            I_bias = -1.3; % -1.3, Khosravi-Hashemi, 2012; -4.3, Khosravi-Hashemi, 2011
        % Stochastic Gaussian White Noise Current Parameters %
            % See Khosravi-Hashemi, 2011; Guo, 2011
            gaussMean = 0;      % Mean value for stochastic process
            gaussVariance = 0.8;  % Variance of process (from McGillivray, 2012)
            N_xi = 1;   % Scaling factor for process (Gaussian Intensity)
        % Synaptic properties
            Ws = 1;     % The overall synaptic weight (scales Isyn after summing)
            sigmaB = 0.5; % Balance factor between E and I cell input strength (1=all E)
    end
    
    properties(GetAccess = public, SetAccess = private)
        alphaSynapses = []; % List of Input alpha synapses
        vi;                 % mV (resting). Always reset.
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%                    PUBLIC NON-STATIC METHODS                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	methods(Static = false, Access = public)
        
    % MAIN VOLTAGE DERIVATIVE FUNCTION %
        function dvdt = voltageDerivative(this,v,n,h,hca,t)
            
            % Depolarizing Sodium Current %
            am = TorusNeuronMod.alphaM(v);
            minf =  am/(am + TorusNeuronMod.betaM(v));
            I_Na = this.gNa*(minf^3)*(0.85-n)*(v-this.ENa);
            
            % Delayed rectifier Potassium Current %
            I_K = this.gK*(n^4)*(v-this.EK);
            
            % HCN Channel-mediated Hyperpolarization Activated Current %
            I_H = (this.g_h*this.f_h*h) * (v-this.Eh);
            
            % Leak Current %
            I_L = this.gL*(v-this.EL);
            
            % T-type Calcium Current %
            sinf = 1/(1+exp(-1*(v+69)/7.8));
            I_Ca = this.gT_Ca*(sinf^3)*hca*(v-this.ECa)*this.f_ca;
            
            % Synaptic Input Current %
            I_syns = zeros(length(this.alphaSynapses),1); 
            as = this.alphaSynapses;
            for i = 1:length(as)
                if strcmpi(as(i).type,'E')
                    I_syns(i) = (2*this.sigmaB) * as(i).I_syn(v,t);
                else
                    I_syns(i) = (2*(1-this.sigmaB)) * as(i).I_syn(v,t);
                end
            end 
            I_syn = this.Ws * sum(I_syns);
            
            % Sum Currents (adds negative sign here)
            I_total = -I_Na - I_K - I_H - I_L - I_Ca - I_syn + this.I_bias;
            
            % Get final dv/dt
            dvdt = I_total/this.Cm;
        end
        
    % % INITIATION HELPER % %
        % Returns the initial (Steady-state) vector describing the system %%
        function iv = getInitialVector(this) 
            this.resetRestingVi();
            v = this.vi;
            n = this.nInf(v);
            h = this.hInf(v);
            hCa = this.hInf_Ca(v);
            iv = [v,n,h,hCa]';
        end
        
    % % HELPER FUNCTIONS FOR SYNAPTIC INPUT % %
        function addAlphaSynapse(this,alphasyn)
            if ~isa(alphasyn,'AlphaSynapse'); 
                error('Adding non-alpha-synapse to neuron!');
            end
            this.alphaSynapses = [this.alphaSynapses alphasyn];
        end
        function as = getAlphaSynapse(this,num)
            as = this.alphaSynapses(num);
        end    
        function pruneSynapses(this)
            this.alphaSynapses = [];
        end
        
	% % CONSTRUCTOR % %
		function t = TorusNeuronMod(varargin) %Permits 0-arg constructor calls
           if ~any([nargin==0,nargin==3])
               error('Improper argument number to neuron constructor')
           end
           if nargin == 3 
                cm = varargin{1};Gs=varargin{2};Es=varargin{3};
                if length(Gs)~=3 || length(Es)~=3
                    error('G and E vector length too short')
                end
                t.Cm = cm; 
                t.gNa = Gs(1); t.gK = Gs(2); t.gL = Gs(3);
                t.ENa = Es(1); t.EK = Es(2); t.EL = Es(3);
                t.resetRestingVi();
           end
        end
        
    % % DEEP COPY METHOD % %
        function newNeuron = deepCopy(this)
            warning off MATLAB:structOnObject % Suppress warnings
            if ~isa(this,'TorusNeuronMod')
                error('Unexpected Variable Type');
            end
            newNeuron = feval(class(this)); % Instantiate class
            p = fieldnames(struct(this));   % Procure all field names (incl. private)
            for i = 1:length(p)
                newNeuron.(p{i}) = this.(p{i}); % Set new neuron properties equivalently
            end
            for j = 1:length(this.alphaSynapses)
                newNeuron.alphaSynapses(j) = this.alphaSynapses(j).deepCopy();
            end
            warning on MATLAB:structOnObject % Turn warnings back on
        end
    end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%                    PRIVATE NON-STATIC METHODS                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    methods(Static = false, Access = private)
  
        % Reset estimated proper initial membrane voltage by running short sim %
        % No noise and no synaptic input. Basic Euler method to solve the
        % DE.
        function resetRestingVi(this) 
            % Ignore synaptic conductances. 
            ss = 0.025; % in ms
            endTime = 150; % in ms
            ts = ss:ss:endTime;
            N = length(ts);
            v0 = -55;          % Initial guess for v
            y0 = [v0;this.nInf(v0);this.hInf(v0);this.hInf_Ca(v0)];
            y = [ y0 zeros(4,N-1) ];
            for i = 2 : N
               % t = ts(i);
                v = y(1,i-1); %vec(1);
                n = y(2,i-1); %vec(2);
                h = y(3,i-1); %vec(3);
                hca = y(4,i-1); %vec(4);
                %vp = [0;0;0;0];
                %%%% Get dv/dt %%%%
                % Depolarizing Sodium Current %
                amh = 0.1*(v+40.7);
                am = (amh)/(1-exp(-amh));
                bm = 4*exp(-0.05*(v+49.7));
                minf =  am/(am + bm);
                I_Na = this.gNa*(minf^3)*(0.85-n)*(v-this.ENa);
                % Delayed rectifier Potassium Current %
                I_K = this.gK*(n^4)*(v-this.EK);
                % HCN Channel-mediated Hyperpolarization Activated Current %
                I_H = (this.g_h*this.f_h*h) * (v-this.Eh);
                % Leak Current %
                I_L = this.gL*(v-this.EL);
                % T-type Calcium Current %
                sinf = 1/(1+exp(-1*(v+69)/7.8));
                I_Ca = this.gT_Ca*(sinf^3)*hca*(v-this.ECa)*this.f_ca;
                % Synaptic Input Current %
                %ind = round(t/ss);
                %sb = this.sigmaB;
                I_syn = 0; %this.Ws * (2*(v-esAmpa))*((asynValsE(ind)*(sb))+(asynValsI(ind)*(1-sb)));
                % Sum Currents (adds negative sign here)
                I_total = -I_Na - I_K - I_H - I_L - I_Ca - I_syn + this.I_bias;
                % Get final dv/dt
                dvdt = I_total/this.Cm;
                %vp(1) = dvdt;
                %%%% End get dv/dt %%%%

                %%%% Get N' %%%%
                an = (0.01*(v+40.7))/(1-exp(-0.1*(v+40.7)));
                bn = 0.125*exp(-0.0125*(v+50.7));
                ninf = an/(an+bn);
                taun = 0.05/(an+bn);
                ndot = (ninf-n)/taun;
                %vp(2) = ndot;
                %%%% End Get N' %%%%

                %%%% Get h' %%%%
                tauh = exp(0.033*(v+75))/( 0.011*( 1 + exp(0.083*(v+75)) ) );
                v_hm = -73; % half maximal voltage in time constrant of h
                hinf = 1/(1+exp(0.151*(v-v_hm)));
                hdot = (hinf-h)/tauh;
                %vp(3) = hdot;
                %%%% End Get h' %%%%

                %%%% Get T' (Ca) %%%%
                q = sqrt(0.25+exp( (v+82)/6.3 ) );
                hinfca = 1/(0.5+q);
                hCadot = 2*(hinfca - hca)/this.tau_hCa;
                %vp(4) = hCadot;
                %%%% End Get T' (Ca) %%%%

                %%%% Combine the vectors %%%%
                vp = [dvdt; ndot; hdot; hCadot];
                
                %%%%%%%%%%%%%%%%%%%%%% Take the step %%%%%%%%%%%%%%%%%%%%%%
                y(:,i) = y(:,i-1) + (ss * vp);
            end
            this.vi = y(1,end);
        end
        
    end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%                    PUBLIC STATIC METHODS                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
	methods (Static = true, Access = public)
        
    % PRIMARY VECTOR OF DERIVATIVES FOR ODE SOLVER TO USE %     
    function vp = vprime(t,vec,this)
        % Preparation %
        %[v,n,h,hca]=deal(vec(1),vec(2),vec(3),vec(4));
        v = vec(1);
        n = vec(2);
        h = vec(3);
        hca = vec(4);
        vp = zeros(4,1);
        
        %%%% Get dv/dt %%%%
        % Depolarizing Sodium Current %
        amh = 0.1*(v+40.7);
        am = (amh)/(1-exp(-amh));
        bm = 4*exp(-0.05*(v+49.7));
        minf =  am/(am + bm);
        I_Na = this.gNa*(minf^3)*(0.85-n)*(v-this.ENa);
        % Delayed rectifier Potassium Current %
        I_K = this.gK*(n^4)*(v-this.EK);
        % HCN Channel-mediated Hyperpolarization Activated Current %
        I_H = (this.g_h*this.f_h*h) * (v-this.Eh);
        % Leak Current %
        I_L = this.gL*(v-this.EL);
        % T-type Calcium Current %
        sinf = 1/(1+exp(-1*(v+69)/7.8));
        I_Ca = this.gT_Ca*(sinf^3)*hca*(v-this.ECa)*this.f_ca;
        % Synaptic Input Current %
       % as = this.alphaSynapses;
        % Get Isyns. ASSUMES ONLY TWO SYNAPSES! AND ASSUMES BOTH ARE SAME AMPAR! %
        Is = 0; %#ok<NASGU>
        if strcmpi(this.alphaSynapses(1).type,'E')
            ind = round(t/this.alphaSynapses(1).ss2);
            sb = this.sigmaB;
            esAmpa = this.alphaSynapses(1).E_syn;
            Is = (2*(v-esAmpa))*(this.alphaSynapses(1).preCalcInterpVals(ind)*(sb)...
               + this.alphaSynapses(2).preCalcInterpVals(ind)*(1-sb));
            %Is = (this.alphaSynapses(1).preCalcInterpVals(ind)*(v-esAmpa)*(2*sb)...
            %    + this.alphaSynapses(2).preCalcInterpVals(ind)*(v-esAmpa)*(2*(1-sb)));
            %I_syns(1) = (2*this.sigmaB) * as(1).I_syn(v,t);
            %I_syns(2) = (2*(1-this.sigmaB)) * as(2).I_syn(v,t);
        else
            ind = round(t/this.alphaSynapses(1).ss2);
            sb = this.sigmaB;
            esAmpa = this.alphaSynapses(1).E_syn;
            Is = this.alphaSynapses(2).preCalcInterpVals(ind)*(v-esAmpa)*(2*sb)...
                + this.alphaSynapses(1).preCalcInterpVals(ind)*(v-esAmpa)*(2*(1-sb));
        end
        I_syn = this.Ws * Is;
        % Sum Currents (adds negative sign here)
        I_total = -I_Na - I_K - I_H - I_L - I_Ca - I_syn + this.I_bias;
        % Get final dv/dt
        dvdt = I_total/this.Cm;
        vp(1) = dvdt;
        %%%% End get dv/dt %%%%
        
        %%%% Get N' %%%%
        an = (0.01*(v+40.7))/(1-exp(-0.1*(v+40.7)));
        bn = 0.125*exp(-0.0125*(v+50.7));
        ninf = an/(an+bn);
        taun = 0.05/(an+bn);
        ndot = (ninf-n)/taun;
        vp(2) = ndot;
        %%%% End Get N' %%%%
    
        %%%% Get h' %%%%
        tauh = exp(0.033*(v+75))/( 0.011*( 1 + exp(0.083*(v+75)) ) );
        v_hm = -73; % half maximal voltage in time constrant of h
        hinf = 1/(1+exp(0.151*(v-v_hm)));
        hdot = (hinf-h)/tauh;
        vp(3) = hdot;
        %%%% End Get h' %%%%
    
        %%%% Get T' (Ca) %%%%
        q = sqrt(0.25+exp( (v+82)/6.3 ) );
        hinfca = 1/(0.5+q);
        hCadot = 2*(hinfca - hca)/this.tau_hCa;
        vp(4) = hCadot;
        %%%% End Get T' (Ca) %%%%
    end
    % END LOW OVERHEAD VECTOR DERIVATIVE EVALUATOR %
    
    % (Older Version) %
    function vp = vprimeOld(t,vec,cell)
        [v,n,h,hca]=deal(vec(1),vec(2),vec(3),vec(4));
        vp = zeros(4,1);
        % Compute vector of derivatives
        vp(1) = cell.voltageDerivative(v,n,h,hca,t);
        vp(2) = cell.nPrime(v,n);
        vp(3) = cell.hPrime(v,h);
        vp(4) = cell.hCaPrime(v,hca,cell);
    end

	% SECONDARY DIFFERENTIAL FUNCTIONS (GATING AND ACTIVATION VARIABLES) %
    % From Khosravi-Hashemi, 2012 and 2011 %%
        function ndot = nPrime(v,n)
            ndot = (TorusNeuronMod.nInf(v)-n)/TorusNeuronMod.tau_n(v);
        end
        function ninf = nInf(v)
            an = TorusNeuronMod.alphaN(v);
            ninf = an/(an+TorusNeuronMod.betaN(v));
        end
        function taun = tau_n(v)
            taun = 0.05/(TorusNeuronMod.alphaN(v)+TorusNeuronMod.betaN(v));
        end
        function alphan = alphaN(v)
            alphan = (0.01*(v+40.7))/(1-exp(-0.1*(v+40.7)));
        end
        function betan = betaN(v)
            betan = 0.125*exp(-0.0125*(v+50.7));
        end
        function minf = mInf(v)
            am = TorusNeuronMod.alphaM(v);
            minf = am/(am+TorusNeuronMod.betaM(v));
        end
        function alpham = alphaM(v)
            alpham = (0.1*(v+40.7))/(1-exp(-0.1*(v+40.7)));
        end
        function betam = betaM(v)
            betam = 4*exp(-0.05*(v+49.7));
        end
    % Largely from Neymotin 2013 (See also Poolos 2002, Migliore 2012) %%
        function hdot = hPrime(v,h) 
            hdot = (TorusNeuronMod.hInf(v)-h)/TorusNeuronMod.tau_h(v);
        end
        function hinf = hInf(v)
            v_hm = -73; % half maximal voltage in time constrant of h
            hinf = 1/(1+exp(0.151*(v-v_hm)));
        end
        function tauh = tau_h(v)
           tauh = exp(0.033*(v+75))/( 0.011*( 1 + exp(0.083*(v+75)) ) ); 
        end
    % Calcium Channel Equations %%
        function hCadot = hCaPrime(v,hCa,neuronInstance)
            hCadot = 2*(TorusNeuronMod.hInf_Ca(v) - hCa)/neuronInstance.tau_hCa;
        end
        function sinf = sInf(v)
            sinf = 1/(1+exp(-1*(v+69)/7.8));
        end
        function hinfca = hInf_Ca(v)
            q = sqrt(0.25+exp( (v+82)/6.3 ) );
            hinfca = 1/(0.5+q);
        end
    end
    %#ok<*DSPS>
end