简化神经元模型7 -- GIF Model

Hindmarsh-Rose 模型

0. 写在前面

前面介绍了:
Hodgkin-Huxley Model
简化神经元模型1 – LIF Model
简化神经元模型2 – QIF Model
简化神经元模型3 – ExpIF Model
简化神经元模型4 – AdEx Model
简化神经元模型5 – Izhikevich Model
简化神经元模型6 – Hindmarsh-Rose Model

下面继续介绍简化神经元模型，本篇博客主要是介绍泛化整合发放模型（Generalized Integrate-and-Fire model，GIF model），并用BrainPy实现模型的定义、模拟和分析。

GIF 模型也叫 Mihala-Niebur 模型，它由 Mihala和 Niebur 于 2009 年提出。人们之所以叫它泛化整合发放模型，是因为该模型具有相当广泛的表征能力，能通过参数的调节产生几十种发放模式。

资料参考： 《神经计算建模实战》吴思等著

1. GIF 模型的定义

GIF 模型的数学表达形式如下：
$\tau \frac{dV}{dt} = - (V - V_{\text{rest}}) + R \sum_{j} I_j + RI \tag{1}$ $\frac{d\Theta}{dt} = a (V - V_{\text{rest}}) - b (\Theta - \Theta_{\infty}) \tag{2}$ $\frac{dI_j}{dt} = -k_j I_j, \quad j = 1, 2, ..., n \tag{3}$ $\text{if } V > \Theta, \quad I_j \leftarrow R_j I_j + A_j, \quad V \leftarrow V_{\text{reset}}, \quad \Theta \leftarrow \max(\Theta_{\text{reset}}, \Theta) \tag{4}$ 该模型一共包含 $n + 2$ 个变量，其中第一个变量为膜电位 $V$ ，第二个变量为膜电位阈值 $\Theta$ ，之后 $n$ 个变量都被称为内部电流（internal currents），也被称为脉冲诱导电流（spike-induced currents）。方程 (3) 表示每个内部电流 $I_j$ 都以速率 $k_j$ 衰减。方程 (1) 在 LIF 模型的膜电位微分方程的基础上增加了内部电流 $\sum_{j} I_j$ 对膜电位的影响，它们在形式上与外部电流 $I$ 相同。方程 (2) 表明，在 GIF 模型中，膜电位阈值不再是一个常数，而成为了一个变量 $\Theta$ 。方程 (2) 右侧的第一项描述了膜电位的高低（以静息膜电位 $V_{\text{rest}}$ 为基准）对 $\Theta$ 的影响，第二项描述了 $\Theta$ 以 $\Theta_{\infty}$ 为基准、以速率 $b$ 衰减，其中 $\Theta_{\infty}$ 是没有外部输入的情况下时间趋于无穷时 $\Theta$ 的稳定值， $\Theta_{\infty}$ 之于 $\Theta$ 可类比 $V_{\text{rest}}$ 之于 $V$ 。

当神经元发放动作电位时，除了膜电位 $V$ 被重置为 $V_{\text{reset}}$ 外，膜电位阈值 $\Theta$ 如果低于设定的重置值 $\Theta_{\text{reset}}$ 则被重置为 $\Theta_{\text{reset}}$ ，内部电流 $I_j$ 则按照 $R_j I_j + A_j$ 的规则更新，其中 $R_j$ 、 $A_j$ 是自由参数。在参数设置上，要求 $\Theta_{\text{reset}} > V_{\text{reset}}$ 。

GIF 模型在 $n = 2$ 时就已经能够模拟几十种发放模式。下面以两个内部电流为例构建 GIF 模型:

import brainpy as bp
import brainpy.math as bmclass GIF(bp.dyn.NeuDyn):def __init__(self, size, V_rest=-70., V_reset=-70., theta_inf=-50., theta_reset=-60.,R=20., tau=20., a=0., b=0.01, k1=0.2, k2=0.02, R1=0., R2=1., A1=0.,A2=0., **kwargs):# 初始化父类时计算了self.num供下文使用super(GIF, self).__init__(size=size, **kwargs)# 初始化参数self.V_rest = V_restself.V_reset = V_resetself.theta_inf = theta_infself.theta_reset = theta_resetself.R = Rself.tau = tauself.a = aself.b = bself.k1 = k1self.k2 = k2self.R1 = R1self.R2 = R2self.A1 = A1self.A2 = A2# 初始化变量self.V = bm.Variable(bm.zeros(self.num) + V_reset)self.theta = bm.Variable(bm.ones(self.num) * theta_inf)self.input = bm.Variable(bm.zeros(self.num))self.spike = bm.Variable(bm.zeros(self.num, dtype=bool))self.I1 = bm.Variable(bm.zeros(self.num))self.I2 = bm.Variable(bm.zeros(self.num))# 定义积分器self.integral = bp.odeint(f=self.derivative, method='exp_auto')def dI1(self, I1, t):return - self.k1 * I1def dI2(self, I2, t):return - self.k2 * I2def dVth(self, V_th, t, V):return self.a * (V - self.V_rest) - self.b * (V_th - self.theta_inf)def dV(self, V, t, I1, I2, Iext):return (- (V - self.V_rest) + self.R * Iext + self.R * I1 + self.R * I2) / self.tau# 将所有微分方程联合为一个，以便同时积分@propertydef derivative(self):return bp.JointEq([self.dI1, self.dI2, self.dVth, self.dV])def update(self,):I1, I2, V_th, V = self.integral(self.I1, self.I2, self.theta, self.V, bp.share['t'],self.input, bp.share['dt'])  # 更新变量I1, I2, Vspike = self.theta <= V  # 将大于阈值的神经元标记为发放了脉冲V = bm.where(spike, self.V_reset, V)  # 将发放了脉冲的神经元V置为V_reset，其余赋值为更新后的VI1 = bm.where(spike, self.R1 * I1 + self.A1, I1)  # 按照公式更新发放了脉冲的神经元的I1I2 = bm.where(spike, self.R2 * I2 + self.A2, I2)  # 按照公式更新发放了脉冲的神经元的I2reset_th = bm.logical_and(V_th < self.theta_reset, spike)  # 判断哪些神经元的V_th需要重置V_th = bm.where(reset_th, self.theta_reset, V_th)  # 将需要重置的神经元V_th重置为theta_reset# 将更新后的结果赋值给self.*self.spike.value = spikeself.I1.value = I1self.I2.value = I2self.theta.value = V_thself.V.value = Vself.input[:] = 0.  # 重置外界输入

图3.30展示了 GIF 模型产生的二十种发放模式（可视化代码放在最后），各个模式下 $V$ 和 $\Theta$ 随时间的变化以及外部电流输入 $I$ 的设置都展现在图中。
在这里插入图片描述
所有发放模式共享的参数包括： $b = 0.01$ , $\tau = 50$ , $k_1 = 0.2$ , $k_2 = 0.02$ , $V_{\text{rest}} = -70$ , $\Theta_{\infty} = -50$ , $R_1 = 0$ , $R_2 = 1$ , $V_{\text{reset}} = -70$ , $\Theta_{\text{reset}} = -60$ 。

各种模式对应的其他可变参数请参看下表 3.3。
在这里插入图片描述

2. GIF 模型的动力学分析

虽然都是以 LIF 模型为基础进行改进，但相比之前提及的其他简化模型，GIF 模型在建模思路上有着很大的不同。之前的模型都力求神经元膜电位的阈下变化更加接近真实情况，并通过在微分方程中加入非线性项达到增加模型表征能力的目的；但是 GIF 模型中每个微分方程都是线性的，这不仅使得所有方程都可以求得解析解，而且还降低了计算的复杂度。而且，这些微分方程构成的方程组的系数矩阵还是一个三角矩阵，这意味着变量之间的依赖存在顺次关系：方程 (3) 中的各个 $I_j$ 可独立求解；而 $V$ 又只依赖于各个 $I_j$ ，于是解得 $I_j$ 后就可以求解 $V$ ；最后 $\Theta$ 也就能随之求解了。我们可以按照特定的顺序依次求解这些方程，从而轻易地获取各个变量随时间变化的函数。

既然 GIF 模型的微分方程都是线性的，那它强大、广泛的表征能力从何而来呢？下面从两个角度分别讨论。

2.1. 角度一：GIF 模型的变量较多

首先，相比其他模型，GIF 模型的变量更多，即使取 $n = 2$ 也存在 4 个变量。线性微分方程决定了神经元的阈下膜电位变化曲线较为单一，于是 GIF 模型将膜电位阈值 $\Theta$ 设置为一个变量，这样在不同状态下，即使神经元的膜电位值相同，它们的行为（是否发放动作电位）也可能不同。例如在适应性发放（spike frequency adaptation, 图3.30C）模式中，神经元发放频率的降低就是因为膜电位阈值 $\Theta$ 在不断升高，膜电位 $V$ 达到阈值的时间也随之增加。在另一种适应（accommodation, 图3.30E）模式中，将电流直接升高到基强度以上可以激发动作电位，但阶梯式地将电流升高到同一值却不能激发动作电位，这也是因为 $\Theta$ 在电流阶梯式增加的过程中不断提高，因此基强度电流也随之改变，原先的电流强度已经小于新的基强度，故不能产生动作电位。

此外，内部电流 $I_j$ 也能在一定程度上增加膜电位变化的多样性。图3.31（a）展示了 GIF 模型产生发放后电位（afterpotentials）模式时 $I_1$ 、 $I_2$ 随时间的变化。其中参数设置为 $a = 0.005$ , $A_1 = 5$ , $A_2 = -0.3$ 。由于 $I_1$ 的更新幅度比 $I_2$ 大，但动力学变化比 $I_2$ 快，两者叠加后会形成一个先正后负再趋于 0 的函数（类似两个高斯函数叠加形成的墨西哥草帽的右侧形状），这也使得膜电位有一个先去极化后超极化再趋于稳定的电位变化模式。类似的原理也会让神经元对外部刺激的频率有特定的偏好（频率偏好，preferred frequency，图3.31（b），参数设置为 $a = 0.005$ , $A_1 = -3$ , $A_2 = 0.5$ ）。第一次脉冲电流的刺激（图中有两次，分别在第 0 ms 和第 400 ms）让神经元发放动作电位，第二次的刺激如果紧接第一次刺激，在 $I_1$ 抑制作用较强的时候施加，则神经元无法产生新的动作电位；如果间隔一段时间后再施加，此时 $I_1$ 已基本衰减为 0，但动力学变化更慢的 $I_2$ 的促进作用仍然存在，因此新的脉冲电流无需第一次那么强就能够让神经元产生新的发放。在这里插入图片描述

2.2. 角度二：GIF 模型的更新规则较为复杂多

GIF 模型能模拟多种发放模式的第二个原因是它在发放动作电位后的更新规则较为复杂。更新规则会对所有变量都产生影响，且除了 $V$ 以外都不是直接被重置为一个恒定值—— $\Theta$ 的重置值取决于当前 $\Theta$ 和 $\Theta_{\text{reset}}$ 的大小，而 $I_j$ 则会经历一个线性变换。从参数数量也能看出，更新规则具有很大的调整空间，相应地，神经元也能具有多种多样的发放模式。

例如，簇发放（bursting/tonic bursting，图3.32（a），参数设置为 $a = 0.005$ , $A_1 = 10$ , $A_2 = -0.6$ ）在 GIF 模型中就是通过内部电流的更新实现的。注意参数设定中 $R_1 = 0$ , $R_2 \neq 0$ ，这表示 $I_1$ 的更新是重置为固定值，而 $I_2$ 的更新是累加的。在集中发放阶段，由于 $A_1$ 很大，每次更新时内部电流 $I_1$ 都会显著增大，系统接受到一个很大的电流输入，因此膜电位快速上升并达到阈值，形成密集的脉冲发放。与此同时，每次更新时 $I_2$ 的（绝对值的）增量不大，但由于 $I_2$ 的时间常数很大，衰减很慢，因此在多次高频的更新中增量不断积累，抑制性电流强度不断增大，最终让神经元停止发放，进入间隔期。进入间隔期后， $I_2$ 逐渐衰减，当衰减得足够小之后，神经元又能再次发放。

除了内部电流的更新规则，膜电位的更新规则也使得 GIF 模型具有多样性。在反弹发放（rebound spike，图3.32（b），参数设置为 $a = 0.005$ , $A_1 = 0$ , $A_2 = 0$ ）中，神经元先接收一个抑制性电流，在抑制性电流撤走后，即使没有正向电流，神经元也能产发放。为模拟这一现象，GIF 模型首先让膜电位阈值 $\Theta$ 在抑制性电流输入阶段随膜电位一起降低，这样当外部电流输入变为 0 后，时间常数更小的膜电位快速上升，达到膜电位阈值，由此产生了脉冲。如何让神经元在发放了一次脉冲后就停止呢？由于当前膜电位阈值 $\Theta$ 小于阈值重置值 $\Theta_{\text{reset}}$ ，更新规则将 $\Theta$ 重置为 $\Theta_{\text{reset}}$ ，膜电位阈值升高到膜电位以上，因此神经元停止发放。如果没有膜电位阈值的更新规则，则重置后 $V$ 将超过 $\Theta$ ，这是生物学不合理的。在这里插入图片描述

综上，我们知道了 GIF 模型广泛的表征能力来源于多变量和复杂的更新规则，它们弥补了线性微分方程在表征能力上的不足。另外，从图3.32中我们可以发现，内部电流 $I_j$ 对应的参数在簇发放中起关键作用，膜电位阈值 $\Theta$ 对应的参数在反弹发放中起关键作用，说明不同的变量及其对应的参数可以分别导致 GIF 模型产生一些发放模式。进一步而言，如果将簇发放中 $I_j$ 相关的参数和反弹发放中 $\Theta$ 相关的参数整合到一起，就可以得到一种新的发放模式：反弹簇发放（rebound bursting，图3.30 O）。通过这样的参数组合，GIF 模型能模拟的发放模式就变得更多了。

虽然 GIF 模型能产生各种各样的发放模式，但我们也应该意识到它的局限性：为了减少计算复杂度，它通过增加变量、复杂化更新规则的方法间接实现了神经元中的一些非线性变化，但这些设定（例如膜电位阈值的更新规则）是否生物合理却有待商榷。即使 GIF 模型能复现许多发放模式，但膜电位随时间的变化曲线却与真实情况并不相似，因此它并不一定能替代其他模型成为一个包揽一切的泛化模型。事实上，GIF 模型在工程上的贡献大于在神经科学上的贡献——模型中的微分方程和更新规则都可以通过简单的逻辑实现，这为电子神经元的构造提供了思路。

BrainPy_span_138">3. 分析过程所用到的一系列 BrainPy 代码

将上述一系列定义模型代码和可视化代码全部给出：

import brainpy as bp
import brainpy.math as bm
import matplotlib.pyplot as pltplt.rcParams.update({"font.size": 15})
plt.rcParams['font.sans-serif'] = ['Times New Roman']class GIF(bp.dyn.NeuDyn):def __init__(self, size, V_rest=-70., V_reset=-70., theta_inf=-50., theta_reset=-60.,R=20., tau=20., a=0., b=0.01, k1=0.2, k2=0.02, R1=0., R2=1., A1=0.,A2=0., **kwargs):# 初始化父类时计算了self.num供下文使用super(GIF, self).__init__(size=size, **kwargs)# 初始化参数self.V_rest = V_restself.V_reset = V_resetself.theta_inf = theta_infself.theta_reset = theta_resetself.R = Rself.tau = tauself.a = aself.b = bself.k1 = k1self.k2 = k2self.R1 = R1self.R2 = R2self.A1 = A1self.A2 = A2# 初始化变量self.V = bm.Variable(bm.zeros(self.num) + V_reset)self.theta = bm.Variable(bm.ones(self.num) * theta_inf)self.input = bm.Variable(bm.zeros(self.num))self.spike = bm.Variable(bm.zeros(self.num, dtype=bool))self.I1 = bm.Variable(bm.zeros(self.num))self.I2 = bm.Variable(bm.zeros(self.num))# 定义积分器self.integral = bp.odeint(f=self.derivative, method='exp_auto')def dI1(self, I1, t):return - self.k1 * I1def dI2(self, I2, t):return - self.k2 * I2def dVth(self, V_th, t, V):return self.a * (V - self.V_rest) - self.b * (V_th - self.theta_inf)def dV(self, V, t, I1, I2, Iext):return (- (V - self.V_rest) + self.R * Iext + self.R * I1 + self.R * I2) / self.tau# 将所有微分方程联合为一个，以便同时积分@propertydef derivative(self):return bp.JointEq([self.dI1, self.dI2, self.dVth, self.dV])def update(self,):I1, I2, V_th, V = self.integral(self.I1, self.I2, self.theta, self.V, bp.share['t'],self.input, bp.share['dt'])  # 更新变量I1, I2, Vspike = self.theta <= V  # 将大于阈值的神经元标记为发放了脉冲V = bm.where(spike, self.V_reset, V)  # 将发放了脉冲的神经元V置为V_reset，其余赋值为更新后的VI1 = bm.where(spike, self.R1 * I1 + self.A1, I1)  # 按照公式更新发放了脉冲的神经元的I1I2 = bm.where(spike, self.R2 * I2 + self.A2, I2)  # 按照公式更新发放了脉冲的神经元的I2reset_th = bm.logical_and(V_th < self.theta_reset, spike)  # 判断哪些神经元的V_th需要重置V_th = bm.where(reset_th, self.theta_reset, V_th)  # 将需要重置的神经元V_th重置为theta_reset# 将更新后的结果赋值给self.*self.spike.value = spikeself.I1.value = I1self.I2.value = I2self.theta.value = V_thself.V.value = Vself.input[:] = 0.  # 重置外界输入def run_GIF():fig, gs = bp.visualize.get_figure(1, 2, 4, 6)# 模拟相位脉冲（phasic spiking）group = GIF(10, a=0.005, A1=0., A2=0.)runner = bp.DSRunner(group, monitors=['V', 'theta'], inputs=('input', 1.5), dt=0.01)runner(500)fig.add_subplot(gs[0, 0])bp.visualize.line_plot(runner.mon.ts, runner.mon.V, legend='V', zorder=10, show=False)bp.visualize.line_plot(runner.mon.ts, runner.mon.theta, legend='theta',title='phasic spiking', show=False)# 模拟超极化爆发（hyperpolarization-induced bursting）group = GIF(10, a=0.03, A1=10., A2=-0.6)runner = bp.DSRunner(group, monitors=['V', 'theta'], inputs=('input', -1), dt=0.01)runner(500)fig.add_subplot(gs[0, 1])bp.visualize.line_plot(runner.mon.ts, runner.mon.V, legend='V', zorder=10, show=False)bp.visualize.line_plot(runner.mon.ts, runner.mon.theta, legend='theta',title='hyperpolarization-induced bursting', show=True)def plot_gallery():def _run(ax1, model, duration, I_ext, title=''):runner = bp.DSRunner(model,inputs=('input', I_ext, 'iter'),monitors=['V', 'theta'])runner.run(duration)ts = runner.mon.tsax1.plot(ts, runner.mon.V[:, 0], label=r'$V$', linestyle='-')ax1.plot(ts, runner.mon.theta[:, 0], label=r'$\theta$', linestyle='--')ax1.set_ylabel('Potential (mV)')ax1.set_xlabel(r'$t$ (ms)')ax1.set_xlim(-0.1, ts[-1] + 0.1)if title: plt.title(title)# plt.legend()ax1.spines['top'].set_visible(False)ax1.spines['right'].set_visible(False)row, col = 7, 3fig, gs = bp.visualize.get_figure(row, col, 3, 5)# Tonic Spikingax = fig.add_subplot(gs[0 // col, 0 % col])Iext, duration = bp.inputs.section_input((1.5,), (200.,), return_length=True)_run(ax, GIF(1), duration, Iext, 'A. Tonic Spiking')# Class 1 Excitabilityax = fig.add_subplot(gs[1 // col, 1 % col])Iext, duration = bp.inputs.section_input([1. + 1e-6], [500.], return_length=True)_run(ax, GIF(1), duration, Iext, 'B. Class 1 Excitability')# Spike Frequency Adaptationax = fig.add_subplot(gs[2 // col, 2 % col])Iext, duration = bp.inputs.section_input([2.], [200.], return_length=True)_run(ax, GIF(1, a=0.005), duration, Iext, 'C. Spike Frequency Adaptation')# Phasic Spikingax = fig.add_subplot(gs[3 // col, 3 % col])Iext, duration = bp.inputs.section_input([1.5], [500.], return_length=True)_run(ax, GIF(1, a=0.005), duration, Iext, 'D. Phasic Spiking')# Accommodationax = fig.add_subplot(gs[4 // col, 4 % col])Iext, duration = bp.inputs.section_input([1.5, 0., 0.5, 1., 1.5, 0.],[100., 500., 100., 100., 100., 100.],return_length=True)_run(ax, GIF(1, a=0.005), duration, Iext, 'E. Accommodation')# Threshold Variabilityax = fig.add_subplot(gs[5 // col, 5 % col])Iext, duration = bp.inputs.section_input([1.5, 0., -1.5, 0., 1.5, 0.],[20., 180., 20., 20., 20., 140.],return_length=True)_run(ax, GIF(1, a=0.005), duration, Iext, 'F. Threshold Variability')# Rebound Spikingax = fig.add_subplot(gs[6 // col, 6 % col])Iext, duration = bp.inputs.section_input([0., -3.5, 0.], [50., 750., 200.], return_length=True)_run(ax, GIF(1, a=0.005), duration, Iext, 'G. Rebound Spiking')# Class 2 Excitabilityax = fig.add_subplot(gs[7 // col, 7 % col])Iext, duration = bp.inputs.section_input([2 * (1. + 1e-6)], [200.], return_length=True)neu = GIF(1, a=0.005)neu.theta[:] = -30._run(ax, neu, duration, Iext, 'H. Class 2 Excitability')# Integratorax = fig.add_subplot(gs[8 // col, 8 % col])Iext, duration = bp.inputs.section_input([1.5, 0., 1.5, 0., 1.5, 0., 1.5, 0.],[20., 10., 20., 250., 20., 30., 20., 30.],return_length=True)_run(ax, GIF(1, a=0.005), duration, Iext, 'I. Integrator')# Input Bistabilityax = fig.add_subplot(gs[9 // col, 9 % col])Iext, duration = bp.inputs.section_input([1.5, 1.7, 1.5, 1.7],[100., 400., 100., 400.],return_length=True)_run(ax, GIF(1, a=0.005), duration, Iext, 'J. Input Bistability')# Hyperpolarization-induced Spikingax = fig.add_subplot(gs[10 // col, 10 % col])Iext, duration = bp.inputs.section_input([-1.], [400.], return_length=True)neu = GIF(1, theta_reset=-60., theta_inf=-120.)neu.theta[:] = -50._run(ax, neu, duration, Iext, 'K. Hyperpolarization-induced Spiking')# Hyperpolarization-induced Burstingax = fig.add_subplot(gs[11 // col, 11 % col])Iext, duration = bp.inputs.section_input([-1.], [400.], return_length=True)neu = GIF(1, theta_reset=-60., theta_inf=-120., A1=10., A2=-0.6)neu.theta[:] = -50._run(ax, neu, duration, Iext, 'L. Hyperpolarization-induced Bursting')# Tonic Burstingax = fig.add_subplot(gs[12 // col, 12 % col])Iext, duration = bp.inputs.section_input([2.], [500.], return_length=True)_run(ax, GIF(1, a=0.005, A1=10., A2=-0.6), duration, Iext, 'M. Tonic Bursting')# Phasic Burstingax = fig.add_subplot(gs[13 // col, 13 % col])Iext, duration = bp.inputs.section_input([1.5], [500.], return_length=True)neu = GIF(1, a=0.005, A1=10., A2=-0.6)_run(ax, neu, duration, Iext, 'N. Phasic Bursting')# Rebound Burstingax = fig.add_subplot(gs[14 // col, 14 % col])Iext, duration = bp.inputs.section_input([0., -3.5, 0.],[100., 500., 400.],return_length=True)_run(ax, GIF(1, a=0.005, A1=10., A2=-0.6), duration, Iext, 'O. Rebound Bursting')# Mixed Modeax = fig.add_subplot(gs[15 // col, 15 % col])Iext, duration = bp.inputs.section_input([2.], [500.], return_length=True)_run(ax, GIF(1, a=0.005, A1=5., A2=-0.3), duration, Iext, 'P. Mixed Mode')# Afterpotentialsax = fig.add_subplot(gs[16 // col, 16 % col])Iext, duration = bp.inputs.section_input((2., 0.), [15., 185.], return_length=True)_run(ax, GIF(1, a=0.005, A1=5., A2=-0.3), duration, Iext, 'Q. Afterpotentials')# Basal Bistabilityax = fig.add_subplot(gs[17 // col, 17 % col])Iext, duration = bp.inputs.section_input([5., 0., 5., 0.],[10., 90., 10., 90.],return_length=True)_run(ax, GIF(1, A1=8., A2=-0.1), duration, Iext, 'R. Basal Bistability')# Preferred Frequencyax = fig.add_subplot(gs[18 // col, 18 % col])Iext, duration = bp.inputs.section_input([5., 0., 4., 0., 5., 0., 4., 0.],[10., 10., 10., 370., 10., 90., 10., 290.],return_length=True)_run(ax, GIF(1, a=0.005, A1=-3., A2=0.5), duration, Iext, 'S. Preferred Frequency')# Spike Latencyax = fig.add_subplot(gs[19 // col, 19 % col])Iext, duration = bp.inputs.section_input([8., 0.], [2., 48.], return_length=True)_run(ax, GIF(1, a=-0.08), duration, Iext, 'T. Spike Latency')# plt.savefig('GIF_gallery.pdf', dpi=500, transparent=True)plt.show()def detailed_running():def _run(model, duration, I_ext, title=''):runner = bp.DSRunner(model,inputs=('input', I_ext, 'iter'),monitors=['V', 'theta', 'I1', 'I2'])runner.run(duration)fig, gs = bp.visualize.get_figure(5, 1, 0.9, 6.)ax1 = fig.add_subplot(gs[:3, 0])ax1.plot(runner.mon.ts, runner.mon.V[:, 0], label=r'$V$', linestyle='-')ax1.plot(runner.mon.ts, runner.mon.theta[:, 0], label=r'$\theta$', linestyle='--')ax1.set_ylabel('Potential (mV)')ax1.set_xlim(-0.1, runner.mon.ts[-1] + 0.1)plt.legend()ax1.spines['top'].set_visible(False)ax1.spines['right'].set_visible(False)ax1 = fig.add_subplot(gs[3:, 0])ax1.plot(runner.mon.ts, runner.mon.I1, label=r'$I_1$', linestyle='-')ax1.plot(runner.mon.ts, runner.mon.I2, label=r'$I_2$', linestyle='--')ax1.plot(runner.mon.ts, bm.as_numpy(I_ext), label=r'$I$', linestyle='dotted')ax1.set_xlim(-0.1, runner.mon.ts[-1] + 0.1)ax1.set_xlabel(r'$t$ (ms)')ax1.set_ylabel('Current')plt.legend()ax1.spines['top'].set_visible(False)ax1.spines['right'].set_visible(False)if title:plt.savefig(f'GIF_{title.replace(" ", "-")}.pdf', transparent=True, dpi=500)Iext, duration = bp.inputs.section_input((2., 0.), [15., 185.], return_length=True)_run(GIF(1, a=0.005, A1=5., A2=-0.3), duration, Iext, 'Afterpotentials')Iext, duration = bp.inputs.section_input([5., 0., 4., 0., 5., 0., 4., 0.], [10., 10., 10., 370., 10., 90., 10., 290.], return_length=True)_run(GIF(1, a=0.005, A1=-3., A2=0.5), duration, Iext, 'Preferred Frequency')Iext, duration = bp.inputs.section_input([2.], [500.], return_length=True)_run(GIF(1, a=0.005, A1=10., A2=-0.6), duration, Iext, 'Tonic Bursting')Iext, duration = bp.inputs.section_input([0., -3.5, 0.], [50., 750., 200.], return_length=True)_run(GIF(1, a=0.005), duration, Iext, 'Rebound Spiking')plt.show()if __name__ == '__main__':run_GIF()plot_gallery()detailed_running()