Brian2_脉冲神经网络_神经元学习记录

脉冲神经网络被称作是第三代神经网络，生物可信度更高，在近年来兴起的类脑科学研究中，SNN也一直占据着核心地位。并且在性能相近的情况下，基于脉冲神经网络制成的芯片相较于人工神经网络功耗更低，稳定性、鲁棒性等也更为优异。
通过pip方式安装了Brian2，并运行了一些官网上提供的程序，也附了一些自己的学习心得。

from brian2 import *
start_scope()
#start_scope()函数确保在调用该函数之前创建的任何 
#Brian对象都不会包含在下一次模拟运行中。
tau = 10*ms
eqs ='dv/dt = (1-v)/tau : 1'
G = NeuronGroup(1, eqs,method='exact')
print('Before v = %s' % G.v[0])
run(100*ms)
print('After v = %s' % G.v[0])
start_scope()
G = NeuronGroup(1, eqs, method='exact')
M = StateMonitor(G, 'v', record=0)
run(30*ms)
plot(M.t/ms, M.v[0], 'C0', label='Brian')
plot(M.t/ms, 1-exp(-M.t/tau), 'C1--',label='Analytic')
xlabel('Time (ms)')
ylabel('v')
legend()
show()

程序中要注意量纲，如将不同单位的量进行相加、等式左右量纲不一致程序都会报错。如在eqs ='dv/dt = (1-v)/tau : 1'中去掉/tau就会出现量纲不一致而报错。

该程序中，创建神经元组函数NeuronGroup（）创建了一个以微分方程dv/dt = (1-v)/tau : 1定义好模型的神经元，我们使用了 StateMonitor函数来监测神经元状态。前两个参数是要记录的组和要记录的变量。我们还指定record=0. 这意味着我们记录神经元 0 的所有值。我们必须指定我们想要记录哪些神经元，因为在具有许多神经元的大型模拟中，它通常会占用太多 RAM 来记录所有神经元的值。

蓝色实线即为该神经元v的变化，橙色虚线为微分方程dv/dt = (1-v)/tau : 1的解析解1-exp(-t/tau)的函数曲线，二者是重合的。

start_scope()
tau = 10*ms
eqs = '''
dv/dt = (1-v)/tau : 1
'''
G = NeuronGroup(1, eqs, threshold='v>0.8', reset='v = 0', method='exact')
M = StateMonitor(G, 'v', record=0)
run(50*ms)
plot(M.t/ms, M.v[0])
xlabel('Time (ms)')
ylabel('v')
show()

我们在NeuronGroup声明中添加了两个新关键字： threshold='v>0.8'和reset='v = 0'。这意味着当我们发射尖峰时，并在尖峰后立即重置。

from brian2 import *
start_scope()
tau = 10*ms
eqs ='dv/dt = (1-v)/tau : 1'
G = NeuronGroup(1, eqs, threshold='v>0.8', reset='v = 0', method='exact')
statemon = StateMonitor(G, 'v', record=0)
spikemon = SpikeMonitor(G)
run(50*ms)
print('Spike times: %s' % spikemon.t[:])
plot(statemon.t/ms, statemon.v[0])
for t in spikemon.t:
    axvline(t/ms, ls='--', c='C1', lw=3)
xlabel('Time (ms)')
ylabel('v')
show()

脉冲时间：Spike times: [16. 32.1 48.2] ms

神经元模型的一个共同特征是不应期。这意味着在神经元发射一个尖峰之后，它会在一段时间内变得顽固，并且在这段时间结束之前不能再发射另一个尖峰。以下是我们在 Brian 中的做法。

from brian2 import *
start_scope()
tau = 10*ms
eqs ='dv/dt = (1-v)/tau : 1'
start_scope()
tau = 10*ms
eqs = '''
dv/dt = (1-v)/tau : 1 (unless refractory)
'''
G = NeuronGroup(1, eqs, threshold='v>0.8', reset='v = 0', refractory=5*ms, method='exact')
statemon = StateMonitor(G, 'v', record=0)
spikemon = SpikeMonitor(G)
run(50*ms)
plot(statemon.t/ms, statemon.v[0])
for t in spikemon.t:
    axvline(t/ms, ls='--', c='C1', lw=3)
xlabel('Time (ms)')
ylabel('v')
show()

虽然我们在NeuronGroup函数中添加了参数refractory=5*ms，但是在微分方程中还是要添加(unless refractory)，否则神经元的行为不会产生不应期。

from brian2 import *
start_scope()
tau = 5*ms
eqs = '''
dv/dt = (1-v)/tau : 1
'''
G = NeuronGroup(1, eqs, threshold='v>0.8', reset='v = 0', refractory=15*ms, method='exact')
statemon = StateMonitor(G, 'v', record=0)
spikemon = SpikeMonitor(G)
run(50*ms)
plot(statemon.t/ms, statemon.v[0])
for t in spikemon.t:
    axvline(t/ms, ls='--', c='C1', lw=3)
axhline(0.8, ls=':', c='C2', lw=3)
xlabel('Time (ms)')
ylabel('v')
print("Spike times: %s" % spikemon.t[:])
show()

那么这里发生了什么？第一个尖峰的行为是相同的： v上升到 0.8，然后神经元在0到8 毫秒的时间内触发尖峰，然后立即重置为 0。由于不应期现在是 15 毫秒，这意味着神经元将无法尖峰再次直到时间 8 + 15 = 23 毫秒。在第一个尖峰之后，v的值立即开始上升，因为我们没有在NeuronGroup的定义中指定(unless refractory)，此时虽然它可以在大约 8 毫秒的时间达到值 0.8（绿色虚线），但因为神经元现在正处于不应期，故到达阈值后不会重置反而继续上升，直到不应期结束电位才会被复位。

下面演示多个神经元。

from brian2 import *
start_scope()
N = 100
tau = 10*ms
eqs = '''
dv/dt = (2-v)/tau : 1
'''
G = NeuronGroup(N, eqs, threshold='v>1', reset='v=0', method='exact')
G.v = 'rand()'
spikemon = SpikeMonitor(G)
run(50*ms)
plot(spikemon.t/ms, spikemon.i, '.k')
xlabel('Time (ms)')
ylabel('Neuron index')
show()

变量N 来确定神经元的数量，G.v = 'rand()'是用 0 到 1 之间不同的随机值初始化每个神经元，这样做只是为了让每个神经元做一些不同的事情。spikemon.t表示神经元尖峰时间，spikemon.i为每个尖峰提供相应的神经元索引值。这是神经科学中使用的标准“光栅图”。

from brian2 import *
start_scope()
N = 100
tau = 10*ms
v0_max = 3.
duration = 1000*ms
eqs = '''
dv/dt = (v0-v)/tau : 1 (unless refractory)
v0 : 1
'''
G = NeuronGroup(N, eqs, threshold='v>1', reset='v=0', refractory=5*ms, method='exact')
M = SpikeMonitor(G)
G.v0 = 'i*v0_max/(N-1)'
run(duration)
figure(figsize=(12,4))
subplot(121)
plot(M.t/ms, M.i, '.k')
xlabel('Time (ms)')
ylabel('Neuron index')
subplot(122)
plot(G.v0, M.count/duration)
xlabel('v0')
ylabel('Firing rate (sp/s)')
show()

在这个例子中，我们将神经元以v0 指数方式驱动到该值，但是当v>1时，它会触发一个尖峰并重置。结果是它发射尖峰的速率将与v0的值相关。因为v0<1它永远不会发射尖峰，而且随着v0变大它会以更高的速率发射尖峰。右图显示了v0的函数的点火率 。这是这个神经元模型的 If 曲线。

请注意，在图中我们使用了SpikeMonitor中的变量count: 这是组中每个神经元被激发的尖峰数量的数组。将其除以运行的持续时间得出点火率。

from brian2 import *
start_scope()
N = 100
tau = 10*ms
v0_max = 3.
duration = 1000*ms
sigma = 0.2
eqs = '''
dv/dt = (v0-v)/tau+sigma*xi*tau**-0.5 : 1 (unless refractory)
v0 : 1
'''
G = NeuronGroup(N, eqs, threshold='v>1', reset='v=0', refractory=5*ms, method='euler')
M = SpikeMonitor(G)
G.v0 = 'i*v0_max/(N-1)'
run(duration)
figure(figsize=(12,4))
subplot(121)
plot(M.t/ms, M.i, '.k')
xlabel('Time (ms)')
ylabel('Neuron index')
subplot(122)
plot(G.v0, M.count/duration)
xlabel('v0')
ylabel('Firing rate (sp/s)')
show()

Often when making models of neurons, we include a random element to model the effect of various forms of neural noise. In Brian, we can do this by using the symbol xi in differential equations. Strictly speaking, this symbol is a “stochastic differential(随机微分)” but you can sort of thinking of it as just a Gaussian random variable with mean 0 and standard deviation 1. We do have to take into account the way stochastic differentials scale with time, which is why we multiply it by tau**-0.5 in the equations below (see a textbook on stochastic differential equations for more details). Note that we also changed the method keyword argument to use 'euler' (which stands for the Euler-Maruyama method); the 'exact' method that we used earlier is not applicable to stochastic differential equations(随机微分方程).

这一段程序我没有看懂，故把原文也附了上来，看些文献后研究，相比于上一块代码，神经网络中增加了一些噪音，使得点火率以S形方式上升，以下为运行结果。