%STDP

%create one postsynaptic neuron with many presynaptic inputs

clear all;
nneurons=100; %number of neurons
nsteps = 200; %number of iteration steps
prob = 0.03; %probability to emit a spike at each step
spike_threshold = 2.5; %threshold for postsynaptic spiking
muu = 0.1; %learning rate
time = 1:nsteps; %time steps for plotting
amin = 0.1;
aplus = 0.01;
gmax = 10;
old_m = 0.0;

%initialize synapses to 0.5

for n=1:nneurons
synapse(n) = 0.5;
p(n) = 0.0;
old_p(n) = 0.0;
end
for step=1:nsteps
output(step) = 0.0;
input(step) = 0.0;
end;
m=0.0;

%create randomized spiking activity in presynaptic neurons
for n=1:nneurons
for step=1:nsteps
if rand < prob
neuron(n,step) = 1.0;
else
neuron(n,step) = 0.0;
end

end
end;

figure ;
subplot (1,2,1);
axis ([-1, nsteps, 0.5, nneurons+20]);
hold;
for i=1:nneurons;
plot (time, i*neuron(i,:), 'green .');
end

%at each time step, calculate weighted input to postsynaptic neurons
%and also calculate output

for step=1:nsteps
input(step) = 0.0;
for n=1:nneurons
input(step) = input(step) + synapse(n)*neuron(n,step);
end;

if input(step) > spike_threshold
output(step) = 1.0;

else
output (step) = 0.0;
end;
%STDP rule
if (output(step) ~= 0.0)
in_m = -amin;
else
in_m = 0.0;
end
m = exp (-1/1)*old_m+(1-exp(-1/1))*in_m;
old_m = m;

for n=1:nneurons
if (neuron(n,step) ~= 0.0)
in_n = aplus;
else
in_n = 0.0;
end
p(n) = exp(-1/10)*old_p(n)+(1-exp(-1/10))*in_n;
old_p(n) = p(n);
end
change synaptic weights proprotionally to their value
postsynaptic spikes
if (output(step) ~= 0.0)
for n=1:nneurons
synapse(n) = synapse(n) + p(n)*gmax;
if (synapse(n) > gmax)
synapse(n) = gmax;
end
end
end
presynaptic spikes
for n=1:nneurons
if (neuron(n,step) ~= 0.0)
synapse(n) = synapse(n) + m*gmax;
if (synapse(n) < 0.0)
synapse(n) = 0.0;
end
end
end

end;
pause (0.05);
subplot (1,2,1);
plot(time, (nneurons+2)*output, 'red .');
plot (time, nneurons+5+2*input, 'blue');
subplot (1,2,2);
axis ([0,1,2.0,nneurons+9]);
barh(synapse);
end;