% Brunel (2000) Dynamics of Sparsely Connected Networks of Excitatory and Inhibitory Spiking Neurons

order = 2500;     % scales size of network (total 5*order neurons) 

N  = 5 * order;
Ne = (8 * N) / 10;  % number of excitatory neurons 
Ni = (2 * N) / 10;  % number of inhibitory neurons 

epsilon = 0.1;     % connectivity probability 

Ce   = epsilon * Ne; % total number of excitatory synapses 
Ci   = epsilon * Ni; % total number of inhibitory synapses 
C    = Ce + Ci;      % total number of internal synapses 
Cext = Ce;           % total number of external synapses 

Np = Cext; % total number of external poisson sources

tau     = 20.0;  % neuron membrane time constant [ms] 
tau_rp  = 2.0;   % refractory time [ms] 
theta   = 20.0;  % spike threshold [mV] 
Vreset  = 0.0;   % reset potential [mV] 

delay = 1.5; % synaptic delay, all connections [ms] 
J     = 0.1; % synaptic weight [mV] 

% initial vector of weights
I0 = zeros(N,1);

% time steps
h = 0.1;

% delay expressed as % time steps 
delaysteps = round (delay/h);
% delay matrix
D  = repmat(I0, [1, delaysteps]);


% case B : fast oscillation, irregular activity 
g    =  6.0; % rel strength, inhibitory synapses 
eta  =  4.0; % nu_ext / nu_thr 


% case A : almost full synchronization 
g    =  3.0; % rel strength, inhibitory synapses 
eta  =  2.0; % nu_ext / nu_thr 


Je = J; % excitatory weights 
Ji = -g  * Je; % inhibitory weights 

Jext = Je; % external weights 
                   
nu_thresh = theta / (Ce * Je * tau); % threshold rate
nu_ext    = eta * nu_thresh; % external rate per synapse
p_rate    = 0.001 * nu_ext;  % Poisson spiking rate [Hz]

% setup connection matrix 
[ei, ej] = find(epsilon >= rand(N,Ne));
[ii, ij] = find(epsilon >= rand(N,Ni));
S=cat(2,
      Je*sparse(ei,ej,1.0,N,Ne), % excitatory weights
      Ji*sparse(ii,ij,1.0,N,Ni) % inhibitory weights
     )

Vs = zeros(N,1);    % Initial values of V
Rs = zeros(N,1);    % refractory neurons
tend = 1200;        % simulation of 1200 ms

firings=cell(tend/h);     % spike timings

for t=1:h:tend            % simulation loop

  t

  % indices of spikes
  fired  = find(Vs >= theta);
  if (size(fired))(1) > 0
    fired'
    i = round(t/h);
    firings{i} = [t, fired'];
  endif

  % reset the spikng neurons to reset potential
  Vs(fired) = Vreset;

  % vector of refractory times (0 if the neuron is not in refractory
  % mode, otherwise the time remaining of its refractory period)
  refractory     = find(Rs > 0);
  Rs(refractory) = Rs(refractory) - h;
  Rs(fired)      = tau_rp;
  nonrefractory  = find(Rs <= 0);
  Rs(nonrefractory) = 0;

  % Poisson input
  I = Jext * randp (1000*h*p_rate, N, Np); 

  % the current synaptic weight vector is the sum of column 1 of the
  % delay matrix, plus the external weights
  W = D(:,1) + sum(I,2);

  % shift the delay matrix and append to the end the weights of the
  % neurons that spiked in the current time step
  D = shift(D,-1,2);
  D(:,delaysteps) = sum(S(:,fired),2);

  % Voltage equation right hand side
  Vrhs = (-Vs (nonrefractory) ./ tau) .+ W (nonrefractory);

  % Euler integration of the voltage equation
  Vs(nonrefractory) = Vs(nonrefractory) .+ (h * Vrhs);

end;

save("-v7", "firings.dat", "firings");

% plot spike times
figure ();
hold on;
for i = 1:size(firings)(1)
  i
  data = firings{i};
  n = size(data);
  if (n(1) > 0)
    plot(data(1,1),data(1,2:n(2)),'b.');
  endif
end;