%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 = 14 ;
% 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.5     -65      6      30 ;...    % tonic spiking PD cell
    0.02      0.5     -65      6       30 ;...   % tonic spiking PD cell
    0.02      0.2     -50      2       30 ;...    % tonic bursting AB Cell
    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      .8      -60     4        0;...      % rebound spike PY cells
    0.1       .8      -60     0        0;...      % rebound spike LP Cell
    0.1       .5      -60     0        0  ;...    % rebound spike VD cell
    0.03      1       -60     4        0  ;...    % rebound spike IC cell
    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 = [3 2 1 14 13 13 13 13 13 13 13 13 15 16] ;
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 0 0 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1.6 -1;...
      0 0 0 -1 -1 -1 -1 -1 -1 -1 -1 -1 0 -1;...
      0 0 0 -1 -1 -1 -1 -1 -1 -1 -1 -1 0 -1;...
      0 -1 -1 0 -1 -1 -1 -1 -1 -1 -1 -1 -1 0;...
      0 0 0 -1 0 0 0 0 0 0 0 0 0 -1;...
      0 0 0 -1 0 0 0 0 0 0 0 0 0 -1;...
      0 0 0 -1 0 0 0 0 0 0 0 0 0 -1;...
      0 0 0 -1 0 0 0 0 0 0 0 0 0 -1;...
      0 0 0 -1 0 0 0 0 0 0 0 0 0 -1;...
      0 0 0 -1 0 0 0 0 0 0 0 0 0 -1;...
      0 0 0 -1 0 0 0 0 0 0 0 0 0 -1;...
      0 0 0 -1 0 0 0 0 0 0 0 0 0 -1;...
      0 0 0 -1 -1 -1 -1 -1 -1 -1 -1 -1 0 -1.2; ...
      0 0 0 0 0 0 0 0 0 0 0 0 -1 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 0 0 0 0 0 0 0 0 0 0;...
      0 1 0 0 0 0 0 0 0 0 0 0 0 0;...
      0 0 1 0 0 0 0 0 0 0 0 0 0 0;...
      0 0 0 1 0 0 0 0 0 0 0 0 0 0;...
      0 0 0 0 1 0 0 0 0 0 0 0 0 0;...
      0 0 0 0 0 1 0 0 0 0 0 0 0 0;...
      0 0 0 0 0 0 1 0 0 0 0 0 0 0;...
      0 0 0 0 0 0 0 1 0 0 0 0 0 0;...
      0 0 0 0 0 0 0 0 1 0 0 0 0 0;...
      0 0 0 0 0 0 0 0 0 1 0 0 0 0;...
      0 0 0 0 0 0 0 0 0 0 1 0 0 0;...
      0 0 0 0 0 0 0 0 0 0 0 1 0 0;...
      0 0 0 0 0 0 0 0 0 0 0 0 1 0;...
      0 0 0 0 0 0 0 0 0 0 0 0 0 1;]  ;
  
%Sense inputs are time vectors on the inputs
SenseInput(:,1) = 0 * (time>50 & time<100);
SenseInput(:,2) = 0 * (time>50 & time<100);
SenseInput(:,3) = 0 * (time>50 & time<100);
SenseInput(:,4) = 0 * (time>50 & time<100);
SenseInput(:,5) = 0 * (time>50 & time<100);
SenseInput(:,6) = 0 * (time>50 & time<100);
SenseInput(:,7) = 0 * (time>50 & time<100);
SenseInput(:,8) = 0 * (time>50 & time<100);
SenseInput(:,9) = 0 * (time>50 & time<100);
SenseInput(:,10) = 0 * (time>50 & time<100);
SenseInput(:,11) = 0 * (time>50 & time<100);
SenseInput(:,12) = 0 * (time>50 & time<100);
SenseInput(:,13)= 0* (time>50 & time<100);
SenseInput(:,14) = 0* (time>50 & time<100);


%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
    Itrans1= .5.* (v(1)-v(2)); % AB-PD1 connection
    Itrans2= .5 .*(v(1)-v(3)); % AB-PD2 connection
    Itrans3= .4 .*(v(2)-v(3)); %PD1-PD2 connection
    
    % coupling parameter between PY's
    G= .5; 
    %PY1's electric connections to other PY's and the one LP 
    Itrans2ONE= G .*(v(5)-v(6));
    Itrans3ONE= G .* (v(5)-v(7));
    Itrans4ONE= G .* (v(5)- v(8));
    Itrans5ONE= G .* (v(5)- v(9));
    Itrans6ONE= G .* (v(5)- v(10));
    Itrans7ONE= G .* (v(5)- v(11));
    Itrans8ONE= G .* (v(5)- v(12));
    ItransLPONE= G .* (v(5)-v(4));
    
    %PY2's electric connections to other PY's and the one LP
    Itrans3TWO= G .* (v(6)- v(7));
    Itrans4TWO= G .* (v(6)- v(8));
    Itrans5TWO= G .* (v(6)- v(9));
    Itrans6TWO= G .* (v(6)- v(10));
    Itrans7TWO= G .* (v(6)- v(11));
    Itrans8TWO= G .* (v(6)- v(12));
    ItransLPTWO= G .* (v(6)-v(4));
    
    %PY3's electric connections to other PY's and the one LP
    Itrans4THREE= G .* (v(7)- v(8));
    Itrans5THREE= G .* (v(7)- v(9));
    Itrans6THREE= G .* (v(7)- v(10));
    Itrans7THREE= G .* (v(7)- v(11));
    Itrans8THREE= G .* (v(7)- v(12));
    ItransLPTHREE= G .* (v(7)- v(4));
    
    %PY4's electric connections to other PY's and the one LP
    Itrans5FOUR= G .* (v(8)- v(9));
    Itrans6FOUR= G .* (v(8)- v(10));
    Itrans7FOUR= G .* (v(8)- v(11));
    Itrans8FOUR= G .* (v(8)- v(12));
    ItransLPFOUR= G .* (v(8)- v(4));
    
    %PY5's electric connections to other PY's and the one LP
    Itrans6FIVE= G .* (v(9)- v(10));
    Itrans7FIVE= G .* (v(9)- v(11));
    Itrans8FIVE= G .* (v(9)- v(12));
    ItransLPFIVE= G .* (v(9)- v(4));
    
    %PY6's electric connections to other Py's and the one LP
    Itrans7SIX= G .* (v(10)- v(11));
    Itrans8SIX= G .* (v(10)- v(12));
    ItransLPSIX= G .* (v(10)- v(4));
    
    %PY7's electric connections to other PY's and the one LP
    Itrans8SEVEN= G .* (v(11)- v(12));
    ItransLPSEVEN= G .* (v(11)- v(4));
    
    %PY8's electric connection to LP
    ItransLPEIGHT= G .* (v(12)- v(4));   
    
    %VD's electric connection to PD and AB
    Itrans1THIRTEEN=.05.*(v(13)-v(1));
    Itrans2THIRTEEN=.05.*(v(13)-v(2));
    Itrans3THIRTEEN= .05.*(v(13)-v(3));
    
    %IC's electric connection to PD cells
    Itrans2FOURTEEN=.05.*(v(14)-v(3));
    Itrans3FOURTEEN=.05.*(v(14)-v(2));
    
    Istate =  s*SynStrength + Istate .* SynDecay ; 
    I = Istate +  Is + SenseInput(tindex,:)*SenseStrength ;

    %update state variables
    v(1) = v(1) + dt * (0.04*v(1).^2+5*v(1)+140-u(1)+I(1) - Itrans1-Itrans2 + Itrans1THIRTEEN);
    v(2) = v(2) + dt * (0.04*v(2).^2+5*v(2)+140-u(2)+I(2)+ Itrans1 -Itrans3 + Itrans2THIRTEEN +Itrans2FOURTEEN);
    v(3) = v(3) + dt * (0.04*v(3).^2+5*v(3)+140-u(3)+I(3) + Itrans2 +Itrans3 + Itrans3THIRTEEN+Itrans2FOURTEEN);
    v(4) = v(4) + dt * (0.04*v(4).^2+5*v(4)+140-u(4)+I(4)+ItransLPONE+ItransLPTWO+ItransLPTHREE+ItransLPFOUR+ItransLPFIVE+ItransLPSIX+ItransLPSEVEN+ItransLPEIGHT);
    v(5) = v(5) + dt * (0.04*v(5).^2+5*v(5)+140-u(5)+I(5)-Itrans2ONE-Itrans3ONE-Itrans4ONE-Itrans5ONE-Itrans6ONE-Itrans7ONE-Itrans8ONE-ItransLPONE);
    v(6) = v(6) + dt * (0.04*v(6).^2+5*v(6)+140-u(6)+I(6)+Itrans2ONE - Itrans3TWO-Itrans4TWO-Itrans5TWO-Itrans6TWO-Itrans7TWO-Itrans8TWO-ItransLPTWO);
    v(7) = v(7) + dt * (0.04*v(7).^2+5*v(7)+140-u(7)+I(7)+Itrans3TWO + Itrans3ONE - Itrans4THREE-Itrans5THREE-Itrans6THREE-Itrans7THREE-Itrans8THREE );
    v(8) = v(8) + dt * (0.04*v(8).^2+5*v(8)+140-u(8)+I(8)+Itrans4ONE+Itrans4TWO+Itrans4THREE- Itrans5FOUR-Itrans6FOUR-Itrans7FOUR-Itrans8FOUR-ItransLPFOUR);
    v(9) = v(9) + dt * (0.04*v(9).^2+5*v(9)+140-u(9)+I(9)+Itrans5ONE+Itrans5TWO+Itrans5THREE+Itrans5FOUR-Itrans6FIVE-Itrans7FIVE-Itrans8FIVE-ItransLPFIVE);
    v(10)= v(10) + dt * (0.04*v(10).^2+5*v(10)+140-u(10)+I(10)+Itrans6ONE+Itrans6TWO+Itrans6THREE+Itrans6FOUR+Itrans6FIVE-Itrans7SIX-Itrans8SIX-ItransLPSIX);
    v(11)= v(11) + dt * (0.04*v(11).^2+5*v(11)+140-u(11)+I(11)+Itrans7ONE+Itrans7TWO+Itrans7THREE+Itrans7FOUR+Itrans7FIVE+Itrans7SIX-Itrans8SEVEN-ItransLPSEVEN);
    v(12)= v(12) + dt * (0.04*v(12).^2+5*v(12)+140-u(12)+I(12)+Itrans8ONE+Itrans8TWO+Itrans8THREE+Itrans8FOUR+Itrans8FIVE+Itrans8SIX+Itrans8SEVEN-ItransLPEIGHT);
    v(13)= v(13) + dt * (0.04*v(13).^2+5*v(13)+140-u(13)+I(13)-Itrans1THIRTEEN-Itrans2THIRTEEN-Itrans3THIRTEEN);
    v(14)= v(14) + dt * (0.04*v(14).^2+5*v(14)+140-u(14)+I(14)-Itrans2FOURTEEN-Itrans3FOURTEEN);
    
    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:14
    if(i~=3 && i~=6 && i~=7 && i~=5 && i~=9 && i~=10 && i~=11 && i~=12)
   
    %plot(time,Iout(:,i),'r')
    %plot(time, SenseInput(:,i)*SenseStrength(i,:), 'g')
        if(i==1)
            subplot(6,1,1)
            plot(time,Vout(:,i))
            hold on
            ylabel('AB')
        elseif(i==2)
            subplot(6,1,2)
            plot(time,Vout(:,i))
            hold on
            ylabel('PD')
        elseif(i==4)
            subplot(6,1,5)
            plot(time,Vout(:,i))
            hold on
            ylabel('LP')
        elseif(i==8)
            subplot(6,1,6)
            plot(time,Vout(:,i))
            hold on
            ylabel('PY')
        elseif(i==13)
            subplot(6,1,3)
            plot(time,Vout(:,i))
            hold on
            ylabel('VD')
        elseif(i==14)
            subplot(6,1,4)
            plot(time,Vout(:,i))
            hold on
            ylabel('IC')                   
        end
        
    set(gca,'ylim', [-100 50])
    end
end
   