适应性指数整合发放（AdEx）模型
目录
0. 写在前面
1. AdEx 模型的定义
2. AdEx 模型的发放模式
3. 用相平面分析法研究 AdEx 模型产生不同发放模式的动力学机制
3.1. 激发锋发放（tonic spiking）
3.2. 适应（adaptation）
3.3. 簇发放（bursting）
3.4. 瞬时锋发放（transient spiking）
3.5. 延迟发放（delayed spiking）
3.6. 起始簇发放（initial bursting）
4. 分析过程所用到的一系列 BrainPy 代码
0. 写在前面
前面介绍了:
Hodgkin-Huxley Model
简化神经元模型1 – LIF Model
简化神经元模型2 – QIF Model
简化神经元模型3 – ExpIF Model
下面继续介绍简化神经元模型，本篇博客主要是介绍适应性指数整合发放（AdEx）模型，并用BrainPy实现模型的定义、模拟和分析。
QIF 模型和 ExpIF 模型通过引入非线性项使得神经元的膜电位变化与真实情况更加相似，但它们仍然存在和 LIF 模型相同的限制：神经元对膜电位变化的历史没有记忆。因此，它们不能模拟出丰富多样的神经元动力学行为。
下面的情景大家一定都经历过或者听到过：
- Sensory adaptation: ‘ 入芝兰之室，久而不闻其香 ’
- 夜晚，你独自来到海边，海风拂面，闻到海水的咸腥味，忍不住深吸一口气（或捂住鼻子）。但过了一段时间，这种味道逐渐减弱，最后甚至闻不到了。
- 久居鲍鱼之肆，不闻其臭
- 久处闹市，不闻其喧
实际上，这并不是产生这种气味的气体分子在逐渐减少，而是我们长时间处于某种环境后，对该环境的感知逐渐减弱或消失，是一种适应性的表现。单神经元尺度上也存在类似的表现：某些神经元面对恒定的外部刺激时，一开始高频发放，随后发放率逐渐降低，最终稳定在一个较小值，这就是神经元的适应（adaptation）行为。适应行为是有好处的，它使得神经元对新的刺激敏感，对持续重复的刺激脱敏，降低了信息编码的能量消耗，提高了编码效率。
资料参考 ： 《神经计算建模实战》 吴思等 著
1. AdEx 模型的定义
对于上述提到的适应性行为，如何去进行模拟呢？之前我们都是仅使用了一个变量（即膜电位 𝑉）和一个微分方程来刻画神经元的动力学性质，这使得它们无法存储膜电位以外的动力学信息。因此我们就会自然而然思考是否需要再引入一些其余的变量呢？Brette 和 Gerstner 于 2005 年提出了适应性指数整合发放模型（Adaptive Exponential Integrate-and-Fire model, AdEx model），在已有的整合发放模型上增加了一个变量 𝑤 来描述神经元的适应性：
AdEx 模型表示为 ：
τ
m
d
V
d
t
=
−
(
V
−
V
rest
)
+
Δ
T
e
V
−
V
T
Δ
T
−
R
w
+
R
I
(
t
)
(1)
\tau_m \frac{dV}{dt} = - (V - V_{\text{rest}}) + \Delta_T e^{\frac{V - V_T}{\Delta_T}} - R_w + R I(t)\tag{1}
τmdtdV=−(V−Vrest)+ΔTeΔTV−VT−Rw+RI(t)(1)
τ
w
d
w
d
t
=
a
(
V
−
V
rest
)
−
w
+
b
τ
w
∑
t
(
f
)
δ
(
t
−
t
(
f
)
)
(2)
\tau_w \frac{dw}{dt} = a (V - V_{\text{rest}}) - w + b \tau_w \sum_{t^{(f)}} \delta(t - t^{(f)})\tag{2}
τwdtdw=a(V−Vrest)−w+bτwt(f)∑δ(t−t(f))(2)
if
V
>
θ
,
V
←
V
reset
for last
t
ref
(3)
\text{if } V > \theta, \quad V \leftarrow V_{\text{reset}} \text{ for last } t_{\text{ref}}\tag{3}
if V>θ,V←Vreset for last tref(3)式中，
w
w
w 被称为适应变量（adaptation variable）。方程 (1) 与 ExpIF 模型的动力学方程 (简化神经元模型3 – ExpIF Model) 相似，但右侧增加了一项
−
R
w
-Rw
−Rw 以表示膜电位受适应变量
w
w
w 的调节。
w
w
w 越大，神经元的适应性越强，膜电位变化率
d
V
/
d
t
dV/dt
dV/dt 越小。方程 (2) 刻画了
w
w
w 的变化，其主要由三部分组成：
- 第一项 a ( V − V rest ) a(V - V_{\text{rest}}) a(V−Vrest) 表示膜电位 V V V 对 w w w 起正向调控，即膜电位越高，神经元适应性增加得越快；
- 第二项 − w -w −w 为衰减项；
- 第三项描述了动作电位对 w w w 的影响，其中 t ( f ) t^{(f)} t(f) 表示神经元动作电位发放的时刻，当 t = t ( f ) t = t^{(f)} t=t(f) 时， δ \delta δ 函数积分为 1，因此 w w w 会瞬时增加 b b b（此项中的 τ w \tau_w τw 和等式左侧的 τ w \tau_w τw 相抵消）。这一项表明，每当神经元产生一个动作电位，神经元的适应性都会增加一个固定值。
AdEx 模型同样也需要通过条件判断式 (3) 来显式地重置发放动作电位后的膜电位。
根据上述定义的 AdEx ，可以使用 BrainPy 进行模型的搭建：
import brainpy as bp
import brainpy.math as bm
import matplotlib.pyplot as plt
import numpy as np
plt.rcParams.update({"font.size": 15})
plt.rcParams['font.sans-serif'] = ['Times New Roman']
# 适应性指数整合发放模型
class AdEx(bp.dyn.NeuDyn):
def __init__(self, size, V_rest=-65., V_reset=-68., V_th=20., V_T=-60., delta_T=1., a=1.,
b=2.5, R=1., tau=10., tau_w=30., tau_ref=0., name=None):
# 初始化父类
super(AdEx, self).__init__(size=size, name=name)
# 初始化参数
self.V_rest = V_rest
self.V_reset = V_reset
self.V_th = V_th
self.V_T = V_T
self.delta_T = delta_T
self.a = a
self.b = b
self.tau = tau
self.tau_w = tau_w
self.R = R
self.tau_ref = tau_ref
# 初始化变量
self.V = bm.Variable(bm.random.randn(self.num) + V_reset)
self.w = bm.Variable(bm.zeros(self.num))
self.input = bm.Variable(bm.zeros(self.num))
self.spike = bm.Variable(bm.zeros(self.num, dtype=bool)) # 脉冲发放状态
self.t_last_spike = bm.Variable(bm.ones(self.num) * -1e7) # 上一次脉冲发放时间
self.refractory = bm.Variable(bm.zeros(self.num, dtype=bool)) # 是否处于不应期
# 定义积分器
self.integral = bp.odeint(f=self.derivative, method='exp_auto')
def dV(self, V, t, w, Iext):
_tmp = self.delta_T * bm.exp((V - self.V_T) / self.delta_T)
dVdt = (- V + self.V_rest + _tmp - self.R * w + self.R * Iext) / self.tau
return dVdt
def dw(self, w, t, V):
dwdt = (self.a * (V - self.V_rest) - w) / self.tau_w
return dwdt
# 将两个微分方程联合为一个，以便同时积分
@property
def derivative(self):
return bp.JointEq([self.dV, self.dw])
def update(self, ):
_t, _dt = bp.share['t'], bp.share['dt']
# 以数组的方式对神经元进行更新
V, w = self.integral(self.V, self.w, _t, self.input, dt=_dt) # 更新膜电位V和权重值w
refractory = (_t - self.t_last_spike) <= self.tau_ref # 判断神经元是否处于不应期
V = bm.where(refractory, self.V, V) # 若处于不应期，则返回原始膜电位self.V，否则返回更新后的膜电位V
w = bm.where(refractory, self.w, w) # w同理
spike = V > self.V_th # 将大于阈值的神经元标记为发放了脉冲
self.spike.value = spike # 更新神经元脉冲发放状态
self.t_last_spike.value = bm.where(spike, _t, self.t_last_spike) # 更新最后一次脉冲发放时间
self.V.value = bm.where(spike, self.V_reset, V) # 将发放了脉冲的神经元膜电位置为V_reset，其余不变
self.w.value = bm.where(spike, w + self.b, w) # 发放了脉冲的神经元 w = w + b
self.refractory.value = bm.logical_or(refractory, spike) # 更新神经元是否处于不应期
self.input[:] = 0. # 重置外界输入
def reset_state(self, batch_size=None):
self.V[:] = self.V_reset
self.w[:] = 0.
self.t_last_spike[:] = -1e7
self.refractory.value = bm.zeros_like(self.refractory)
self.input[:] = 0.
self.spike[:] = False
# self.spike.value=False 好像在最新的BrainPy版本中需要写成这样
在定义好上述模型之后，编写如下简单的可视化代码怪观察其在恒定刺激的情况下产生的行为：
def compare_Adex_models():
def simple_run_Adex():
neu = AdEx(1)
runner = bp.DSRunner(
neu, monitors=['V', 'w', 'spike'], inputs=('input', 9.), dt=0.01)
runner(500)
runner.mon.V = np.where(runner.mon.spike, 20., runner.mon.V)
return runner.mon.ts, runner.mon.V, runner.mon.w
def simple_run_Adex_adaptive():
# 调整参数以模拟发放频率逐渐降低，最后停止发放
neu = AdEx(1, a=1.6, b=0.4, tau_w=110.0) # 这里进行调参
runner = bp.DSRunner(
neu, monitors=['V', 'w', 'spike'], inputs=('input', 9.), dt=0.01)
runner(500)
runner.mon.V = np.where(runner.mon.spike, 20., runner.mon.V)
return runner.mon.ts, runner.mon.V, runner.mon.w
ts1, V1, w1 = simple_run_Adex()
ts2, V2, w2 = simple_run_Adex_adaptive()
fig, axes = plt.subplots(1, 2, figsize=(12, 8))
# 左边的子图
axes[0].plot(ts1, V1, label='V')
axes[0].plot(ts1, w1, label='w')
axes[0].set_ylabel('$V$ (mV)')
axes[0].set_xlabel('$t$ (ms)')
axes[0].set_title('Original AdEx')
axes[0].spines['top'].set_visible(False)
axes[0].spines['right'].set_visible(False)
axes[1].plot(ts2, V2, label='V')
axes[1].plot(ts2, w2, label='w')
axes[1].set_ylabel('$w$ (mA)')
axes[1].set_xlabel('$t$ (ms)')
axes[1].spines['top'].set_visible(False)
axes[1].spines['right'].set_visible(False)
plt.tight_layout()
plt.show()
上图左侧我们给定了一个恒定的刺激，发现神经元在恒定刺激下发放频率逐渐降低，最终趋于一个恒定值；右侧我们改变了a=1.6, b=0.4, tau_w=110.0。发现神经元发放频率逐渐降低，最终停止发放，不再兴奋。
实际上，AdEx模型不仅能模拟神经元的适应行为，还能模拟其他多种发放模式。下面具体讨论一下各类发放模式以及形成它们的动力学原因。
2. AdEx 模型的发放模式
我们知道：不同脑区、不同类型的神经元具有不同的发放模式。根据神经元在恒定刺激下产生动作电位的时间间隔的不同，可以将神经元的发放分为四种模式：
- 激发锋发放（tonic spiking）：也叫规律发放（regular spiking），动作电位之间的时间间隔相同；
- 适应性发放（adapting）：动作电位之间的时间间隔逐渐变大，最后趋于一个恒定值；
- 簇发放（bursting）：神经元在短时间内连续产生多个脉冲，然后在较长时间间隔内不再发放，并如此周期性重复；
- 不规律发放（irregular spiking）：动作电位之间的时间间隔不具有明显的规律。
上述分类是基于神经元稳定发放的状态而决定的。有的神经元在刚接受刺激时具有独特的发放模式，这是因为神经元的起始状态和接受刺激后的稳定状态差别较大，从起始状态转变为稳定状态需要一定时间。神经元起始阶段的发放模式又可以被分为三类：
- 典型发放（classic/tonic spiking）：起始阶段和稳定阶段的发放没有差别；
- 起始簇发放（initial burst）：起始阶段出现短时间的连续发放，发放频率远高于稳定状态；
- 延迟发放（delayed spiking）：神经元在经过一段较长时间的阈下变化后才开始发放。
理论上，神经元起始状态的分类和稳定状态的分类可以两两组合，形成 12 种发放模式。所有发放模式的组合都能在电生理实验中被观察到。有趣的是，这些发放模式都能够被 AdEx 模型重现。表3.1中列出了几种典型的发放模式对应的参数，其余参数为： V rest = − 70 V_{\text{rest}} = -70 Vrest=−70， θ = 0 \theta = 0 θ=0， V T = − 50 V_T = -50 VT=−50， Δ T = 2 \Delta_T = 2 ΔT=2， R = 0.5 R = 0.5 R=0.5。它们的模拟结果展示在图3.15中。
3. 用相平面分析法研究 AdEx 模型产生不同发放模式的动力学机制
为什么 AdEx 模型能够模拟出如此多样的发放模式？我们可以通过相平面分析的方法来考察不同参数对d𝑉/d𝑡 和 d𝑤/d𝑡 的影响。
def _ppa2d(group, title, v_range=None, w_range=None, Iext=65.,
duration=400, num_text_sp=0, sp_text=None, extra_fun=None):
v_range = [-70., -40.] if not v_range else v_range
w_range = [-10., 50.] if not w_range else w_range
if sp_text is None:
sp_text = dict()
fig, gs = bp.visualize.get_figure(1, 1, 4.5, 6)
ax = fig.add_subplot(gs[0, 0])
# 使用BrainPy中的相平面分析工具
phase_plane_analyzer = bp.analysis.PhasePlane2D(
model=group,
target_vars={'V': v_range, 'w': w_range, }, # 待分析变量
pars_update={'Iext': Iext}, # 需要更新的变量
resolutions=0.05 # 精度范围
)
# 画出V, w的零增长曲线
phase_plane_analyzer.plot_nullcline()
# 画出奇点
phase_plane_analyzer.plot_fixed_point()
# 画出向量场
phase_plane_analyzer.plot_vector_field(with_return=True)
# 分段画出V, w的变化轨迹
group.V[:] = group.V_reset
runner = bp.DSRunner(
group, monitors=['V', 'w', 'spike'], inputs=('input', Iext))
runner(duration)
spike = runner.mon.spike.squeeze()
s_idx = np.where(spike)[0] # 找到所有发放动作电位对应的index
s_idx = np.concatenate(([0], s_idx, [len(spike) - 1])) # 加上起始点和终止点的index
for i in range(len(s_idx) - 1):
vs = runner.mon.V[s_idx[i]: s_idx[i + 1]]
ws = runner.mon.w[s_idx[i]: s_idx[i + 1]]
plt.plot(vs, ws, color='darkslateblue')
if i < num_text_sp:
plt.text(group.V_reset - 1, ws[0] - 0.5, sp_text.get(i, str(i)))
# 画出虚线 x = V_reset
plt.plot([group.V_reset, group.V_reset],
w_range, '--', color='grey', zorder=-1)
if extra_fun:
extra_fun()
ax.get_legend().remove()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.show()
定义好分析方法后，依次带入每种发放模式对应的参数，画出想平面分析结果。
ppa2d(AdEx(1, tau=20., a=0., tau_w=30., b=60., V_reset=-55.),
title='Tonic Spiking', w_range=[-5, 75.])
ppa2d(AdEx(1, tau=20., a=0., tau_w=100., b=5., V_reset=-55.),
title='Adaptation', w_range=[-5, 45.])
ppa2d(AdEx(1, tau=5., a=0.5, tau_w=100., b=7., V_reset=-51.),
title='Initial Bursting', w_range=[-5, 50.])
ppa2d(AdEx(1, tau=5., a=-0.5, tau_w=100., b=7., V_reset=-47.),
title='Bursting', w_range=[-5, 60.])
ppa2d(AdEx(1, tau=10., a=1., tau_w=100., b=10., V_reset=-60.),
title='Transient Spiking', w_range=[-5, 60.], Iext=55.)
ppa2d(AdEx(1, tau=5., a=-1., tau_w=100., b=5., V_reset=-60.),
title='Delayed Spiking', w_range=[-30, 20.], Iext=25.)
下面依次分析每种发放模式背后的动力学原因。笔者将膜电位随时间变化的 𝑉/𝑤 −𝑡 图和相平面分析图放在一起，以便读者对比观察。同时，为保证每次模拟的结果相同，笔者将 AdEx 模型的起始膜电位均设置为 0。
3.1. 激发锋发放（tonic spiking）
我们首先研究激发锋发放（图3.16）。图3.16（b）标注了 V V V 和 w w w 的零增长等值线（nullcline），在零增长等值线上变量对时间的导数为 0。在相平面图中我们人为标记了 V − w V-w V−w 轨迹的次序，起始点记为 0，第一次发放动作电位后的重置点记为 1，以此类推。第一次发放后，神经元膜电位重置为 V reset V_{\text{reset}} Vreset（图中虚线表示 V = V reset V = V_{\text{reset}} V=Vreset），同时变量 w w w 升高了 b b b。轨迹点沿着向量场移动，完成了第二次发放，然后从 2 开始进行第三次发放。从图3.16（右）中可以发现，第三次发放后的重置点 3 与 2 重合，因此 2 → 3 形成了一个周期，轨迹点会沿着相同的路径周期性运动，从而形成了有规律的激发锋发放模式。
从相平面图中我们还可以发现，每次 V V V 和 w w w 被重置后，轨迹点都要经历一个迂回的路线才能发放动作电位，这在 V / w − t V/w - t V/w−t 图中表现为膜电位先下降后上升。这是因为激发锋发放模式中的 b b b 取值较大（本例中 b = 60 b = 60 b=60），因此发放动作电位后的重置点能够落在 V V V-零增长等值线以上。由于在 V V V-零增长等值线以上 d V / d t < 0 dV/dt < 0 dV/dt<0（向量场箭头向左），在 V V V-零增长等值线以下 d V / d t > 0 dV/dt > 0 dV/dt>0（向量场箭头向右），所以轨迹点在 V V V-零增长等值线以上时朝膜电位减小的方向走，越过 V V V-零增长等值线后膜电位才开始增长，这使得每次发放之间的时间间隔较大，神经元不能在短时间多次发放。在后文中我们将会看到，并不是所有轨迹都会经历这样的迂回。
3.2. 适应（adaptation）
神经元适应性发放的 V / w − t V/w - t V/w−t 图及相平面分析图如图3.17所示。每经历一次发放后， w w w 的重置点都比上一次重置点的位置高，轨迹点从起始到发放经历的时间会增加，因此表现出发放的时间间隔增大。向量场使得越靠上的轨迹在 w w w 轴上下降得越多，当某一次重置的增量 b b b 和上一次轨迹中 w w w 的减少量相同，即本次重置点和上一次重合时， w w w 不再增加，神经元的发放间隔变为一个恒定值。这就是 AdEx 模型产生适应行为的动力学原因。
相比激发锋发放，适应模式最大的不同在于轨迹点不会越过 V V V 的零增长等值线，这是因为模型中的参数 b b b 取值很小（ b = 5 b = 5 b=5），重置时 w w w 的增量很小。因此，适应行为的神经元不会在阈下经历 V V V 先减小后增大的过程，而是直接从重置点 V reset V_{\text{reset}} Vreset 上升直至产生动作电位。
3.3. 簇发放（bursting）
我们跳过起始簇发放，先考察簇发放的动力学特性。根据图3.18（a）中 V V V 和 w w w 随时间的变化，除了第一次簇发放以外，之后每次簇发放的形状都比较相似，因此我们从第一次簇发放的结尾点开始标记 V − w V-w V−w 轨迹的次序（图3.18（b））。从 0 开始，膜电位在经历了一次迂回后达到发放阈值，随后轨迹点被重置到点 1，经历第二次发放后被重置到点 2。由于点 1 和点 2 都在 V V V-零增长等值线的下方， d V d t > 0 \frac{dV}{dt} > 0 dtdV>0，因此膜电位变化和适应性发放一样直接上升，不经历迂回；又因为 V > V T V > V_T V>VT 时微分方程中的指数项很大，膜电位上升极快，所以膜电位多次快速达到阈值并重置，在 V / w − t V/w - t V/w−t 图中表现为短时间内的集中发放。而第三次发放后的重置点 3 在 V V V-零增长等值线的上方，所以膜电位又会经历一次迂回，这就形成了两个发放簇之间较长的时间间隔。
形成簇发放的必要条件是每次经历迂回路径的时候 w w w 的下降值 Δ w \Delta w Δw 需满足 Δ w ≥ 2 b \Delta w \geq 2b Δw≥2b，这样才能保证两次迂回路径（即长时间间隔）之间能够有至少两次连续的快速发放。本例中 Δ w = 3 b \Delta w = 3b Δw=3b，所以每个发放簇有三个锋。
3.4. 瞬时锋发放（transient spiking）
瞬时锋发放的特点是在恒定刺激下，神经元先快速发放多个动作电位，然后停止发放，膜电位趋于一个定值。膜电位存在稳定值要求相平面图中存在稳定点。相平面图（图3.19（b））也证明该动力学系统确实存在一个稳定焦点。此外，系统还存在一个鞍点，这使得轨迹点不会总是被稳定焦点所吸引。最开始时 w w w 较小，起始点和稳定焦点的距离较远，因此轨迹不会被稳定焦点吸引，神经元的膜电位会不断升高并产生动作电位。如此反复两次后，重置点（图中 2）掉入了稳定焦点的吸引域，轨迹被稳定焦点所吸引，因此 V V V、 w w w 变量都将趋于稳定，神经元将不再产生动作电位。
根据 AdEx 模型的微分方程 (1)，令左侧
d
V
d
t
=
0
\frac{dV}{dt} = 0
dtdV=0，可以求得
V
V
V-零增长等值线为：
f
(
V
)
=
−
R
w
+
R
I
(
t
)
(4)
f(V) = -Rw + RI(t) \tag{4}
f(V)=−Rw+RI(t)(4)其中
f
(
V
)
=
V
−
V
rest
−
Δ
T
exp
(
V
−
V
T
Δ
T
)
f(V) = V - V_{\text{rest}} - \Delta_T \exp\left(\frac{V - V_T}{\Delta_T}\right)
f(V)=V−Vrest−ΔTexp(ΔTV−VT)。这是一个隐式方程，我们不必求出
V
V
V 关于
w
w
w 的显函数就能够发现，若要
V
V
V 保持不变，
I
I
I 增大会使得
w
w
w 也增大，这在相平面图中表现为外部电流增大会使
V
V
V-零增长等值线上移。因此，我们可以通过调节外部输入
I
I
I 的大小来改变
V
V
V-零增长等值线的高低，从而影响神经元的发放模式。
I
I
I 太大时两条零增长等值线不相交，即不存在稳定焦点；
I
I
I 太小时稳定焦点和鞍点距离太远，稳定焦点的吸引域在相空间覆盖面积过大，轨迹点可能直接被稳定焦点吸引而不产生最初的动作电位发放，也不容易形成瞬时锋发放。只有当
I
I
I 取值适中时才能够形成瞬时锋的发放模式。
3.5. 延迟发放（delayed spiking）
延迟发放的特点是膜电位在第一个动作电位之前会经历一个较长时间的阈下变化（图3.20（a））。从相平面分析图（图3.20（b））中我们可以发现，第 0 条轨迹的中间部分和 V V V-零增长等值线非常靠近，这对应着 V / w − t V/w - t V/w−t 图中第一次发放前膜电位变化的中间部分在较长时间内增长缓慢。为什么靠近 V V V-零增长等值线的轨迹运动缓慢呢？回顾 V V V-零增长等值线的定义，线上各点 d V d t = 0 \frac{dV}{dt} = 0 dtdV=0，即当轨迹经过 V V V-零增长等值线时不会有水平方向而只会有竖直方向的运动；靠近 V V V-零增长等值线的轨迹点存在水平运动，但速度十分缓慢。图3.20（右）中第 0 条轨迹有一大段都与 V V V-零增长等值线十分靠近，水平方向的速度分量很小决定了整体速度也很小，由此形成了第一次发放前较长时间的延迟。
要形成延迟发放，不仅要求起始点比之后的重置点距离 V V V-零增长等值线更近，还要求 V V V-零增长等值线附近向量场的竖直分量较小，否则轨迹点将不会紧贴 V V V-零增长等值线运动。若要向量的竖直分量较小，即 d w d t \frac{dw}{dt} dtdw 较小，可以让 τ w / τ \tau_w / \tau τw/τ 足够大，即让 w w w 的动力学变化速率显著小于 V V V 的动力学变化速率。对比激发锋发放的参数设置和相平面图（图3.16（b）），可以印证 τ w / τ \tau_w / \tau τw/τ 对向量场竖直分量的影响。
3.6. 起始簇发放（initial bursting）
对于起始簇发放（initial bursting），分析方法和前面的3.3. 簇发放（bursting）分析一致，这里仅给出了它的发放模式图和相平面分析图，但不再详细分析产生它的动力学原因。读者可以结合两幅图分析每一条轨迹对应图中哪一段曲线，它们为何形成了起始簇发放，以及在参数上需要满足什么条件。
4. 分析过程所用到的一系列 BrainPy 代码
在上述过程中，我们主要是以分析 AdEx Model 为主，相关的代码给的比较的分散，不便于读者研读代码，因此，笔者在此处将上述所有用到的代码放置在下面，供读者研读，部分结果可能与上述分析结果不一致，但是这些可以通过调参的方式进行实现，感兴趣的读者可以在研读代码后进行调参。代码部分也有一些扩展的内容，读者也可以自行学习，以便更好的理解 AdEx Model 。
完整代码如下：
import brainpy as bp
import brainpy.math as bm
import matplotlib.pyplot as plt
import numpy as np
plt.rcParams.update({"font.size": 15})
plt.rcParams['font.sans-serif'] = ['Times New Roman']
bp.math.enable_x64()
# 适应性指数整合发放模型
class AdEx(bp.dyn.NeuDyn):
def __init__(self, size, V_rest=-65., V_reset=-68., V_th=20., V_T=-60., delta_T=1., a=1.,
b=2.5, R=1., tau=10., tau_w=30., tau_ref=0., name=None):
# 初始化父类
super(AdEx, self).__init__(size=size, name=name)
# 初始化参数
self.V_rest = V_rest
self.V_reset = V_reset
self.V_th = V_th
self.V_T = V_T
self.delta_T = delta_T
self.a = a
self.b = b
self.tau = tau
self.tau_w = tau_w
self.R = R
self.tau_ref = tau_ref
# 初始化变量
self.V = bm.Variable(bm.random.randn(self.num) + V_reset)
self.w = bm.Variable(bm.zeros(self.num))
self.input = bm.Variable(bm.zeros(self.num))
self.spike = bm.Variable(bm.zeros(self.num, dtype=bool)) # 脉冲发放状态
self.t_last_spike = bm.Variable(bm.ones(self.num) * -1e7) # 上一次脉冲发放时间
self.refractory = bm.Variable(
bm.zeros(self.num, dtype=bool)) # 是否处于不应期
# 定义积分器
self.integral = bp.odeint(f=self.derivative, method='exp_auto')
def dV(self, V, t, w, Iext):
_tmp = self.delta_T * bm.exp((V - self.V_T) / self.delta_T)
dVdt = (- V + self.V_rest + _tmp - self.R *
w + self.R * Iext) / self.tau
return dVdt
def dw(self, w, t, V):
dwdt = (self.a * (V - self.V_rest) - w) / self.tau_w
return dwdt
# 将两个微分方程联合为一个，以便同时积分
@property
def derivative(self):
return bp.JointEq([self.dV, self.dw])
def update(self, ):
_t, _dt = bp.share['t'], bp.share['dt']
# 以数组的方式对神经元进行更新
V, w = self.integral(self.V, self.w, _t,
self.input, dt=_dt) # 更新膜电位V和权重值w
refractory = (_t - self.t_last_spike) <= self.tau_ref # 判断神经元是否处于不应期
# 若处于不应期，则返回原始膜电位self.V，否则返回更新后的膜电位V
V = bm.where(refractory, self.V, V)
w = bm.where(refractory, self.w, w) # w同理
spike = V > self.V_th # 将大于阈值的神经元标记为发放了脉冲
self.spike.value = spike # 更新神经元脉冲发放状态
self.t_last_spike.value = bm.where(
spike, _t, self.t_last_spike) # 更新最后一次脉冲发放时间
# 将发放了脉冲的神经元膜电位置为V_reset，其余不变
self.V.value = bm.where(spike, self.V_reset, V)
self.w.value = bm.where(spike, w + self.b, w) # 发放了脉冲的神经元 w = w + b
self.refractory.value = bm.logical_or(
refractory, spike) # 更新神经元是否处于不应期
self.input[:] = 0. # 重置外界输入
def reset_state(self, batch_size=None):
self.V[:] = self.V_reset
self.w[:] = 0.
self.t_last_spike[:] = -1e7
self.refractory.value = bm.zeros_like(self.refractory)
self.input[:] = 0.
self.spike[:] = False
# self.spike.value=False 好像在最新的BrainPy版本中需要写成这样
def compare_Adex_models():
def simple_run_Adex():
neu = AdEx(1)
runner = bp.DSRunner(
neu, monitors=['V', 'w', 'spike'], inputs=('input', 9.), dt=0.01)
runner(500)
runner.mon.V = np.where(runner.mon.spike, 20., runner.mon.V)
return runner.mon.ts, runner.mon.V, runner.mon.w
def simple_run_Adex_adaptive():
# 调整参数以模拟发放频率逐渐降低，最后停止发放
neu = AdEx(1, a=1.6, b=0.4, tau_w=110.0) # 这里进行调参
runner = bp.DSRunner(
neu, monitors=['V', 'w', 'spike'], inputs=('input', 9.), dt=0.01)
runner(500)
runner.mon.V = np.where(runner.mon.spike, 20., runner.mon.V)
return runner.mon.ts, runner.mon.V, runner.mon.w
ts1, V1, w1 = simple_run_Adex()
ts2, V2, w2 = simple_run_Adex_adaptive()
fig, axes = plt.subplots(1, 2, figsize=(12, 8))
# 左边的子图
axes[0].plot(ts1, V1, label='V')
axes[0].plot(ts1, w1, label='w')
axes[0].set_ylabel('$V$ (mV)')
axes[0].set_xlabel('$t$ (ms)')
axes[0].set_title('Original AdEx')
axes[0].spines['top'].set_visible(False)
axes[0].spines['right'].set_visible(False)
axes[1].plot(ts2, V2, label='V')
axes[1].plot(ts2, w2, label='w')
axes[1].set_ylabel('$w$ (mA)')
axes[1].set_xlabel('$t$ (ms)')
axes[1].spines['top'].set_visible(False)
axes[1].spines['right'].set_visible(False)
plt.tight_layout()
plt.show()
def simple_run_Adex():
neu = AdEx(200)
runner = bp.DSRunner(
neu, monitors=['V', 'w', 'spike'], inputs=('input', 9.), dt=0.01)
runner(500)
runner.mon.V = np.where(runner.mon.spike, 20., runner.mon.V)
plt.plot(runner.mon.ts, runner.mon.V, label='V')
plt.plot(runner.mon.ts, runner.mon.w, label='w')
plt.xlabel('$t$ (ms)')
plt.ylabel('$V$ (mV)')
plt.show()
def run():
for epoch in range(1):
print(f'Epoch {epoch}')
if epoch > 98:
simple_run_Adex()
def run_AdEx_model():
# 运行AdEx模型
neu = AdEx(2)
neu.V[:] = bm.asarray([-68.79061, -66.51926])
runner = bp.DSRunner(
neu, monitors=['V', 'w', 'spike'], inputs=('input', 9.), dt=0.01)
runner(400)
# 可视化V和w的变化
runner.mon.V = np.where(runner.mon.spike, 20., runner.mon.V)
for i in range(2):
fig, gs = bp.visualize.get_figure(2, 1, 2.25, 6)
ax = fig.add_subplot(gs[0, 0])
plt.plot(runner.mon.ts, runner.mon.V[:, i], label='V')
plt.ylabel(r'$V$ (mV)')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax = fig.add_subplot(gs[1, 0])
plt.plot(runner.mon.ts, runner.mon.w[:, i], label='w')
plt.ylabel(r'$w$ (mA)')
plt.xlabel(r'$t$ (ms)')
plt.xlim(-30, 430)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.show()
def AdEx_patterns():
group = AdEx(
size=6,
a=bm.asarray([0., 0., 0.5, -0.5, 1., -1.]),
b=bm.asarray([60., 5., 7., 7., 10., 5.]),
tau=bm.asarray([20., 20., 5., 5., 10., 5.]),
tau_w=bm.asarray([30., 100., 100., 100., 100., 100.]),
V_reset=bm.asarray([-55., -55., -51., -47., -60., -60.]),
R=.5, delta_T=2., V_rest=-70, V_th=-30, V_T=-50
)
group.V.value = group.V_reset
par_I = bm.asarray([65., 65., 65., 65., 55., 25.])
runner = bp.DSRunner(
group, monitors=['V', 'w', 'spike'], inputs=('input', par_I))
runner.run(500.)
runner.mon.V = np.where(runner.mon.spike, 20., runner.mon.V)
names = ['Tonic', 'Adapting', 'Init Bursting',
'Bursting', 'Transient', 'Delayed']
for i_col in range(2):
for i_row in range(3):
i = i_col * 3 + i_row
fig, gs = bp.visualize.get_figure(2, 1, 2.25, 6)
ax = fig.add_subplot(gs[0, 0])
plt.plot(runner.mon.ts, runner.mon.V[:, i], label='V')
plt.ylabel(r'$V$ (mV)')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax = fig.add_subplot(gs[1, 0])
plt.plot(runner.mon.ts, runner.mon.w[:, i], label='w')
plt.ylabel(r'$w$ (mA)')
plt.xlabel(r'$t$ (ms)')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.show()
def _ppa2d(group, title, v_range=None, w_range=None, Iext=65.,
duration=400, num_text_sp=0, sp_text=None, extra_fun=None):
v_range = [-70., -40.] if not v_range else v_range
w_range = [-10., 50.] if not w_range else w_range
if sp_text is None:
sp_text = dict()
fig, gs = bp.visualize.get_figure(1, 1, 4.5, 6)
ax = fig.add_subplot(gs[0, 0])
# 使用BrainPy中的相平面分析工具
phase_plane_analyzer = bp.analysis.PhasePlane2D(
model=group,
target_vars={'V': v_range, 'w': w_range, }, # 待分析变量
pars_update={'Iext': Iext}, # 需要更新的变量
resolutions=0.05 # 精度范围
)
# 画出V, w的零增长曲线
phase_plane_analyzer.plot_nullcline()
# 画出奇点
phase_plane_analyzer.plot_fixed_point()
# 画出向量场
phase_plane_analyzer.plot_vector_field(with_return=True)
# 分段画出V, w的变化轨迹
group.V[:] = group.V_reset
runner = bp.DSRunner(
group, monitors=['V', 'w', 'spike'], inputs=('input', Iext))
runner(duration)
spike = runner.mon.spike.squeeze()
s_idx = np.where(spike)[0] # 找到所有发放动作电位对应的index
s_idx = np.concatenate(([0], s_idx, [len(spike) - 1])) # 加上起始点和终止点的index
for i in range(len(s_idx) - 1):
vs = runner.mon.V[s_idx[i]: s_idx[i + 1]]
ws = runner.mon.w[s_idx[i]: s_idx[i + 1]]
plt.plot(vs, ws, color='darkslateblue')
if i < num_text_sp:
plt.text(group.V_reset - 1, ws[0] - 0.5, sp_text.get(i, str(i)))
# 画出虚线 x = V_reset
plt.plot([group.V_reset, group.V_reset],
w_range, '--', color='grey', zorder=-1)
if extra_fun:
extra_fun()
ax.get_legend().remove()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.show()
def _vt_plot(neu, title, input=('input', 65.), duration=400):
runner = bp.DSRunner(neu, monitors=['V', 'w', 'spike'], inputs=input)
runner(duration)
fig, gs = bp.visualize.get_figure(2, 1, 2.25, 6)
ax = fig.add_subplot(gs[0, 0])
runner.mon.V = np.where(runner.mon.spike, 0., runner.mon.V)
ax.plot(runner.mon.ts, runner.mon.V, label='V')
plt.ylabel(r'$V$ (mV)')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax = fig.add_subplot(gs[1, 0])
ax.plot(runner.mon.ts, runner.mon.w, label='w', color=u'#ff7f0e')
ax.set_xlabel(r'$t$ (ms)')
plt.ylabel(r'$w$ (mA)')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
# plt.show()
fig.align_ylabels()
def AdEx_phase_plane_analysis():
bp.math.enable_x64()
# Tonic Spiking
model = AdEx(1, tau=20., a=0., tau_w=30., b=60., V_reset=-55.,
V_rest=-70., V_th=0., V_T=-50., R=0.5, delta_T=2.)
def f():
plt.plot(np.linspace(-70, -40, 500), np.zeros(500),
'.', color='lightcoral', alpha=.7)
plt.text(-54.5, 0, '0')
plt.text(-54.5, 60, '1')
plt.text(-54.5, 69.5, '2, 3, ...')
plt.text(-68., 62.8, 'V nullcline')
plt.text(-68., 2, 'w nullcline')
_ppa2d(model, title='Tonic Spiking', w_range=[-5, 75.], extra_fun=f)
model.reset_state()
_vt_plot(model, title='Tonic Spiking')
# Adaptation
model = AdEx(1, tau=20., a=0., tau_w=100., b=5., V_reset=-55.,
V_rest=-70., V_th=0., V_T=-50., R=0.5, delta_T=2.)
def f():
plt.plot(np.linspace(-70, -40, 500), np.zeros(500),
'.', color='lightcoral', alpha=.7)
plt.text(-47.7, 30.5, 'V nullcline')
plt.text(-68., 2., 'w nullcline')
plt.text(-50.7, 19., 'Trajectory')
_ppa2d(model, title='Adaptation',
w_range=[-5, 45.], extra_fun=f, num_text_sp=5)
model.reset_state()
_vt_plot(model, title='Adaptation')
# Initial Bursting
model = AdEx(1, tau=5., a=0.5, tau_w=100., b=7., V_reset=-51.,
V_rest=-70., V_th=0., V_T=-50., R=0.5, delta_T=2.)
def f():
plt.text(-59, 45.3, 'V nullcline')
plt.text(-68., 7., 'w nullcline')
plt.text(-48.1, 27.3, 'Trajectory')
_ppa2d(model, title='Initial Bursting',
w_range=[-5, 50.], extra_fun=f, num_text_sp=6)
model.reset_state()
_vt_plot(model, title='Initial Bursting')
# Bursting
model = AdEx(1, tau=5., a=-0.5, tau_w=100., b=7., V_reset=-47.,
V_rest=-70., V_th=0., V_T=-50., R=0.5, delta_T=2.)
def f():
plt.text(-48, 29.5, '1')
plt.text(-48, 35.5, '2')
plt.text(-46.8, 38.8, '0')
plt.text(-46.8, 42, '3')
plt.text(-60.8, 48.83, 'V nullcline')
plt.text(-68., 1.2, 'w nullcline')
plt.text(-45.88, 36.8, 'Trajectory')
_ppa2d(model, title='Bursting', w_range=[-5, 60.], extra_fun=f)
model.reset_state()
_vt_plot(model, title='Bursting')
# Transient Spiking
model = AdEx(1, tau=10., a=1., tau_w=100., b=10., V_reset=-60.,
V_rest=-70., V_th=0., V_T=-50., R=0.5, delta_T=2.)
def f():
plt.text(-65., 47.5, 'V nullcline')
plt.text(-69., 8., 'w nullcline')
plt.text(-58.65, 2., 'Trajectory')
plt.annotate('stable focus', xy=(-50.750613238954806, 19.25038786580243),
xytext=(-54., 28.44), arrowprops=dict(arrowstyle="->"))
plt.annotate('saddle node', xy=(-47.9494123367877, 22.050461841990344),
xytext=(-46.4, 16.5), arrowprops=dict(arrowstyle="->"))
_ppa2d(model, title='Transient Spiking',
w_range=[-5, 60.], Iext=55., extra_fun=f, num_text_sp=3)
model.reset_state()
_vt_plot(model, title='Transient Spiking', input=('input', 55.))
# Delayed Spiking
model = AdEx(1, tau=5., a=-1., tau_w=100., b=5., V_reset=-60.,
V_rest=-70., V_th=0., V_T=-50., R=0.5, delta_T=2.)
def f():
plt.text(-62., 11., 'V nullcline')
plt.text(-67., -12.6, 'w nullcline')
plt.text(-48.9, -16.3, 'Trajectory')
_ppa2d(model, title='Delayed Spiking', w_range=[-30, 20.], Iext=25., extra_fun=f,
num_text_sp=2, sp_text={1: '1,2,..'})
model.reset_state()
_vt_plot(model, title='Delayed Spiking', input=('input', 25.))
plt.show()
def AdEx_bifurcation_analysis():
bp.math.enable_x64()
bif = bp.analysis.Bifurcation2D(
model=AdEx(1, tau=10, tau_w=100, a=1, b=10, V_reset=-
60, V_rest=-70, V_T=-50, delta_T=2., R=0.5),
target_vars={'V': [-70., -40.], 'w': [0., 40.]}, # 设置变量的分析范围
target_pars={'Iext': [20., 70.]}, # 设置参数的范围
resolutions={'Iext': 0.1} # 设置分辨率
)
bif.plot_bifurcation(show=True)
def AdEx_ppa2d_transient_spiking():
bp.math.enable_x64()
group = AdEx(1, tau=10., a=1., tau_w=100., b=10., V_reset=-60.,
V_rest=-70., V_th=0., V_T=-50., R=0.5, delta_T=2.)
v_range = [-70., -40.]
w_range = [-5, 60.]
Iext = 56.
fig, gs = bp.visualize.get_figure(1, 1, 4.5, 6)
ax = fig.add_subplot(gs[0, 0])
# 使用BrainPy中的相平面分析工具
phase_plane_analyzer = bp.analysis.PhasePlane2D(
model=group,
target_vars={'V': v_range, 'w': w_range, }, # 待分析变量
pars_update={'Iext': Iext}, # 需要更新的变量
resolutions=0.05
)
# 画出V, w的零增长曲线
phase_plane_analyzer.plot_nullcline()
# 画出奇点
phase_plane_analyzer.plot_fixed_point()
# 画出向量场
phase_plane_analyzer.plot_vector_field(with_return=True)
# 分段画出V, w的变化轨迹
group.V[:] = group.V_reset
runner = bp.DSRunner(
group, monitors=['V', 'w', 'spike'], inputs=('input', Iext))
runner(500)
spike = runner.mon.spike.squeeze()
s_idx = np.where(spike)[0] # 找到所有发放动作电位对应的index
s_idx = np.concatenate(([0], s_idx, [len(spike) - 1])) # 加上起始点和终止点的index
for i in range(len(s_idx) - 1):
vs = runner.mon.V[s_idx[i]: s_idx[i + 1]]
ws = runner.mon.w[s_idx[i]: s_idx[i + 1]]
plt.plot(vs, ws, color='darkslateblue')
# 画出虚线 x = V_reset
plt.plot([group.V_reset, group.V_reset],
w_range, '--', color='grey', zorder=-1)
plt.show()
def ppa2d(group, title, v_range=None, w_range=None, Iext=65., duration=400):
v_range = [-70., -40.] if not v_range else v_range
w_range = [-10., 50.] if not w_range else w_range
phase_plane_analyzer = bp.analysis.PhasePlane2D(
model=group,
target_vars={'V': v_range, 'w': w_range},
pars_update={'Iext': Iext},
resolutions=0.05
)
phase_plane_analyzer.plot_nullcline()
phase_plane_analyzer.plot_fixed_point()
phase_plane_analyzer.plot_vector_field()
runner = bp.DSRunner(
group, monitors=['V', 'w', 'spike'], inputs=('input', Iext))
runner(duration)
spike = runner.mon.spike.squeeze()
s_idx = np.where(spike)[0]
s_idx = np.concatenate(([0], s_idx, [len(spike) - 1]))
for i in range(len(s_idx) - 1):
vs = runner.mon.V[s_idx[i]: s_idx[i + 1]]
ws = runner.mon.w[s_idx[i]: s_idx[i + 1]]
plt.plot(vs, ws, color='darkslateblue')
plt.plot([group.V_reset, group.V_reset],
w_range, '--', color='grey', zorder=-1)
plt.xlim(v_range)
plt.ylim(w_range)
plt.title(title)
plt.show()
if __name__ == '__main__':
ppa2d(AdEx(1, tau=20., a=0., tau_w=30., b=60., V_reset=-55.),
title='Tonic Spiking', w_range=[-5, 75.])
ppa2d(AdEx(1, tau=20., a=0., tau_w=100., b=5., V_reset=-55.),
title='Adaptation', w_range=[-5, 45.])
ppa2d(AdEx(1, tau=5., a=0.5, tau_w=100., b=7., V_reset=-51.),
title='Initial Bursting', w_range=[-5, 50.])
ppa2d(AdEx(1, tau=5., a=-0.5, tau_w=100., b=7., V_reset=-47.),
title='Bursting', w_range=[-5, 60.])
ppa2d(AdEx(1, tau=10., a=1., tau_w=100., b=10., V_reset=-60.),
title='Transient Spiking', w_range=[-5, 60.], Iext=55.)
ppa2d(AdEx(1, tau=5., a=-1., tau_w=100., b=5., V_reset=-60.),
title='Delayed Spiking', w_range=[-30, 20.], Iext=25.)
# compare_Adex_models()
# run()
# simple_run_Adex()
# run_AdEx_model()
# AdEx_patterns()
# AdEx_ppa2d_transient_spiking()
# AdEx_bifurcation_analysis()
扫一扫
1090
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
