FitzHugh-Nagumo analysis

The FitzHugh-Nagumo model is given by:

\[\frac {dV} {dt} = V(1 - \frac {V^2} 3) - w + I_{ext}\]
\[\tau \frac {dw} {dt} = V + a - b w\]
[1]:
from collections import OrderedDict
import brainpy as bp

bp.backend.set(dt=0.02)
[2]:
class FitzHughNagumo(bp.NeuGroup):
    target_backend = 'general'

    def __init__(self, size, a=0.7, b=0.8, tau=12.5, Vth=1.9, **kwargs):
        self.a = a
        self.b = b
        self.tau = tau
        self.Vth = Vth

        self.V = bp.ops.zeros(size)
        self.w = bp.ops.zeros(size)
        self.spike = bp.ops.zeros(size)
        self.input = bp.ops.zeros(size)

        super(FitzHughNagumo, self).__init__(size=size, **kwargs)

    @staticmethod
    @bp.odeint(method='rk4')
    def integral(V, w, t, Iext, a, b, tau):
        dw = (V + a - b * w) / tau
        dV = V - V * V * V / 3 - w + Iext
        return dV, dw

    def update(self, _t):
        V, self.w = self.integral(self.V, self.w, _t, self.input, self.a, self.b, self.tau)
        self.spike = (V >= self.Vth) * (self.V < self.Vth)
        self.V = V
        self.input[:] = 0.

Simulation

[3]:
FNs = FitzHughNagumo(2, monitors=['V'])

FNs.run(duration=300., inputs=('input', 1.), report=True)
bp.visualize.line_plot(FNs.mon.ts, FNs.mon.V, show=True)

FNs.run(duration=(300., 600.), inputs=('input', 0.6), report=True)
bp.visualize.line_plot(FNs.mon.ts, FNs.mon.V, show=True)
Compilation used 0.0000 s.
Start running ...
Run 10.0% used 0.145 s.
Run 20.0% used 0.309 s.
Run 30.0% used 0.465 s.
Run 40.0% used 0.597 s.
Run 50.0% used 0.737 s.
Run 60.0% used 0.881 s.
Run 70.0% used 1.027 s.
Run 80.0% used 1.172 s.
Run 90.0% used 1.322 s.
Run 100.0% used 1.466 s.
Simulation is done in 1.466 s.

../../_images/tutorials_dynamics_analysis_FitzHugh_Nagumo_analysis_5_1.png
Compilation used 0.0010 s.
Start running ...
Run 10.0% used 0.149 s.
Run 20.0% used 0.303 s.
Run 30.0% used 0.460 s.
Run 40.0% used 0.610 s.
Run 50.0% used 0.758 s.
Run 60.0% used 0.922 s.
Run 70.0% used 1.134 s.
Run 80.0% used 1.310 s.
Run 90.0% used 1.451 s.
Run 100.0% used 1.606 s.
Simulation is done in 1.606 s.

../../_images/tutorials_dynamics_analysis_FitzHugh_Nagumo_analysis_5_3.png

Phase plane analysis

[4]:
phase = bp.analysis.PhasePlane(FNs.integral,
                               target_vars={'V': [-3, 2], 'w': [-2, 2]},
                               fixed_vars=None,
                               pars_update={'Iext': 1., "a": 0.7, 'b': 0.8, 'tau': 12.5})
phase.plot_nullcline()
phase.plot_fixed_point()
# phase.plot_trajectory(initials={'V': -1, 'w': 1}, duration=100.)
phase.plot_limit_cycle_by_sim(initials={'V': -1, 'w': 1}, duration=100.)
phase.plot_vector_field(show=True)
plot nullcline ...
SymPy solve "integral(V, w) = 0" to "w = f(V, )", success.
SymPy solve "integral(V, w) = 0" to "w = f(V, )", success.
plot fixed point ...
SymPy solve derivative of "integral(V, w)" by "V", success.
SymPy solve derivative of "integral(V, w)" by "w", success.
SymPy solve derivative of "integral(V, w)" by "V", success.
SymPy solve derivative of "integral(V, w)" by "w", success.
Fixed point #1 at V=0.4088658369434122, w=1.3860822961792651 is a unstable node.
plot limit cycle ...
plot vector field ...
../../_images/tutorials_dynamics_analysis_FitzHugh_Nagumo_analysis_7_1.png

Codimension 1 bifurcation analysis

[5]:
bifurcation = bp.analysis.Bifurcation(FNs.integral,
                                      target_pars={'Iext': [-1, 1]},
                                      target_vars={'V': [-3, 2], 'w': [-2, 2]},
                                      fixed_vars=None,
                                      pars_update={'a': 0.7, 'b': 0.8, 'tau': 12.5},
                                      numerical_resolution=0.01)
bifurcation.plot_bifurcation(show=True)
plot bifurcation ...
SymPy solve "integral(V, w, Iext) = 0" to "w = f(V, Iext)", success.
SymPy solve derivative of "integral(V, w, Iext)" by "V", success.
SymPy solve derivative of "integral(V, w, Iext)" by "w", success.
SymPy solve derivative of "integral(V, w, Iext)" by "V", success.
SymPy solve derivative of "integral(V, w, Iext)" by "w", success.
../../_images/tutorials_dynamics_analysis_FitzHugh_Nagumo_analysis_9_1.png
../../_images/tutorials_dynamics_analysis_FitzHugh_Nagumo_analysis_9_2.png

Codimension 2 bifurcation analysis

[6]:
bifurcation = bp.analysis.Bifurcation(FNs.integral,
                                      target_pars=OrderedDict(a=[0.5, 1.], Iext=[0., 1.]),
                                      target_vars=OrderedDict(V=[-3, 3], w=[-3., 3.]),
                                      fixed_vars=None,
                                      pars_update={'b': 0.8, 'tau': 12.5},
                                      numerical_resolution=0.01)
bifurcation.plot_bifurcation(show=True)
plot bifurcation ...
SymPy solve "integral(V, w, a, Iext) = 0" to "w = f(V, a,Iext)", success.
SymPy solve derivative of "integral(V, w, a, Iext)" by "V", success.
SymPy solve derivative of "integral(V, w, a, Iext)" by "w", success.
SymPy solve derivative of "integral(V, w, a, Iext)" by "V", success.
SymPy solve derivative of "integral(V, w, a, Iext)" by "w", success.
../../_images/tutorials_dynamics_analysis_FitzHugh_Nagumo_analysis_11_1.png
../../_images/tutorials_dynamics_analysis_FitzHugh_Nagumo_analysis_11_2.png