Quick example to create a cluster of interconnected spiking neurons controlled by a leaky integrator for the action potential. A simple input is provided characterised by the same numerosity of the neural cluster. All parameters can be changes to see the effect of different types of connectivity on the overall activity.
NB the simulation relies on a “step” function. See below for code and instructions.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 |
clear all close all SECS=0.5; %simulated time (in seconds) dt=0.5; T= SECS*1000/dt; dim1=100; %number of units per cluster A.potential = ones(T,dim1)*(-70); %initialization for all the units in the cluster A.activity = ones(T,dim1)*0; input= zeros(T,dim1); %initialization of the input %matrix of connections in the cluster perc_GABA=0.5; %0-1 scale perc_conn=0.25; %0-1 scale: 1=fully intercconnected (independently of value) intercon=1.00; %0-1 scale: strength of connections in the cluster self_mat=(rand(dim1)>(1-perc_conn)).*(((rand(dim1) - perc_GABA)*intercon)); %matrix of connections between input and cluster perc_GABA=0.0; %0-1 scale perc_conn=0.15; %0-1 scale: 1=fully intercconnected (independently of value) intercon=1.00; %0-1 scale: strength of connections in the cluster inp_cl_mat=(rand(dim1)>(1-perc_conn)).*(((rand(dim1) - perc_GABA)*intercon)); decay=3; baseline=-70; threshold=baseline+15; noise=0; %if noise=1 -> bl=bl+-0.1 for t=2:T A=cluster_step(t, dt, A, self_mat, input, inp_cl_mat, decay, baseline, threshold, noise); %%INPUT if t>10 && t<100 input(t,:)= (rand(1,dim1)<0.12)*40; %percentage of active units in the input, per step (0.5 milliseconds) elseif t>100 && t<=300 input(t,:)= (rand(1,dim1)<0.03)*40; elseif t>300 && t<=500 input(t,:)= (rand(1,dim1)<0.06)*40; elseif t>500 && t<=700 input(t,:)= (rand(1,dim1)<0.09)*40; elseif t>700 && t<=900 input(t,:)= (rand(1,dim1)<0.03)*40; elseif t>900 && t<=1000 input(t,:)= (rand(1,dim1)<0.15)*40; end if t>7 for ii=1:dim1 %control to avoid multiple spikes in refractory period if sum(input((t-6):(t-1),ii))>0 input(t,ii)=0; end end end end h=figure; screen_size=get(0,'ScreenSize'); screen_size(1:2)=screen_size(3:4)*1/10; screen_size(3:4)=[1600 1000]; set(h,'Position',screen_size) colormap(gray) subplot(3,1,1) imagesc (input') xlabel ('Milliseconds', 'FontSize', 14) title('INPUT', 'FontSize', 18) subplot(3,1,2) plot(A.potential) hold on plot ([0 T],[-55 -55], ':k') xlabel ('Milliseconds', 'FontSize', 14) axis([0 T (baseline-15) (threshold+10)]) title('Cluster potential', 'FontSize', 18) subplot(3,1,3) imagesc(A.activity') title('Cluster activity', 'FontSize', 18) xlabel ('Milliseconds', 'FontSize', 14) %average for cluster for ind=10:length(A.activity(:,1)) avg_cl(ind)=mean(sum(A.activity(ind-9:ind,:))); inp_cl(ind)=mean(sum(input(ind-9:ind,:))); end g=figure; screen_size=get(0,'ScreenSize'); screen_size(1:2)=screen_size(3:4)*1/10; screen_size(3:4)=[1600 1000]; set(g,'Position',screen_size) subplot(3,1,1) plot(inp_cl) title('Input (average for 5 milliseconds)', 'FontSize', 18) xlabel ('Milliseconds', 'FontSize', 14) subplot(3,1,3) plot(avg_cl) title('Cluster activity (average for 5 milliseconds)', 'FontSize', 18) xlabel ('Milliseconds', 'FontSize', 14) |
Step function to update the activity of the units in the cluster. Save this code as a separate file with the name “cluster_step” in the same folder as the main file above.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
function output = cluster_step(t, dt, structure, self_w, input, input_mat, decay, baseline, th_out, noise_switch) dim=length(structure.activity(t,:)); baseline=baseline + ((rand(1, dim)*20)-10) .*noise_switch; %action potential structure.potential(t,:) = structure.potential(t,:) + (dt./decay).*(-structure.potential(t-1,:) + baseline +... input(t-1,:)*input_mat + (structure.activity(t-1,:)*self_w)); for in=1:dim if t>6 && sum(structure.activity(t-6:t,in))>=1 structure.potential(t,in) = -70; %disp ('Check') end end %transfer function structure.activity(t,:) = 40*(structure.potential(t,:)>th_out); output = structure; |
