echo off;
%----------------------------------------------------------------------
% Physiology 317 - Methods in Computational Neurobiology
%
% hw2p3.m - Homework #2, Problem #3
%
% Integrate & Fire neuron simulation -- SINUSOIDAL STIMULUS
% Spikes are blocked (simulated TTX)
%
% 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

%%%%%%%%%%%%%%%%
% Intialization
%%%%%%%%%%%%%%%%

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

VTHR = 999;                     % simulate TTX block by setting VTHR high
VSPK = 70; % mV

TMAX = 1000; %msec   
DT = 1;                                 % integration time step
t =( 0:DT:TMAX)';                       % time vector
v = zeros(size(t));                     % pre-allocate voltage vector

FREQ_LIST = [1, 2, 5, 10, 20, 50, 100]; % test frequencies


%%%%%%%%%%%%
% FREQ LOOP
%%%%%%%%%%%%

for itrial=1:length(FREQ_LIST)          % BEGIN FREQ LOOP
    
    nspikes(itrial) = 0;                % initialize count for this trial
    
    FREQ = FREQ_LIST(itrial);           % get stimulus frequency
    I = sin(2*pi*FREQ*t/1000);          % build current injection vector
    
    %%%%%%%%%%%%%
    % TIME LOOP
    %%%%%%%%%%%%%
    v(1) = 0;
    for i=2:length(t);                  % BEGIN TIME 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;                    % generate a spike
            nspikes(itrial) = nspikes(itrial) + 1;
        end
    end                                 % END TIME LOOP
    
    maxv(itrial) = max(v);              % save maximum voltage for this trial
    
                                        % Plot current vs. time
    subplot(2,2,1)
    offset = 5*(itrial-1);
    plot(t,I+offset,'r-'); hold on;
    
                                        % Plot voltage vs. time
    subplot(2,2,2);
    offset = 100*(itrial-1);
    plot(t,v+offset,'b-'); hold on;
    
end                                     % END FREQ LOOP

                                        % Set axis limits, labels, etc.
subplot(221)
xlabel('time (msec)');
ylabel('current (nA)');
axis([0 TMAX -5 35]);
title('STIMULUS');
grid on;

subplot(222)
xlabel('time (msec)');   
ylabel('voltage (mV)');   
axis([0 TMAX -100 700]);            
title('RESPONSE');
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT maxv versus test frequency
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(212)
semilogx(FREQ_LIST, maxv, 'b*-');
xlabel('stimulus frequency (Hz)')
ylabel('peak depolarization (mV)')
title('PEAK DEPOLARIZATION vs. STIMULUS FREQUENCY');
grid on;

% PUT BANNER TITLE AT TOP OF PAGE
axes('position',[.05 .05 0.9 0.9]); axis off;      
title('HW 2.3: I&F MODEL WITH SPIKES BLOCKED -- SINUSOIDAL STIMULUS',...
    'fontw','bold');
