echo off;
%----------------------------------------------------------------------
% Physiology 317 - Methods in Computational Neurobiology
%
% hw4p3.m - Homework #4, Problem #3
%
% Two neuron central pattern generator based on
% reciprocal inhibition and adaptive threshold
%
% UNITS: R in MOhm, C in nF, T in msec, V in mV, I in nA
%
% M. Nelson, 2004 
%----------------------------------------------------------------------
clear;                          % clear MATLAB variables
figure(1); clf;                 % clear MATLAB graphcis
orient portrait;                   

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

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

VTHR = 5;  % mV
VSPK = 70; % mV
EINH = -15; % mV

TAU_SYN = 15; % msec
GPEAK = 0.1; % uS
GNORM = GPEAK ./(TAU_SYN*exp(-1));

TAU_THRESH = 50; %msec

TMAX = 1500; %msec   
DT = 1;  % msec                  
t=(0:DT:TMAX)';                

% pre-allocate state variables
nNeurons = 2;
v = zeros(length(t), nNeurons);         
thresh = zeros(length(t), nNeurons);
z = zeros(length(t), nNeurons);
g = zeros(length(t), nNeurons); 
Isyn = zeros(length(t), nNeurons);

% define network connections
preSyn(1) = 2;
preSyn(2) = 1;

% define injection current
Inject(1) = 1.1;
Inject(2) = 0.9;

%%%%%%%%%%%%%%%
% MAIN LOOP
%%%%%%%%%%%%%%%
for i=2:length(t);                 
    for n = 1:nNeurons
        
        % update the synaptic conductance
        u = (v(i-1,preSyn(n)) == VSPK);
        z(i,n) = z(i-1,n) + ...
            (-z(i-1,n)/TAU_SYN + GNORM*u)*DT;
        g(i,n) = g(i-1, n) + ...
            (-g(i-1,n)/TAU_SYN + z(i-1,n))*DT;
        
        % update synaptic current
        Isyn(i,n) = g(i-1,n)*(v(i-1,n) - EINH);
        
        % update threshold
        dthresh_dt = (-thresh(i-1,n) + v(i-1,n))/TAU_THRESH;
        thresh(i,n) = thresh(i-1,n) + dthresh_dt*DT;

        % update membrane voltage
        if(v(i-1, n) == VSPK)                  
            v(i, n) = EINH; 
        else                           
            dvdt = (-v(i-1, n)/R + Inject(n) - Isyn(i,n))/C; 
            v(i,n) = v(i-1,n) + dvdt*DT; 
        end
       
        % check for threshold crossing
        if(v(i,n) >= thresh(i,n))             
            v(i,n) = VSPK; 
        end
        
    end  
end   

%%%%%%%%%%%%%%%
% PLOT RESULTS
%%%%%%%%%%%%%%%
plot(t,v(:,1),'b-'); hold on;
plot(t,thresh(:,1),'r-');
plot(t,100+v(:,2),'g-'); hold on;
plot(t,100+thresh(:,2),'r-');
title('Homework 4 Problem 3');
xlabel('time (msec)');
grid on;