close all;
clear;

% The number of excitatory neurons in the network. The mammalian cortex has
% about 4 times as many excitatory neurons as inhibitory ones.
Ne=800;
%The number of inhibitory neurons in the network.
Ni=200;

%% Initializing the neuronal network 

%Random numbers
re=rand(Ne, 1);
ri=rand(Ni, 1);
%This will set the value of a for all excitatory neurons to 0.02 and the value of a
%for inhibitory neurons to a random number between 0.02 and 0.1
a=[0.02*ones(Ne, 1); 0.02+0.08*ri];
%This will allow b to range from 0.2 to –0.25
b=[0.2*ones(Ne, 1); 0.25-0.05*ri];
%This will allow the spike reset membrane potential to range between -65 and -50
c=[-65+15*re.^2; -65*ones(Ni,1)];
%This will allow the recovery reset value to range between 2 and 8
d=[8-6*re.^2; 2*ones(Ni, 1)];

%%  Adjust these parameters
% Connectivity matrix S = [S_exc, S_inh]
% make sure max synaptic weight for S_inh is negative
S_exc=.8;
S_inh=-.4;

%Create some random input external to the network with mean defined here
excitatory_mean=0;
inhibitory_mean=0;

% Uncomment the connetion you prefer: all-to-all or sparse
% All-to-all connection for neurons
S=[S_exc*rand(Ne+Ni, Ne), S_inh*rand(Ne+Ni, Ni)];

% Sparse connectivity of excitatory synapses, full connectivity of inhibitory
% synapses
 S=[S_exc*sprand(Ne+Ni, Ne, 0.5), S_inh*rand(Ne+Ni, Ni)];


%% Numerical code for solving differential equation begins

%The initial values for v and u
v=-65*ones(Ne+Ni,1);
u=b.*v;
%Firings will be a two-column matrix. The first column will indicate the time that a
%neuron’s membrane potential crossed 30, and the second column will be a number
%between 1 and Ne+Ni that identifies which neuron fired at that time.
firings=[];
for t=1:2000

    I=[excitatory_mean+4*randn(Ne, 1); inhibitory_mean+2*randn(Ni,1)];
    %Determine which neurons crossed threshold at the current time step t.
    fired=find(v>=30);
    %Add the times of firing and the neuron number to firings.
    if isempty(firings)
        firings=horzcat(t*ones(length(fired),1), fired);
    else
        firings=[firings; t*ones(length(fired),1), fired];
    end
    %Reset the neurons that fired to the spike reset membrane potential and
    %recovery variable.
    v(fired)=c(fired);
    u(fired)=u(fired)+d(fired);
    %Add to the input, I, for each neuron a value equal to the sum of the synaptic
    %strengths of all other neurons that fired in the last time step connected to that
    %neuron.
    I=I+sum(S(:,fired), 2);
    %Move the simulation forward using Euler’s method
    % Slightly different integrate a fire model...involves two variables 
    v=v+0.5*(0.04*v.^2+5*v+140-u+I);
    v=v+0.5*(0.04*v.^2+5*v+140-u+I);
    u=u+a.*(b.*v-u);
end;
%Plot the raster plot of the network activity.
plot(firings(:,1), firings(:,2),'.');
