%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Neural Network Simulator
% Using models from Izhikevich
% http://nsi.edu/users/izhikevich/publications/spikes.htm
% Bruce Land -- BioNB 222, March 2008
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clf
clear all

% Number of neurons
nNeuron = 3 ;
% Number of Sensory cells
% In this version they must be equal
nSense = nNeuron;

%timing
dt = 1/16 ; %millisecond -- time step
tEnd = 200; %maximum simulation time
%build the time vector
time = 0:dt:tEnd ;

%cell parameters from Izhikevich web page. See:
% http://nsi.edu/users/izhikevich/publications/spikes.htm
%   a           b       c       d       current
pars=[0.02      0.2     -65      6      14 ;...    % tonic spiking
    0.02      0.25    -65      6       0.5 ;...   % phasic spiking
    0.02      0.2     -50      2       15 ;...    % tonic bursting
    0.02      0.25    -55     0.05     0.6 ;...   % phasic bursting
    0.02      0.2     -55     4        10 ;...    % mixed mode
    0.01      0.2     -65     8        30 ;...    % spike frequency adaptation
    0.02      -0.1    -55     6        0  ;...    % Class 1
    0.2       0.26    -65     0        0  ;...    % Class 2
    0.02      0.2     -65     6        7  ;...    % spike latency
    0.05      0.26    -60     0        0  ;...    % subthreshold oscillations
    0.1       0.26    -60     -1       0  ;...    % resonator
    0.02      -0.1    -55     6        0  ;...    % integrator
    0.03      0.25    -60     4        0;...      % rebound spike
    0.03      0.25    -52     0        0;...      % rebound burst
    0.03      0.25    -60     4        0  ;...    % threshold variability
    1         1.5     -60     0      -65  ;...    % bistability
    1       0.2     -60     -21      0  ;...    % DAP
    0.02      1       -55     4        0  ;...    % accomodation
    -0.02      -1      -60     8        80 ;...    % inhibition-induced spiking
    -0.026     -1      -45     0        80];       % inhibition-induced bursting

%extract the cell type for each neuron
nType = [8 8 2] ;
a=pars(nType,1)';
b=pars(nType,2)';
c=pars(nType,3)';
d=pars(nType,4)';
Is = pars(nType,5)'; %steady current

%synaptic connections
%The weight matrix -- (+)excitatory: (-)inhibitory
%each column represents the input weights to one neuron
%e.g. element(1,2) is the weight connecting:
% (neuron 1 output spike) to (neuron 2 input synapse)
SynStrength = ...
10 * [0   -.95   3 ;...
      -1   0    3  ;...
       0   0   0  ] ;

%each neuron also potentially has several sensory inputs
%each column represents the input weights to one neuron
%from a sensory unit.
%e.g. element(1,2) is the weight connecting:
% (sensory unit 1 output current) to (neuron 2 input)
SenseStrength = ...
     [1   0   0 ;...
      0   1   0  ;...
      0   0   0  ] ;
  
%Sense inputs are time vectors on the inputs
SenseInput(:,1) = 10 * (time>50 & time<100);
SenseInput(:,2) = 15 * (time>75 & time<100);
SenseInput(:,3) = 0 * time ;

%synaptic current rises instantly, then decays exponentially
%set the decay rate of synaptic current to about 3 mSec
SynRate = .3*ones(1,nNeuron);
SynDecay = 1 - SynRate*dt;

%init state variables
v = c; % reset voltage;
s = zeros(1,nNeuron); % spike occurance
u = zeros(1,nNeuron); % revocery variable

%initial input currents
Istate = 0;  %state variable used for synaptic currents decay
I = Is ;

%step the time
tindex = 1; %pointer for output arrays

%init plotting variables
Vout = zeros(length(time),nNeuron) ;
Aout = zeros(length(time),nNeuron) ;

for t=time

    %add values to plotting variables
    Vout(tindex,:) = v ;
    Iout(tindex,:) = I ;
    
    %input current to cells
    %decaying synaptic currents
    Istate =  s*SynStrength + Istate .* SynDecay ; 
    I = Istate +  Is + SenseInput(tindex,:)*SenseStrength ;

    %update state variables
    v = v + dt * (0.04*v.^2+5*v+140-u+I);
    u = u + dt * a .* (b.*v-u);
    
    s = zeros(1,nNeuron) ;
    fired=find(v>=30); % indices of cells that spiked
    if ~isempty(fired)
        v(fired) = c(fired) ;
        u(fired) = u(fired)+d(fired) ;
        s(fired) = 1 ; %used in the next time step for synaptic input
    end;
   
    %update time index
    tindex = tindex+1;

end %main time step for-loop

figure(1)
clf
np=nNeuron ;
for i=1:nNeuron
    subplot(np,1,i)
    plot(time,Vout(:,i))
    hold on
    plot(time,Iout(:,i),'r')
    plot(time, SenseInput(:,i)*SenseStrength(i,:), 'g')
    ylabel(num2str(i))
    set(gca,'ylim', [-100 50])
end