Simulation of the response of the target population to stimuli

This script has been used to generate Figures 5A and 5B for the MS

Blessing of Dimensionality in Spiking Neural Networks: The by-chance functional learning by Makarov & Lobov

version 2 by V. Makarov (vmakarov@ucm.es). 13-05-2025 Some adjustments and performance improvements

version 1 by V. Makarov (vmakarov@ucm.es). 30-12-2024 Working code to prepare the last figure of the ms.

Contents

Set parameters, load gray-scale images

clear
addpath('./MatFunc')
load('./Images_LastFig/ImagesForFig4.mat') % Load images

nProc = 8;                % number of stimuli to process

P.Tstim = 1;              % pattern (stimulus) duration [s]
P.Fmax = 1000;            % max frequency [Hz]
P.h = 0.05;               % integration step [ms]
P.tau = 10;               % synaptic time scale [ms]
P.gamma = 10;             % magnitude of the inhibitory current
P.eta = 0.25;             % relative learning rate
P.Theta = 2.6;            % firing threshold
P.mu = 1.2;               % learning rate
P.beta = 5.05;            % excitation rate

P.alpha = 1/(P.eta*P.tau*1e-3*P.Fmax);
P.lambda = 4*(P.Theta + P.gamma)/(P.beta*P.tau*1e-3*P.Fmax);
P0 = P; P0.mu = 0;        % No learning

Prepare stimulus vectors from the images and generate input spike trains

rng(4)
n = size(ImgGr{1},1)*size(ImgGr{1},2); % stimulus dimension
x = zeros(n,nProc); % Allocate memory
tij = cell(nProc,1);
aux = linspace(0,1,n);
for k = 1:nProc
    Imin = min(ImgGr{k}(:)); Imax = max(ImgGr{k}(:));
    x(:,k) = 1 - (ImgGr{k}(:)-Imin)/(Imax - Imin); % scale to [0,1] and invert
    [~,id] = sort(x(:,k));
    x(id,k) = aux; % change the pixel distribution to U([0,1])
    tij{k} = GenerateFrequencyPatterns(x(:,k), P.Tstim, P.Fmax); % Generate spike trains
end

Examples of the normalization of images (Fig. 5A)

figure('Color','w','Position',[100 100 600 300])
PlotFig5A(ImgGr, x)

Simulate the population response to input stimuli

L = 500;         % number of neurons
w0 = rand(n,L);  % initial synaptic weights

% Response before learning
disp('Before learning')
[v1, firings1] = ProcessStimuli(nProc, tij, w0, P0); % no learning

% Learning
disp('Learning')
[v2, firings2, wlrn] = ProcessStimuli(nProc, tij, w0, P);

% Response after learning
disp('Testing')
[v3, firings3] = ProcessStimuli(nProc, tij, wlrn, P0);

Before learning
1 of 8 done.
2 of 8 done.
3 of 8 done.
4 of 8 done.
5 of 8 done.
6 of 8 done.
7 of 8 done.
8 of 8 done.
Learning
1 of 8 done.
2 of 8 done.
3 of 8 done.
4 of 8 done.
5 of 8 done.
6 of 8 done.
7 of 8 done.
8 of 8 done.
Testing
1 of 8 done.
2 of 8 done.
3 of 8 done.
4 of 8 done.
5 of 8 done.
6 of 8 done.
7 of 8 done.
8 of 8 done.

Plot Figure 5B

neurons = unique([firings3; firings1; firings2],'stable');% fired neurons

figure('color','w','Position',[100 100 900 700])
PlotFig5B(x, neurons, v1, v2, v3)

AUXILARY FUNCTIONS

function PlotAux(v, firings, Fs)
nf = unique(firings,'stable');
chttl = cell(length(nf),1);
for k = 1:length(nf)
    chttl{k} = num2str(nf(k));
end
if nargin < 3, Fs = 2000; end
PlotEEG(v(nf,:) + 65,Fs,chttl,100)
end

function [v, firings, w] = ProcessStimuli(nProc, tij, w, P)
v = [];
firings = [];
for k = 1:nProc
    [dum, w, aux] = IntegrateIzhikevich_L(P.h, P.Tstim, tij{k}, w, P);
    v = [v aux(:,1:10:end)]; % downsampling
    if ~isempty(dum)
        dum = unique(dum(:,2),'stable');
        firings = [firings; dum];
    end
    disp([num2str(k) ' of ' num2str(nProc) ' done.'])
    pause(0.1)
end
end

function PlotFig5A(ImgGr, x)
for k = 1:2
    subplot(2,4,1 + 4*(k-1))
    image(ImgGr{k+1})
    axis square
    axis off
    hold on
    if k == 1, title('image before'); end

    subplot(2,4,3 + 4*(k-1))
    image(255*reshape(1 - x(:,k+1),[32 32]))
    axis square
    axis off
    if k == 1, title('image after'); end

    subplot(2,4,2 + 4*(k-1))
    histogram(ImgGr{k}(:)/255,0:0.1:1,'Normalization','probability');
    axis square
    axis([0 1 0 0.5])
    if k == 1, title('pxs distribution'); end

    subplot(2,4,4 + 4*(k-1))
    histogram(1 - x(:,k),'Normalization','probability');
    axis square
    axis([0 1 0 0.5])
    if k == 1, title('pxs distribution'); end
end
colormap('gray')
end

function PlotFig5B(x, neurons, v1, v2, v3)
nProc = size(x,2);
for k = 1:nProc
    subplot(4,nProc,k)
    image(255*reshape(1-x(:,k),[32 32]))
    axis square
    axis off
end
colormap('gray')

subplot(4,1,2)
PlotAux(v1, neurons, 2000)
title('before learning')

subplot(4,1,3)
PlotAux(v2, neurons, 2000)
title('learning')

subplot(4,1,4)
PlotAux(v3, neurons, 2000)
title('after learning')
end

Published with MATLAB® R2024a