function backward_shift(dt)
close all

cells = 12;
ee = ones(cells-1, 1);
connect = diag(ee, 1) + diag(ee, -1);
connect(1, 12) = 1;
connect(12, 1) = 1;
figure(1)
spy(connect)
ylabel('post')
xlabel('pre')
drawnow
%no plasticity modeled
%input: excitatory cells in MEC with independent Poisson spike trains
%input given to excitatory cells within a network of CA3 E and I cells

%stdp4 Input
% T: total time of simulation (sec)
% dt: timestep (ms)

laptime = 120; %T is ms
laps = 20;

simtime = laps*laptime;


Ne_M = 270;  %number of MEC excitatory cells

%Generate input spike trains
rate_M = 10;    %rate of MEC input spikes


%parameters
tauM_E = 20;    %exc. membrane voltage decay constant (ms)

Vrest_E = -70;  %exc. resting potential (ms)

Vth_E = -54;    %exc. threshold potential (mV)

Vreset_E = -60; %exc. reset potential (mV)


tau_gE = 5;

wMax_net = 1;
W = .1*wMax_net*connect;
Winp = 2;

%Initiate vars
v_E = Vrest_E*ones(cells,1);  %voltage for CA3 E cells

gI_E = 0*v_E;  %inhibitory conductance of CA3 E cells
% gI_I = 0; %no I-I connections
gE_E = gI_E;   %excitatory conductance of CA3 E cells


j=0;

%store excitatory MEC spikes
%assumption: no more than one spike per MEC neuron per timestep
% ct = 1;
% ind = 0;
% spiked1_M = zeros(ceil(1.5*rate_M*dt*Ne_M/1000),1);
% while Tinj(ct,1) < j*dt
%     ind = ind + 1;
%     spikedNeur = Tinj(ct,2);
%     spiked1_M(ind) = spikedNeur;
%     ct = ct + 1;
% end
% spiked_M = spiked1_M(1:ind);

%Parameters for plasticity from CA3 E to CA3 E cells
tauPlus_net = 10;
tauMin_net = 10;
APlus_net = wMax_net/10;
AMin_net = APlus_net*1.05;

%Variables for plasticity 
M_net = zeros(cells,1);  %related to timings of post spikes (CA3-CA3)
P_net = zeros(1,cells);    %related to timings of pre spikes (CA3-CA3)

%constants to be used in while loop
gEbot = tau_gE + dt;
quot_E = tauM_E/dt;

Pbot_net = tauPlus_net + dt;
Mbot_net = tauMin_net + dt;

%Assume all CA3 cells are initially at rest
spikedEx_C = [];

trefe = 5;  % refractory period (ms)
tref  = trefe*ones(cells,1);
refracCells = zeros(cells,1);   % initialize refractory cell vector
Wstore = spalloc(cells*cells, simtime, 24*simtime);
vstore = zeros(cells, simtime/dt);
gstore = vstore;
mstore = vstore;
pstore = vstore;


NumSpikedCA3 = 0;
figure(8)
fprintf('entering loop\n')
delta_t = laptime/12;
outputStep = ceil(simtime/100);
while j*dt < simtime
    %conductances decay
%     if mod(j*dt,outputStep)==0
%         fprintf('%g%% complete\n', (j*dt/simtime)*100)
%     end
    gEtop = tau_gE*gE_E;
    gE = gEtop/gEbot;
    gE_E = gE(1:cells);
    
    mstore(:, j+1) = M_net;
    pstore(:, j+1) = P_net;

       %P_net and M_net decay
    Ptop_net = tauPlus_net*P_net;
    P_net = Ptop_net/Pbot_net;
    
    Mtop_net = tauMin_net*M_net;
    M_net = Mtop_net/Mbot_net;    
    %adding mec input
    tdiff = mod(j*dt, laptime);
    if tdiff == 0
        pause
    end
    if mod(j*dt, delta_t)==0
        
        inp = tdiff/delta_t;
        if inp==0
            inp=12;
        end
        gE_E(inp) = gE_E(inp) + Winp;
    end
    vstore(:, j+1) = v_E;
    gstore(:, j+1) = gE_E;
    
    %Update voltages
    vtop_E = quot_E*v_E + Vrest_E; % +gEx*Eex
    vbot_E = quot_E + 1 + gE_E;
    v = vtop_E./vbot_E;
    
    refracIndex = find(refracCells > j*dt); % get cells that are still refractory
    v(refracIndex) = v_E(refracIndex);    % hold refractory cells at rest

    %find CA3 spiked cells
    spikedEx_C = find(v >= Vth_E);
    gE_E = gE_E + sum(W(:,spikedEx_C),2);

    if ~isempty(spikedEx_C)
        v(spikedEx_C) = Vreset_E;
        refracCells(spikedEx_C) = j*dt + tref(spikedEx_C);    % set new refractory times
        NumSpikedCA3 = NumSpikedCA3 + length(spikedEx_C);
        NumSpikedEx = length(spikedEx_C);
        
        %net post to net
        M_net(spikedEx_C) = M_net(spikedEx_C) - AMin_net;        

        Pbig_net = repmat(P_net*wMax_net,NumSpikedEx,1); 


        
        W_hot = W(spikedEx_C,:);
                
        W_hot = min(W_hot + Pbig_net, wMax_net);
        W(spikedEx_C,:) = W_hot;
        
        %net pre to net
        P_net(spikedEx_C) = P_net(spikedEx_C) + APlus_net;
         
        Mbig_net = repmat(M_net*wMax_net,1,NumSpikedEx);        
        
        W_hot = W(:,spikedEx_C);
        W_hot = max(W_hot + Mbig_net, 0);
        W(:,spikedEx_C) = W_hot;
        
        W(spikedEx_C,spikedEx_C) = W(spikedEx_C,spikedEx_C) - spdiags(diag(W(spikedEx_C,spikedEx_C)),0,NumSpikedEx,NumSpikedEx);   %no cell creates a synapse onto itself

    
    end

    v_E = v;
    %Set up next step
    j = j+1;


        Wlong = reshape(W', cells*cells, 1);
        Wstore(:, ceil(j*dt)) = Wlong; 
%         if ~isempty(spikedEx_C)
%             figure(2)
%             ff = ones(length(spikedEx_C), 1);
%             plot(ceil(j*dt*ff), spikedEx_C, '.r')
%             hold on
%             drawnow
%         end      
        
%         t = linspace(0.001, 2*pi, 100);
    if ~isempty(spikedEx_C)
        
        pos = (j*dt)/laptime;
        lap = ceil(pos);
        if j*dt <= laptime
            polar(pos*2*pi, laps,'.r')
            hold on
        else
            polar(pos*2*pi, laps-lap,'.b')
        end
        drawnow
    end
end


%plot voltages
figure(4)
    for k=1:cells
        subplot(3, 4,k)
        plot(dt*(1:(simtime)),vstore(k,:),'b')
        title(['v cell ' num2str(k)])
    end
    figure(5)
    for k=1:cells
        subplot(3, 4,k)
        plot(dt*(1:(simtime)),gstore(k,:),'b')
        title(['gE cell ' num2str(k)])

    end
    figure(6)
    
    for k=1:cells
        subplot(3,4,k)
        plot(dt*(1:simtime), mstore(k,:), 'b')
        title(['M_net ' num2str(k)])
    end
        figure(7)
    for k=1:cells
        subplot(3,4,k)
        plot(dt*(1:simtime), pstore(k,:), 'b')
        title(['P_net ' num2str(k)])
    end

    [row, col, a] = find(connect == 1);
    
    n = length(a);
    
    if ceil(sqrt(n))*floor(sqrt(n)) >= n
        subdimx = floor(sqrt(n));
        subdimy = ceil(sqrt(n));
    else
        subdimx = ceil(sqrt(n));
        subdimy = ceil(sqrt(n));
    end

    figure(3)
    for k=1:n

        subplot(subdimx,subdimy,k)
        plot(dt*(1:simtime),Wstore((row(k)-1)*12+col(k),:))
        title(['CA3  ' num2str(col(k)) ' to CA3 ' num2str(row(k))])
    end

return