% discriminaion.m
% file for making memory based decisions
% cells 1 and 2 get opposite inputs
% cells 3 and 4 are integrators with input 
% from cells 1 and 2 respectively
% cells 5 and 6 are the decision-making cells 
% with cross-inhibition and input from cells 1 and 2.
% decision-making cells must be suppressed until the 2nd stimulus

Ncells = 6;
dt = 0.0005;                    % time step
delay = 3.0;                    % time between 2 stimuli
duration = 0.5;                 % duration of each stimulus
tmax = 1.0 + delay +2.0*duration;             % maximum time
t = 0:dt:tmax;                  % set up vector of time
r=zeros(length(t),Ncells);        % set up vector of rates, initially zero
rinf = zeros(1,Ncells);
taufast = 0.010;                % fast AMPA time constant for the network
tauslow = 0.050;                % slow NMDA time constant
tau = taufast*ones(1,Ncells);   % by default set all cells to have taufast
%tau(3) = tauslow;               % give memory cells slower timeconstant
%tau(4) = tauslow;                   

W = zeros(Ncells,Ncells);       
W(3,3) = 1.0;                   % recurrent excitatory feedback strength
W(4,4) = 1.0;                   % is 1 for the integrators (memory cells)
W(1,3) = 0.2;                   % feedforward to integrator
W(2,4) = 0.2;                   % feedforward to integrator
W(3,1) = -0.2;                  % Feedback inhibition
W(4,2) = -0.2;                  % Feedback inhibition
 
W(1,5) = 0.2;                   % to decision-making circuit
W(2,6) = 0.2;                   % to decision-making circuit

W(5,5) = 1.1;                  % decision-making cells have 
W(6,6) = 1.1;                  % strong self-excitation
W(5,6) = -1.2;                  % and stronger
W(6,5) = -1.2;                  % cross-inhibition

Ith =zeros(1,Ncells);             % threshold current for firing

Iapp1 = 10;                      % value of applied current
Ion1 = 1.0;                       % time to start applied current
Ioff1 = Ion1+duration;                    % time to stop applied current
Iapp2 = 14;
Ion2 = Ioff1+delay;                     % time to start 2nd applied current
Ioff2 = Ion2+duration;                      % time to stop 2ndapplied current
Imax = 20;

sigma = 0.2;                    % level of noise in current

Iapp = zeros(length(t),Ncells);          % initialize vector of applied currents as zero
for i = 1:length(t)                 % run through all the elements
    if t(i) > Ion1 && t(i) <= Ioff1   % if the time is between the start as stop times for applied current
        Iapp(i,1) = Iapp1;            % make the applied current value non-zero
        Iapp(i,2) = Imax-Iapp1;
    end
    if t(i) > Ion2 && t(i) <= Ioff2   % if the time is between the start as stop times for applied current
        Iapp(i,1) = Iapp2;            % make the applied current value non-zero
        Iapp(i,2) = Imax-Iapp2;
    end
    if t(i) < Ion2
        Iapp(i,5) = -100;              % ensure decision-making cells 
        Iapp(i,6) = -100;              % do not fire before 2nd stimulus
    end
end

for i = 2:length(t)                             % run through time steps for integration
    Itot = Iapp(i,:) + r(i-1,:)*W + ...             % total current is applied current + feedback
        sigma*randn(1,Ncells)/sqrt(dt);                % + noise  

        
    rinf = max(Itot - Ith,0.0);    % steady state firing rate depends on current and must be at least zero             
    rinf = min(rinf,100);           % set maximum firing rate
    
    r(i,:) = rinf + (r(i-1,:)-rinf).*exp(-dt./tau);   % integrate the firing rate towards stedy state
end

plot(t,r)                        % plot firing rate against time
xlabel('time,t')                   
ylabel('rate,r')