clearvars -global
clear all
%This script sets the parameters and launches a training script
%Choices: Hardcoded weight matrix (synfire chain) or learned weight matrix
%Run generateSequence script to simulate spontaneous dynamics
%ABCBA supervisor is hardcoded in this script, supervisor of Fig 6 is loaded from data folder
%different lengths
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulation time parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tStart = 0;
tEnd = 12000; %Simulation in milli-seconds
tStep = 0.1; %0.1 millisecond time step
time = [tStart:tStep:tEnd];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% WEIGHT MATRIX OF TEMPORAL BACKBONE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hardcoded = false;
thirty_clusters = true;
if hardcoded
synfire = true;
if synfire
mult = 0.15;
numClusters = 120;
EneuronNum = 120;
IneuronNum = 30;
neuronNum = 150;
createWeightMatrixEIF
weightsEE = diag(800*ones(EneuronNum-1,1),1);
weightsEE(end,1) = 800;
else
mult = 3;
EneuronNum = mult*800; %Number of excitatory neurons in the network
numClusters = mult*10; %Number of clusters
IneuronNum = round(0.25*EneuronNum); %Number of inhibitory neurons in the network
neuronNum = EneuronNum + IneuronNum; %Total number of neurons
createWeightMatrixEIF %manually set the connectivities
end
else
if thirty_clusters
%matrix60b has 30 exc clusters
A = load('FF60b_2400_nc80.mat'); %import weight matrix
A = A.A;
else
%This matrix has 80 exc clusters, very large simulation: may run
%into memory errors on local pc. Cluster recommended.
A = load('FF180_6400b.mat');
A = A.A;
end
neuronNum = size(A,1);
EneuronNum = 0.80*size(A,1);
IneuronNum = 0.20*size(A,1);
numClusters = EneuronNum/80; %80 neurons/cluster
weightsEE = A(1:EneuronNum,1:EneuronNum);
weightsEI = A(1:EneuronNum,EneuronNum+1:neuronNum);
weightsIE = A(EneuronNum+1:neuronNum,1:EneuronNum);
weightsII = A(EneuronNum+1:neuronNum,EneuronNum+1:neuronNum);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUPERVISOR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if thirty_clusters
actionNeuronNum = 3; %read-out neurons
superNeuronNum = 3; %supervisor neurons
interNeuronNum = 3; %interneurons
%75 ms stimulations
superLength = 400; %[ms]
supervisor = ones(actionNeuronNum,superLength/tStep+1); %baseline input
supervisor(1,25/tStep:100/tStep) = 10.;
supervisor(1,325/tStep+1:400/tStep) = 10.;
supervisor(2,100/tStep+1:175/tStep) = 10.;
supervisor(2,250/tStep+1:325/tStep) = 10.;
supervisor(3,175/tStep+1:250/tStep) = 10.;
else
supervisor = load('supervisor_600.mat'); %600 ms part of a bird song
supervisor = supervisor.supervisor;
actionNeuronNum = size(supervisor,1);
superNeuronNum = actionNeuronNum;
interNeuronNum = actionNeuronNum;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% WEIGHT MATRIX TO ACTION READ OUT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%initially all zero read-out weights, all-to-all connected
W_AE = zeros(actionNeuronNum,EneuronNum);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parameters for both E and I Neurons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Vreset = -60; %Reset for both exc and inh neurons
C = 300; %capacitance
tau_abs = 5; %refractory period
tau_absA = 1; %refractory period for read-out neurons
tau_absS = 1; %refractory period for supervisor neurons
tau_absH = 1; %refractory period for interneurons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parameters for the E-Neurons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Vthres = 20; %Spiking threshold for exc neurons
tau_E = 20; %Membrane time constant
V_E = -70; %resting potential
DET = 2; %slope of exponential
E_E = 0; %reversal potential
V_T = -52; %threshold potential (the spiking threshold for inh neurons)
A_T = 10; %post spike threshold potential increase
tau_T = 30; %adaptive threshold time scale
EVthreshold = V_T*ones(1,EneuronNum); %neuronal threshold vector for all E neurons
EVthresholdA = V_T*ones(1,actionNeuronNum); %neuronal threshold vector for all read-out neurons
EVthresholdS = V_T*ones(1,superNeuronNum); %neuronal threshold vector for all supervisor neurons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parameters for the I-Neurons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tau_I = 20; %Membrane time constant
V_I = -62; %resting potential
E_I = -75; %reversal potential
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parameters for the adaptation (Exc only)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tau_w = 100; %adaptation time constant in the recurrent network
tau_wA = 10; %adaptation time constant in the read-out
tau_wS = 10; %adaptation time constant in the supervisor
a = 0; %adaptation slope in the recurrent network
aA = 0; %adaptation slope of read-out neurons
aS = 0; %adaptation slope of supervisor neurons
b = 1000; %adaptation amplitude in the recurrent network
bA = 0; %adaptation amplitude of read-out neurons
bS = 0; %adaptation amplitude of supervisor neurons
w = a*(Vreset-V_E)*ones(1,EneuronNum); %initial adaptation vector for all recurrent network exc neurons
wA = aA*(Vreset-V_E)*ones(1,actionNeuronNum); %initial adaptation vector for all read-out neurons
wS = aS*(Vreset-V_E)*ones(1,superNeuronNum); %initial adaptation vector for all supervisor neurons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parameters for the synapses and neuronal conductances
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tauedecay = 6; %decay time for e-synapses
tauerise = 1; %rise time of e-synapses
tauidecay = 2; %decay time for i-synapses
tauirise = 0.5; %rise time of i-synapses
xedecay = zeros(1,neuronNum);
xerise = zeros(1,neuronNum);
xidecay = zeros(1,neuronNum);
xirise = zeros(1,neuronNum);
xedecayA = zeros(1,actionNeuronNum);
xeriseA = zeros(1,actionNeuronNum);
xidecayA = zeros(1,actionNeuronNum);
xiriseA = zeros(1,actionNeuronNum);
xedecayS = zeros(1,superNeuronNum);
xeriseS = zeros(1,superNeuronNum);
xedecayH = zeros(1,interNeuronNum);
xeriseH = zeros(1,interNeuronNum);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% External input to both E and I clockneurons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rex = 4.50; %external rate to E-neurons, for synfire chain: 2.75
rix = 2.25; %external rate to I-neurons
Jeex = 1.6; %weights for ee external input
Jiex = 1.52; %weights for ie external input
nextx = zeros(1,neuronNum); %vector containing the next external input spike times for recurrent network neurons
nextx(1,1:EneuronNum) = exprnd(1,1,EneuronNum)/rex;
nextx(1,1+EneuronNum:end) = exprnd(1,1,IneuronNum)/rix;
rx = zeros(1,neuronNum);
rx(1,1:EneuronNum) = rex;
rx(1,EneuronNum+1:end) = rix;
forwardInputsEPrev = zeros(1,neuronNum);
forwardInputsIPrev = zeros(1,neuronNum);
rhx = 1; %external rate to H-neurons
jHX = 1.78; %external weight
jAH = 200; %H to A weight (non-plastic)
jHA = 200; %A to H weight (non-plastic)
nextxH = exprnd(1,1,interNeuronNum)/rhx; %initialize vector containing the next external input spike times for interneurons
forwardInputsHPrev = zeros(1,interNeuronNum);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Input to read-out neurons (from supervisor neurons and E-neurons)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Jas = 200; %weight strenght for connections from the supervisor neurons to read-out neurons
forwardInputsAEPrev = zeros(1,actionNeuronNum); %vector containing excitatory input to the read-outs (rec. netw. + supervisor input)
forwardInputsAIPrev = zeros(1,actionNeuronNum); %vector containing inhibitory input to the read-outs (interneuron input)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% External input to supervisor neurons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Jsex = 1.78; %weights for se external target sequence input
nextxS = exprnd(1,1,superNeuronNum); %initialize next incoming external spike to supervisor
forwardInputsSPrev = zeros(1,superNeuronNum); %vector containing excitatory input to the supervisor neuron (external only)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parameters for the synaptic plasticity (vSTDP)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
th_LTP = -49; %LTP threshold constant
th_LTD = -70; %LTD threshold constant
w_max = 25; %maximal read-out weight strength
A_LTP = 0.0008; %LTP amplitude constant
A_LTD = 0.0014; %LTD amplitude constant
tau_x = 5; %time constant of presynaptic low pass filtered spike train
tau_u = 10; %time constant of postsynaptic low pass filtered membrane voltage (LTD)
tau_vs = 7; %time constant of postsynaptic low pass filtered membrane voltage (LTP)
x = zeros(1,EneuronNum); %low pass filtered presynaptic spike train
u = Vreset+(V_T-Vreset)*rand(actionNeuronNum,(tEnd - tStart)/tStep + 1); %low pass filtered postsynaptic membrane voltage (LTD)
vs = Vreset+(V_T-Vreset)*rand(actionNeuronNum,(tEnd - tStart)/tStep + 1);%low pass filtered postsynaptic membrane voltage (LTP)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Training
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
trainer
%plot of read-out weights
figure;
imagesc(W_AE(end:-1:1,:))
xlabel('Recurrent network neuron index', 'FontSize',25)
ylabel('Read-out neuron', 'FontSize',25)
xticks([1 400 800 1200 1600 2000 2400])
yticks([1 2 3])
yticklabels({'C','B','A'})
colorbar
%generateSequence