【论文】Polychronization: Computation with Spikes

背景回顾：
Izhikevich模型是怎麼回事：http://tieba.baidu.com/p/1967164816
传统ANN（频率编码）的理论基础：http://tieba.baidu.com/p/2060148850
目前所有工程上应用的神经网络都是ANN，它认为单个神经元所携带的一切信息都包含在动作电位的频率之中，每个动作电位在何时产生并不重要。ANN确实解决了一些工程问题，比如最近的Deep learning，但是ANN显然并不足以对人脑中的大量事实做出解释：
最明显的比如
1. 突触权重的改变，遵循STDP原则，对单个spike的时间高度敏感。
2. 视皮质表征的topology的作用，无法用任何频率编码的理论解释，而可以解释为操纵传输速度的需求。
3. 固定频率的脑电波的成因及作用
Izh神经元属于SNN(Spike NN)，它能用很少的计算代价模拟神经元的每个spike在何时发生，而非像ANN一样仅用频率赋值。据估算2017年以后的计算能力已经足够Izh模型人脑规模的神经网络。
Izh模型的缺陷：单纯izh模型虽然比ANN有巨大进步，但仍然是一种简化，它不包含对树突的计算，而树突很可能起到逻辑门的作用，会对神经元输出造成很大影响。

本论文中的神经网络与传统ANN的区别：
1. 采用Izh神经元模型而非ANN神经元模型，ANN只负责描述神经元在何时输出何种频率，而Izh模型描述神经元在每个时间点上是否输出spike(动作电位)
2. 突触的区别，ANN中突触只有一个变量即权重；而本文中突触有两个变量，权重和传输时间，传输时间的生理学解释为，动作电位从轴丘发生，经过轴突-突触-树突，对后一个神经元的轴丘施加潜在影响的时间。

这两天干别的事去了，刚刚才开始看代码
izhikevich模型：
V=V+0.04*V^2+5*V+140-u+I;
u=u+a*(b*V-u);
if V>30
V=c;
u=u+d;
end
填上参数就是了，V是膜电位；u是隐变量；I是输入电流；abcd是神经元参数
izhikevich模型演示：http://www.izhikevich.org/publications/izhikevich.m

加注释版
% spnet.m: Spiking network with axonal conduction delays and STDP% Created by Eugene M.Izhikevich. February 3, 2004M=100; % number of synapses per neuronD=20; % maximal conduction delay% excitatory neurons % inhibitory neurons % total numberNe=800; Ni=200; N=Ne+Ni;a=[0.02*ones(Ne,1); 0.1*ones(Ni,1)];%izh神经元参数a初始化d=[ 8*ones(Ne,1); 2*ones(Ni,1)];%izh神经元参数d初始化 sm=10; % maximal synaptic strengthpost=ceil([N*rand(Ne,M);Ne*rand(Ni,M)]);%1000*100矩阵，每个神经元的突触前100个神经元序号，%抑制性只连兴奋性，s=[6*ones(Ne,M);-5*ones(Ni,M)]; % synaptic weightssd=zeros(N,M); % their derivativesfor i=1:Nif i<=Nefor j=1:Ddelays{i,j}=M/D*(j-1)+(1:M/D);%第i个神经元延迟j时间的突触序号end;elsedelays{i,1}=1:M;%抑制性神经元都是1ms延时end;pre{i}=find(post==i&s>0); % pre excitatory neurons % 第i神经元之前的兴奋性神经元aux{i}=N*(D-1-ceil(ceil(pre{i}/N)/(M/D)))+1+mod(pre{i}-1,N);%an auxiliary table of indices%needed to speed up STDP implementation.end;STDP = zeros(N,1001+D);v = -65*ones(N,1); % initial valuesu = 0.2.*v; % initial valuesfirings=[-D 0]; % spike timingsfor sec=1:60*60*24 % simulation of 1 dayfor t=1:1000 % simulation of 1 secI=zeros(N,1);I(ceil(N*rand))=20; % random thalamic inputfired = find(v>=30); % indices of fired neuronsv(fired)=-65;%izhikevich规则u(fired)=u(fired)+d(fired);%izhikevich规则STDP(fired,t+D)=0.1;for k=1:length(fired)sd(pre{fired(k)})=sd(pre{fired(k)})+STDP(N*t+aux{fired(k)});%synapse dirivatives adjustmentend;firings=[firings;t*ones(length(fired),1),fired];%何时何神经元激活过的记录k=size(firings,1);%截止到t，总神经元激活次数while firings(k,1)>t-Ddel=delays{firings(k,2),t-firings(k,1)+1};%firings列表第k个神经元，应在t时间输出动作电位的突触的序号%k初始为列表长度，随后逐渐-1，直到firings(k,1)<=t-Dind = post(firings(k,2),del);%此突触后的神经元序号I(ind)=I(ind)+s(firings(k,2), del)';%突触后神经元的输入电流+权重（*单位电流）sd(firings(k,2),del)=sd(firings(k,2),del)-1.2*STDP(ind,t+D)';k=k-1;end;v=v+0.5*((0.04*v+5).*v+140-u+I); % for numericalv=v+0.5*((0.04*v+5).*v+140-u+I); % stability timeu=u+a.*(0.2*v-u); % step is 0.5 msSTDP(:,t+D+1)=0.95*STDP(:,t+D); % tau = 20 msend;plot(firings(:,1),firings(:,2),'.');axis([0 1000 0 N]); drawnow;STDP(:,1:D+1)=STDP(:,1001:1001+D);%下一秒的STDP初，更改为上一秒的STDP末ind = find(firings(:,1) > 1001-D);%对下一秒初可能施加影响的神经元firings=[-D 0;firings(ind,1)-1000,firings(ind,2)];%重置firings列表s(1:Ne,:)=max(0,min(sm,0.01+s(1:Ne,:)+sd(1:Ne,:)));%s<--s+0.1+sdsd=0.9*sd;clcfprintf([num2str(sec/24/60/60),'%%\n'])end;

关于stdp有个问题，如果synapse权重的改变量与两个spike的相对时间又如上曲线关系，那么应该是E（Δw）=∑p(t)*pdf(t)才对
p(t)是间隔时间t的概率，pdf(t)是如上曲线的方程，E()是期望。那么岂不是如果E(Δw)>0，权重就会趋向于无限？以上曲线应该会随w变动吧

%spnet.m: Spiking network with axonal conduction delays and STDP% Created by Eugene M.Izhikevich. February 3, 2004
M=100; % number of synapses per neuron
D=20; % maximal conduction delay
%excitatory neurons % inhibitory neurons % total number
Ne=800;
Ni=200;
N=Ne+Ni;
a=[0.02*ones(Ne,1); 0.1*ones(Ni,1)];%izh神经元参数a初始化
d=[ 8*ones(Ne,1); 2*ones(Ni,1)];%izh神经元参数d初始化
sm=10; % maximal synaptic strength
post=ceil([N*rand(Ne,M);Ne*rand(Ni,M)]);%1000*100矩阵，每个神经元的突触前100个神经元序号，%抑制性只连兴奋性，
s=[6*ones(Ne,M);-5*ones(Ni,M)]; % synaptic weights
sd=zeros(N,M); % their derivatives
for i=1:N
if i<=Ne
for j=1:D
delays{i,j}=M/D*(j-1)+(1:M/D);%第i个神经元延迟j时间的突触序号,即(1:5)+5＊（j-1）
%比如delay{i,2}=[6,7,8,9,10],delay{i,3}=[11,12,13,14,15]
end;
else
delays{i,1}=1:M;%抑制性神经元都是1ms延时
end;
pre{i}=find(post==i&s>0); % 第i个神经元的pre excitatory synapse 在post列表中的序号，约有80个元素，最大值为100*1000
aux{i}=N*(D-1-ceil(ceil(pre{i}/N)/(M/D)))+1+mod(pre{i}-1,N);
% ceil(pre{i}/N) 是这些synapse在post列表中的列数，即此神经元的第几个突触
% ceil(突触序号／（M/D）)即是其延迟时间，因为延迟j毫秒的突触为M/D*(j-1)+(1:M/D)见delays(i,j)的设置
% mod(pre{i}-1,N) 定位 其在post 列表中的行数，即突触前神经元序号
% an auxiliary table of indices
% needed to speed up STDP implementation.
end;
STDP = zeros(N,1001+D);
v = -65*ones(N,1); % initial values
u = 0.2.*v; % initial values
firings=[-D 0]; % spike timings
for sec=1:60*60*24 % simulation of 1 day
for t=1:1000 % simulation of 1 sec
I=zeros(N,1);
I(ceil(N*rand))=20; % random thalamic input
fired = find(v>=30); % indices of fired neurons
v(fired)=-65;%izhikevich规则
u(fired)=u(fired)+d(fired);%izhikevich规则
STDP(fired,t+D)=0.1;
for k=1:length(fired)
sd(pre{fired(k)})=sd(pre{fired(k)})+STDP(N*t+aux{fired(k)});%synapse dirivatives adjustment
end;
firings=[firings;t*ones(length(fired),1),fired];%何时何神经元激活过的记录
k=size(firings,1);%截止到t，总神经元激活次数
while firings(k,1)>t-D
del=delays{firings(k,2),t-firings(k,1)+1};%firings列表第k个神经元，应在t时间输出动作电位的突触的序号%k初始为列表长度，随后逐渐-1，直到firings(k,1)<=t-D
ind = post(firings(k,2),del);%此突触后的神经元序号
I(ind)=I(ind)+s(firings(k,2), del)';%突触后神经元的输入电流+权重（*单位电流）
sd(firings(k,2),del)=sd(firings(k,2),del)-1.2*STDP(ind,t+D)';
k=k-1;
end;
v=v+0.5*((0.04*v+5).*v+140-u+I); % for numerical
v=v+0.5*((0.04*v+5).*v+140-u+I); % stability time
u=u+a.*(0.2*v-u); % step is 0.5 ms
STDP(:,t+D+1)=0.95*STDP(:,t+D); % tau = 20 ms
end;
plot(firings(:,1),firings(:,2),'.');
axis([0 1000 0 N]);
drawnow;
STDP(:,1:D+1)=STDP(:,1001:1001+D);%下一秒的STDP初，更改为上一秒的STDP末
ind = find(firings(:,1) > 1001-D);%对下一秒初可能施加影响的神经元
firings=[-D 0;firings(ind,1)-1000,firings(ind,2)];%重置firings列表
s(1:Ne,:)=max(0,min(sm,0.01+s(1:Ne,:)+sd(1:Ne,:)));%s<--s+0.1+sd
sd=0.9*sd;
clc
fprintf([num2str(sec/24/60/60),'%%\n'])
end;
t=-20:20;
dw=0.1*exp(-t/20);
plot(t,dw)