目录
一、基础参考内容
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神经形态硬件仿真器
二、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]*ms
S = 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_spacing
G = 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)*ms
W = 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，运行方法和上述相同，就不再列出来了。
扫一扫
2123
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
