echo off;
%----------------------------------------------------------------------
% Physiology 317 - Methods in Computational Neurobiology
%
% hw2p1.m - Homework #2, Problem #1
%
% Integrate & Fire neuron simulation -- STEP STIMULUS
% UNITS: R in MOhm, C in nF, T in msec, V in mV, I in nA
%
% M. Nelson, 1995, 1997, 2002, 2004 
%----------------------------------------------------------------------
clear;                          % clear MATLAB variables
figure(1); clf;                 % clear MATLAB graphcis
orient tall;                    % use the full page

%%%%%%%%%%%%%%%%%%
% INITIALIZATION
%%%%%%%%%%%%%%%%%%

DO_EXTRA_CREDIT = 1;            % extra-credit control flag

C = 1;   % nF
R = 10; % MOhm

VTHR = 5;  % mV
VSPK = 70; % mV

TMAX = 100; %msec   
DT = 1;                         % integration time step
t=(0:DT:TMAX)';                 % time vector
STIM_ON = 10;                   % stimulus onset
STIM_OFF = 60;                  % stimulus offset
ISTEP = 1.0;                    % current injection magnitude
I = ISTEP*(t>STIM_ON & t<STIM_OFF);      % current injection vector 

v = zeros(length(t),1);         % pre-allocate voltage array

%%%%%%%%%%%%%%%
% MAIN LOOP
%%%%%%%%%%%%%%%
for i=2:length(t);                 
    
    if(v(i-1) == VSPK)                  % if we just spiked, 
        v(i) = 0;                       %   reset the voltage to 0
    else                                % otherwise,
        dvdt = (-v(i-1)/R + I(i))/C;    %   calculate dv/dt using
        v(i) = v(i-1) + dvdt*DT;        %   Euler integration
    end
    
    if(v(i) >= VTHR)                    % check for threshold crossing
        v(i) = VSPK;                    % generate a spike
    end
    
end   

%%%%%%%%%%%%%%%
% PLOT RESULTS
%%%%%%%%%%%%%%%
subplot(2,2,1);                     % upper left corner of page
plot(t,I,'r-');                     % current vs. time
xlabel('time (msec)');              % x-axis label
ylabel('current (nA)');             % y-axis label
axis([0 TMAX -2 10]);
title('STIMULUS');
grid on;

subplot(2,2,2);                     % upper right corner of page
plot(t,v,'b-');                     % voltage vs. time
xlabel('time (msec)');
ylabel('voltage (mV)');
axis([0 TMAX -20 100]);            
title('RESPONSE');
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now plot firing frequency vs. injected current
% increase TMAX to 1000 msec, so we get more spikes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TMAX = 1000; %msec   
DT = 1;                                 % integration time step
t=(0:DT:TMAX)';                         % time vector
v = zeros(size(t));                     % pre-allocate

ISTEP = [0:.1:1.8 2:0.5:10];            % current injection for each trial

for itrial=1:length(ISTEP)
    
    nspikes(itrial) = 0;                % spike counter for this trial
    I = ISTEP(itrial)*ones(size(t));    % current injection vector
    
    %%%%%%%%%%%%%%%
    % MAIN LOOP
    %%%%%%%%%%%%%%%
    v(1) = 0;
    for i=2:length(t);
        
        if(v(i-1) == VSPK)                  % if we just spiked, 
            v(i) = 0;                       %   reset the voltage to 0
        else                                % otherwise,
            dvdt = (-v(i-1)/R + I(i-1))/C;  %   calculate dv/dt using
            v(i) = v(i-1) + dvdt*DT;        %   Euler integration
        end
        
        if(v(i) >= VTHR)         % check for threshold crossing
            v(i) = VSPK;
            nspikes(itrial) = nspikes(itrial) + 1;
        end
    end   
end

subplot(212)
stim_duration = TMAX*.001;   % stim duration in SECONDS 
firing_rate = nspikes/stim_duration;
plot(ISTEP, firing_rate, 'b.-');
xlabel('injected current (nA)')
ylabel('firing frequency (Hz)')
axis([min(ISTEP) max(ISTEP) 0 600]);
title('FIRING FREQUENCY vs. INJECTED CURRENT');
grid on;

% PUT BANNER TITLE AT TOP OF PAGE
axes('position',[.05 .05 0.9 0.9]); axis off;
title('HW 2.1: INTEGRATE-AND-FIRE NEURON MODEL -- STEP STIMULUS',...
    'fontw','bold');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EXTRA CREDIT:
% add Gaussian noise to the injection current;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(DO_EXTRA_CREDIT)
    
    ISTEP = [-3:.1:1.8 2:0.5:10];           % current injection for each trial
    
    for itrial=1:length(ISTEP)
        
        nspikes(itrial) = 0;                % spike counter for this trial
        I = ISTEP(itrial)*ones(size(t));    % current injection vector
        Inoise = 2*randn(size(t));          % random noise
        I = I + Inoise;
        
        v(1) = 0;
        for i=2:length(t);                      % MAIN LOOP
            
            if(v(i-1) == VSPK)                  % if we just spiked, 
                v(i) = 0;                       %   reset the voltage to 0
            else                                % otherwise,
                dvdt = (-v(i-1)/R + I(i-1))/C;  %   calculate dv/dt using
                v(i) = v(i-1) + dvdt*DT;        %   Euler integration
            end
            
            if(v(i) >= VTHR)                    % check for threshold crossing
                v(i) = VSPK;
                nspikes(itrial) = nspikes(itrial) + 1;
            end
        end   
    end
    
    subplot(212)
    hold on;
    stim_duration = TMAX*.001;   % stim duration in SECONDS 
    firing_rate = nspikes/stim_duration;
    plot(ISTEP, firing_rate, 'm.:');
    axis([min(ISTEP) max(ISTEP) 0 600]);
    legend('noise-free', '2 nA noise', 4);
end
