泄漏整合发放（LIF）模型
目录
0. HH 模型的优势与不足
1. LIF 模型的定义
2. LIF 模型的动力学性质
3. 基于 BrainPy 的 LIF 模型
4. LIF 模型的优点和不足
4.1. LIF 模型的优点
4.2. LIF 模型的不足
5. 思考与总结
0. HH 模型的优势与不足
关于HH Model ，可以查看我写的这篇博客:Hodgkin-Huxley Model
霍奇金-赫胥黎模型能够精确模拟神经元动作电位的产生，但由于需要小数值积分步长和多个变量，计算成本高。在研究神经网络时，我们更关注神经元之间的相互作用而非膜电位变化的精确细节，因此出现了简化模型，以降低计算开销并保留神经元脉冲发放的主要特征，在生物真实性和计算效率之间取得平衡。
发展至今，神经科学界提出了一系列流行简化模型。这些模型包括：泄露整合发放（Leaky Integrate-and-Fire, LIF）模型、二次整合发放（Quadratic Integrate-and-Fire, QIF）模型、指数整合发放（Exponential Integrate-and-Fire, ExpIF）模型、适应性指数整合发放（Adaptive Exponential Integrateand-Fire, AdEx）模型、Izhikevich 模型、Hindmarsh-Rose（HR）模型和泛化整合发放 (GeneralizedIntegrate-and-Fire, GIF) 模型。下面将依次介绍这些模型，并用BrainPy实现模型的定义、模拟和分析。
1. LIF 模型的定义
我们首先介绍一个最简单的神经元模型泄露整合发放模型(以下简称 LIF)，LIF是 Lapicque 于 1907 年提出的。具体而言：LIF 模型由一个线性微分方程和一个条件判断组成：
τ
d
V
d
t
=
−
(
V
−
V
rest
)
+
R
I
(
t
)
(1)
\tau \frac{dV}{dt} = -(V - V_{\text{rest}}) + R I(t) \tag{1}
τdtdV=−(V−Vrest)+RI(t)(1)
if
V
>
V
th
,
V
←
V
reset
(2)
\text{if } V > V_{\text{th}}, \quad V \leftarrow V_{\text{reset}} \tag{2}
if V>Vth,V←Vreset(2)式中，
τ
=
R
C
\tau = RC
τ=RC 是LIF模型的时间常数（
τ
\tau
τ 越大，模型的动力学变化就越慢），
V
rest
V_{\text{rest}}
Vrest 表示神经元的静息电位，
V
th
V_{\text{th}}
Vth 表示神经元脉冲发放的膜电位阈值。膜电位一旦超过了
V
th
V_{\text{th}}
Vth，则视为神经元发放了脉冲；
V
reset
V_{\text{reset}}
Vreset 表示神经元产生脉冲后的重置膜电位值，
R
R
R 表示细胞膜电阻。具体来说，方程 (1) 表示神经元膜电位
V
V
V 在阈值下接收输入
I
(
t
)
I(t)
I(t) 时不断积分的过程,条件判断式 (2) 表示
V
V
V 达到膜电位阈值后神经元将发放动作电位，同时膜电位将被重置为
V
reset
V_{\text{reset}}
Vreset。
调研发现，对于原始的LIF，其实是没有神经元的不应期的，即神经元在发放动作电位后的一段时间内不再兴奋。要模拟不应期，必须再补充一个条件判断：如果当前时刻距离上次发放的时间小于不应期时长 t ref t_{\text{ref}} tref，则神经元处于不应期，膜电位 𝑉 不会随时间变化。此时条件判断式 (2) 可以改写为： if V > V th , V ← V reset last t ref (3) \text{if } V > V_{\text{th}}, \quad V \leftarrow V_{\text{reset}} \text{ last } t_{\text{ref}} \tag{3} if V>Vth,V←Vreset last tref(3)
2. LIF 模型的动力学性质
与 HH 模型相比，LIF 模型对细胞膜进行了简化处理，将其视为单一通道的电路。具体来说，在膜电位未达到阈值之前，神经元的膜电位主要受泄漏通道调控，该通道的电阻被视为一个固定值，并在模型中表示为常数电阻 R。钠离子和钾离子通道在动作电位发放过程中起关键作用，它们导致膜电位发生非线性变化。然而，在LIF模型中，这一复杂过程被简化为膜电位达到阈值后直接重置为某个设定值，因此不再单独建模这些离子通道的行为。
LIF模型的等效电路图如下，它展示了这种简化的电路结构，其中膜电位的变化主要通过一个恒定电阻和一个电容来描述，考虑膜电位的泄漏效应和重置机制：
根据基尔霍夫电流定律，可以列出等式如下：
I
(
t
)
=
C
d
V
d
t
+
V
−
V
rest
R
(4)
I(t) = C \frac{dV}{dt} + \frac{V - V_{\text{rest}}}{R} \tag{4}
I(t)=CdtdV+RV−Vrest(4)对等式 (4) 进行移项、化简后可以得到：
R
C
d
V
d
t
=
−
(
V
−
V
rest
)
+
R
I
(
t
)
(5)
RC \frac{dV}{dt} = -(V - V_{\text{rest}}) + RI(t) \tag{5}
RCdtdV=−(V−Vrest)+RI(t)(5)代入
τ
=
R
C
\tau = RC
τ=RC 即可得到 LIF 模型的微分方程 (1)。
由于LIF模型形式简单，因此可以对微分方程 (1) 进行理论求解。假设LIF神经元接受恒定输入，即
I
(
t
)
=
I
c
I(t) = I_c
I(t)=Ic，且初始条件为
V
(
t
0
)
=
V
rest
V(t_0) = V_{\text{rest}}
V(t0)=Vrest 时，通过一阶非齐次线性微分方程的常数变易法或分离变量法，可以得到方程 (1) 的解：
V
(
t
)
=
V
reset
+
R
I
c
(
1
−
e
−
t
−
t
0
τ
)
(6)
V(t) = V_{\text{reset}} + RI_c \left(1 - e^{-\frac{t - t_0}{\tau}}\right) \tag{6}
V(t)=Vreset+RIc(1−e−τt−t0)(6)
上述详细的解法过程可参考下图：
上述是一个一阶线性微分方程，关于其求解公式可以参考这篇文章：高等数学学习笔记（1）——微分方程解法公式中的2.1节部分。
由此可见，在给定常值输入
I
c
I_c
Ic 下，神经元从静息电位
V
rest
V_{\text{rest}}
Vrest 积分到阈值发放电位
V
th
V_{\text{th}}
Vth 所需的时间
T
T
T 为：
T
=
−
τ
ln
(
1
−
V
th
−
V
rest
R
I
c
)
(7)
T = -\tau \ln \left(1 - \frac{V_{\text{th}} - V_{\text{rest}}}{RI_c}\right) \tag{7}
T=−τln(1−RIcVth−Vrest)(7)
若不考虑不应期，LIF神经元在恒定电流输入
I
c
I_c
Ic 下的发放率为：
f
=
1
T
=
1
−
τ
ln
(
1
−
V
th
−
V
rest
R
I
c
)
(8)
f = \frac{1}{T} = \frac{1}{-\tau \ln \left(1 - \frac{V_{\text{th}} - V_{\text{rest}}}{RI_c}\right)} \tag{8}
f=T1=−τln(1−RIcVth−Vrest)1(8)
若考虑不应期
t
ref
t_{\text{ref}}
tref，则LIF神经元的发放率为：
f
=
1
T
+
t
ref
=
1
t
ref
−
τ
ln
(
1
−
V
th
−
V
rest
R
I
c
)
(9)
f = \frac{1}{T + t_{\text{ref}}} = \frac{1}{t_{\text{ref}} - \tau \ln \left(1 - \frac{V_{\text{th}} - V_{\text{rest}}}{RI_c}\right)} \tag{9}
f=T+tref1=tref−τln(1−RIcVth−Vrest)1(9)
另一方面，要引发LIF神经元的动作电位，外部电流输入
I
c
I_c
Ic 的值也需要高于某一阈值
I
θ
I_\theta
Iθ，这一阈值被称为最小电流（minimal current），或基强度电流（rheobase current）。根据式 (6)，我们可以容易地得出，当
t
→
∞
t \to \infty
t→∞ 时，
V
(
t
)
→
V
reset
+
R
I
c
V(t) \to V_{\text{reset}} + RI_c
V(t)→Vreset+RIc。因此，引发神经元膜电位到达阈值电位
V
th
V_{\text{th}}
Vth 的基强度电流为：
I
θ
=
V
th
−
V
reset
R
(10)
I_\theta = \frac{V_{\text{th}} - V_{\text{reset}}}{R} \tag{10}
Iθ=RVth−Vreset(10)
只有当
I
>
I
θ
I > I_\theta
I>Iθ 时，
V
V
V 才能越过
V
th
V_{\text{th}}
Vth 产生动作电位。
3. 基于 BrainPy 的 LIF 模型
下面利用 BrainPy 实现了一个包含不应期的 LIF 模型，代码如下：
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 LIF(bp.dyn.NeuDyn):
def __init__(self, size, V_rest=0., V_reset=-5., V_th=20., R=1., tau=10., t_ref=5., **kwargs):
# 初始化父类，继承本身
super(LIF, self).__init__(size=size, **kwargs)
# 初始化参数
self.V_rest = V_rest # 静息膜电位
self.V_reset = V_reset
self.V_th = V_th
self.R = R
self.tau = tau
self.t_ref = t_ref # 不应期时长
# 初始化变量
self.V = bm.Variable(bm.random.randn(self.num) + V_reset)
self.input = bm.Variable(bm.zeros(self.num))
# 相对于HH模型多出的三个变量，用于判断脉冲时间，不应期等
self.t_last_spike = bm.Variable(bm.ones(self.num) * -1e7) # 上一次脉冲发放时间
self.refractory = bm.Variable(
bm.zeros(self.num, dtype=bool)) # 是否处于不应期
self.spike = bm.Variable(bm.zeros(self.num, dtype=bool)) # 脉冲发放状态
# 使用指数欧拉方法进行积分
self.integral = bp.odeint(
f=self.derivative, method='exponential_euler')
# 定义膜电位关于时间变化的微分方程
def derivative(self, V, t, Iext):
dvdt = (-V + self.V_rest + self.R * Iext) / self.tau
return dvdt
def update(self):
_t, _dt = bp.share['t'], bp.share['dt']
# 以数组的方式对神经元进行更新
refractory = (_t - self.t_last_spike) <= self.t_ref # 判断神经元是否处于不应期
V = self.integral(self.V, _t, self.input, dt=_dt) # 根据时间步长更新膜电位
# 若处于不应期，则返回原始膜电位self.V，否则返回更新后的膜电位V
V = bm.where(refractory, self.V, V)
spike = V > self.V_th # 将大于阈值的神经元标记为发放了脉冲
self.spike[:] = spike # 更新神经元脉冲发放状态
self.t_last_spike[:] = bm.where(
spike, _t, self.t_last_spike) # 更新最后一次脉冲发放时间
# 将发放了脉冲的神经元膜电位置为V_reset，其余不变
self.V[:] = bm.where(spike, self.V_reset, V)
self.refractory[:] = bm.logical_or(refractory, spike) # 更新神经元是否处于不应期
self.input[:] = 0. # 重置外界输入
上述类的定义和HH Model类似，相关的微分方程按照公式写入即可。
在定义好 LIF 模型后，可以初始化一个 LIF 模型并运行数值模拟：
def run_LIF():
# 运行LIF模型
group = LIF(1)
runner = bp.DSRunner(group, monitors=['V'], inputs=('input', 22.))
runner(200) # 运行时长为200ms
# 结果可视化
fig, gs = bp.visualize.get_figure(1, 1, 4.5, 6)
# 创建一个1行1列的网格，宽为4.5，高为6
ax = fig.add_subplot(gs[0, 0])
# 将子图添加到图形中
plt.plot(runner.mon.ts, runner.mon.V)
plt.xlabel(r'$t$ (ms)')
plt.ylabel(r'$V$ (mV)')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
# plt.savefig('LIF_output.pdf', transparent=True, dpi=500)
plt.show()
得到的图像如下：
上述图形确实模拟了动作电位的产生，但是并不是模拟的真实的动作电位的情况。在发放动作电位前，LIF神经元膜电位的增长速度是逐渐降低的，并不是像真实的神经元先缓慢增长，到达一定值后再迅速增长。
对于上面的式(9)，可以通过数值模拟验证理论分析的结果。当 I c → ∞ I_c \to \infty Ic→∞ 时，发放间隔 T → 0 T \to 0 T→0，而发放率 f → 1 t ref f \to \frac{1}{t_{\text{ref}}} f→tref1。比如，当 t ref = 5 ms t_{\text{ref}} = 5 \, \text{ms} tref=5ms 时，发放率会随 I c I_c Ic 的升高逐渐趋于 1 5 ms = 200 Hz \frac{1}{5 \, \text{ms}} = 200 \, \text{Hz} 5ms1=200Hz。
具体来说：
- 当输入电流 I c I_c Ic 趋向于无穷大时，神经元的发放间隔 T T T 趋向于零。
- 发放率 f f f 将趋近于 1 t ref \frac{1}{t_{\text{ref}}} tref1。
- 例如，若不应期 t ref t_{\text{ref}} tref 为 5 毫秒，则发放率将逐渐接近 1 5 ms = 200 Hz \frac{1}{5 \, \text{ms}} = 200 \, \text{Hz} 5ms1=200Hz。
编写如下代码进行实验仿真：
def fi_curve():
duration = 1000 # 设定仿真时长
I_cur = np.arange(0, 600, 2) # 定义大小为10, 20, ..., 1000的100束电流
neu = LIF(len(I_cur), tau=5., t_ref=5.) # 定义神经元群
runner = bp.DSRunner(neu, monitors=['spike'], inputs=(
'input', I_cur), dt=0.01) # 设置运行器，其中每一个神经元接收I_cur的一个恒定电流
runner(duration=duration) # 运行神经元模型
f_list = runner.mon.spike.sum(axis=0) / (duration / 1000) # 计算每个神经元的脉冲发放次数
fig, gs = bp.visualize.get_figure(1, 1, 4.5, 6)
ax = fig.add_subplot(gs[0, 0])
plt.plot(I_cur, f_list)
plt.xlabel('Input current')
plt.ylabel('spiking frequency')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
# plt.savefig('LIF_fi_curve.pdf', transparent=True, dpi=500)
plt.show()
可以发现，在恒定电流输入的条件下，LIF 神经元的发放频率随电流的增加而增加，最终趋于一个恒定值。前期Input current 在小于一定的值时，是没有发放频率的，该值被称为基强度电流。
对于式(10)，可以假设
V
th
=
20
V_{\text{th}} = 20
Vth=20，
V
reset
=
0
V_{\text{reset}} = 0
Vreset=0，
R
=
1
R = 1
R=1，则计算得输入的阈值
I
θ
=
20
I_\theta = 20
Iθ=20。也就是说，当外部恒定电流的输入值大于 20 时，神经元才能产生动作电位。
如下代码模拟了上述情况:
def threshold_1():
fig, gs = bp.visualize.get_figure(1, 1, 4.5, 6)
ax = fig.add_subplot(gs[0, 0])
duration = 200
neu1 = LIF(1, t_ref=5.)
neu1.V[:] = bm.array([-5.]) # 设置V的初始值
runner = bp.DSRunner(neu1, monitors=['V'], inputs=(
'input', 20)) # 给neu1一个大小为20的恒定电流
runner(duration)
plt.plot(runner.mon.ts, runner.mon.V, label='$I = 20$')
plt.text(153, 21, r'$I = 20$')
neu2 = LIF(1, t_ref=5.)
neu2.V[:] = bm.array([-5.]) # 设置V的初始值
runner = bp.DSRunner(neu2, monitors=['V'], inputs=(
'input', 21)) # 给neu2一个大小为21的恒定电流
runner(duration)
plt.plot(runner.mon.ts, runner.mon.V, label='$I = 21$')
plt.text(153, -2, r'$I = 21$')
plt.xlabel(r'$t$ (ms)')
plt.ylabel(r'$V$ (mV)')
plt.ylim(-6, 23)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.show()
上图中，强度为 20 的电流不能激发动作电位，因为此时神经元的膜电位随时间无限逼近但并不能达到
V
th
V_{\text{th}}
Vth。将输入电流的强度稍微升高至 21 即可以激发神经元的发放。
为了更好的观察 LIF 模型的发放规律，我们在上述图片的基础上引入了 I = 10 I=10 I=10和 I = 22 I=22 I=22,具体实现代码如下：
def threshold_2():
fig, gs = bp.visualize.get_figure(1, 1, 4.5, 6)
ax = fig.add_subplot(gs[0, 0])
duration = 200
neu1 = LIF(1, t_ref=5.)
neu1.V[:] = bm.array([-5.]) # 设置V的初始值
runner = bp.DSRunner(neu1, monitors=['V'], inputs=(
'input', 20)) # 给neu1一个大小为20的恒定电流
runner(duration)
plt.annotate('$I = 20$', xy=(130, 20), xytext=(153, 21),
arrowprops=dict(facecolor='black', shrink=0.005, width=1, headwidth=5, headlength=5))
plt.plot(runner.mon.ts, runner.mon.V, label='$I = 20$')
neu2 = LIF(1, t_ref=5.)
neu2.V[:] = bm.array([-5.]) # 设置V的初始值
runner = bp.DSRunner(neu2, monitors=['V'], inputs=(
'input', 21)) # 给neu2一个大小为21的恒定电流
runner(duration)
plt.annotate('$I = 21$', xy=(146, -5), xytext=(153, -2),
arrowprops=dict(facecolor='black', shrink=0.05, width=1, headwidth=5, headlength=5))
plt.plot(runner.mon.ts, runner.mon.V, label='$I = 21$')
neu3 = LIF(1, t_ref=5.)
neu3.V[:] = bm.array([-5.]) # 设置V的初始值
runner = bp.DSRunner(neu3, monitors=['V'], inputs=(
'input', 22)) # 给neu3一个大小为22的恒定电流
runner(duration)
plt.annotate('$I = 22$', xy=(88, 18), xytext=(100, 21),
arrowprops=dict(facecolor='black', shrink=0.005, width=1, headwidth=5, headlength=5))
plt.plot(runner.mon.ts, runner.mon.V, label='$I = 22$')
neu4 = LIF(1, t_ref=5.)
neu4.V[:] = bm.array([-5.]) # 设置V的初始值
runner = bp.DSRunner(neu4, monitors=['V'], inputs=(
'input', 10)) # 给neu4一个大小为10的恒定电流
runner(duration)
plt.annotate('$I = 10$', xy=(50, 10), xytext=(10, 21),
arrowprops=dict(facecolor='black', shrink=0.005, width=1, headwidth=5, headlength=5))
plt.plot(runner.mon.ts, runner.mon.V, label='$I = 10$')
plt.xlabel(r'$t$ (ms)')
plt.ylabel(r'$V$ (mV)')
plt.ylim(-6, 23)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.show()
结果展示如下:
可以发现，当输入电流
I
I
I小于基强度电流时，在给定的参数条件下，膜电位会逐渐增大并稳定在某一特定值，且输入电流越大，稳定的膜电位值也越高；当输入电流超过基强度电流时，神经元会激发动作电位，且激发频率随着输入电流的增大而加快。
if __name__ == '__main__':
run_LIF()
fi_curve()
threshold_1()
threshold_2()
笔者通过理论分析，认为在LIF模型中，当输入电流足够大时，神经元的发放频率会受到不应期（refractory period）的限制，从而导致频率达到一个理论上限。然而，在上述的模拟过程中，笔者并未观察到这一现象。随着输入电流的不断增大，神经元发放频率虽然逐渐增加，但并未达到一个明确的上限，而是呈现出无限接近某一值的趋势。笔者这种现象可能与以下原因有关：
- 不应期设置不正确：不应期的设置需合理，并且在代码中正确应用。
- 模拟时间不足：模拟时间需足够长，以便观察到稳定的发放频率。
- 电流输入设置问题：电流输入需设置正确，并且足够大以达到预期的效果。
至于上述现象出现的具体原因，限于笔者的能力，不能找出具体的原因。有兴趣的读者可以自行 C o d i n g Coding Coding，并可以和我进一步交流。
4. LIF 模型的优点和不足
4.1. LIF 模型的优点
LIF模型是一个高度简化的神经元模型，但同时也具有不可忽视的优点。它的优点包括：
- 模型简洁、直观：微分方程 (1) 和条件判断式 (2)（或 (3)）中每个参数都有明确的意义，便于理解；
- 计算效率高：仅由一个线性微分方程组成，在数值模拟时具有很高的计算效率；
- 较好的精确性：在模拟阈下膜电位的变化时仍然具有较好的精确性。也就是说，如果神经元不发放脉冲，它的膜电位变化是可以用线性微分方程捕捉的。
基于以上的优点，LIF模型至今仍然在很多研究中被广泛使用。
4.2. LIF 模型的不足
我们也应该意识到LIF模型存在的缺陷：
- 过度简化的微分方程：使它不能模拟真实的神经元发放过程。在大多数情况下，接受外部刺激的神经元的膜电位会先缓慢上升，一旦超过一个阈值就变为快速上升，再快速下降，由此形成了动作电位。LIF模型只建模了阈下膜电位的变化，而把膜电位快速上升和下降的脉冲发放过程简化为了膜电位的重置，这使得模型缺少了对非线性过程的表征。
- 缺乏记忆功能：不保有对之前脉冲活动的任何记忆。换言之，对于每次重置后的时刻，LIF神经元的状态都是相同的，它受到的调控有且仅有外部电流，因此LIF模型更像是定义了一个从外部电流输入 I I I 到膜电位 V V V 的映射，其膜电位变化与之前的活动历史没有任何关联。
- 无法模拟复杂活动：不能模拟更为复杂的神经元活动。正因为LIF神经元的膜电位变化只与外部电流输入有关，这决定了LIF模型能够模拟的神经元发放模式非常单一，一些复杂的发放模式，如适应性发放和簇发放（后面会写博客进行介绍），都不能被LIF模型重现。
为弥补LIF模型的缺陷，计算神经科学家们在LIF模型的基础上陆续提出了新的简化模型，试图从不同的角度增加模型的动力学表征能力。
5. 思考与总结
LIF模型（Leaky Integrate-and-Fire Model）作为计算神经科学领域中最经典且广泛使用的神经元模型之一，具有重要的历史地位和实用价值。它以简洁的数学形式捕捉了神经元膜电位的基本动态特性，尽管忽略了动作电位产生的详细生物物理机制，但其高效的计算性能和直观的物理意义使其成为研究大规模神经网络和神经动力学的理想工具。
从历史角度来看，LIF模型的提出为神经科学和计算科学之间的桥梁奠定了重要基础。它不仅是早期神经元建模的核心方法之一，也为后续更复杂的模型的发展提供了重要的理论框架和启发。LIF模型的简洁性和可扩展性使其在理论研究和工程应用中得到了广泛应用，例如在神经形态计算、脑机接口和人工智能领域。
我对LIF模型的态度是高度肯定的。尽管它在生物细节上存在一定的简化，但其在描述神经元整体行为方面的能力令人印象深刻。它不仅是教学和理论分析中的有力工具，也为实际应用提供了高效的解决方案。LIF模型的存在，让我们知道，科学建模的核心在于在复杂性和实用性之间找到平衡，而这一点正是LIF模型的成功之处。
总之，LIF模型作为神经元建模的基石之一，不仅在历史上具有重要意义，也在当今的神经科学和工程应用中持续发挥着不可替代的作用。它的简洁性和普适性使其成为理解神经元行为和设计神经形态系统的关键工具。
扫一扫
5011
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
