QIF vs LIF模型实战对比：用BrainPy可视化神经元放电的5个关键差异

QIF vs LIF模型实战对比：用BrainPy可视化神经元放电的5个关键差异

在神经计算建模的实践中，选择哪种简化神经元模型，往往取决于你究竟想“看”到什么。是追求计算效率，还是希望捕捉动作电位那微妙的上扬曲线？对于已经熟悉Hodgkin-Huxley这类生物物理模型的研究者来说，LIF（Leaky Integrate-and-Fire）和QIF（Quadratic Integrate-and-Fire）是迈向大规模网络模拟时绕不开的两个基石。但它们的差异绝非仅仅是数学公式上的一个平方项那么简单。

这篇文章不是对单一模型的复述，而是一次并排运行的实战对比。我们将直接使用BrainPy，像在实验室里并排放置两台示波器一样，同时运行LIF和QIF模型，从动态的电压轨迹、精确的频率-电流响应，到核心参数如何实时塑造神经元行为，进行逐项拆解。你会发现，a_0系数的一个微小滑动，就能在QIF模型中引发从平滑整合到爆发式放电的戏剧性转变，而这种直观的调控在LIF中是无法实现的。我们关注的是在交叉验证中才能凸显的细节：阈值附近的动态、对输入信号的编码保真度，以及它们各自在大规模仿真中可能带来的陷阱与惊喜。文中的所有代码都经过精心设计，力求清晰、可复用，你可以直接将其移植到Colab笔记本中，调整参数，亲眼见证这些差异是如何从方程中“生长”出来的。

1. 从电路到方程：理解两种模型的根本出发点

在深入代码对比之前，我们必须厘清LIF和QIF模型在设计哲学上的分岔点。两者都是对真实神经元的高度抽象，但抽象的侧重点不同，这直接决定了它们的行为边界。

LIF模型的核心隐喻是一个漏电的电容器。它将神经元膜视为一个简单的RC电路：膜电容C并联一个泄漏电阻R。输入电流 I(t) 一方面为电容充电，使膜电位 V(t) 上升；另一方面，电荷通过泄漏电阻不断流失，试图将电位拉回静息电位 V_rest。其微分方程简洁明了：

τ * dV/dt = -(V - V_rest) + R * I(t)

这里的 τ = R*C 是膜时间常数，决定了电位变化的“惯性”。当 V 达到设定的阈值 V_th 时，模型触发一个“脉冲”事件，V 瞬间被重置为 V_reset，并可能进入一个短暂的不应期。关键在于，这个充电过程是指数趋近的，电位上升速率会随着接近阈值而逐渐放缓，这与真实动作电位在阈值附近加速上升的特性相悖。

QIF模型正是为了修补这一缺陷而生。它在LIF的线性恢复项基础上，引入了一个二次项，其方程形式为：

τ * dV/dt = a_0 * (V - V_rest) * (V - V_c) + R * I(t)

这里多出了一个关键参数 V_c（临界电压或转折电压），以及系数 a_0。这个二次项 (V - V_rest)(V - V_c) 在 V 介于 V_rest 和 V_c 之间时为负（抑制），当 V > V_c 时为正（兴奋），且其绝对值随 V 偏离零点而增大。这使得膜电位的动力学在 V_c 附近发生质变：一旦跨过这个点，去极化会自我强化，产生一个急剧上升的“锋电位”雏形，从而更逼真地模拟了钠通道激活带来的正反馈过程。

为了更直观地看到这种数学结构上的差异如何转化为动力学的不同，我们可以对比它们的相图（nullcline）和流场。下面的表格概括了这种核心区别：

提示：V_c 在QIF模型中不是一个发放阈值，而是一个动力学分岔点。它决定了模型从“整合”模式切换到“发放”模式的临界电压，通常设置在 V_rest 和 V_th 之间。

理解了这个根本区别，我们就能明白，为何在同样的输入电流下，两个模型会讲出不同的“电生理故事”。接下来，我们就用BrainPy让这两个故事同时上演。

2. 并排实现：构建可对比的BrainPy模型类

为了进行公平的对比，我们需要构建结构清晰的模型类，确保除了核心动力学方程外，其他部分（如状态更新逻辑、不应期处理、脉冲重置）尽可能一致。这能让我们将观察到的任何差异都归因于模型本身，而非代码实现上的偶然性。

首先，我们导入必要的库，并定义一个基础的神经元组类，将共享的逻辑封装起来。

import brainpy as bp
import brainpy.math as bm
import numpy as np
import matplotlib.pyplot as plt

class BaseSpikingNeuron(bp.dyn.NeuGroup):
    """一个包含脉冲重置和不应期处理的基类"""
    def __init__(self, size, V_rest, V_reset, V_th, t_ref, name=None):
        super().__init__(size=size, name=name)
        # 共享的参数
        self.V_rest = V_rest
        self.V_reset = V_reset
        self.V_th = V_th
        self.t_ref = t_ref  # 不应期时长

        # 共享的状态变量
        self.V = bm.Variable(bm.ones(self.num) * V_reset)
        self.input = bm.Variable(bm.zeros(self.num))
        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))

    def update_shared(self, t, dt):
        """处理不应期、脉冲检测和重置的共享逻辑"""
        # 判断是否处于不应期
        is_refractory = (t - self.t_last_spike) <= self.t_ref

        # 更新膜电位（由子类的integral具体实现）
        V_new = self.integrate_V(t, dt)

        # 不应期内电位保持不变
        V_new = bm.where(is_refractory, self.V, V_new)

        # 检测脉冲
        spike_new = V_new >= self.V_th

        # 更新状态
        self.spike.value = spike_new
        self.t_last_spike.value = bm.where(spike_new, t, self.t_last_spike)
        self.V.value = bm.where(spike_new, self.V_reset, V_new)
        self.refractory.value = bm.logical_or(is_refractory, spike_new)

        # 清空输入缓存，为下一步准备
        self.input[:] = 0.

基于这个基类，我们可以非常简洁地实现LIF和QIF模型，只需重写膜电位积分的部分。

class LIFNeuron(BaseSpikingNeuron):
    def __init__(self, size, V_rest=-65., V_reset=-68., V_th=-50.,
                 R=1., tau=10., t_ref=2., name=None):
        super().__init__(size, V_rest, V_reset, V_th, t_ref, name)
        self.R = R  # 膜电阻
        self.tau = tau  # 膜时间常数
        # 定义LIF特有的微分方程和积分器
        def dV(V, t, I):
            return (-(V - self.V_rest) + self.R * I) / self.tau
        self.integral = bp.odeint(dV, method='exp_auto')

    def integrate_V(self, t, dt):
        # 调用积分器计算新的膜电位
        return self.integral(self.V, t, self.input, dt=dt)

class QIFNeuron(BaseSpikingNeuron):
    def __init__(self, size, V_rest=-65., V_reset=-68., V_th=-50.,
                 V_c=-55., a_0=0.07, R=1., tau=10., t_ref=2., name=None):
        super().__init__(size, V_rest, V_reset, V_th, t_ref, name)
        self.V_c = V_c  # 临界电压
        self.a_0 = a_0   # 二次项系数
        self.R = R
        self.tau = tau
        # 定义QIF特有的微分方程和积分器
        def dV(V, t, I):
            return (self.a_0 * (V - self.V_rest) * (V - self.V_c) + self.R * I) / self.tau
        self.integral = bp.odeint(dV, method='rk4')  # 对于非线性更强的方程，RK4可能更稳定

    def integrate_V(self, t, dt):
        return self.integral(self.V, t, self.input, dt=dt)

注意：在QIF实现中，我们选择了rk4积分方法。这是因为二次项在电位接近V_c和V_th时变化剧烈，使用exp_auto（针对线性系统优化）有时可能导致数值不稳定或精度下降。在实际对比中，确保两个模型都使用稳定、精度足够的方法是关键。

有了这两个类，我们就拥有了并排对比的“实验仪器”。接下来，我们接通电流，观察第一个也是最直观的差异：动作电位的形态。

3. 关键差异一：动作电位上升相的“表情”

给神经元一个恒定的去极化电流，观察其产生的动作电位（更准确地说，是达到阈值触发重置的电压轨迹），是检验模型生物逼真度的第一关。我们设置相同的仿真环境和刺激电流，让LIF和QIF同时运行。

# 仿真参数
duration = 150  # 毫秒
I_stim = 12.0   # 纳安 (nA)

# 创建神经元实例（使用相近的默认参数）
lif_neuron = LIFNeuron(1, V_rest=-65., V_reset=-68., V_th=-50., tau=10., t_ref=2.)
qif_neuron = QIFNeuron(1, V_rest=-65., V_reset=-68., V_th=-50., V_c=-55., a_0=0.07, tau=10., t_ref=2.)

# 创建运行器并施加恒定电流
runner_lif = bp.dyn.DSRunner(lif_neuron,
                             monitors=['V', 'spike'],
                             inputs=[('input', I_stim)],
                             dt=0.01)
runner_qif = bp.dyn.DSRunner(qif_neuron,
                             monitors=['V', 'spike'],
                             inputs=[('input', I_stim)],
                             dt=0.01)

print("正在运行LIF模型...")
runner_lif.run(duration)
print("正在运行QIF模型...")
runner_qif.run(duration)

# 可视化对比
fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True)

axes[0].plot(runner_lif.mon.ts, runner_lif.mon.V[:, 0], 'b-', label='LIF', linewidth=2)
axes[0].axhline(y=lif_neuron.V_th, color='b', linestyle='--', alpha=0.5, label='LIF Threshold')
axes[0].set_ylabel('Membrane Potential (mV)')
axes[0].set_title('LIF Model - Voltage Trace')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(runner_qif.mon.ts, runner_qif.mon.V[:, 0], 'r-', label='QIF', linewidth=2)
axes[1].axhline(y=qif_neuron.V_th, color='r', linestyle='--', alpha=0.5, label='QIF Threshold')
axes[1].axhline(y=qif_neuron.V_c, color='orange', linestyle=':', alpha=0.7, label='QIF V_c (Critical)')
axes[1].set_xlabel('Time (ms)')
axes[1].set_ylabel('Membrane Potential (mV)')
axes[1].set_title('QIF Model - Voltage Trace')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

运行这段代码，你会立刻看到鲜明的对比。LIF的电压上升是一条光滑的、斜率逐渐减小的指数曲线，它“温柔”地触及阈值。而QIF的曲线则生动得多：在初始阶段，上升缓慢，类似于LIF；但当电位跨过那条橙色的虚线（V_c）后，上升速度陡然增加，画出一条近乎垂直的上升支，然后被阈值截断并重置。这个“弯道加速”的特性，正是QIF模型能更好近似真实神经元动作电位上升支的原因。它捕捉到了电压门控钠通道一旦被大量激活，便势不可挡的正反馈过程。

这种形态差异不仅仅是视觉上的。它意味着QIF神经元对达到阈值前最后一段时间的输入电流更为敏感，因为此时增益（导数）很大。而LIF神经元在接近阈值时反而变得“迟钝”。这直接影响了下一点要对比的内容：它们对输入强度的频率编码特性。

4. 关键差异二：频率-电流响应曲线的斜率与线性度

神经元的频率-电流（f-I）曲线描述了其输出脉冲频率如何随输入电流强度变化，这是神经元作为信息编码器的核心传输函数。我们通过扫描一系列输入电流，测量每个电流下神经元的平均发放频率，来绘制这条曲线。

def calculate_fI_curve(neuron_class, neuron_params, I_range, duration, dt=0.01):
    """计算给定神经元类和参数下的f-I曲线"""
    frequencies = []
    for I in I_range:
        # 为每个电流强度创建新的神经元实例，避免状态残留
        neu = neuron_class(1, **neuron_params)
        runner = bp.dyn.DSRunner(neu,
                                 monitors=['spike'],
                                 inputs=[('input', I)],
                                 dt=dt)
        runner.run(duration)
        # 计算平均发放频率 (Hz)。忽略初始瞬态，例如前50ms
        start_idx = int(50 / dt)
        spike_times = runner.mon.ts[start_idx:][runner.mon.spike[start_idx:, 0]]
        if len(spike_times) > 1:
            avg_isi = np.mean(np.diff(spike_times))  # 平均脉冲间隔 (ms)
            freq = 1000.0 / avg_isi  # 转换为 Hz
        else:
            freq = 0.0
        frequencies.append(freq)
    return np.array(frequencies)

# 定义参数和电流扫描范围
base_params = {
    'V_rest': -65., 'V_reset': -68., 'V_th': -50., 'tau': 10., 't_ref': 2.
}
qif_extra_params = {'V_c': -55., 'a_0': 0.07}
lif_params = {**base_params, 'R': 1.0}
qif_params = {**base_params, **qif_extra_params, 'R': 1.0}

I_values = np.linspace(0, 20, 41)  # 从0到20 nA，41个点
sim_duration = 500  # 每个仿真运行500ms以获得稳定频率

print("正在扫描LIF的f-I曲线...")
lif_freqs = calculate_fI_curve(LIFNeuron, lif_params, I_values, sim_duration)
print("正在扫描QIF的f-I曲线...")
qif_freqs = calculate_fI_curve(QIFNeuron, qif_params, I_values, sim_duration)

# 可视化对比
plt.figure(figsize=(8, 5))
plt.plot(I_values, lif_freqs, 'bo-', label='LIF', markersize=4, linewidth=1.5)
plt.plot(I_values, qif_freqs, 'rs-', label='QIF (a_0=0.07)', markersize=4, linewidth=1.5)
plt.xlabel('Input Current I (nA)')
plt.ylabel('Firing Frequency f (Hz)')
plt.title('Frequency-Current (f-I) Curves Comparison')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

得到的图表会揭示几个重要信息：

1. 阈值电流：LIF模型有一个非常尖锐的阈值。电流低于某个临界值（约10nA，取决于参数）时，频率严格为零；一旦超过，频率从零跳变到一个有限值。而QIF模型的阈值区域更“柔和”，在阈值附近，频率是连续地从零开始增加的。
2. 曲线形状：在远高于阈值后，LIF的f-I曲线通常接近线性（f ∝ I）。QIF的曲线则可能保持线性，也可能呈现轻微的上凸或下凹，这高度依赖于 a_0 和 V_c 的取值。QIF提供了额外的自由度来塑造神经元的输入-输出关系。
3. 增益（斜率）：在相同参数下，QIF曲线在阈值以上的初始斜率（即增益）往往高于LIF。这意味着对于小幅超阈值电流，QIF神经元能产生更高的频率响应，编码更敏感。

为了更精确地量化阈值行为，我们可以编写一个函数来寻找“基强度电流”（rheobase），即能引发神经元发放的最小恒定电流。

def find_rheobase(neuron_class, neuron_params, I_low, I_high, tolerance=0.01, duration=200):
    """使用二分法寻找基强度电流"""
    low, high = I_low, I_high
    while (high - low) > tolerance:
        mid = (low + high) / 2
        neu = neuron_class(1, **neuron_params)
        runner = bp.dyn.DSRunner(neu, monitors=['spike'], inputs=[('input', mid)])
        runner.run(duration)
        # 检查后一半仿真时间内是否有脉冲
        if runner.mon.spike[int(duration/runner.dt/2):, 0].any():
            high = mid  # 能发放，电流过高
        else:
            low = mid   # 不能发放，电流过低
    return (low + high) / 2

lif_rheo = find_rheobase(LIFNeuron, lif_params, 8, 12, tolerance=0.05)
qif_rheo = find_rheobase(QIFNeuron, qif_params, 3, 7, tolerance=0.05)
print(f"LIF模型的基强度电流（理论近似）: {lif_rheo:.2f} nA")
print(f"QIF模型的基强度电流（理论近似）: {qif_rheo:.2f} nA")

你会发现，在给定的参数集下，QIF的基强度电流显著低于LIF。这再次印证了QIF在阈值附近更高的敏感性。这种差异直接影响到网络动力学：一个由QIF神经元构成的网络，其整体活跃度可能对平均输入强度的微小变化反应更剧烈。

5. 关键差异三：参数 a_0 的动态调控与分岔行为

在LIF模型中，一旦 τ、V_th 等参数固定，其行为模式就基本确定了。而QIF模型中的 a_0 参数是一个强大的“旋钮”，它能动态且连续地改变模型的兴奋性，甚至引发动力学的分岔（Bifurcation）。这使得QIF不仅能模拟常规的整合发放，还能模拟Izhikevich模型中某些复杂的放电模式。

我们可以设计一个实验：在单次仿真中，让 a_0 随时间线性变化，观察神经元放电频率和模式的实时响应。

class QIFWithDynamicA0(QIFNeuron):
    """扩展QIF类，使其a_0可以随时间变化"""
    def __init__(self, size, a_0_func=None, **kwargs):
        super().__init__(size, **kwargs)
        self.a_0_func = a_0_func if a_0_func else (lambda t: self.a_0)  # 默认是常数
        # 需要重写积分函数以使用动态的a_0
        def dV_dynamic(V, t, I):
            a0_now = self.a_0_func(t)
            return (a0_now * (V - self.V_rest) * (V - self.V_c) + self.R * I) / self.tau
        self.integral = bp.odeint(dV_dynamic, method='rk4')

# 定义a_0随时间变化的函数：从0.02线性增加到0.12
def a0_ramp(t):
    a0_start, a0_end = 0.02, 0.12
    ramp_duration = 400  # 毫秒
    if t <= ramp_duration:
        return a0_start + (a0_end - a0_start) * (t / ramp_duration)
    else:
        return a0_end

# 创建动态神经元并仿真
dynamic_neuron = QIFWithDynamicA0(1, a_0_func=a0_ramp, V_rest=-65., V_reset=-68.,
                                   V_th=-50., V_c=-55., a_0=0.07, R=1., tau=10., t_ref=2.)
I_stim_dynamic = 5.0  # 一个固定的中等强度电流
runner_dynamic = bp.dyn.DSRunner(dynamic_neuron,
                                  monitors=['V', 'spike', ('a_0', lambda t: a0_ramp(t))],
                                  inputs=[('input', I_stim_dynamic)],
                                  dt=0.01)
runner_dynamic.run(500)

# 可视化结果
fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True, gridspec_kw={'height_ratios': [1, 2, 1]})

# 子图1：a_0参数的变化
axes[0].plot(runner_dynamic.mon.ts, runner_dynamic.mon.a_0, 'g-', linewidth=2)
axes[0].set_ylabel('$a_0$ parameter')
axes[0].set_title('Dynamic Modulation of $a_0$')
axes[0].grid(True, alpha=0.3)

# 子图2：膜电位轨迹
axes[1].plot(runner_dynamic.mon.ts, runner_dynamic.mon.V[:, 0], 'k-', linewidth=1)
axes[1].axhline(y=dynamic_neuron.V_th, color='r', linestyle='--', alpha=0.7, label='Threshold')
axes[1].axhline(y=dynamic_neuron.V_c, color='orange', linestyle=':', alpha=0.7, label='V_c')
axes[1].set_ylabel('Membrane Potential (mV)')
axes[1].legend(loc='upper right')
axes[1].grid(True, alpha=0.3)

# 子图3：脉冲发放（点图）
spike_times = runner_dynamic.mon.ts[runner_dynamic.mon.spike[:, 0]]
axes[2].eventplot(spike_times, colors='black', lineoffsets=0, linelengths=0.8)
axes[2].set_xlabel('Time (ms)')
axes[2].set_ylabel('Spikes')
axes[2].set_yticks([])
axes[2].grid(True, alpha=0.3, axis='x')

plt.tight_layout()
plt.show()

# 计算并打印不同a_0区间的平均频率
print("放电模式分析:")
spike_mask = runner_dynamic.mon.spike[:, 0]
if spike_mask.any():
    all_spike_times = runner_dynamic.mon.ts[spike_mask]
    # 划分时间段：a_0 < 0.05, 0.05 <= a_0 < 0.09, a_0 >= 0.09
    a0_vals = runner_dynamic.mon.a_0[spike_mask]
    intervals = [(0, 0.05), (0.05, 0.09), (0.09, 0.13)]
    for low, high in intervals:
        idx = (a0_vals >= low) & (a0_vals < high)
        if idx.any():
            seg_times = all_spike_times[idx]
            if len(seg_times) > 1:
                avg_freq = 1000.0 / np.mean(np.diff(seg_times))
                print(f"  a_0在[{low:.2f}, {high:.2f})区间: 平均频率 ≈ {avg_freq:.1f} Hz")
            else:
                print(f"  a_0在[{low:.2f}, {high:.2f})区间: 脉冲太少，无法计算稳定频率")
else:
    print("  在整个仿真过程中未发放脉冲。")

运行这个动态实验，你会观察到随着 a_0 从低到高增加，神经元的放电行为可能经历几个阶段：从静息，到低频不规则发放，再到高频规律发放，甚至可能出现频率的饱和或下降（取决于其他参数）。a_0 本质上控制了二次项“势阱”的深度和形状，从而改变了膜电位的有效驱动力和恢复力之间的平衡。

• 低 a_0：模型行为接近LIF，二次项影响弱，需要更强的电流才能达到阈值，发放频率较低。
• 中等 a_0：二次项效应显著，在 V_c 附近产生强烈的正反馈，导致快速去极化和较高的发放频率。
• 高 a_0：恢复力项 (V - V_rest)(V - V_c) 在亚阈值区也变得很大，可能使神经元更难去极化，反而可能抑制发放，或者导致发放后复位到更深的超极化状态。

这种通过单一参数连续调节放电模式的能力，让QIF模型在模拟具有不同兴奋性类型的神经元（如规则发放型、初始爆发型等）时更具灵活性。相比之下，LIF模型要模拟这些差异，通常需要引入额外的自适应电流或更复杂的结构。

6. 关键差异四：在脉冲序列编码与噪声鲁棒性上的表现

当我们不仅仅给神经元施加恒定电流，而是输入一个随时间变化的信号或叠加了噪声的电流时，LIF和QIF在编码信息方面的差异会更加明显。一个常见的测试是给它们输入一个振荡电流（如正弦波），观察输出脉冲序列是如何对输入相位进行编码的。

# 定义振荡输入电流
def oscillatory_current(t, freq=5, amplitude=3, bias=8):
    """正弦振荡电流，单位nA"""
    return bias + amplitude * np.sin(2 * np.pi * freq * t / 1000)  # 频率freq Hz

duration_osc = 1000  # 仿真1000ms
time_points = np.arange(0, duration_osc, 0.01)
I_osc = oscillatory_current(time_points, freq=5, amplitude=4, bias=10)

# 运行两个模型
lif_osc = LIFNeuron(1, **lif_params)
qif_osc = QIFNeuron(1, **qif_params)

# 使用函数式输入
runner_lif_osc = bp.dyn.DSRunner(lif_osc,
                                 monitors=['V', 'spike'],
                                 inputs=[('input', oscillatory_current, 't')],
                                 dt=0.01)
runner_qif_osc = bp.dyn.DSRunner(qif_osc,
                                 monitors=['V', 'spike'],
                                 inputs=[('input', oscillatory_current, 't')],
                                 dt=0.01)
runner_lif_osc.run(duration_osc)
runner_qif_osc.run(duration_osc)

# 绘制输入电流和脉冲响应
fig, axes = plt.subplots(3, 1, figsize=(12, 8), sharex=True, gridspec_kw={'height_ratios': [1, 2, 1]})

# 子图1：输入电流
axes[0].plot(runner_lif_osc.mon.ts, I_osc, 'gray', linewidth=1.5, alpha=0.8)
axes[0].set_ylabel('Input Current (nA)')
axes[0].set_title('Oscillatory Input and Neural Responses')
axes[0].grid(True, alpha=0.3)

# 子图2：膜电位对比（叠加显示）
axes[1].plot(runner_lif_osc.mon.ts, runner_lif_osc.mon.V[:, 0], 'b-', alpha=0.7, label='LIF V', linewidth=1)
axes[1].plot(runner_qif_osc.mon.ts, runner_qif_osc.mon.V[:, 0], 'r-', alpha=0.7, label='QIF V', linewidth=1)
axes[1].set_ylabel('Membrane Potential (mV)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# 子图3：脉冲序列对比（点图）
lif_spike_t = runner_lif_osc.mon.ts[runner_lif_osc.mon.spike[:, 0]]
qif_spike_t = runner_qif_osc.mon.ts[runner_qif_osc.mon.spike[:, 0]]
axes[2].eventplot([lif_spike_t, qif_spike_t], colors=['blue', 'red'], lineoffsets=[0.9, 0.1], linelengths=0.7)
axes[2].set_xlabel('Time (ms)')
axes[2].set_ylabel('Spikes\n(Blue:LIF, Red:QIF)')
axes[2].set_yticks([0.9, 0.1])
axes[2].set_yticklabels(['LIF', 'QIF'])
axes[2].grid(True, alpha=0.3, axis='x')

plt.tight_layout()
plt.show()

# 分析脉冲相位锁定
def compute_phase_locking(spike_times, stim_freq, total_duration, ignore_transient=200):
    """计算脉冲相位相对于刺激正弦波的分布"""
    spike_times = spike_times[spike_times > ignore_transient]
    if len(spike_times) == 0:
        return None
    # 将发放时间转换为刺激周期内的相位 (0到2π)
    phases = (spike_times / 1000 * stim_freq) % 1.0 * 2 * np.pi
    return phases

lif_phases = compute_phase_locking(lif_spike_t, 5, duration_osc)
qif_phases = compute_phase_locking(qif_spike_t, 5, duration_osc)

if lif_phases is not None and qif_phases is not None:
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), subplot_kw=dict(projection='polar'))
    # LIF相位分布
    ax1.hist(lif_phases, bins=20, color='blue', alpha=0.6, density=True)
    ax1.set_title('LIF Phase Locking', pad=20)
    # QIF相位分布
    ax2.hist(qif_phases, bins=20, color='red', alpha=0.6, density=True)
    ax2.set_title('QIF Phase Locking', pad=20)
    plt.suptitle('Spike Phase Distribution Relative to Input Oscillation (5 Hz)')
    plt.tight_layout()
    plt.show()

在这个振荡输入场景下，你可能会发现QIF神经元的脉冲发放往往更紧密地锁定在输入正弦波的特定相位上（例如上升相），其相位分布直方图峰值更尖锐。而LIF神经元的发放可能更分散，或者锁定在另一个相位。这是因为QIF在阈值附近更高的增益使其对输入电流的瞬时值更敏感，从而能更精确地“捕捉”到输入波形的上升沿。

此外，当我们在输入电流中加入高斯白噪声时，QIF模型由于其非线性的阈值机制，有时会表现出与LIF不同的噪声鲁棒性或随机共振现象。例如，在次阈值电流附近，适量的噪声反而能帮助QIF神经元更规律地发放，而LIF可能只是产生完全随机的泊松式发放。探索这些差异需要更系统的噪声分析，但了解这一点有助于你在构建对噪声具有特定响应特性的网络时做出模型选择。

7. 关键差异五：计算开销与大规模网络仿真中的权衡

最后，我们无法回避工程现实：计算效率。LIF模型因其线性本质，计算极其高效。其微分方程甚至可以在没有数值积分器的情况下，利用指数函数的解析更新公式来求解，这在超大规模神经网络仿真（如百万神经元级别）中是至关重要的优势。

# LIF模型的解析更新步骤（伪代码风格，展示原理）
def lif_analytical_update(V_prev, I_input, dt, V_rest, R, tau, V_th, V_reset):
    """
    使用解析解更新LIF神经元。
    假设在时间步长dt内输入电流I_input恒定。
    """
    # 膜电位在没有阈值情况下的解析解
    V_inf = V_rest + R * I_input  # 稳态电位
    V_next = V_inf + (V_prev - V_inf) * np.exp(-dt / tau)

    # 检查并处理脉冲
    spike = V_next >= V_th
    V_next = np.where(spike, V_reset, V_next)
    return V_next, spike

这种解析更新比任何数值积分方法都快几个数量级。而QIF模型的非线性使其必须依赖数值积分器（如欧拉法、RK4），每一步都需要计算二次项，这无疑增加了计算负担。虽然对于现代计算机和GPU加速库（如BrainPy利用JAX后端）来说，仿真几千个QIF神经元仍然很快，但当规模扩大到数十万时，与LIF的效率差距就会变得显著。

因此，在选择模型时，你需要做一个根本性的权衡：

• 选择LIF，如果你：

  • 仿真的网络规模极大（>10^5神经元）。
  • 关注的是脉冲时序或群体平均频率，而非单个神经元的精确波形。
  • 输入电流动态相对平缓，或者噪声是主要驱动因素。
  • 需要极快的仿真速度进行参数扫描或学习算法迭代。

• 选择QIF，如果你：

  • 网络规模中等（<10^4神经元），可以承受适度的计算开销。
  • 研究内容涉及动作电位形态、阈值动力学或频率编码的细节。
  • 希望模型能通过参数（如a_0）方便地调节兴奋性类型。
  • 输入的动态范围较宽，需要神经元对快速变化的信号有更真实的响应。

在实际项目中，我经常采用混合策略：在网络的大部分区域使用LIF模型以保证效率，而在需要精细模拟的关键节点（如特定的输入层或读出层）使用QIF模型。BrainPy的灵活性和面向对象设计让这种混合建模变得非常直接。你可以让不同类型的神经元在同一个网络中无缝交互，只需确保它们都继承自相同的通信接口（如通过spike变量传递脉冲）。

通过这五个维度的并排对比，从电压轨迹、f-I曲线、参数动态、编码特性到计算考量，LIF和QIF不再是教科书上两个并列的公式，而是拥有了清晰应用边界和独特“性格”的工具。下次当你开始一个神经计算建模项目时，不妨先问自己：在这个项目中，是效率优先，还是逼真度优先？答案会自然地指引你做出选择。所有的对比代码都已准备好，调整参数，运行它们，亲眼看看这些差异如何在你的屏幕上展开，这才是建模工作真正开始的地方。