NSUtils.m

classdef NSUtils 
	methods (Static = true, Access = public)
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        % % % % % % % % % % SIMULATION UTILITIES % % % % % % % % %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        % Master calling function for calculating the CSI and VPDs of the
        % just-run simulation. Note that it should compute VPD_full.
        % The RMSE will also be computed by xcorr and circshift
        function stats = computePostSimulationStats...
                (psths,nChirps,ss,stsCa,q,numTrials,alpha,siglev)
            [allcsias,allcsifs,avgFRc] = NSUtils.computeAvgCSIS(psths,nChirps,ss,siglev); 
            avgvpd = NSUtils.computeTotalAvgVPD(stsCa,q,nChirps,numTrials);
            stats.csias = allcsias; % CSIavgs
            stats.csifs = allcsifs; % CSImaxfrs
            stats.vpd_mt = avgvpd;  % Average cross-trial VPD
            stats.avgFRc = avgFRc;  % Average FR during the window after the chirp
            % Compute SI score
            SI = mean(allcsifs) - (alpha * avgvpd);
            if SI < 0; SI = 0; end
            stats.SI = SI;
            % Store params
            stats.params.q = q;
            stats.params.alphaSI = alpha;
            % Compute pairwise RMSE average
            rmse = 0; counter = 0;
            for c1 = 1 : nChirps
                for c2 = (c1+1) : nChirps
                    p1 = psths{c1};
                    p2 = psths{c2};
                    rmse = NSUtils.calcOptRMSE(p1,p2);
                    counter = counter + 1;
                end
            end
            stats.rmse_avgpw = rmse / counter;
        end
        
        % Input: Two PSTHs
        % Output: the RMSE between them, after optimal circ-shifting
        function rmse = calcOptRMSE(psth1,psth2)
            N = length(psth1);
            % Compute cross-correlation
            [cc,lags] = xcorr(psth1,psth2);
            shift = lags(cc==max(cc));
            % Circular shift the data (psth2) to match
            newpsth2 = circshift(psth2,shift);
            % Compute the RMSE now
            rmse = sqrt( sum( (psth1(:)-newpsth2(:)).^2 ) / N );
        end   

        % Compute the cross-trial chirp selectivities.
        % Computes the CSIAs, which use the average FRs across the beat and
        % chirp, and the CSIFs, which use the max FRs across the beat and
        % chirp. 
        function [allcsias,allcsifs,cFRavg] = computeAvgCSIS(psths,numChirps,ss,siglev)
            n = numChirps;
            allcsias = zeros(n,1);
            allcsifs = zeros(n,1);
            cFRavg = 0;
            offset = 5; % Because inexactness of the spike time (i.e.
               % to before the chirp time), ignore offset ms before it,
               % else the timing detection may be off as well,
               % especially for the BC. Spike bleeding post-smoothing also
               % requires this.
            %%% Prepare input variables %%%
            % % % % Fixed Input variables % % % % 
            TimeWindow = 100;  % ms (for a 5hz stimului), after the chirp onset
            endTime = 1600;   % The time length of the stimulus in ms
            % % % % Compute necessary variables % % % %
            ChirpTimeOnset = endTime/2 - offset; % How many ms into the stim the chirp occurs
            % % % % Compute actual indices for chirp window % % % %
            t = ss:ss:endTime;
            ind1 = find(abs(t-ChirpTimeOnset)<0.0001);
            ind2 = ind1+(TimeWindow/ss); % Time of ind1 + window (in time points)
            for i = 1 : n
                % Get PSTH %
                psth = psths{i};
                
                %%% Update the FRavg during the post-chirp window %%%
                cFRavg = ( cFRavg*(i-1) + mean(psth(ind1:ind2)) )/i;
                
                %%% Compute individual CSIA %%%
                % % % % Get R_beat % % % %
                % The average FR of the neuron across the whole PSTH, except
                % for the window after the chirp onset (i.e. for TimeWindow ms
                % after ChirpTimeOnset) in two segments (pre/post chirp).
                % Using averages weighted by the number of points in the two
                % intervals, so as not to overwight mean2, albeit slightly.
                wmean1 = mean(psth(1:ind1-1))*length(psth(1:ind1-1));
                wmean2 = mean(psth(ind2+1:end))*length(psth(ind2+1:end));
                rb = (wmean1+wmean2)/(length(psth)-length(ind1:ind2));
                % % % % Get R_Chirp % % % %
                rc = mean(psth(ind1:ind2));
                % % % % Get CSI and Store It % % % %
                if ( (rc+rb)==0 )
                    csia = 0; % No firing at occurred, so there is no selectivity
                else
                    csia = ( (rc-rb) / (rc+rb) );
                end
                allcsias(i,1) = csia;
                
                %%% Compute individual CSIF %%%
                % % % % Get R_beat % % % %
                % The max FR of the neuron anywhere outside the TimeWindow,
                % looking both before and after the chirp:chirp+TW period 
                mfr1 = max(psth(1:ind1-1));
                mfr2 = max(psth(ind2+1:end));
                rb = max(mfr1,mfr2);
                % % % % Get R_Chirp % % % %
                rc = max(psth(ind1:ind2));
               % fprintf('rb=%f,rc=%f\n',rb,rc);
                % % % % Get CSI and Store It % % % %
                if ( (rc+rb)==0 )
                    csif = 0; % No firing at occurred, so there is no selectivity
                else
                    csif = ( (rc-rb) / (rc+rb) );
                end
                
                % % Check if the csi passes our significance requirements % %
                if rc < siglev; csif = 0; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                
                allcsifs(i,1) = csif;
            end
        end
        
        % Computes the all-trial average Victor-Purpura distance.
        %   In detail, it calculates:
        %     vpdavg = (1/n) sum_(forall chirps) sum_(forall trials) d( m(i,j),m(k,p) )
        %   where m(u,v) is the spike-time vector for chirp u & trial v, d
        %   is the VPD function, and n is the number of combinations checked.
        % In essence, we compute the VPD between all possible pairs of
        %   chirps and trials, sum them, and take the average.
        % In case the binaries hold responses to non-chirp
        %   input as well, the function asks for the number of chirps and
        %   assumes they are the first numChirps inputs (i.e. binaries).
        % Input:
        %   sts = spike times, as a |nchirps|x|ntrials| cell array (stsCa)
        %   q = the VPD q parameter
        %   numchirps = the total number of chirps to look across (in case
        %       the sts contains more than just the stims we want to analyze)
        %   numTrials = the number of trials for which each stim was run
        function avg = computeTotalAvgVPD(sts,q,numChirps,numTrials)
            n = numChirps; m = numTrials;
            sum = 0; counter = 0;
            for i = 1 : n % across chirps
                for j = 1 : m % across trials
                    sts1 = sts{i,j};
                    for k = i : n % across all other chirps
                        Lij = 1;
                        if i==k; Lij = j+1; end
                        for p = Lij : m
                            sts2 = sts{k,p};
                            d = NSUtils.spkd(sts1,sts2,q);
                            sum = sum + d;
                            counter = counter + 1;
                        end
                    end
                end
            end
            avg = sum / counter;
        end
        
        % Computes the total VPD (then averages it) but only considers
        % matched trials (i.e. computes d across all chirps c_i, but only
        % compares to the same-number trial, i.e:
        %   d_avg = (1/n) sum_{c in chirps} sum_{trial j} sum_{other
        %       chirps} d( M(ci,j), M(ci,j) )
        % Necessary when computational time is problematic.
        function avg = computeMatchedTrialAvgVPD(sts,cost,numChirps,numTrials)
            n = numChirps; m = numTrials;
            sum = 0; counter = 0;
            for i = 1 : n % across chirps
                for j = 1 : m % across trials
                    tli = sts{i,j};
                    for k = (i+1) : n % across all other chirps
                        tlj = sts{k,j};
                        %%% RUN VPD CALCULATOR %%%
                        %d = NSUtils.spkd(sts1,sts2,q);
                        nspi=length(tli);
                        nspj=length(tlj);
                        if cost==0
                           d=abs(nspi-nspj);
                           sum = sum + d;
                           counter = counter + 1;
                           continue;
                           %return
                        elseif cost==Inf
                           d=nspi+nspj;
                           sum = sum + d;
                           counter = counter + 1;
                           continue;
                           %return
                        end
                        scr=zeros(nspi+1,nspj+1);
                        %INITIALIZE MARGINS WITH COST OF ADDING A SPIKE
                        scr(:,1)=(0:nspi)';
                        scr(1,:)=(0:nspj);
                        if nspi && nspj 
                           for ii=2:nspi+1
                              for jj=2:nspj+1
                                 scr(ii,jj)=min([scr(ii-1,jj)+1 scr(ii,jj-1)+1 ...
                                     scr(ii-1,jj-1)+cost*abs(tli(ii-1)-tlj(jj-1))]);
                              end
                           end
                        end
                        d=scr(nspi+1,nspj+1); 
                        %%% END VPD CALCULATOR %%%
                        sum = sum + d;
                        counter = counter + 1;
                    end
                end
            end
            avg = sum / counter;            
        end
        
        %%%%%% Begin External Code %%%%%%
        % The Victor-Purpura Distance between a spike train pair
        % Procured on 5/26/14 at:
        % http://www-users.med.cornell.edu/~jdvicto/spkdm.html
        % Note: this is not the most computationally efficient method
        function d=spkd(tli,tlj,cost)
            % d=spkd(tli,tlj,cost) calculates the "spike time" distance
            % (Victor & Purpura 1996) for a single cost
            % tli: vector of spike times for first spike train
            % tlj: vector of spike times for second spike train
            % cost: cost per unit time to move a spike
            %  Copyright (c) 1999 by Daniel Reich and Jonathan Victor.
            %  Translated to Matlab by Daniel Reich from FORTRAN code by Jonathan Victor.
            nspi=length(tli);
            nspj=length(tlj);
            if cost==0
               d=abs(nspi-nspj);
               return
            elseif cost==Inf
               d=nspi+nspj;
               return
            end
            scr=zeros(nspi+1,nspj+1);
            %INITIALIZE MARGINS WITH COST OF ADDING A SPIKE
            scr(:,1)=(0:nspi)';
            scr(1,:)=(0:nspj);
            if nspi & nspj %#ok<AND2>
               for i=2:nspi+1
                  for j=2:nspj+1
                     scr(i,j)=min([scr(i-1,j)+1 scr(i,j-1)+1 ...
                         scr(i-1,j-1)+cost*abs(tli(i-1)-tlj(j-1))]);
                  end
               end
            end
            d=scr(nspi+1,nspj+1);        
        end
        %%%%%% End External Code %%%%%%        
        

    % Simple method to streamline access to the psth, spike times,
        % and binaries after running a simulation.
        function [psth,bins,sts]=getPsthStsAndBinaries(vals,thresh,ss)
            vals = vals';
            vs = vals(:,1);
            mt = length(vs)*ss;
            times = ss:ss:mt;
            bins = NSUtils.transformHHrespToBinaries(vs,thresh);
            psth = NSUtils.constructPsthFromBinaries(bins);
            sts = NSUtils.convBinariesToSpikeTimes(bins,times);
        end
        
    % Converts a series of binary spike times (implicit time) to actual 
    % spike times in a vector. 
    function sts = convBinariesToSpikeTimes(binaries,timeSeries)
        sts = timeSeries(binaries==1); % Logical vector indexing
    end

    % Applies a simple thresholding to find spikes
        % Returns a series of 1s and 0s, where each 1 represents a spike
        % at the time of its index. 
        function b = transformHHrespToBinaries(voltageTimeSeries,threshold)
            b = zeros(length(voltageTimeSeries),1);
            belowThresh = 1;
            for i = 1:length(voltageTimeSeries)
                v = voltageTimeSeries(i);
                if (v > threshold) & (belowThresh==1) %#ok<AND2>
                    belowThresh = 0;
                    b(i) = 1;
                elseif (v < threshold) & (belowThresh==0) %#ok<AND2>
                    belowThresh = 1;
                end
            end    
        end    
        
        % Computes the PSTH given by a binary array
        % Currently using fixed bin size. 
        % Based on Bretschneider 2006 description
        % Gives an estimated instantaneous firing rate per millisecond
        %   for every point in the binaries
        function psth = constructPsthFromBinaries(binaries)
            delta = 400; % 10 ms -> binaries are in 0.025ms steps -> 400 points
            ss = 0.025; % Step size in ms (note: points per ms = 10)
            N = ceil(length(binaries)/delta); % Number of bins
            spikesPerBin = zeros(N,1);
            i = 1; currBin = 1;
            while i < length(binaries)
                endNum = i + delta - 1;
                if ( (i + delta) > length(binaries) ) % Last bin may be partial delta
                    endNum = length(binaries);
                end
                numSpikes = sum(binaries(i:endNum));
                spikesPerBin(currBin) = numSpikes;
                i = i + delta; 
                currBin = currBin + 1;
            end
            x = (spikesPerBin / (delta*ss))*1000; % Normalize to spikes per s
            r = repmat(x,1,delta)';
            temp = r(:);
            psth = temp(1:length(binaries));
        end
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        % % % % % % % % GENERAL PLOTTING UTILITIES % % % % % % % % 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % Plots across multi-trial output of the 7 current stims:
        %   5sc, 1bc, 2n. 
        function plotMultiTrialOutput(cell,ss,synapses,stsCa,psths,graphNoise)
            figure; % Plot the chirps and the response
            numStims = 5;
            numPlts = 5;
            k = size(stsCa);
            numTrials = k(2); % Get the number of trials from the sts cell array matrix
            tmax = round(length(synapses(1).Esyn.conved)/10); % Assumes 1st one is chirp!          
            tsR = 0.1:0.1:tmax;
            tsS = ss:ss:tmax;
            mn = min(tsR);
            mx = max(tsR);
            ssr = tsR(2)-tsR(1);
            for i = 1:numStims
                % Plot chirp
                chirp = synapses(i).Esyn.getChirp();
                subplot(numPlts,numStims,i);
                plot(chirp); xlim([ssr,mx/ssr]);
                % Plot Convolved synaptic activity (unscaled)
                s1e = synapses(i).Esyn; 
                s1i = synapses(i).Isyn;
                eNo = s1e.conved; %.preCalcInterpVals; %s1e.gSynArray(tsR);
                iNo = s1i.conved; %.preCalcInterpVals; %s1i.gSynArray(tsR);
                subplot(numPlts,numStims,i+numStims);
                plot(eNo); xlim([mn,mx/ssr]);
                subplot(numPlts,numStims,i+2*numStims);
                plot(iNo); xlim([mn,mx/ssr]);
                % Plot TS neuron raster
                subplot(numPlts,numStims,i+3*numStims);
                for j = 1:numTrials
                    sts = stsCa{i,j};
                    hold on;
                    plot(sts,ones(length(sts),1)*j,'*');
                    xlim([mn/10,mx/ssr/10]); % Its in real time
                    ylim([0,numTrials+1]); % Little extra room
                end
                % Plot TS neuron psth
                subplot(numPlts,numStims,i+4*numStims);
                plot(psths{i}); xlim([min(tsS),max(tsS)/(tsS(2)-tsS(1))]);
            end
            if ~graphNoise; return; end % If GN is not requested, end now
            figure; % Plot the noise and the response
            numNoises = 2;
            numPlts = 4;
            counter = 1;
            for i = (numStims+1):length(synapses)
                % Get times
                tmax = round(length(synapses(1).Esyn.conved)/10);
                tsR = 0.1:0.1:tmax;
                tsS = ss:ss:tmax;  
                ss1 = tsR(2)-tsR(1);
                % Plot noise
                subplot(numPlts,numNoises,counter);
                chirp = synapses(i).Esyn.getChirp();
                plot(chirp); xlim([ss1,tmax]);
                % Plot synapses
                s1e = synapses(i).Esyn; 
                s1i = synapses(i).Isyn;
                eNo = s1e.gSynArray(tsR);
                iNo = s1i.gSynArray(tsR);
                subplot(numPlts,numNoises,counter+numNoises);
                plot(eNo); xlim([ss1,tmax]);
                subplot(numPlts,numNoises,counter+2*numNoises);
                plot(iNo); xlim([ss1,tmax]);
                % Plot TS neuron psth
                subplot(numPlts,numNoises,counter+3*numNoises);
                plot(psths{i}); xlim([min(tsS),max(tsS)/(tsS(2)-tsS(1))]);
                counter = counter + 1;
            end    
            %title
%             ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',[0 ...
%                 1],'Box','off','Visible','off','Units','normalized',...
%                 'clipping' , 'off');
            s = sprintf('f_h:%.2f,f_ca:%.2f,I_bias:%.2f,Ws:%.2f,sigmaB:%.2f,N_xi:%.2f',...
                cell.f_h,cell.f_ca,cell.I_bias,cell.Ws,cell.sigmaB,cell.N_xi);
            text(0.5, 1,sprintf('\bf %s',s),...
                'HorizontalAlignment','center','VerticalAlignment', 'top');
        end
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        % % % % % % % % % MATHEMATICAL UTILITIES % % % % % % % % %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        % Simple Numerical integration scheme for stochastic differential
        %   equations based on the Euler-Maruyama scheme with a Gaussian
        %   white noise term in the voltage.
        % General Refs: Maruyama 1955; Kloeden & Platen 1992
        % Neuroscience Methods Refs: Manwani & Koch 1999
        % Neuroscience Use Refs: McGillivray, 2012; Guo, 2011; 
        function [ts,y] = FastNeuronEulerMaruyama(neuron,ss,tmax) 
            
            % % % % % % % % % CHECK INPUT % % % % % % % % %
                if ~isa(neuron,'TorusNeuronMod')
                    error('Attempted EulerMaruyama without proper cell type.');
                end
                if ~( (ss > 0) && (ss < 1) && (length(ss)==1) )
                    error('Unacceptable Step Size');
                end
                if ~( (length(tmax)==1) && (tmax > 1) )
                    error('Unacceptable max time');
                end    
            % % % % % % % % % PARAMETERS % % % % % % % % %
                % Retrieve parameters from neuron
                gMu = neuron.gaussMean;
                gSigma = neuron.gaussVariance;
                D = neuron.N_xi;   
                y0 = neuron.getInitialVector();
                Cm = neuron.Cm;
                %dvecByDt = @(t,v) neuron.vprime(t,v,neuron);
                % Fixed parameters
                ts = ss:ss:tmax;     % times of evaluation
                
            % % % % % % % % % ALGORITHM % % % % % % % % %
                % Generate and Precalculate variables
                N = length(ts);
                xi = D * normrnd(gMu,gSigma,N,1);
                y = [ y0 zeros(4,N-1) ];
                sqrtdt = sqrt(ss);
                if strcmpi(neuron.alphaSynapses(1).type,'E')
                    asynValsE = neuron.alphaSynapses(1).preCalcInterpVals;
                    asynValsI = neuron.alphaSynapses(2).preCalcInterpVals;
                else
                    asynValsE = neuron.alphaSynapses(2).preCalcInterpVals;
                    asynValsI = neuron.alphaSynapses(1).preCalcInterpVals;                     
                end    
                esAmpa = neuron.alphaSynapses(1).E_syn;
                
                % Run algorithm
                for i = 2:N
                    normNoise = [sqrtdt * (xi(i)/Cm); 0; 0; 0]; 
                    %step = ss * neuron.vprime( ts(i), y(:,i-1), neuron );
                    
                    %%%%%%%%%%%%%% START: Compute Step %%%%%%%%%%%%%%
                    t = ts(i);
                    %vec = y(:,i-1);
                    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 = neuron.gNa*(minf^3)*(0.85-n)*(v-neuron.ENa);
                    % Delayed rectifier Potassium Current %
                    I_K = neuron.gK*(n^4)*(v-neuron.EK);
                    % HCN Channel-mediated Hyperpolarization Activated Current %
                    I_H = (neuron.g_h*neuron.f_h*h) * (v-neuron.Eh);
                    % Leak Current %
                    I_L = neuron.gL*(v-neuron.EL);
                    % T-type Calcium Current %
                    sinf = 1/(1+exp(-1*(v+69)/7.8));
                    I_Ca = neuron.gT_Ca*(sinf^3)*hca*(v-neuron.ECa)*neuron.f_ca;
                    % Synaptic Input Current %
                    ind = round(t/ss);
                    sb = neuron.sigmaB;
                    I_syn = neuron.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 + neuron.I_bias;
                    % Get final dv/dt
                    dvdt = I_total/neuron.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)/neuron.tau_hCa;
                    %vp(4) = hCadot;
                    %%%% End Get T' (Ca) %%%%
                    
                    %%%% Combine the vectors %%%%
                    vp = [dvdt; ndot; hdot; hCadot];
                    
                    %%%%%%%%%%%%%% END: Compute Step %%%%%%%%%%%%%%
                    step = ss * vp;
                    y(:,i) = y(:,i-1) + step + normNoise;
                end
        end
        
         % Simple Numerical integration scheme for stochastic differential
        %   equations based on the Euler-Maruyama scheme with a Gaussian
        %   white noise term in the voltage.
        % General Refs: Maruyama 1955; Kloeden & Platen 1992
        % Neuroscience Methods Refs: Manwani & Koch 1999
        % Neuroscience Use Refs: McGillivray, 2012; Guo, 2011; 
        function [ts,y,Is] = FastNeuronEulerMaruyamaRetIs(neuron,ss,tmax) 
            
            % % % % % % % % % CHECK INPUT % % % % % % % % %
                if ~isa(neuron,'TorusNeuronMod')
                    error('Attempted EulerMaruyama without proper cell type.');
                end
                if ~( (ss > 0) && (ss < 1) && (length(ss)==1) )
                    error('Unacceptable Step Size');
                end
                if ~( (length(tmax)==1) && (tmax > 1) )
                    error('Unacceptable max time');
                end    
            % % % % % % % % % PARAMETERS % % % % % % % % %
                % Retrieve parameters from neuron
                gMu = neuron.gaussMean;
                gSigma = neuron.gaussVariance;
                D = neuron.N_xi;   
                y0 = neuron.getInitialVector();
                Cm = neuron.Cm;

                % Fixed parameters
                ts = ss:ss:tmax;     % times of evaluation
                
            % % % % % % % % % ALGORITHM % % % % % % % % %
                % Generate and Precalculate variables
                N = length(ts);
                xi = D * normrnd(gMu,gSigma,N,1);
                y = [ y0 zeros(4,N-1) ];
                sqrtdt = sqrt(ss);
                if strcmpi(neuron.alphaSynapses(1).type,'E')
                    asynValsE = neuron.alphaSynapses(1).preCalcInterpVals;
                    asynValsI = neuron.alphaSynapses(2).preCalcInterpVals;
                else
                    asynValsE = neuron.alphaSynapses(2).preCalcInterpVals;
                    asynValsI = neuron.alphaSynapses(1).preCalcInterpVals;                     
                end    
                esAmpa = neuron.alphaSynapses(1).E_syn; % Nernst potential
                
                % Prepare to store copies of the currents as well
                Is.I_syn = zeros(1,N);
                Is.I_Na = zeros(1,N);
                Is.I_K = zeros(1,N);
                Is.I_H = zeros(1,N);
                Is.I_L = zeros(1,N);
                Is.I_Ca = zeros(1,N);
                
                % Run algorithm
                for i = 2:N
                    normNoise = [sqrtdt * (xi(i)/Cm); 0; 0; 0]; 
                    
                    %%%%%%%%%%%%%% START: Compute Step %%%%%%%%%%%%%%
                    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);

                    %%%% 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 = neuron.gNa*(minf^3)*(0.85-n)*(v-neuron.ENa);
                    % Delayed rectifier Potassium Current %
                    I_K = neuron.gK*(n^4)*(v-neuron.EK);
                    % HCN Channel-mediated Hyperpolarization Activated Current %
                    I_H = (neuron.g_h*neuron.f_h*h) * (v-neuron.Eh);
                    % Leak Current %
                    I_L = neuron.gL*(v-neuron.EL);
                    % T-type Calcium Current %
                    sinf = 1/(1+exp(-1*(v+69)/7.8));
                    I_Ca = neuron.gT_Ca*(sinf^3)*hca*(v-neuron.ECa)*neuron.f_ca;
                    % Synaptic Input Current %
                    ind = round(t/ss);
                    sb = neuron.sigmaB;
                    I_syn = neuron.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 + neuron.I_bias;
                    % Get final dv/dt
                    dvdt = I_total/neuron.Cm;
                    %%%% 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;
                    %%%% 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 constant of h
                    hinf = 1/(1+exp(0.151*(v-v_hm)));
                    hdot = (hinf-h)/tauh;
                    %%%% End Get h' %%%%

                    %%%% Get T' (Ca) %%%%
                    q = sqrt(0.25+exp( (v+82)/6.3 ) );
                    hinfca = 1/(0.5+q);
                    hCadot = 2*(hinfca - hca)/neuron.tau_hCa;
                    %%%% End Get T' (Ca) %%%%
                    
                    %%%% Store the currents %%%%
                    Is.I_syn(1,i) = I_syn;
                    Is.I_Na(1,i) = I_Na;
                    Is.I_K(1,i) = I_K;
                    Is.I_H(1,i) = I_H;
                    Is.I_L(1,i) = I_L;
                    Is.I_Ca(1,i) = I_Ca;
                    
                    %%%% Combine the vectors %%%%
                    vp = [dvdt; ndot; hdot; hCadot];
                    
                    %%%%%%%%%%%%%% END: Compute Step %%%%%%%%%%%%%%
                    
                    step = ss * vp;
                    y(:,i) = y(:,i-1) + step + normNoise;
                end
        end
        
        % Simple Numerical integration scheme for stochastic differential
        %   equations based on the Euler-Maruyama scheme with a Gaussian
        %   white noise term in the voltage.
        % General Refs: Maruyama 1955; Kloeden & Platen 1992
        % Neuroscience Methods Refs: Manwani & Koch 1999
        % Neuroscience Use Refs: McGillivray, 2012; Guo, 2011; 
        function [ts,y,gs] = FastNeuronEulerMaruyamaRetGs(neuron,ss,tmax) 
            
            % % % % % % % % % CHECK INPUT % % % % % % % % %
                if ~isa(neuron,'TorusNeuronMod')
                    error('Attempted EulerMaruyama without proper cell type.');
                end
                if ~( (ss > 0) && (ss < 1) && (length(ss)==1) )
                    error('Unacceptable Step Size');
                end
                if ~( (length(tmax)==1) && (tmax > 1) )
                    error('Unacceptable max time');
                end    
            % % % % % % % % % PARAMETERS % % % % % % % % %
                % Retrieve parameters from neuron
                gMu = neuron.gaussMean;
                gSigma = neuron.gaussVariance;
                D = neuron.N_xi;   
                y0 = neuron.getInitialVector();
                Cm = neuron.Cm;

                % Fixed parameters
                ts = ss:ss:tmax;     % times of evaluation
                
            % % % % % % % % % ALGORITHM % % % % % % % % %
                % Generate and Precalculate variables
                N = length(ts);
                xi = D * normrnd(gMu,gSigma,N,1);
                y = [ y0 zeros(4,N-1) ];
                sqrtdt = sqrt(ss);
                if strcmpi(neuron.alphaSynapses(1).type,'E')
                    asynValsE = neuron.alphaSynapses(1).preCalcInterpVals;
                    asynValsI = neuron.alphaSynapses(2).preCalcInterpVals;
                else
                    asynValsE = neuron.alphaSynapses(2).preCalcInterpVals;
                    asynValsI = neuron.alphaSynapses(1).preCalcInterpVals;                     
                end    
                esAmpa = neuron.alphaSynapses(1).E_syn; % Nernst potential
                
                % Prepare to store copies of the currents as well
                gs.g_syn = zeros(1,N);
                gs.g_Na = zeros(1,N);
                gs.g_K = zeros(1,N);
                gs.g_H = zeros(1,N);
                gs.g_L = zeros(1,N);
                gs.g_Ca = zeros(1,N);
                
                % Run algorithm
                for i = 2:N
                    normNoise = [sqrtdt * (xi(i)/Cm); 0; 0; 0]; 
                    
                    %%%%%%%%%%%%%% START: Compute Step %%%%%%%%%%%%%%
                    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);

                    %%%% 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);
                    GNA = neuron.gNa*(minf^3)*(0.85-n);
                    I_Na = GNA*(v-neuron.ENa);
                    % Delayed rectifier Potassium Current %
                    GK = neuron.gK*(n^4);
                    I_K = GK*(v-neuron.EK);
                    % HCN Channel-mediated Hyperpolarization Activated Current %
                    GH = (neuron.g_h*neuron.f_h*h);
                    I_H = GH * (v-neuron.Eh);
                    % Leak Current %
                    I_L = neuron.gL*(v-neuron.EL);
                    % T-type Calcium Current %
                    sinf = 1/(1+exp(-1*(v+69)/7.8));
                    GTCA = neuron.gT_Ca*(sinf^3)*hca*neuron.f_ca;
                    I_Ca = GTCA*(v-neuron.ECa);
                    % Synaptic Input Current %
                    ind = round(t/ss);
                    sb = neuron.sigmaB;
                    GSYN = 2 * neuron.Ws * ( (asynValsE(ind)*(sb)) + (asynValsI(ind)*(1-sb)) );
                    %I_syn = neuron.Ws * (2*(v-esAmpa))*((asynValsE(ind)*(sb))+(asynValsI(ind)*(1-sb)));
                    I_syn = (v-esAmpa) * GSYN;
                    % Sum Currents (adds negative sign here)
                    I_total = -I_Na - I_K - I_H - I_L - I_Ca - I_syn + neuron.I_bias;
                    % Get final dv/dt
                    dvdt = I_total/neuron.Cm;
                    %%%% 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;
                    %%%% 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 constant of h
                    hinf = 1/(1+exp(0.151*(v-v_hm)));
                    hdot = (hinf-h)/tauh;
                    %%%% End Get h' %%%%

                    %%%% Get T' (Ca) %%%%
                    q = sqrt(0.25+exp( (v+82)/6.3 ) );
                    hinfca = 1/(0.5+q);
                    hCadot = 2*(hinfca - hca)/neuron.tau_hCa;
                    %%%% End Get T' (Ca) %%%%
                    
                    %%%% Store the conductances %%%%
                    gs.g_syn(1,i) = GSYN;
                    gs.g_Na(1,i) = GNA;
                    gs.g_K(1,i) = GK;
                    gs.g_H(1,i) = GH;
                    gs.g_Ca(1,i) = GTCA;
                    
                    %%%% Combine the vectors %%%%
                    vp = [dvdt; ndot; hdot; hCadot];
                    
                    %%%%%%%%%%%%%% END: Compute Step %%%%%%%%%%%%%%
                    
                    step = ss * vp;
                    y(:,i) = y(:,i-1) + step + normNoise;
                end
        end       
        
        % Takes a folder of data structures in any format
        % Returns a column vector of strings of target filenames
        function list = getFilesInDir(folder)
            list = {};
            pat = '.mat';
            fs = dir(folder);
            for i = 1 : length(fs)
                targ = fs(i);
                if targ.isdir==0
                    s = targ.name;
                    %  Determines whether string s ends with pattern pat
                    sl = length(s);
                    pl = length(pat);
                    b = (sl >= pl && strcmp(s(sl-pl+1:sl), pat)) || isempty(pat);
                    % End external code
                    if b==1
                        list{end+1} = [folder,'/',s]; %#ok<AGROW>
                    end
                end
            end
        end
        
        
    end
end