%
% Tsonet.m 
%
% A simple 5 cell network (4 pryamidal(ex), 1 theta(in)) of
% iaf cells based on Tsodyks et al. 1996 Hippocampus 6:271-280
%
% 
%
%  example: Tsonet   
%

T = 100; %ms total time
dt = 0.01;
Nex = 4; %Number of excitatory cells
Nin = 1; %Number of inhibitory cells - change code above 1?

Vreset = 0.85;
Vth = 1;
tau = 20; %ms
tauex = 6;
tauin = 4;
J1 = 0.015; % ex syn scale
J2 = 0.02;  % in syn scale
pex = 0.2;  % prob of ex syn firing
pin = 0.7;  % prob of in syn firing
el = 0.15;  % extent of spatial connections
x = [0.1;0.2;0.3;0.4]; % location of place cell sensitivity - meters?
delt = 0.1; % dirac(t-t1) = 1 if (t in tim +/- delt)
I0 = 0.2;   % max ext stim
lame = 0.03; % amp of spatial comp of ext stim
lami = 0.02; % amp of Theta wave comp of ext stim
Th = 180; %ms time between Theta wave peeks (100-200)

y = 0.85*ones(3*Nex+Nin,1);%5 v's, 4 Iex's, 4 Iin's

sf = 0.8; %sigma factor added to 1 (1.8) above diagonal

% Syn weight factor weighted toward the forward cells
sig = ones(Nex)+triu(sf*ones(Nex))-diag((1+sf)*ones(1,Nex));

% Syn weight matrix (ex)
for i = 1:Nex
    for j = 1:Nex
        Jex(i,j) = sig(i,j)*J1*exp(-abs(x(i)-x(j))/el);
    end
end

% Syn weight matrix (in)
Jin = J2*ones(Nex,1);

%Diagonals for update matrix
d0v = [-1/tau*ones(1,5),-1/(tau*tauex)*ones(1,4),-1/(tau*tauin)*ones(1,4)];
d5v = [ones(1,4),zeros(1,4)];
d9v = -ones(1,4);

A = diag(d0v,0)+diag(d5v,5)+ diag(d9v,9);

B = inv(eye(3*Nex+Nin)-dt*A);

t = 0;
j = 1;

% Initiation of spike timer
tim = -100*ones(Nex+Nin,1);

while t < T

    j = j + 1;
    t = (j-1)*dt;

    % Delta functions
    diracex = ((t>(tim(1:Nex)-delt)).*(t<(tim(1:Nex)+delt)));
    diracin = ((t>(tim(Nex+1:Nex+Nin)-delt)).*(t<(tim(Nex+1:Nex+Nin)+delt)));

    % RHS of Iex and Iin
    fex = (rand(Nex,1)>pex).*(Jex*diracex);
    fin = (rand(Nex,1)>pin).*(Jin*diracin);
    
    x0 = t/4000; %Rat position - linear in time
    
    % External stimuli
    Iext1 = I0*(1+lame*exp(-abs(x-x0)));
    Iext2 = I0*(1+lami*cos(2*pi*t/Th))*ones(Nin,1);
    
    % Update via Backward Euler
    y = B*(y+dt*[Iext1;Iext2;fex;fin]);
    
    v(j-1,:) = y(1:5);
    tvec(j-1) = t;

    sind = find(y(1:5)>Vth);	% indices of spiking cells

    tim(sind) = t;	% spike times of spiking cells
    ras(j-1,:) = tim;
    
    y(sind) = Vreset; % reset spiked cells 
    
%    figure(1)
%     plot(t*ones(size(sind)),sind,'*')	% raster plot of spike times
%    hold on

%    figure(2)
%    plot(t,v,'*')
%    hold on
end

figure
plot(tvec,v)
xlabel('time (ms)')
ylabel('Potential')

figure

for i = 1:Nex+Nin
  spk=find(diff(ras(:,i))>0);
  plot(tvec(spk),i,'b*')
  hold on
end

axis([0 T 0 Nex++Nin+1])
xlabel('t (ms)','fontsize',16)
ylabel('cell','fontsize',16)
head = ['dt =',num2str(dt) ];
title(head,'fontsize',16)

hold off