% cerebellar.m
% model of cerebellar conditioning with a bunch of granule cells 
% excited by the CS, that excite one Purkinje cell which also 
% receives inputs from a climbing fibre is excited by the US.
% Coincident granule cell and climbing fibre firing leads to 
% LTD.

clear
dt = 0.005;
tmax = 1;
t = 0:dt:tmax;

Ntrials = 150;

NGcells = 40;               % No. of granule cells
NIcells = NGcells;          % No. of inhibitory (Golgi) cells

% Granule cell f-I curves:
IthGC = 10*ones(1,NGcells)+5*rand(1,NGcells);   % Randomness in thresholds of granule cells
IwidthGC = 10*ones(1,NGcells);                  % Width of f-I curve
rmaxGC = 100*ones(1,NGcells);                   % Max. firing rate
tauGC = 0.02;                                   % membrane time constant

taus = zeros(1,NGcells);                        % synaptic time constant
tauf = zeros(1,NGcells);                        % time constant for facilitation
taud = zeros(1,NGcells);                        % time constant for depression

taus0 = 0.075;                       % base value for synaptic time constant 
taud0 = 0.15;                                    % base value for depression time constant
tauf0 = 0.3;                                    % base value for facilitation time constant

for cell = 1:NGcells
    taus(cell) = taus0;                         % fixed synaptic time consants for all cells
    tauf(cell) = tauf0*(1.0+rand(1));       % heterogeneity in time constant for facilitation
    taud(cell) = taud0*(1.0+rand(1));       % heterogeneity in time constant for deression
end

% Inhibitory (Golgi) cell f-I curves
IthIC = 13;                 % fixed threshold
IwidthIC = 8;               % fixed width of sigmoid
rmaxIC = 300;               % fixed maximum rate 
tauIC = 0.010;              % fixed time constant

% Climbing fiber f-I curve
IthCF = 15;                 % fixed threshold
IwidthCF = 5;               % fixed width of sigmoid
rmaxCF = 300;               % fixed maximum rate 
tauCF = 0.002;              % fixed time constant

IthPC = 15;                 % fixed threshold
IwidthPC = 5;               % fixed width of sigmoid
rmaxPC = 30;                % fixed maximum rate 
tauPC = 0.010;              % fixed time constant

IthNI = 10;                 % fixed threshold
IwidthNI = 3;               % fixed width of sigmoid
rmaxNI = 100;               % fixed maximum rate 
tauNI = 0.005;              % fixed time constant

% Now set all the connection strengths only WGP will vary

WGP0 = 30.0/NGcells;        % Base granule cell to Purkinje Cell connectin strength
WGP = WGP0*ones(NGcells,1); % Vector of individual connections
dW0 = 0.025;                % Maximum fractional change in weight per trial

WGN = 4.0;                  % Connection from granule cell to nucleus
WPN = 10.0;                 % Inhibitory strength from Purkinje Cll to nucleus

WGG = zeros(NGcells);           % Granule cells to granule cells
WGI = zeros(NGcells,NIcells);   % Granule cells to interneurons
WIG = zeros(NIcells,NGcells);   % Interneurons to granule cells

WGG0 = 1.5;                 % Base connection strengths
WGI0 = 0.8;
WIG0 = 16/NIcells;

sigmaW = 1;
connectprob = 0.15;          % connection probability between granule cells

for cell1 = 1:NGcells
    for cell2 = 1:NGcells
        if rand(1) > connectprob    % only connect if random number is less than probability
            WGG(cell1,cell2) = 0.0; % otherwise no connection
        else                        % connection strength has structure and randomness
            WGG(cell1,cell2) = WGG0*(NGcells-mod(cell2-cell1+NGcells,NGcells))/NGcells ...
                * (1.0+5.0*rand(1))/connectprob;
        end
    end
    WGG(cell1,cell1) = 0.25*WGG0;   % self-connection is fixed
    WGI(cell1,cell1) = WGI0;        % each granule cell connects to one interneuron
end

for cell1 = 1:NIcells
    for cell2 = 1:NIcells           % all to all connections from interneuron to granule cells
        WIG(cell1,cell2) = WIG0*(0.5+rand(1)); % but with randomness in strengths
    end
end

tCSon = 0.2;                % time of CS onset
lengthCS = tmax;            % length of CS
ICS0 = 100.0;               % steady state CS applied current
ICS1 = 100.0;               % intitial peak CS current
tauCS = 0.05;               % time for decay of peak CS current

tUSon = 0.5;                % time for US onset
lengthUS = 0.025;           % length of US
IUS0 = 100.0;               % applied current for US

ICS = zeros(size(t));       % set up input current vectors, one for CS
IUS = zeros(size(t));       % and one for US
Iapp = zeros(length(t),NGcells);    % each granule cell can get a different applied current
for i = max(1,floor(tCSon/dt)+1):min(floor((tCSon+lengthCS)/dt)+1,length(ICS))
    ICS(i) = ICS0+ICS1*exp(-(t(i)-tCSon)/tauCS);
    for cell = 1:NGcells
        Iapp(i,cell) = ICS(i)*(NGcells-cell)/(NGcells-1);
    end
end

for trial = 1:Ntrials
    trial
    rGC = zeros(length(t),NGcells);
    D = ones(length(t),NGcells);
    F = ones(length(t),NGcells);
    S = zeros(length(t),NGcells);
    rI = zeros(length(t),NIcells);
    rPC = zeros(length(t),1);
    rCF = zeros(length(t),1);
    rNI = zeros(length(t),1);
    
    for i = floor(tUSon/dt)+1:floor((tUSon+lengthUS)/dt)+1
        IUS(i) = IUS0;
    end
    ItotPC = zeros(size(t));
    
    for i = 2:length(t)
        ItotGC = Iapp(i-1,:) + S(i-1,:)*WGG - rI(i-1,:)*WIG;
        
        rinfGC = rmaxGC./(1+exp(-(ItotGC-IthGC)./IwidthGC)).^2;    % rinf is a vector of the steady state firing rate
        % given the total
        % current vector to
        % each cell

        rGC(i,:) = rinfGC + (rGC(i-1,:)-rinfGC)*exp(-dt/tauGC);   % update r from r(i-1) for all cells

        Finf = rGC(i,:).*tauf./(rGC(i,:).*tauf + 1.0);
        F(i,:) = Finf + (F(i-1,:)-Finf).*exp(-dt./tauf);
        
        Dinf = 1.0./(rGC(i,:).*taud.*F(i,:) + 1.0);
        D(i,:) = Dinf + (D(i-1,:)-Dinf).*exp(-dt./taud);        
       
        Sinf = rGC(i,:).*D(i,:).*F(i,:);
        S(i,:) = Sinf + (S(i-1,:)-Sinf).*exp(-dt./taus);
        
        ItotIC = rGC(i-1,:)*WGI;
        
        rinfI = rmaxIC./(1+exp(-(ItotIC-IthIC)/IwidthIC)).^2;    % rinf is a vector of the steady state firing rate
        % given the total
        % current vector to
        % each cell

        rI(i,:) = rinfI + (rI(i-1,:)-rinfI)*exp(-dt/tauIC);   % update r from r(i-1) for all cells
        
        % Now climbing fiber
        ItotCF = IUS(i);
        rinfCF = rmaxCF./(1+exp(-(ItotCF-IthCF)/IwidthCF)).^2;
        rCF(i) = rinfCF + (rCF(i-1)-rinfCF)*exp(-dt/tauCF);
        
        ItotPC(i) = S(i-1,:)*WGP;
        rinfPC = rmaxPC./(1+exp(-(ItotPC(i)-IthPC)/IwidthPC)).^2;
        rPC(i) = rinfPC + (rPC(i-1)-rinfPC)*exp(-dt/tauPC);
%        if t(i) > tUSon && t(i) < tUSon+lengthUS
%            rPC(i) = 0.0;
%        end
        
        ItotNI = mean(S(i-1,:))*WGN - rPC(i-1)*WPN;
        rinfNI = rmaxNI./(1+exp(-(ItotNI-IthNI)/IwidthNI)).^2;
        rNI(i) = rinfNI + (rNI(i-1)-rinfNI)*exp(-dt/tauNI);


    end
    figure(1)
    plot(t,rNI);
    hold on
    plot(t,rPC,'r');
    xlabel('time')
    ylabel('Firing rate (red=Purkinje; blue=Nucleus)')
    
    % Now to adjust the synapses from GC to PC 
    ratecorr = zeros(NGcells,1);
    for cell = 1:NGcells
        ratecorr(cell) = sum(S(:,cell).*rCF(:))*dt; % Correlation between
                                                    % CF rate and GC to PC
                                                    % synapses
    end
   
    ms = mean(S,1);         % Mean (time-averaged) synaptic output of each cell
    
    % Calculate change in set of Granule Cell to Purkinje Cell weights 
    % depending on correlation between their firing and Climbing Fiber
    dW = dW0*(ms'.*ratecorr-20)/(mean(ratecorr)*mean(ms));
    dW = max(-0.05,dW);             % Each weight change is no more than 
    dW = min(0.05,dW);              % plus or minus 5%
    WGP = WGP.*(1.0-dW);            % Adjust all synapses
    WGP = max(WGP,0.0);             % Each synapse must be greater than zero
    WGP = min(WGP,3*WGP0);          % and less than some maximum of 3 times base

    WGP = WGP*WGP0/mean(WGP);       % This line scales all weights to keep the 
                                    % mean the same
    
    figure(2)
    plot(WGP)
    xlabel('Granule Cell Index')
    ylabel('Connection Strength, WGP')
    hold on
    figure(3) 
    plot(t,ItotPC)
    xlabel('time')
    ylabel('Purkinje Cell Input')
    
    if ( trial == 1 ) 
        figure(4)
        plot(t,rGC)
    end
end