Hindmarsh-Rose 模型

目录

0. 写在前面
1. Hindmarsh-Rose 模型的定义
2. Hindmarsh-Rose 模型簇发放的动力学机制
3. Hindmarsh-Rose 模型的其他发放模式
4. 分析过程所用到的一系列 BrainPy 代码

0. 写在前面

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

下面继续介绍简化神经元模型，本篇博客主要是介绍 Hindmarsh-Rose 模型，并用BrainPy实现模型的定义、模拟和分析。

之前我们讨论了 AdEx 模型和 Izhikevich 模型，它们在原先单变量模型的基础上增加了一个变量，使得模型的动力学表征能力大大增强，能够模拟多种发放模式。不过，以上所有简化模型都需要在神经元发放动作电位后显式地将膜电位 $V$ 重置，这破坏了变量的连续性。接下来我们将介绍一种神经元模型，它通过三个变量刻画神经元的动力学特征，且不再需要显示地重置膜电位。

HindMarsh 和 Rose 于 1984 年提出了 Hindmarsh-Rose 模型。他们在前人研究的双变量模型基础上增加了一个变量来模拟缓慢流入细胞的电流，以丰富神经元的动力学性质。Hindmarsh 和 Rose 最初希望利用该模型来解释簇发放现象，但他们发现这一模型也能产生适应、反弹簇发放等多种模式，因此可以将其视为一个较为通用的神经元模型。

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

1. Hindmarsh-Rose 模型的定义

Hindmarsh-Rose 模型存在三个变量，分别为 $x$、$y$、$z$。它的数学表达形式如下：$$\frac{dx}{dt}=y-a{x}^3+b{x}^2-z+I\text{\hspace{1em}(1)}$$$$\frac{dy}{dt}=c-d{x}^2-y\text{\hspace{1em}(2)}$$$$\frac{dz}{dt}=r[s(x-V_{rest})-z]\text{\hspace{1em}(3)}$$式中，变量 $x$ 表示膜电位，$y$ 和 $z$ 是两个门控变量，$y$ 为快变量，$z$ 为慢变量，$z$ 的时间常数被 $r$ 调节。参数 $a$、$b$、$c$、$d$、$r$、$s$ 均为正数。直观上看，Hindmarsh-Rose 模型方程中各项的意义不甚明确，但我们大致能看出 $y$ 对 $x$ 起促进作用，$x$ 对 $y$ 起抑制作用，从而形成了一个负反馈调节；此外，$x$、$z$ 之间也存在类似的反馈行为。在下文我们将详细讨论模型各个参数对神经元发放模式的影响。

模型定义好后，我们首先利用 BrainPy 构建一个 Hindmarsh-Rose 模型：

import brainpy as bp
import brainpy.math as bmclass HindmarshRose(bp.dyn.NeuDyn):def __init__(self, size, a=1., b=3., c=1., d=5., r=0.001, s=4., x_r=-1.6,theta=1.0, **kwargs):# 初始化父类super(HindmarshRose, self).__init__(size=size, **kwargs)# 初始化参数self.a = aself.b = bself.c = cself.d = dself.r = rself.s = sself.theta = thetaself.x_r = x_r# 初始化变量self.x = bm.Variable(bm.random.randn(self.num) + x_r)self.y = bm.Variable(bm.ones(self.num) * -10.)self.z = bm.Variable(bm.ones(self.num) * 1.7)self.input = bm.Variable(bm.zeros(self.num))self.spike = bm.Variable(bm.zeros(self.num, dtype=bool))  # 脉冲发放状态# 定义积分器self.integral = bp.odeint(f=self.derivative, method='exp_auto')def dx(self, x, t, y, z, Iext):return y - self.a * x * x * x + self.b * x * x - z + Iextdef dy(self, y, t, x):return self.c - self.d * x * x - ydef dz(self, z, t, x):return self.r * (self.s * (x - self.x_r) - z)# 将两个微分方程联合为一个，以便同时积分@propertydef derivative(self):return bp.JointEq([self.dx, self.dy, self.dz])def update(self,):t = bp.share['t']dt = bp.share['dt']x, y, z = self.integral(self.x, self.y, self.z, t, self.input, dt)  # 更新变量x, y, zself.spike.value = bm.logical_and(x >= self.theta, self.x < self.theta)  # 判断神经元是否发放脉冲self.x.value = xself.y.value = yself.z.value = zself.input[:] = 0.  # 重置外界输入

在构建好上述模型之后，我们首先来简单的可视化一下上述模型：

def run_HindmarshRose():group = HindmarshRose(10)runner = bp.DSRunner(group, monitors=['x', 'y', 'z'], inputs=('input', 2.), dt=0.01)runner(20)runner(1000)  # 再运行100msplt.figure(figsize=(6, 4))bp.visualize.line_plot(runner.mon.ts, runner.mon.x, legend='x', show=False)bp.visualize.line_plot(runner.mon.ts, runner.mon.y, legend='y', show=False)bp.visualize.line_plot(runner.mon.ts, runner.mon.z, legend='z', show=True)

可视化结果如下：

图3.25不仅展示了代码的运行结果（图3.25（a）），还将神经元动作电位集中发放阶段（20-200 ms）各个变量的变化情况放大显示了出来（图3.25（b））。

在簇发放模式下，神经元在短时间内快速多次发放动作电位，然后在较长时间内停止发放，并如此周期性重复。从图3.25中可以发现，在一个发放簇中，$x$ 和 $y$ 都随时间快速大幅度变化；相比之下，$z$ 的时间常数 $1/r$ 很大，变化也要平缓很多。正是由于快变量 $y$ 和慢变量 $z$ 的存在，神经元的膜电位才能在某一段时间内快速变化，在另一段时间内缓慢变化。我们将在下节通过动力学分析来探究 Hindmarsh-Rose 模型产生簇发放的原因。

2. Hindmarsh-Rose 模型簇发放的动力学机制

在探究 Hindmarsh-Rose 模型产生簇发放的原因之前，我们先讨论另一个问题：为什么 Hindmarsh-Rose 模型不需要显示地重置膜电位，膜电位自己就能快速上升和下降？

相平面分析无疑是一个很好的用于分析神经元发放模式的工具。当我们想观察两个变量在向量场中的运动时，我们可以把它们可视化到二维相平面图中。但 Hindmarsh-Rose 模型存在三个变量，可视化为三维相平面图并不够直观。有没有办法利用二维相平面来分析 Hindmarsh-Rose 模型呢？

在图3.25中我们已经发现，相对于 $x$ 和 $y$，$z$ 的变化非常平缓。即使在一个发放簇中，$z$ 也没有明显的起伏。这让我们想到，在集中发放阶段，可以将 $z$ 近似为一个恒定值，这样三变量模型就退化为双变量模型，我们就能进行二维相平面分析了。这种分析方法被称为快慢变量分离法，即根据变量的时间常数将它们分为快变量和慢变量，然后将慢变量作为分岔参数，进而研究慢变量的不同取值如何影响快系统的分岔行为。

固定 $z=1.8$，我们可以画出 $x$ 和 $y$ 的相平面分析图：

  # 定义分析器model = HindmarshRose(1)phase_plane_analyzer = bp.analysis.PhasePlane2D(model=model,target_vars={'x': [-2, 3], 'y': [-13., 2.]},  # 待分析变量fixed_vars={'z': 1.8},  # 固定变量pars_update={'Iext': 2.},  # 需要更新的变量resolutions=0.01)# 画出向量场phase_plane_analyzer.plot_vector_field(plot_style=dict(color='lightgrey'))# 画出V, y的零增长曲线phase_plane_analyzer.plot_nullcline()# 画出固定点phase_plane_analyzer.plot_fixed_point()# 画出V, y的变化轨迹phase_plane_analyzer.plot_trajectory({'x': [1.], 'y': [0.]}, duration=100.,color='darkslateblue', linewidth=2, alpha=0.9)

相平面分析图如图3.26（b）所示，为后续分析，我们还额外添加了一些数字标记。图中，变量 $x$ 和 $y$ 构成的轨迹在向量场中形成了一个封闭的圆环，轨迹点沿顺时针方向在环上运动。这种环被称为“极限环”（limit cycle），它表明变量在发生周期性的变化。对应到 Hindmarsh-Rose 模型，$x$、$y$ 在一个发放簇中反复升高和降低。和 AdEx 模型的相平面图相比，Hindmarsh-Rose 模型中膜电位的下降是极限环的一部分，而并非人为的重置。

再来分析在一个极限环中 $x$、$y$ 随时间的变化为何能够产生尖峰。在 $x$、$y$ 随时间的变化图（图3.26（a））中，我们将一次发放分为三个阶段，起始点分别用“1”、“2”、“3”标记，并用虚线标示了 $x$、$y$ 在何处对齐。同时，我们也可以根据 $x$、$y$ 的值将这些标记点对应到相平面图上（图3.26（b））。在 1 → 2 阶段，$x$、$y$ 的轨迹十分靠近 $x$ 和 $y$ 的零增长等值线，这使得轨迹点运动缓慢，因此 $x$ 和 $y$ 的变化也较为平缓。当轨迹到达 2 时，$y$ 达到最大值，形成尖峰。同时由于方程 (1) 中 $y$ 对 $x$ 有促进作用，这在相平面图中表现为向量场的水平分量变大，$x$ 也开始快速上升。此后向量场不断增强（图中表现为灰色箭头变粗），轨迹点移动速度变快，$x$ 急速上升并达到最大值；与此同时，$y$ 迅速下降，几乎在同时达到最小值，即图中的标记点 3。$x$ 的增加使得方程 (1) 中 $-a{x}^3+b{x}^2$ 的三次方项越发显著，$dx/dt$ 由正转负，相平面图中膜电位快速下降至 1，同时伴随着 $y$ 有一个小幅度的急速上升。由此，神经元完成了一个完整的发放。

值得注意的是，轨迹在 1 处有一个较大的转折，这不仅表现为方向的突变，也表现为速度大小的突变：轨迹点在越过近 $x$-零增长等值线前速度很大，越过之后速度迅速变小。这样一个由急到缓的突变构成了神经元在发放动作电位后的快速复极化和缓慢去极化的转折，在不人为重置的情况下非常好地模拟了真实神经元的膜电位变化。

在了解了 Hindmarsh-Rose 模型产生快速连续发放的原因之后，我们再来研究其在两个发放簇之间产生间隔的原因。根据图3.25（a），在变量 $z$ 增加到一定值后，神经元就会停止快速发放，进入间隔期。$z$ 的增加会通过方程 (1) 影响 $x$ 的微分，而不会影响 $y$ 的微分。具体而言，$x$-零增长等值线的表达式为：$$-a{x}^3+b{x}^2+y-z+I=0.\text{\hspace{1em}(4)}$$当 $x$ 保持不变时，$z$ 增大会使 $y$ 也增大时，对应 $x$-$y$ 相平面图中 $x$-零增长等值线上移。我们可以画出 $z=2.05$ 的相平面图，并对比 $z=1.80$ 的相平面图，观察 $z$ 的上升会对系统的动力学性质带来怎样的影响。

图3.27显示，$x$-零增长等值线上移后，两条等值线由一个交点变为了三个，即系统出现了三个奇点。其中，除去原有的不稳定焦点 $U$ 之外，系统多了一个鞍点 $S$ 和稳定结点 $E$。当 $z=1.8$ 时，轨迹点在越过 $x$-零增长等值线时并没有越过 $y$-零增长等值线，因此轨迹点速度的竖直分量仍然向上，轨迹能够围绕 $U$ 形成极限环。但当鞍点 $S$ 出现后，如果轨迹从 $S$ 下方经过，则其越过了过 $y$-零增长等值线，速度的竖直分量变为向下，于是轨迹离开鞍点 $S$ 并被稳定结点 $E$ 吸引。这时神经元膜电位不再上升，也不会发放动作电位，于是形成了两个发放簇之间的间隔。如果 $z$ 保持不变，则 $(x,y)$ 将在点 $E$ 趋于稳定；但由于间隔期 $z$ 在不断下降，当 $z$ 下降到足够低时，鞍点和稳定结点消失，轨迹点速度的竖直分量又变为朝上，极限环再次形成。只不过此时轨迹点的运动速度很小，所以需要一个较长的时间才能进入极限环并产生动作电位。所以，两次簇发放之间会有较长时间的间隔期。

3. Hindmarsh-Rose 模型的其他发放模式

上节我们提到，Hindmarsh-Rose 模型最初被用来解释簇发放现象，但它也能模拟出其他的发放模式，如适应性发放，反弹簇发放等等。除此以外，Hindmarsh-Rose 还能产生一种独特的发放模式：触发簇发放（triggered repetitive firing）。该发放模式是指神经元在接受到一个脉冲电流之后，即使没有外部输入也能够重复多次发放的模式。它的 $x,y,z-t$ 图如图3.28所示。
在理解了 Hindmarsh-Rose 簇发放的动力学原因后，触发簇发放也就不难理解了。在脉冲电流前后，外部电流输入都为 0，且作为慢变量的 $z$ 在接收脉冲电流的这一段很短时间内的变化可以忽略不计，因此 $x-y$ 相平面图并不会发生改变。脉冲电流改变的是变量 $x$ 和 $y$ 的值，我们可以理解为轨迹的起始点 $(x_0,y_0)$ 在相平面图中发生了平移（图3.29（a）和（b））。由于不稳定焦点、鞍点和稳定结点同时存在，轨迹既可能落入极限环，也可能被稳定结点吸引。在接收脉冲电流之前，$(x_0,y_0)$ 位于鞍点下方，落在稳定结点的吸引域，因此会朝稳定结点运动，最终稳定在稳定结点处。虽然图3.29（a）中轨迹的起始点的取值并不等于图3.28中 $t<40ms$ 前 $x$、$y$ 的值，但轨迹仍朝向稳定结点运动，这说明较小的脉冲电流不足以触发动作电位的产生。而当脉冲电流足够大时，$(x_0,y_0)$ 移动到了鞍点上方（图3.29（b）中轨迹的起始点），轨迹能像簇发放一样形成一个极限环，使得神经元重复发放动作电位。
当然，这里的相平面分析没有考虑 $z$ 的变化。随着 $z$ 的增大，$x$-零增长等值线上移，鞍点也将缓慢上移。某一次极限环的过程中，鞍点上升到了轨迹转折点的上方，则轨迹将在越过 $y$-零增长等值线后向下朝着稳定结点运动，神经元不再发放动作电位，簇发放停止。

4. 分析过程所用到的一系列 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 HindmarshRose(bp.dyn.NeuDyn):def __init__(self, size, a=1., b=3., c=1., d=5., r=0.001, s=4., x_r=-1.6,theta=1.0, **kwargs):# 初始化父类super(HindmarshRose, self).__init__(size=size, **kwargs)# 初始化参数self.a = aself.b = bself.c = cself.d = dself.r = rself.s = sself.theta = thetaself.x_r = x_r# 初始化变量self.x = bm.Variable(bm.random.randn(self.num) + x_r)self.y = bm.Variable(bm.ones(self.num) * -10.)self.z = bm.Variable(bm.ones(self.num) * 1.7)self.input = bm.Variable(bm.zeros(self.num))self.spike = bm.Variable(bm.zeros(self.num, dtype=bool))  # 脉冲发放状态# 定义积分器self.integral = bp.odeint(f=self.derivative, method='exp_auto')def dx(self, x, t, y, z, Iext):return y - self.a * x * x * x + self.b * x * x - z + Iextdef dy(self, y, t, x):return self.c - self.d * x * x - ydef dz(self, z, t, x):return self.r * (self.s * (x - self.x_r) - z)# 将两个微分方程联合为一个，以便同时积分@propertydef derivative(self):return bp.JointEq([self.dx, self.dy, self.dz])def update(self,):t = bp.share['t']dt = bp.share['dt']x, y, z = self.integral(self.x, self.y, self.z, t, self.input, dt)  # 更新变量x, y, zself.spike.value = bm.logical_and(x >= self.theta, self.x < self.theta)  # 判断神经元是否发放脉冲self.x.value = xself.y.value = yself.z.value = zself.input[:] = 0.  # 重置外界输入def run_HindmarshRose():group = HindmarshRose(10)runner = bp.DSRunner(group, monitors=['x', 'y', 'z'], inputs=('input', 2.), dt=0.01)runner(20)runner(1000)  # 再运行100msplt.figure(figsize=(6, 4))bp.visualize.line_plot(runner.mon.ts, runner.mon.x, legend='x', show=False)bp.visualize.line_plot(runner.mon.ts, runner.mon.y, legend='y', show=False)bp.visualize.line_plot(runner.mon.ts, runner.mon.z, legend='z', show=True)def bursting_firing():duration = 1200group = HindmarshRose(1)runner = bp.DSRunner(group, monitors=['x', 'y', 'z'], inputs=('input', 2.), dt=0.01)runner(duration)  # 再运行100msdef visualize(mon, duration, xim, text_pos):fig, gs = bp.visualize.get_figure(3, 1, 1.5, 6)ax = fig.add_subplot(gs[0, 0])plt.plot(mon.ts, mon.x, label='x')ax.spines['top'].set_visible(False)ax.spines['right'].set_visible(False)ax.set_xticks([])plt.text(text_pos, (mon.x.max() + mon.x.min()) / 2, r'$x$')ax.set_xlim(*xim)ax = fig.add_subplot(gs[1, 0])plt.plot(mon.ts, mon.y, label='y')ax.spines['top'].set_visible(False)ax.spines['right'].set_visible(False)ax.set_xticks([])plt.text(text_pos, (mon.y.max() + mon.y.min()) / 2, r'$y$')ax.set_xlim(*xim)ax = fig.add_subplot(gs[2, 0])plt.plot(mon.ts, mon.z, label='z')ax.spines['top'].set_visible(False)ax.spines['right'].set_visible(False)plt.text(text_pos, (mon.z.max() + mon.z.min()) / 2, r'$z$')ax.set_xlim(*xim)ax.set_xlabel(r'$t$ (ms)')visualize(runner.mon, duration=duration, xim=(-1, duration + 50), text_pos=duration)plt.savefig('HindmarshRoseModel_output1.pdf', transparent=True, dpi=500)from jax.tree_util import tree_mapvisualize(tree_map(lambda a: a[2000:20000], runner.mon), duration=200,xim=(19, 200 + 8), text_pos=25)plt.savefig('HindmarshRoseModel_output2.pdf', transparent=True, dpi=500)plt.show()def phase_plane_analysis():bp.math.enable_x64()group = HindmarshRose(1)group.x[:] = 1.runner = bp.DSRunner(group, monitors=['x', 'y', 'z'], inputs=('input', 2.), dt=0.01)runner(80)runner(80)  # 再运行100msfig, gs = bp.visualize.get_figure(2, 1, 2.25, 6)ax = fig.add_subplot(gs[0, 0])ax.set_ylabel(r'$x$')ax.set_xlim(80, 160)plt.plot(runner.mon.ts, runner.mon.x, color='tab:blue')plt.axvline(117.8, linestyle='--', color='gray')plt.axvline(135.05, linestyle='--', color='gray')plt.axvline(136.07, linestyle='--', color='gray')plt.text(118.1, -0.85, '1')plt.text(136.57, 1.608, '3')plt.text(137.34, -0.85, "1'")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.y, color='tab:orange')ax.set_ylabel(r'$y$')ax.set_xlim(80, 160)plt.axvline(117.8, linestyle='--', color='gray')plt.axvline(135.05, linestyle='--', color='gray')plt.axvline(136.07, linestyle='--', color='gray')plt.text(118.1, -4.32, '1')plt.text(132.1, 0.22, '2')plt.text(136.72, -6.21, '3')plt.text(137.34, -4.32, "1'")ax.spines['top'].set_visible(False)ax.spines['right'].set_visible(False)plt.savefig('HindmarshRose_output123.pdf', transparent=True, dpi=500)fig, gs = bp.visualize.get_figure(1, 1, 4.5, 6)ax = fig.add_subplot(gs[0, 0])# 定义分析器model = HindmarshRose(1)phase_plane_analyzer = bp.analysis.PhasePlane2D(model=model,target_vars={'x': [-2, 3], 'y': [-13., 2.]},  # 待分析变量fixed_vars={'z': 1.8},  # 固定变量pars_update={'Iext': 2.},  # 需要更新的变量resolutions=0.01)# 画出向量场phase_plane_analyzer.plot_vector_field(plot_style=dict(color='lightgrey'))# 画出V, y的零增长曲线phase_plane_analyzer.plot_nullcline()# 画出固定点phase_plane_analyzer.plot_fixed_point()# 画出V, y的变化轨迹phase_plane_analyzer.plot_trajectory({'x': [1.], 'y': [0.]}, duration=100.,color='darkslateblue', linewidth=2, alpha=0.9)ax.get_legend().remove()plt.text(-1.185, -3.96, '1')plt.text(0.135, 0.92, '2')plt.text(1.836, -6.48, '3')plt.text(1.846, -0.78, 'x nullcline')plt.text(1.676, -10.9, 'y nullcline')plt.text(-0.805, -6.07, 'Trajectory')plt.annotate('unstable focus', xy=(0.6703567413201327, -1.246890802266034),xytext=(-0.415, -3), arrowprops=dict(arrowstyle="->"))ax.spines['top'].set_visible(False)ax.spines['right'].set_visible(False)plt.savefig('HindmarshRose_ppa123.pdf', transparent=True, dpi=500)plt.show()def phase_plane_analysis_v2():bp.math.enable_x64()def analysis(x_range, y_range, z, extra_f=None, name=None):fig, gs = bp.visualize.get_figure(1, 1, 4.5, 6)ax = fig.add_subplot(gs[0, 0])# 定义分析器model = HindmarshRose(1)phase_plane_analyzer = bp.analysis.PhasePlane2D(model=model,target_vars={'x': x_range, 'y': y_range},  # 待分析变量fixed_vars={'z': z},  # 固定变量pars_update={'Iext': 2.},  # 需要更新的变量resolutions=0.01)# 画出向量场phase_plane_analyzer.plot_vector_field(plot_style=dict(color='lightgrey'))# 画出V, y的零增长曲线phase_plane_analyzer.plot_nullcline()# 画出固定点phase_plane_analyzer.plot_fixed_point()# 画出V, y的变化轨迹phase_plane_analyzer.plot_trajectory({'x': [1.], 'y': [0.]}, duration=100.,color='darkslateblue', linewidth=2, alpha=0.9)ax.get_legend().remove()if extra_f: extra_f()ax.spines['top'].set_visible(False)ax.spines['right'].set_visible(False)if name:plt.savefig(name, transparent=True, dpi=500)def f():plt.text(1.846, -0.78, 'x nullcline')plt.text(1.676, -10.9, 'y nullcline')plt.text(-0.805, -6.07, 'Trajectory')plt.annotate('unstable focus', xy=(0.6703567413201327, -1.246890802266034),xytext=(-0.415, -3), arrowprops=dict(arrowstyle="->"))analysis(x_range=[-2, 3], y_range=[-15., 2.], z=1.8, extra_f=f)plt.show()def f():plt.text(-1.0334, -3.441, 'y nullcline')plt.text(-1.05, -4.863, 'x nullcline')plt.text(-0.9539, -4.355, 'Trajectory')analysis(x_range=[-1.2, -0.8], y_range=[-5., -3.], z=1.8, extra_f=f,name='HindmarshRose_I=1.8-v2.pdf')plt.show()def f():plt.text(1.846, -0.78, 'x nullcline')plt.text(1.676, -10.9, 'y nullcline')plt.text(-0.805, -6.07, 'Trajectory')plt.annotate('unstable focus', xy=(0.6040054794020179, -0.8241129286685903),xytext=(-0.415, -3), arrowprops=dict(arrowstyle="->"))plt.annotate('saddle node', xy=(-0.9521776294506809, -3.5332112027188076),xytext=(-1.856, -2.07), arrowprops=dict(arrowstyle="->"))plt.annotate('stable node', xy=(-1.6518281033863245, -12.642680671164262),xytext=(-1.125, -11.58), arrowprops=dict(arrowstyle="->"))analysis(x_range=[-2, 3], y_range=[-15., 2.], z=2.05, extra_f=f,name='HindmarshRose_I=2.05-v1.pdf')plt.show()def f():plt.text(-0.9831, -3.209, 'y nullcline')plt.text(-1.178, -4.863, 'x nullcline')plt.text(-0.9530, -4.095, 'Trajectory')plt.annotate('saddle node', xy=(-0.9521957298446758, -3.5333742767571485),xytext=(-0.9179, -3.755), arrowprops=dict(arrowstyle="->"))analysis(x_range=[-1.2, -0.8], y_range=[-5., -3.], z=2.05, extra_f=f,name='HindmarshRose_I=2.05-v2.pdf')plt.show()def bifurcation_analysis():bp.math.enable_x64()model = HindmarshRose(1)bif = bp.analysis.Bifurcation2D(model=model,target_vars={'x': [-2, 2], 'y': [-20, 5]},fixed_vars={'z': 1.8},target_pars={'Iext': [0., 2.5]},resolutions={'Iext': 0.01})bif.plot_bifurcation(show=True)def bursting_analysis():group = HindmarshRose(1)group.x[:] = -1.6046805group.y[:] = -11.875047group.z[:] = -0.01800228inputs, duration = bp.inputs.section_input([0., 2., 0.], [40, 10, 200], return_length=True)runner = bp.DSRunner(group, monitors=['x', 'y', 'z'], inputs=('input', inputs, 'iter'), dt=0.01)runner.run(duration)fig, gs = bp.visualize.get_figure(3, 1, 1.5, 6)ax = fig.add_subplot(gs[0, 0])plt.plot(runner.mon.ts, runner.mon.x, color='tab:blue')ax.set_ylabel(r'$x$')ax = fig.add_subplot(gs[1, 0])plt.plot(runner.mon.ts, runner.mon.y, color='tab:orange')ax.set_ylabel(r'$y$')ax = fig.add_subplot(gs[2, 0])plt.plot(runner.mon.ts, runner.mon.z, color='tab:brown')ax.set_ylabel(r'$z$')plt.show()def triggered_bursting():I, duration = bp.inputs.section_input([0, 2., 0], [40, 5, 400], return_length=True, dt=0.01)model = HindmarshRose(1, s=0.5)model.x[:] = model.x_rmodel.y[:] = -12.model.z[:] = 0.runner = bp.DSRunner(model, inputs=(model.input, I, 'iter'), monitors=list('xyz'), dt=0.01)runner.run(duration)fig, gs = bp.visualize.get_figure(3, 1, 1.5, 8)ax = fig.add_subplot(gs[0:2, 0])plt.plot(runner.mon.ts, runner.mon.x, label='x')plt.plot(runner.mon.ts, runner.mon.y, linestyle='dashed', label='y')plt.plot(runner.mon.ts, runner.mon.z, linestyle='dotted', label='z')ax.spines['top'].set_visible(False)ax.spines['right'].set_visible(False)plt.legend()plt.xlim(-1, duration + 1)ax = fig.add_subplot(gs[2, 0])plt.plot(runner.mon.ts, bm.as_numpy(I))plt.xlim(-1, duration + 1)plt.ylabel(r'$I$')plt.xlabel(r'$t$ (ms)')ax.spines['top'].set_visible(False)ax.spines['right'].set_visible(False)plt.savefig('HindmarshRose_triggered_bursting.pdf', transparent=True, dpi=500)plt.show()def triggered_bursting_ppa():bp.math.enable_x64()def analysis(x_range, y_range, z, init_xy=None, extra_f=None, name=None):fig, gs = bp.visualize.get_figure(1, 1, 4.5, 6)ax = fig.add_subplot(gs[0, 0])# 定义分析器model = HindmarshRose(1, s=0.5)phase_plane_analyzer = bp.analysis.PhasePlane2D(model=model,target_vars={'x': x_range, 'y': y_range},  # 待分析变量fixed_vars={'z': z},  # 固定变量pars_update={'Iext': 0.},  # 需要更新的变量resolutions=0.01)# 画出向量场phase_plane_analyzer.plot_vector_field(plot_style=dict(color='lightgrey'))# 画出V, y的零增长曲线phase_plane_analyzer.plot_nullcline()# 画出固定点phase_plane_analyzer.plot_fixed_point()# 画出V, y的变化轨迹if init_xy is not None:phase_plane_analyzer.plot_trajectory({'x': [init_xy[0]], 'y': [init_xy[1]]},duration=100., color='darkslateblue',linewidth=2, alpha=0.9, dt=0.01)ax.get_legend().remove()if extra_f:extra_f()ax.spines['top'].set_visible(False)ax.spines['right'].set_visible(False)if name:plt.savefig(name, transparent=True, dpi=500)def f():plt.text(1.5, -9.8, 'y nullcline')plt.text(1.13, -2, 'x nullcline')plt.text(-0.575, -7.76, 'Trajectory')plt.annotate('stable node', xy=(-1.6320958395541862, -12.318684344410208),xytext=(-1, -12.318684344410208), arrowprops=dict(arrowstyle="->"))plt.annotate('saddle node', xy=(-0.9803837299517475, -3.8057584936177298),xytext=(-0.485, -2.5), arrowprops=dict(arrowstyle="->"))plt.annotate('unstable focus', xy=(0.6124732604697042, -0.8756173672185038),xytext=(0.8, 0.6), arrowprops=dict(arrowstyle="->"))analysis([-2, 3], [-13., 2.], z=0.02, extra_f=f, init_xy=(1.4, -10.),name='HindmarshRose_triggered_bursting_ppa1.pdf')def f():plt.text(1.5, -9.8, 'y nullcline')plt.text(1.73, -2, 'x nullcline')plt.text(-0.145, -3.5, 'Trajectory')plt.annotate('stable node', xy=(-1.6320958395541862, -12.318684344410208),xytext=(-1, -12.318684344410208), arrowprops=dict(arrowstyle="->"))plt.annotate('saddle node', xy=(-0.9803837299517475, -3.8057584936177298),xytext=(-0.685, -6.3), arrowprops=dict(arrowstyle="->"))plt.annotate('unstable focus', xy=(0.6124732604697042, -0.8756173672185038),xytext=(0.8, 0.6), arrowprops=dict(arrowstyle="->"))analysis([-2, 3], [-13., 2.], z=0.02, extra_f=f, init_xy=(-0.7, -3),name='HindmarshRose_triggered_bursting_ppa2.pdf')plt.show()if __name__ == '__main__':run_HindmarshRose()bursting_firing()phase_plane_analysis()phase_plane_analysis_v2()bifurcation_analysis()bursting_analysis()triggered_bursting()triggered_bursting_ppa()