.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/fhn.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_fhn.py: FitzHugh-Nagumo with the Analytical Jacobian ============================================ The `FitzHugh-Nagumo (FHN) model `_ is a textbook 2D reduction of the Hodgkin-Huxley equations. The system has a single cubic non-linearity and is a clean target for demonstrating the *analytical-Jacobian* path in `PyCoBi`: PyRates symbolically differentiates the vector field and writes the DFDU / DFDP entries straight into the generated Fortran, and auto-07p uses the analytical Jacobian instead of falling back to finite differences. The model equations read: .. math:: \dot v &= v - \frac{v^3}{3} - w + I_\mathrm{ext}, \dot w &= \varepsilon (v + a - b w), where :math:`v` is the fast membrane-potential variable and :math:`w` is a slow recovery variable. With the standard parameter regime (:math:`a=0.7,\, b=0.8,\, \varepsilon=0.08`) the system has a stable fixed point for :math:`I_\mathrm{ext} \lesssim 0.34` and :math:`I_\mathrm{ext} \gtrsim 1.42` and a stable limit cycle in between, bounded by two Hopf bifurcations. In what follows we will: 1. Load the FHN model from a co-located YAML file and confirm that PyRates emitted analytical DFDU / DFDP entries into the generated Fortran. 2. Continue the steady state in :math:`I_\mathrm{ext}` and locate the two Hopfs. 3. Branch-switch at HB1 to trace the periodic-solution branch (the limit cycle) and overlay it on the equilibrium diagram. 4. Compare the analytical and finite-difference Jacobian paths on the limit-cycle continuation — where the Newton step solves a larger BVP system than on the equilibrium branch, so the analytical path has more room to pay off than it does on a plain equilibrium scan. References ^^^^^^^^^^ .. [1] R. FitzHugh (1961) *Impulses and physiological states in theoretical models of nerve membrane.* Biophysical Journal 1 (6): 445-466, https://doi.org/10.1016/S0006-3495(61)86902-6. .. GENERATED FROM PYTHON SOURCE LINES 42-48 Step 1: load the model and inspect the generated Fortran ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ We ship a small ``fhn.yaml`` next to this script so the example is self-contained. PyRates' slash-notation path takes the form ``//``; we anchor on ``__file__`` so the load works regardless of cwd. .. GENERATED FROM PYTHON SOURCE LINES 48-72 .. code-block:: Python import re import time from pathlib import Path from pycobi import ODESystem import matplotlib.pyplot as plt here = Path(__file__).resolve().parent yaml_path = str(here / 'fhn' / 'fhn') # `analytical_jacobian=True` is the default; spelling it out for emphasis. We also request both an `ivp` # (for the IVP that converges us to the steady state) and an `eq` scenario so the second `run()` call can # load `c.eq` directly. `init_cont=True` opts into the legacy automatic-IVP behaviour — the default since # PyCoBi 0.10.0 is `False`. ode = ODESystem.from_yaml( yaml_path, auto_dir="~/PycharmProjects/auto-07p", init_cont=True, analytical_jacobian=True, auto_constants=('ivp', 'eq', 'lc'), NMX=20000, NPR=20000, ) .. GENERATED FROM PYTHON SOURCE LINES 73-76 Let's check that PyRates actually wrote the analytical Jacobian into the generated Fortran by grepping for ``dfdu(`` and ``dfdp(`` assignments in the source file. We'd also expect ``JAC = 1`` in the generated ``c.eq`` / ``c.lc`` files since the JAC flag is wired up to the analytical-Jacobian emission. .. GENERATED FROM PYTHON SOURCE LINES 76-88 .. code-block:: Python src_path = Path(ode.dir) / 'system_equations.f90' src = src_path.read_text() n_dfdu = len(re.findall(r'\bdfdu\(\s*\d', src, re.IGNORECASE)) n_dfdp = len(re.findall(r'\bdfdp\(\s*\d', src, re.IGNORECASE)) print(f"DFDU entries in {src_path.name}: {n_dfdu}") print(f"DFDP entries in {src_path.name}: {n_dfdp}") ceq = (Path(ode.dir) / 'c.eq').read_text() clc = (Path(ode.dir) / 'c.lc').read_text() print("JAC = 1 in c.eq:", any('JAC = 1' in ln for ln in ceq.splitlines())) print("JAC = 1 in c.lc:", any('JAC = 1' in ln for ln in clc.splitlines())) .. GENERATED FROM PYTHON SOURCE LINES 89-94 You should see four ``DFDU`` entries (the 2x2 Jacobian is dense for FHN — all four entries are non-zero) and four ``DFDP`` entries (one row of partials per state variable, across the parameters the right-hand side actually depends on), all written by PyRates' symbolic differentiation. The cubic ``v^3/3`` term in the RHS becomes ``1 - y(1)**2`` in the DFDU(1,1) slot — that's exactly :math:`\partial \dot v / \partial v` for the FHN equations. .. GENERATED FROM PYTHON SOURCE LINES 96-102 Step 2: continue the steady state in :math:`I_\mathrm{ext}` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Standard 1D equilibrium continuation. We start from ``'EP2'``, which is the converged final point of the IVP that the constructor ran (``'EP1'`` is the :math:`t=0` initial condition). `bidirectional=True` walks the branch both ways so we capture the full S-curve. .. GENERATED FROM PYTHON SOURCE LINES 102-116 .. code-block:: Python eq_sols, eq_cont = ode.run( starting_point='EP2', name='I_branch', c='eq', ICP='Iext', bidirectional=True, NMX=500, NPR=10, DS=0.01, DSMAX=0.05, ) print("\nEquilibrium-branch bifurcation counts:") print(eq_sols['bifurcation'].value_counts()) print(eq_sols.loc[eq_sols['bifurcation'] == 'HB', ['bifurcation', 'Iext']]) ode.plot_continuation('Iext', 'v', cont='I_branch') plt.title('FHN: equilibria as a function of $I_\\mathrm{ext}$') plt.show() .. GENERATED FROM PYTHON SOURCE LINES 117-121 The bifurcation diagram shows the steady state's :math:`v` value as a function of :math:`I_\mathrm{ext}`. Solid line: stable; dotted gray: unstable. The green circles mark the two Hopf bifurcations (:math:`\mathrm{HB}_1 \approx 0.331` and :math:`\mathrm{HB}_2 \approx 1.419`) where the equilibrium loses / regains stability via a complex-conjugate pair of eigenvalues crossing the imaginary axis. .. GENERATED FROM PYTHON SOURCE LINES 123-136 Step 3: branch-switch at HB1 and continue the limit cycle ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ At a Hopf bifurcation a branch of periodic solutions emerges from the equilibrium. We switch onto that branch by passing the Hopf label as ``starting_point`` and setting auto-07p's switches: ``IPS=2`` (continue periodic solutions), ``ISW=-1`` (start from HB / PD / BP), and ``c='lc'`` which loads the limit-cycle preset PyRates wrote for us. ``ICP=['Iext', 11]`` declares Iext as the continuation parameter plus ``PAR(11)`` for the period (auto's convention for LC continuations). For FHN with these parameters the LC born at HB1 is *sub-critical*: a small-amplitude unstable cycle emerges, undergoes a fold-of-cycle (LP) at :math:`I_\mathrm{ext} \approx 0.324`, becomes stable, and then the stable branch sweeps across the oscillatory regime up to HB2. We use a larger ``NMX`` than the equilibrium scan to cover that whole sweep. .. GENERATED FROM PYTHON SOURCE LINES 136-146 .. code-block:: Python lc_sols, lc_cont = ode.run( origin='I_branch', starting_point='HB1', name='lc', c='lc', ICP=['Iext', 11], IPS=2, ISP=2, ISW=-1, NMX=5000, NPR=500, DS=0.01, DSMAX=0.5, get_period=True, ) print("\nLimit-cycle branch bifurcation counts:") print(lc_sols['bifurcation'].value_counts()) .. GENERATED FROM PYTHON SOURCE LINES 147-150 Overlay the LC branch on the equilibrium diagram. For a limit cycle, ``plot_continuation`` packs each period into a (min, max) envelope, so the LC branch appears as a pair of curves bounding the oscillation amplitude in :math:`v`. .. GENERATED FROM PYTHON SOURCE LINES 150-157 .. code-block:: Python fig, ax = plt.subplots() ode.plot_continuation('Iext', 'v', cont='I_branch', ax=ax) ode.plot_continuation('Iext', 'v', cont='lc', ax=ax, ignore=['UZ', 'BP']) ax.set_title('FHN: equilibrium + limit cycle in $I_\\mathrm{ext}$') plt.show() .. GENERATED FROM PYTHON SOURCE LINES 158-167 Step 4: analytical vs. finite-difference Jacobian ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The ``JAC=`` override is per-call, so we can re-run the same LC continuation with the FD fallback without re-instantiating the system. We time both on the limit-cycle branch because that's where the Newton step solves a much bigger BVP system than on the equilibrium scan (auto's BVP system has dimension ``NTST × NCOL × NDIM`` per Newton iteration — for the c.lc defaults of ``NTST=50``, ``NCOL=4``, and ``NDIM=2`` that's a 400-dim Jacobian per iteration, versus a 2x2 on the equilibrium branch). .. GENERATED FROM PYTHON SOURCE LINES 167-190 .. code-block:: Python t0 = time.perf_counter() ode.run( origin='I_branch', starting_point='HB1', name='lc_jac', c='lc', ICP=['Iext', 11], IPS=2, ISP=2, ISW=-1, NMX=5000, NPR=5000, DS=0.01, DSMAX=0.5, JAC=1, # analytical (the default since we passed analytical_jacobian=True) ) t_jac = time.perf_counter() - t0 t0 = time.perf_counter() ode.run( origin='I_branch', starting_point='HB1', name='lc_fd', c='lc', ICP=['Iext', 11], IPS=2, ISP=2, ISW=-1, NMX=5000, NPR=5000, DS=0.01, DSMAX=0.5, JAC=0, # force finite-difference fallback ) t_fd = time.perf_counter() - t0 print(f"\nLC continuation, analytical Jacobian: {t_jac:5.2f}s") print(f"LC continuation, finite-difference : {t_fd:5.2f}s") print(f"speed-up : {t_fd / t_jac:5.2f}x") .. GENERATED FROM PYTHON SOURCE LINES 191-199 For a 2D system the speed-up is modest (the symbolic 2x2 Jacobian only saves us a handful of RHS evaluations per Newton iteration), but the analytical path scales much better with state-vector dimension. The finite-difference Jacobian needs ``NDIM`` RHS evaluations to fill each Newton column — :math:`O(\mathrm{NDIM}^2)` total — while the analytical path fills the matrix in a single symbolic sweep regardless of NDIM. For the medium-sized biophysical models PyCoBi was built for (NDIM in the 10–50 range), the same comparison typically cuts continuation time by an order of magnitude and improves Newton convergence near fold-of-cycle / period-doubling points where the finite-difference Jacobian becomes noisy. .. GENERATED FROM PYTHON SOURCE LINES 201-203 Step 5: clean up ^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 203-206 .. code-block:: Python ode.close_session(clear_files=True) .. GENERATED FROM PYTHON SOURCE LINES 207-209 As usual, ``close_session(clear_files=True)`` removes the temporary Fortran sources and ``c.*`` files. Use ``clear_files=False`` to keep them for inspection. .. _sphx_glr_download_auto_examples_fhn.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: fhn.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: fhn.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: fhn.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_