clear

dt=0.005;
tmax = 1.;
t=0:dt:tmax;                % Create time vector

Ncells = 30;                % Number of cells around the ring

r=zeros(length(t),Ncells);  % Rate matrix for all timesteps and all cells
h=zeros(1,Ncells);          % Applied thalamic input to all cells
Itot=zeros(1,Ncells);       % Total current to all cells


A = 40.;                    % Strength of input
epsilon = 0.1;              % Variation of input across cells
J0 = -2;                    % Lateral inhibition
J2 =3;                      % Local excitation
tau = 0.05;                 % Time constant for rate model

cuecell = Ncells/2;         % Cell with peak of the input

Ith = 10;                   % Threshold current needed for half-maximal firing
Iwidth = 10;                % How close to threshold to get change in firing rate
rmax = 100;                 % Maximum firing rate

hstart = 0.2;               % Time to begin input
hend = 0.5;                 % Time to end input

c = 0.5                     % Contrast affects level of input

for cell = 1:Ncells         % This loop sets input current that varies across cells
    h(cell) = A*c*(1-epsilon + epsilon*cos(2*pi*(cell-cuecell)/Ncells));
end

for i = 2:length(t)         % Negin time integration

    for cell = 1:Ncells     % Need to update every cell

        if ( t(i) >= hstart ) && (t(i) < hend )
            Itot(cell) = h(cell);       % Add input for a short time period
        else    
            Itot(cell) = 0.0;           % Otherwise no input current
        end

        for cell2 = 1:Ncells            % Calculate feedback to cell from all cells (cell2)
            Itot(cell) = Itot(cell) +  ... % Next line depends on the rate of another cell and the 
                r(i-1,cell2)*( J0+J2*cos(2*pi*(cell2-cell)/Ncells) ) / Ncells;  % difference in their 
                                                                                % preferred angle
                                                                                % divide by Ncells 
                                                                                % to average
        end
    end

    rinf = rmax./(1+exp(-(Itot-Ith)/Iwidth)).^2;    % rinf is a vector of the steady state firing rate 
                                                    % given the total
                                                    % current vector to
                                                    % each cell

    r(i,:) = rinf + (r(i-1,:)-rinf)*exp(-dt/tau);   % update r from r(i-1) for all cells


end

subplot(2,1,1)                  
plot(r(floor(hstart/dt),:))             % plot rate of all cells before input
hold on

plot(r(ceil(hend/dt),:),'r')            % plot rate of all cells at end of input

plot(r(length(t),:),'g')                % plot rate of all cells at end of simulation

subplot(2,1,2)
plot(t,r(:,cuecell))                    % plot rate of cuecell as a function of time
hold on
plot(t,r(:,mod(cuecell-1+round(Ncells/4),Ncells)+1),'r')    % plot rate of cell separated by 1/4 circle
plot(t,r(:,mod(cuecell-1-round(Ncells/4),Ncells)+1),'g')    % ahead then behind the cue cell

plotcell = cuecell+round(Ncells/2)                          

plot(t,r(:,mod(cuecell-1+round(Ncells/2),Ncells)+1),'c')    % plot rate of cell with opposite tuning to cuecell