brainpy实战解析：构建动态平衡的兴奋-抑制神经网络模型

1. 从零开始：理解兴奋-抑制平衡网络

如果你对大脑如何工作感到好奇，想知道成千上万的神经元如何在混乱中产生有序的思维，那么兴奋-抑制平衡网络（E-I平衡网络）绝对是一个迷人的起点。这可不是什么遥不可及的纯理论，而是现代计算神经科学用来解释大脑皮层活动“背景噪音”的核心模型之一。简单来说，它试图回答一个关键问题：为什么大脑神经元在安静状态下也在噼里啪啦地随机放电，而不是一片死寂？答案就藏在“平衡”二字里。

想象一下一个热闹的鸡尾酒会。房间里既有大声分享趣事的人（兴奋性信号），也有劝大家冷静、别太吵的人（抑制性信号）。如果这两种声音势均力敌，整个房间就会维持在一个既活跃又不至于失控的嘈杂背景音中。E-I平衡网络就是大脑的“鸡尾酒会”，其中兴奋性神经元（E）和抑制性神经元（I）相互制衡，使得每个神经元接收到的总输入在零附近剧烈波动。这种动态平衡的结果，就是神经元表现出高度不规则、看似随机的脉冲发放模式，而这恰恰是大脑皮层在无明确任务时的典型状态。这种背景活动并非无用，它为信息处理提供了灵活的基底，让大脑能够对外部变化做出快速响应。

为什么要用BrainPy来构建它呢？因为BrainPy作为一个纯Python的神经动力学模拟框架，它的语法非常直观，把微分方程、神经元模型、网络连接用面向对象的方式包装得明明白白。你不需要先去啃一堆C++或CUDA代码，就能快速搭建并运行一个像E-I平衡网络这样中等规模的网络模型，亲眼看到动态平衡是如何涌现的。这对于研究者、学生，甚至是感兴趣的开发者来说，门槛降低了很多。接下来，我们就手把手，从原理到代码，构建一个活生生的动态平衡网络。

2. 模型基石：选择神经元与突触

构建网络的第一步，是选定用来搭建“砖块”的基本模型。在BrainPy里，这就像从丰富的工具箱里挑选合适的零件。对于E-I平衡网络，最经典、最常用的组合是漏电积分发放（LIF）神经元模型和指数型化学突触（Exponential Synapse）。这个组合在计算效率和生物合理性之间取得了很好的平衡。

LIF神经元可以说是计算神经科学的“Hello World”。它大大简化了真实的神经元，只保留最核心的特性：膜电位会像漏电的电容器一样衰减（“漏电”），接收的输入电流会累积（“积分”），当电位超过某个阈值时，神经元就会产生一个脉冲（“发放”），然后电位重置并进入短暂的不应期。虽然简单，但LIF模型能捕捉到神经元响应的基本时间特性和脉冲生成机制。在BrainPy中，调用 bp.neurons.LIF 就能轻松创建一个LIF神经元群体。你需要关心的关键参数包括：静息电位 V_rest、发放阈值 V_th、重置电位 V_reset、膜时间常数 tau 以及不应期 tau_ref。时间常数 tau 尤其重要，它决定了神经元对输入反应的快慢，值越大，反应越“慢”，记忆性越强。

指数型化学突触则模拟了神经递质释放后，在突触后膜引起的电流变化。它的特点是电流会瞬间上升到一个峰值，然后按指数方式衰减。BrainPy中的 bp.synapses.Exponential 类实现了这个模型。这里有一个关键概念叫突触后电流类型，我们通常使用 COBA（Conductance-Based） 模型。这意味着突触输入是以改变膜电导（ conductance，可以理解为离子通道的开放程度）的形式来影响神经元，其产生的电流大小依赖于当前膜电位与反转电位之间的差值。兴奋性突触的反转电位通常设为0 mV左右（如0 mV），而抑制性突触的反转电位则接近静息电位（如-80 mV）。这种设定更贴近生物实际，因为抑制性输入在膜电位高时（去极化）能产生更强的外向电流，效果更显著。

选择好模型后，网络的结构就清晰了：我们将创建两个LIF神经元群，一个代表兴奋性群体（E），一个代表抑制性群体（I）。它们之间通过四种连接交织在一起：E到E（E2E）、E到I（E2I）、I到E（I2E）、I到I（I2I）。正是这四类连接强度的微妙配比，最终决定了网络能否进入并维持那个有趣的动态平衡态。

3. 动手搭建：用BrainPy实现E-I网络

理论说再多，不如跑行代码。让我们打开编辑器，一步步把网络搭起来。首先，我们需要导入必要的库。BrainPy是核心，brainpy.math 是其后端数学计算模块（支持JAX、NumPy等），我们还会用到Matplotlib进行可视化。

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

接下来，我们定义网络类。在BrainPy中，复杂的动力学系统推荐通过继承 bp.dyn.Network 类来构建，这有助于组织和管理内部的各个组件。

class BalancedEINet(bp.dyn.Network):
    def __init__(self, num_E=3200, num_I=800, method='exp_auto', **kwargs):
        super().__init__(**kwargs)

        # 1. 定义神经元参数
        neuron_pars = dict(
            V_rest=-60.,   # 静息电位 (mV)
            V_th=-50.,     # 发放阈值 (mV)
            V_reset=-60.,  # 重置电位 (mV)
            tau=20.,       # 膜时间常数 (ms)
            tau_ref=5.,    # 绝对不应期 (ms)
            method=method
        )

        # 创建兴奋性和抑制性神经元群
        self.E = bp.neurons.LIF(num_E, **neuron_pars)
        self.I = bp.neurons.LIF(num_I, **neuron_pars)

        # 随机初始化膜电位，增加初始多样性
        self.E.V[:] = bm.random.randn(num_E) * 4. - 60.
        self.I.V[:] = bm.random.randn(num_I) * 4. - 60.

这里我设置了兴奋性神经元数量 num_E=3200，抑制性 num_I=800，比例是经典的4:1，这符合大脑皮层中E/I神经元的大致比例。膜电位用正态分布随机初始化，让网络启动时就有一些差异。

神经元建好了，现在用突触把它们连起来。连接的关键在于稀疏性和强度。我们使用 bp.conn.FixedProb 来建立随机连接，每个连接存在的概率（比如2%）很低，这模拟了大脑连接的稀疏特性。同时，抑制性连接的强度（g_max）通常要设得比兴奋性连接强得多，这是实现平衡的关键。

        # 2. 定义突触参数
        # 兴奋性突触参数：峰值电导较小，衰减较快
        E_syn_pars = dict(g_max=0.3, tau=5., method=method)
        # 抑制性突触参数：峰值电导很大，衰减较慢
        I_syn_pars = dict(g_max=3.7, tau=10., method=method)

        # 创建四种连接
        # E -> E: 兴奋到兴奋，使用兴奋性反转电位0mV
        self.E2E = bp.synapses.Exponential(
            self.E, self.E,
            conn=bp.conn.FixedProb(prob=0.02),
            output=bp.synouts.COBA(E=0.),  # 兴奋性反转电位
            **E_syn_pars
        )
        # E -> I: 兴奋到抑制
        self.E2I = bp.synapses.Exponential(
            self.E, self.I,
            conn=bp.conn.FixedProb(prob=0.02),
            output=bp.synouts.COBA(E=0.),
            **E_syn_pars
        )
        # I -> I: 抑制到抑制，使用抑制性反转电位-80mV
        self.I2I = bp.synapses.Exponential(
            self.I, self.I,
            conn=bp.conn.FixedProb(prob=0.02),
            output=bp.synouts.COBA(E=-80.),
            **I_syn_pars
        )
        # I -> E: 抑制到兴奋
        self.I2E = bp.synapses.Exponential(
            self.I, self.E,
            conn=bp.conn.FixedProb(prob=0.02),
            output=bp.synouts.COBA(E=-80.),
            **I_syn_pars
        )

注意看 g_max 的差异：兴奋性突触是0.3，而抑制性突触高达3.7。这是因为在平衡态下，虽然抑制性神经元数量少（只有兴奋性的1/4），但每个抑制性神经元需要发出更强的信号，才能抵消大量兴奋性神经元带来的影响，维持全局平衡。COBA 输出器中的 E 参数指定了反转电位，这是区分兴奋和抑制效果的核心。

4. 运行与可视化：观察动态平衡的涌现

网络类定义完毕，现在让我们运行它，看看会发生什么。我们使用 bp.dyn.DSRunner 这个模拟器，它可以方便地运行动力学系统并记录我们感兴趣的数据。

# 实例化网络
net = BalancedEINet(3200, 800)

# 创建运行器，监控脉冲、输入和部分神经元的膜电位
runner = bp.dyn.DSRunner(
    net,
    monitors=['E.spike', 'I.spike', 'E.input', 'I.input'], # 监控脉冲和总输入电流
    inputs=[('E.input', 12.), ('I.input', 12.)], # 给所有E和I神经元一个恒定的外部输入电流
    dt=0.1  # 积分步长，单位ms
)
# 运行200毫秒
runner.run(200.)

我们给所有神经元施加了一个12 pA的恒定外部电流。运行结束后，runner.mon 里就存储了所有监控数据。接下来，我们画图看看结果。通常我们会看两种图：脉冲 raster 图和平均发放率图。

# 定义一个绘制脉冲点阵图的函数
def plot_spike_raster(spikes, title, ax):
    """spikes是一个二维数组 (time_steps, num_neurons)"""
    # 找出所有发放脉冲的时间和对应的神经元索引
    t_idx, neuron_idx = np.where(spikes)
    # 将时间索引转换为实际时间（毫秒）
    times = t_idx * runner.dt
    ax.scatter(times, neuron_idx, s=0.5, c='k', alpha=0.6)
    ax.set_title(title)
    ax.set_ylabel('Neuron Index')
    ax.set_xlim(0, runner.total_time)

# 定义一个绘制平均发放率的函数
def plot_firing_rate(t, spikes, window_ms=5., ax):
    """计算滑动窗口内的平均发放率"""
    # 使用BrainPy内置函数计算发放率，窗口宽度为5ms
    rate = bp.measure.firing_rate(spikes, window=window_ms, dt=runner.dt)
    ax.plot(t, rate, lw=1.5)
    ax.set_ylabel('Firing Rate (Hz)')
    ax.set_xlabel('Time (ms)')
    ax.grid(True, linestyle=':', alpha=0.6)

# 创建画布
fig, axs = plt.subplots(2, 2, figsize=(14, 10),
                        gridspec_kw={'height_ratios': [3, 1]},
                        sharex='col')

# 绘制兴奋性神经元脉冲点阵图
plot_spike_raster(runner.mon['E.spike'], 'Excitatory Neurons Spike Raster', axs[0, 0])
# 绘制抑制性神经元脉冲点阵图
plot_spike_raster(runner.mon['I.spike'], 'Inhibitory Neurons Spike Raster', axs[0, 1])
# 绘制兴奋性神经元平均发放率
plot_firing_rate(runner.mon.ts, runner.mon['E.spike'], ax=axs[1, 0])
# 绘制抑制性神经元平均发放率
plot_firing_rate(runner.mon.ts, runner.mon['I.spike'], ax=axs[1, 1])

plt.tight_layout()
plt.show()

当你运行这段代码后，会看到四张子图。上面的两张点阵图里，每一个点代表一个神经元在某个时刻发放了一个脉冲。如果网络处于平衡态，你会看到这些点随机而均匀地分布在整个时间和神经元索引空间，没有明显的同步振荡或爆发式活动。下面的两张曲线图展示了整个群体平均发放率随时间的变化。在模拟开始的一小段时间（约20-50ms），由于所有神经元从相似的初始条件开始，并受到相同的外部电流冲击，你可能会看到一个初始的同步峰。但很快，网络内部复杂的E-I相互作用开始占据主导，发放率会下降并稳定在一个基线水平上下波动。这个波动的、非零的基线发放率，就是动态平衡态的标志！兴奋性和抑制性群体的发放率会大致相当，或者抑制性略高，这正是因为强抑制在压制着兴奋，达成了平衡。

5. 关键参数调优：让网络“平衡”起来

第一次运行很可能得不到完美的平衡态，可能看到神经元全部沉默（发放率为0）或者全部疯狂同步爆发。别担心，这很正常。E-I平衡网络的“平衡”非常依赖于几组关键参数的设置，调参是构建这个模型的核心实践环节。我们需要像一个调音师一样，仔细调整这些“旋钮”。

第一组关键参数：神经元参数。 其中最重要的是膜时间常数 tau 和发放阈值 V_th。tau 越大，神经元对输入的记忆越长，反应越慢，更容易整合输入产生发放；tau 越小，则越“健忘”，需要更强烈或更持续的输入才能触发脉冲。V_th 则直接决定了发放的难易程度。通常，为了让网络活跃起来，我们需要确保在给定的外部输入和连接强度下，神经元有合理的发放概率。如果网络总是沉默，可以尝试略微降低 V_th（例如从-50mV降到-52mV）或增加外部输入电流。

第二组，也是最核心的参数：突触连接强度 g_max。 这是我们调节平衡的主要手段。回顾一下代码，我们有四个 g_max：E2E.g_max, E2I.g_max, I2E.g_max, I2I.g_max。在经典设定中，我们通常让 I2E 和 I2I 的强度远大于 E2E 和 E2I（就像我们之前设置的3.7 vs 0.3）。其背后的原理是，为了实现每个神经元接收的净输入均值为零，需要满足一个近似条件：J_E * N_E * f_E ≈ J_I * N_I * f_I，其中J是连接强度，N是神经元数量，f是发放率。由于N_E通常是N_I的4倍，为了抵消更多的兴奋性输入，每个抑制性连接的强度J_I就必须设计得比J_E大得多（通常4-10倍）。如果抑制不够强，网络会陷入兴奋性爆发；如果抑制过强，网络会被完全压制而沉默。我个人的经验是，先固定兴奋性强度（如0.3），然后以2-4倍的比例调整抑制性强度，观察发放率的变化。

第三组参数：外部输入电流。 这个电流是驱动整个网络的“能量源”。没有它，网络会趋于静息。输入电流的大小直接决定了网络平衡态的平均发放率水平。你可以把它看作一个控制网络总体活跃度的旋钮。在我们的测试中，12 pA是一个常见的起始值。

第四组参数：连接概率和网络规模。 连接概率（我们设为0.02）影响了网络的稀疏性和平均每个神经元接收的连接数。概率太低，网络可能无法形成有效的相互作用；概率太高，计算量增大，且可能引入过高的相关性。网络规模（总神经元数）也影响显著。规模越大，统计特性越稳定，越容易观察到清晰的平衡态，但计算成本也越高。3200+800是一个在个人电脑上可以接受的中等规模。

调参时，建议你写一个简单的循环或参数扫描脚本，批量运行不同参数组合，并自动计算网络稳定后的平均发放率和变异系数（CV of ISI，脉冲间隔的变异系数，接近1表示发放接近泊松随机过程，是平衡态的一个特征）。通过系统性的参数扫描，你能更深刻地理解每个参数如何影响网络的宏观行为。

6. 探索网络特性：线性响应与快速跟踪

一个调好的E-I平衡网络不仅仅是产生随机脉冲，它还具有一些非常有趣且可能对计算有用的特性。我们可以设计实验来验证这些特性。

特性一：对外部刺激的线性响应。 在平衡态下，网络的平均发放率会随着外部输入电流的增加而近似线性地增加。我们可以设计一个阶梯状变化的输入电流，让网络在不同强度的电流下各运行一段时间，然后测量每个阶段的平均发放率。

import sklearn.linear_model

# 定义一系列输入电流强度和持续时间
current_levels = np.array([5., 10., 15., 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80.])
duration_per_level = 5000.  # 每个电流水平持续5000ms

# 构建分段恒定输入
input_current, total_duration = bp.inputs.constant_input(
    [(I, duration_per_level) for I in current_levels]
)

# 重新实例化网络并运行
net = BalancedEINet(3200, 800)
runner = bp.dyn.DSRunner(
    net,
    monitors=['E.spike', 'I.spike'],
    inputs=[('E.input', input_current, 'iter'), ('I.input', input_current, 'iter')],
    dt=0.1
)
runner.run(total_duration)

# 分析函数：计算并拟合发放率-电流关系
def analyze_response(spike_monitor, neuron_type, color):
    firing_rates = []
    # 对每个电流水平，取稳定后的时段计算平均发放率（避开切换的瞬态）
    for i, I in enumerate(current_levels):
        start_step = int((i * duration_per_level + 1000) / runner.dt)  # 跳过前1000ms瞬态
        end_step = start_step + int(3000 / runner.dt)  # 取3000ms稳定期
        spikes_in_epoch = spike_monitor[start_step:end_step]
        # 计算该时段内所有神经元的平均发放率 (Hz)
        mean_fr = np.mean(spikes_in_epoch) * (1000. / runner.dt)  # 将脉冲数/毫秒转换为Hz
        firing_rates.append(mean_fr)

    firing_rates = np.array(firing_rates)
    plt.scatter(current_levels, firing_rates, color=color, alpha=0.7, label=f'{neuron_type} data')

    # 线性回归
    regressor = sklearn.linear_model.LinearRegression()
    regressor.fit(current_levels.reshape(-1, 1), firing_rates.reshape(-1, 1))
    slope = regressor.coef_[0][0]
    intercept = regressor.intercept_[0]
    fit_x = np.array([0, 85])
    fit_y = slope * fit_x + intercept
    plt.plot(fit_x, fit_y, color=color, lw=2, label=f'{neuron_type} fit (slope={slope:.4f})')
    return slope, intercept

plt.figure(figsize=(10, 6))
slope_E, intercept_E = analyze_response(runner.mon['E.spike'], 'Excitatory', 'red')
slope_I, intercept_I = analyze_response(runner.mon['I.spike'], 'Inhibitory', 'blue')
plt.xlabel('External Input Current (pA)')
plt.ylabel('Population Firing Rate (Hz)')
plt.legend()
plt.grid(True, linestyle=':', alpha=0.5)
plt.title('Linear Response of E-I Balanced Network')
plt.show()

运行这段代码，你应该能得到两条漂亮的、近似直线的散点拟合图。这表明网络将输入强度线性地编码在了群体发放率中。这种线性响应特性使得网络可以作为信息传递的可靠载体。

特性二：快速跟踪动态刺激。 与单个LIF神经元相比，E-I平衡网络能更快地响应输入的变化。单个神经元由于膜电容的存在，其膜电位变化有惯性，对快速变化的信号响应会滞后和滤波。但在平衡网络中，由于大量神经元处于阈值附近，任何微小的输入变化都能立即导致一部分神经元提前或推迟发放，从而将输入变化迅速地“广播”到整个网络。你可以尝试输入一个方波或正弦波电流，分别观察单个神经元和整个网络的发放率随时间的变化曲线，会发现网络的反应更“敏锐”，跟随性更好。这个特性被认为是大脑能够快速处理感官信息的基础之一。

7. 常见问题与调试心得

在构建和调试E-I平衡网络的过程中，我踩过不少坑，也总结出一些经验。如果你运行代码时遇到了问题，可以对照下面几点检查。

问题一：网络完全沉默（发放率为0）。 这是最常见的问题。首先，检查外部输入电流是否足够大。12 pA只是一个参考值，取决于你的神经元参数（尤其是阈值V_th）。尝试逐步增大输入电流。其次，检查抑制性连接强度 I2E.g_max 和 I2I.g_max 是否过大。过强的抑制会把所有活动都压死。可以尝试将其调小。最后，确认连接概率和网络规模不是太小，要确保神经元之间有足够的相互作用。

问题二：网络同步爆发（Epileptic bursting）。 表现为所有神经元几乎在同一时间点集体放电，在raster图上形成垂直的条纹，发放率曲线出现周期性的尖峰。这通常意味着兴奋性过强或抑制性过弱。你需要增强抑制性连接的强度（增大 I2E.g_max 和 I2I.g_max），或者减弱兴奋性连接的强度（减小 E2E.g_max 和 E2I.g_max）。另一个可能的原因是神经元的不应期 tau_ref 设得太短，导致神经元在不应期结束后立即再次发放，加剧了同步。可以适当增加不应期。

问题三：发放率不稳定，漂移严重。 网络在运行一段时间后，平均发放率持续上升或下降，而不是围绕一个均值波动。这可能是因为网络参数没有精确满足平衡条件，存在微小的净兴奋或净抑制。你需要更精细地调整E和I的连接强度比例。此外，确保模拟时间足够长，网络已经度过了初始瞬态，进入了稳定状态。计算发放率时，要避开模拟开始的那段瞬态时间。

问题四：计算速度太慢。 BrainPy默认使用JAX后端，在CPU上运行大规模网络确实会慢。有几种加速方法：1) 使用GPU。如果你有NVIDIA GPU，确保安装了JAX的GPU版本，BrainPy会自动利用。速度会有数量级的提升。2) 减小网络规模进行原型调试，参数调好后再放大。3) 适当增大积分步长 dt（比如从0.1ms增加到0.2ms），但这可能会影响数值精度和动力学细节。

调试时，我习惯把关键变量的时间序列画出来看，比如 runner.mon[‘E.input’] 和 runner.mon[‘I.input’]，观察单个神经元接收的总输入电流是否在零附近大幅波动。这才是动态平衡最直接的体现。多试几次参数，观察网络行为如何随之变化，这个过程本身就能让你对平衡网络动力学的理解加深不少。记住，找到一个稳定的平衡态参数集，本身就是一项重要的成果。