SNN学习（3）：Brain2Loihi，基于brain2实现的Loihi模拟器

一、基础参考内容

二、Brian2安装与测试

1、brian2安装

（1）安装软件包

（2）安装brian2

（3）测试

三、简单入门

1、神经元Neurons

（1）单神经元

（2）多神经元

（3）初始电压和发送频率的关系

（4）随机神经元

（5）最后的作业

2、突触Synapses

（1）简单突触

（2）为突触增加权重

（3）引入突触传输延迟

（4）复杂连接

（5）不同连接概率

（6）一对一连接

（7）权重图

（8）复杂突触模型STDP

一、基础参考内容

Brain2Loihi，基于brain2实现的Loihi模拟器

C. Michaelis, A. B. Lehr, W. Oed, and C. Tetzlaff, “Brian2loihi: An emulator for the neuromorphic chip loihi using the spiking neural network simulator brian,” Frontiers in Neuroinformatics, vol. 16, 2022.

GitHub - sagacitysite/brian2_loihi: A Loihi emulator based on Brian2

Brian2loihi论文中的三个实验实例：

https://github.com/sagacitysite/brian2_loihi_utils/tree/main

To demonstrate that the Loihi emulator works as expected, we provide three examples covering a single neuron, a recurrently connected spiking neural network, and the application of a learning rule. All three examples are available as Jupyter notebooks

其他可能用到的参考仿真器：dynapse-simulator，Dynap-SE1神经形态硬件仿真器

GitHub - YigitDemirag/dynapse-simulator: A minimal and interpretable Brian2 based DYNAP neuromorphic processor simulator for educational purposes.

二、Brian2安装与测试

官方安装教程参考：https://brian2.readthedocs.io/en/stable/introduction/install.html

1、brian2安装

我的测试实验在服务器的docker中完成的，其中操作系统为ubuntu20.02（纯净版）。

（1）安装软件包

sudo apt-get update
sudo apt-get install vim
sudo apt-get install git
sudo apt-get install python3
sudo apt-get install python3-pip
pip3 install matplotlib

（2）安装brian2

sudo apt update
sudo apt install python3-brian

（3）测试

新建一个python文件LIF.py，运行python3 LIF.py执行程序。

from brian2 import *start_scope()tau = 10*ms
eqs = '''
dv/dt = (1-v)/tau : 1
'''# Change to Euler method because exact integrator doesn't work here
G = NeuronGroup(1, eqs, threshold='v>0.8', reset='v = 0.5', method='exact')
M = StateMonitor(G, 'v', record=0)G.v = 1 # initial value
spikemon = SpikeMonitor(G)statemon = StateMonitor(G, 'v', record=0)
print('Before v = %s' % G.v[0])
run(50*ms)
print('After v = %s' % G.v[0])
print('Spike times: %s' % spikemon.t[:])plot(M.t/ms, M.v[0])
xlabel('Time (ms)')
ylabel('v')
show()
plt.savefig("LIF.png")

（4）运行结果

同目录下生成一个LIF.png，如下图，表示已经安装成功。。

三、简单入门

1、神经元Neurons

（1）单神经元

from brian2 import *start_scope()tau = 10*ms #从非平衡状态到平衡状态，即达到稳态的时间刻度
eqs = '''
dv/dt = (1-v)/tau : 1 (unless refractory)   #如果需要加不应期时间，refractory，就需要加上 (unless refractory)
'''
#一阶微分方程# 1是神经元个数，eqs是定义的方程，threshold是阈值，reset是重置值，refractory为不应期
#method='exact'是使用精确积分器，method='euler'使用欧拉方法
G = NeuronGroup(1, eqs, threshold='v>0.8', reset='v = 0', refractory=5*ms, method='exact')# 创建状态监视器记录神经元组G的电压
# record=True记录每个神经元，record=0记录神经元0
statemon = StateMonitor(G, 'v', record=True)# 创建脉冲寄存器记录神经元组G的脉冲事件
spikemon = SpikeMonitor(G)# 初始电压，如果不设置则为默认值
G.v = 5 # 模拟开始前的电压
print('Before v = %s' % G.v[0])# 运行模拟时间
run(100*ms)# 模拟完成后的电压
print('After v = %s' % G.v[0])# 打印脉冲时间时间
print('Spike times: %s' % spikemon.t[:])# 为每个脉冲事件绘制一条虚线
for t in spikemon.t:axvline(t/ms, ls='--', c='C1', lw=3)#spikemon.t/msspikemon.t是一个数组，包含了所有脉冲的时间戳（以秒为单位），除以ms（即1e-3）将其转换为毫秒。
plot(statemon.t/ms, statemon.v[0])
xlabel('Time (ms)')
ylabel('v')
show()
plt.savefig("tutorial.png")

运行结果：

（2）多神经元

A. 简单的多神经元

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()' #给每个神经元一个0-1范围的随机电压值spikemon = SpikeMonitor(G)run(100*ms)# .k是绘图样式，“.”表示原点，k表示黑色
plot(spikemon.t/ms, spikemon.i, '.k')
xlabel('Time (ms)')
ylabel('Neuron index')
show()
plt.savefig("mul_neuron.png")

这个图理解下来，先看横轴时间，例如在20ms的刻度上，再往y轴看去，如果在（20,99）（20,50）等都有黑色点，说明50、99在20ms都产生了脉冲事件。【纵着看】

（3）初始电压和发送频率的关系

from brian2 import *start_scope()N = 100
tau = 10*ms
v0_max = 3.
duration = 1000*ms# v0 : 1 声明一个无量纲参数v0，其中的1表示无量纲
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')spikemon = SpikeMonitor(G)# 初始化每个神经元的初始电压
G.v0 ='i*v0_max/(N-1)' run(duration)
figure(figsize=(12,4))
subplot(121)
plot(spikemon.t/ms, spikemon.i, '.k')
xlabel('Time (ms)')
ylabel('Neuron index')
subplot(122)# 其中的count是spikemon是spikemon的数组，其中每个元素表示一个神经元
# （sp/s）是要计算单位时间内的脉冲数量
# 初始电压和发放率之间的关系
plot(G.v0, spikemon.count/duration)
xlabel('v0')
ylabel('Firing rate (sp/s)')
show()
plt.savefig("mul_neuron.png")

运行结果如下：

如果v0小于1，则不会发生脉冲行为。如果v0越大，则发送脉冲的速度越快。

之所以能得到上述结论，在G.v0 ='i*v0_max/(N-1)'时就将其中的0-99神经元的初始电压v0控制在（0,3）v的区间范围内。换个说法就是将神经元的初始电压按照神经元编号的大小增加的，在左图中可以看出，编号越大的神经元脉冲发生的越频繁，所以可以看出初始电压和发送频率的关系。

（4）随机神经元

与之前的代码不同的是，这次在微分方程中加入了随机项 sigma*xi*tau**-0.5，以模拟神经元受到的随机噪声影响。

（1）通过在微分方程中使用符号 xi 来实现引入噪声，xi 可以被视为一个均值为0、标准差为1的高斯（正态）随机变量。这意味着 xi 在每个时间点都会随机地产生一个值，这个值遵循高斯分布。由于随机微分方程需要考虑时间的缩放，因此在方程中将 xi 乘以 tau**-0.5。这是随机微分方程的一个特性，确保了随机过程的强度与时间步长相对应。

（2）在处理随机微分方程时，需要使用特定的数值方法。在这里，使用了 'euler' 方法，即欧拉-马里亚马方法，这是一种适用于随机微分方程的数值积分方法。而之前使用的 'exact' 方法不适用于随机微分方程。

from brian2 import *start_scope()N = 100
tau = 10*ms
v0_max = 3.
duration = 1000*ms
sigma = 0.2# v0 : 1 声明一个无量纲参数v0，其中的1表示无量纲
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')
spikemon = SpikeMonitor(G)# 初始化每个神经元的初始电压
G.v0 ='i*v0_max/(N-1)' run(duration)
figure(figsize=(12,4))
subplot(121)
plot(spikemon.t/ms, spikemon.i, '.k')
xlabel('Time (ms)')
ylabel('Neuron index')
subplot(122)# 其中的count是spikemon是spikemon的数组，其中每个元素表示一个神经元
plot(G.v0, spikemon.count/duration)
xlabel('v0')
ylabel('Firing rate (sp/s)')
show()
plt.savefig("mul_neuron.png")

结果如下：

在加入了随机噪声后，神经元的发放率与初始电压之间的关系图发生了变化。原本在没有噪声时，发放率从0突变到一个正值，而现在则呈现出S形（sigmoidal）增长。

由于随机噪声的加入，原本尖锐的跳变被平滑的S形曲线所取代。这表明随机噪声在神经元模型中起到了重要作用，它使得神经元的发放行为更加多样化和不可预测。

（5）最后的作业

作业要求：尝试添加一个 StateMonitor 来记录单个神经元的变量值。看瞬时放电率。

完整代码：

from brian2 import *# 设置初始参数
start_scope()
N = 1000  # 神经元数量
tau = 10*ms  # 膜时间常数
vr = -70*mV  # 重置电位
vt0 = -50*mV  # 初始阈值
tau_t = 100*ms  # 阈值恢复时间常数
sigma = 0.5 * (vt0 - vr)  # 噪声强度
v_drive = 2 * (vt0 - vr)  # 驱动电压
duration = 100*ms  # 模拟时长
delta_vt0 = 2*mV  # 放电后阈值的增加量# 神经元动态方程
eqs = '''
dv/dt = (v_drive + vr - v)/tau + sigma*xi*tau**-0.5 : volt
dvt/dt = (vt0 - vt)/tau_t : volt
'''# 重置条件和动作
reset = '''
v = vr
vt += delta_vt0
'''# 创建神经元组
G = NeuronGroup(N, eqs, threshold='v > vt', reset=reset, refractory=5*ms, method='euler')# 初始化神经元的状态
G.v = 'rand() * (vt0 - vr) + vr'
G.vt = vt0# 监测放电事件
spikemon = SpikeMonitor(G)# 添加 StateMonitor 监测一个神经元的状态变量
statemon = StateMonitor(G, ['v', 'vt'], record=0)  # 记录第一个神经元# 运行模拟
run(duration)# 可视化瞬时放电率
_ = hist(spikemon.t/ms, 100, histtype='stepfilled', facecolor='k', weights=list(ones(len(spikemon))/(N*defaultclock.dt)))
xlabel('Time (ms)')
ylabel('Instantaneous firing rate (sp/s)')
show()# 可视化单个神经元的状态变化
plot(statemon.t/ms, statemon.v[0]/mV, label='v (Neuron 0)')
plot(statemon.t/ms, statemon.vt[0]/mV, label='vt (Neuron 0)')
xlabel('Time (ms)')
ylabel('Instantaneous firing rate (sp/s)')
legend()
show()# 保存放电率图像
plt.savefig("mul_neuron.png")

瞬时放电率直方图结果如下：

瞬时放电率可以看作是一个时间窗口内，神经元群体活动的整体强度。
如果多个神经元在某个时间点集中放电，瞬时放电率会高。
通过观察瞬时放电率，我们可以分析神经元群体的同步化、振荡模式等集体动态行为。

加上神经元0的膜电位v和阈值vt

2、突触Synapses

（1）简单突触

建立两个神经元之间的简单突触连接。

from brian2 import *start_scope()#I : 1  无量纲的数值，冒号后的1表示参数单位是1，即没有具体单位
#tau : second   表示时间常数单位是秒
eqs = '''
dv/dt = (I-v)/tau : 1
I : 1
tau : second
'''# 创建了两个神经元，具有相同的微分方程，但是参数I和tau的值不同
G = NeuronGroup(2, eqs, threshold='v>1', reset='v = 0', method='exact')# 其中神经元1的I=0，则表示没有突触影响就不会发放电位
G.I = [2, 0]
G.tau = [10, 100]*ms# Synapses(source, target, action) 
# 源和目的都是神经元组G，突触前神经元产生脉冲会对突触后产生即时的影响，即v_post += 0.2
S = Synapses(G, G, on_pre='v_post += 0.2')# 创建了从神经元0到神经元1的突触，则0发送电位时1的v值增加0.2
S.connect(i=0, j=1)M = StateMonitor(G, 'v', record=True)
run(100*ms)plot(M.t/ms, M.v[0], label='Neuron 0')
plot(M.t/ms, M.v[1], label='Neuron 1')
xlabel('Time (ms)')
ylabel('v')show()plt.savefig("sim_synapse.png")

运行结果：

（2）为突触增加权重

from brian2 import *start_scope()#简单的IF神经元模型，动态方程意味着电压的变化率等于输入电流减去当前电压，除以时间常数
eqs = '''
dv/dt = (I-v)/tau : 1
I : 1
tau : second
'''
G = NeuronGroup(3, eqs, threshold='v>1', reset='v = 0', method='exact')
G.I = [2, 0, 0]
G.tau = [10, 100, 100]*ms# Comment these two lines out to see what happens without Synapses
S = Synapses(G, G, 'w : 1', on_pre='v_post += w')
S.connect(i=0, j=[1, 2])# 权重，其中j是神经元索引，则权重根据索引变化不同
S.w = 'j*0.2'M = StateMonitor(G, 'v', record=True)run(50*ms)plot(M.t/ms, M.v[0], label='Neuron 0')
plot(M.t/ms, M.v[1], label='Neuron 1')
plot(M.t/ms, M.v[2], label='Neuron 2')
xlabel('Time (ms)')
ylabel('v')
legend()
show()plt.savefig("sim_synapse.png")

在突触中定义一个无量纲变量w，用来表示突触连接权重。

其中神经元1的权重是0.2，神经元2的权重是0.4。

（3）引入突触传输延迟

这里的延迟不是不应期，而是传输过程中的脉冲延迟。

from brian2 import *start_scope()#简单的IF神经元模型，动态方程意味着电压的变化率等于输入电流减去当前电压，除以时间常数
eqs = '''
dv/dt = (I-v)/tau : 1
I : 1
tau : second
'''
G = NeuronGroup(3, eqs, threshold='v>1', reset='v = 0', method='exact')
G.I = [2, 0, 0]
G.tau = [10, 100, 100]*msS = Synapses(G, G, 'w : 1', on_pre='v_post += w')
S.connect(i=0, j=[1, 2])
S.delay = 'j*2*ms'# 权重，其中j是神经元索引，则权重根据索引变化不同
S.w = 'j*0.2'M = StateMonitor(G, 'v', record=True)run(50*ms)plot(M.t/ms, M.v[0], label='Neuron 0')
plot(M.t/ms, M.v[1], label='Neuron 1')
plot(M.t/ms, M.v[2], label='Neuron 2')
xlabel('Time (ms)')
ylabel('v')
legend()
show()plt.savefig("sim_synapse.png")

结果如下：

（4）复杂连接

左图是连接图，右图中（x,y）坐标表示神经元x与神经元y之间有连接。

from brian2 import *start_scope()
N = 10
G = NeuronGroup(N, 'v:1')# 突触是自连接，每个神经元可以与其他神经元连接
S = Synapses(G, G)# 确保不会创建自连接，任意两个神经元的连接概率为20%
S.connect(condition='i!=j', p=0.2)#可视化突触组S的函数
def visualise_connectivity(S):Ns = len(S.source)  # 源神经元数量Nt = len(S.target)  # 目标神经元数量figure(figsize=(10, 4))  # 设置画布大小# 左侧子图：以连线形式可视化subplot(121)#  zeros(Ns)：生成长度为Ns的数组, arange(Ns)生成从0到Ns-1的数组，表示源神经元索引# 'ok'表示绘制的样式为黑色圆形，ms=10表示点的大小plot(zeros(Ns), arange(Ns), 'ok', ms=10)  # 源神经元位置plot(ones(Nt), arange(Nt), 'ok', ms=10)  # 目标神经元位置for i, j in zip(S.i, S.j):  # 遍历所有突触连接plot([0, 1], [i, j], '-k')  # 画出[i,j]的连接线xticks([0, 1], ['Source', 'Target'])ylabel('Neuron index')  # 标注神经元索引xlim(-0.1, 1.1)ylim(-1, max(Ns, Nt))subplot(122)plot(S.i, S.j, 'ok')xlim(-1, Ns) #x轴的显示范围ylim(-1, Nt)xlabel('Source neuron index')ylabel('Target neuron index')visualise_connectivity(S)show()plt.savefig("sim_synapse.png")

右侧为左侧的连接关系图。

（5）不同连接概率

from brian2 import *# 定义可视化函数
def visualise_connectivity(S):Ns = len(S.source) Nt = len(S.target) figure(figsize=(10, 4))# 左侧子图：以连线形式可视化subplot(121)# zeros(Ns)：生成长度为Ns的数组, arange(Ns)生成从0到Ns-1的数组，表示源神经元索引# 'ok'表示绘制的样式为黑色圆形，ms=10表示点的大小plot(zeros(Ns), arange(Ns), 'ok', ms=10)  # 源神经元位置plot(ones(Nt), arange(Nt), 'ok', ms=10)  # 目标神经元位置for i, j in zip(S.i, S.j):  # 遍历所有突触连接plot([0, 1], [i, j], '-k')  # 画出[i,j]的连接线xticks([0, 1], ['Source', 'Target'])ylabel('Neuron index')  # 标注神经元索引xlim(-0.1, 1.1)ylim(-1, max(Ns, Nt))subplot(122)plot(S.i, S.j, 'ok')xlim(-1, Ns) #x轴的显示范围ylim(-1, Nt)xlabel('Source neuron index')ylabel('Target neuron index')# 初始化神经元组
start_scope()N = 10
G = NeuronGroup(N, 'v:1')for p in [1.0]:S = Synapses(G, G)S.connect(condition='i!=j', p=p)visualise_connectivity(S)suptitle('p = '+str(p))show()
plt.savefig("sim_synapse.png")

（6）一对一连接

from brian2 import *# 定义可视化函数
def visualise_connectivity(S):Ns = len(S.source) Nt = len(S.target) figure(figsize=(10, 4))# 左侧子图：以连线形式可视化subplot(121)# zeros(Ns)：生成长度为Ns的数组, arange(Ns)生成从0到Ns-1的数组，表示源神经元索引# 'ok'表示绘制的样式为黑色圆形，ms=10表示点的大小plot(zeros(Ns), arange(Ns), 'ok', ms=10)  # 源神经元位置plot(ones(Nt), arange(Nt), 'ok', ms=10)  # 目标神经元位置for i, j in zip(S.i, S.j):  # 遍历所有突触连接plot([0, 1], [i, j], '-k')  # 画出[i,j]的连接线xticks([0, 1], ['Source', 'Target'])ylabel('Neuron index')  # 标注神经元索引xlim(-0.1, 1.1)ylim(-1, max(Ns, Nt))subplot(122)plot(S.i, S.j, 'ok')xlim(-1, Ns) #x轴的显示范围ylim(-1, Nt)xlabel('Source neuron index')ylabel('Target neuron index')# 初始化神经元组
start_scope()N = 10
G = NeuronGroup(N, 'v:1')S = Synapses(G, G)
S.connect(j='i')
visualise_connectivity(S)show()
plt.savefig("sim_synapse.png")

（7）权重图

from brian2 import *start_scope()N = 30
neuron_spacing = 50*umetre #神经元间距
width = N/4.0*neuron_spacingG = NeuronGroup(N, 'x : metre')     #神经元组的每个神经元都有一个状态变量
G.x = 'i*neuron_spacing'    #为每个神经元设置空间位置S = Synapses(G, G, 'w : 1')
S.connect(condition='i!=j')# 设置突触权重 w 的表达式，这里使用高斯函数来计算权重
# 权重随着源神经元和目标神经元之间的空间距离的增加而减小
S.w = 'exp(-(x_pre-x_post)**2/(2*width**2))'# 散点图，S.w*20是调整散点的大小
scatter(S.x_pre/um, S.x_post/um, S.w*20)
xlabel('Source neuron position (um)')
ylabel('Target neuron position (um)')show()
plt.savefig("sim_synapse.png")

其中，索引越相近，权重越大，则散点图大小越大。

（8）复杂突触模型STDP

参考资料：STDP学习 — spikingjelly alpha 文档

突触权重受到突触连接的前神经元(pre)和后神经元(post)的脉冲发放的影响，具体而言是：

如果pre神经元先发放脉冲，post神经元后发放脉冲，则突触的权重会增大；
如果pre神经元后发放脉冲，post神经元先发放脉冲，则突触的权重会减小。

from brian2 import *start_scope()tau_pre = tau_post = 20*ms# 定义前突触可塑性系数，是突触权重增加的幅度
A_pre = 0.01# 定义后突触可塑性系数，是突触权重减少的幅度
A_post = -A_pre*1.05# 时间差数组delta_t，范围从-50ms到50ms，共100个点
delta_t = linspace(-50, 50, 100)*msW = where(delta_t>0, A_pre*exp(-delta_t/tau_pre), A_post*exp(delta_t/tau_post))
plot(delta_t/ms, W)
xlabel(r'$\Delta t$ (ms)')
ylabel('W')
axhline(0, ls='-', c='k')show()
plt.savefig("sim_synapse.png")

后面还有一部分是simulation，运行方法和上述相同，就不再列出来了。