Automated Codim-2 Search and Period-Doubling Cascades

PyCoBi ships two convenience helpers for the bookkeeping-heavy parts of a bifurcation analysis:

• codim2_search() — given a list of codim-1 starting points (folds, Hopfs, period-doublings) on an existing continuation, runs a 2-parameter continuation of each, walks the resulting curves for codim-2 points, and recursively continues the codim-1 bifurcations that emerge from them. Supports recursive handling of zero-Hopf (ZH), generalised-Hopf (GH / Bautin), and Bogdanov-Takens (BT) points.

• continue_period_doubling_bf() — chases a cascade of period-doubling bifurcations in 2 parameters, recursing on every new PD encountered. Useful for tracing the boundaries of period-doubling routes to chaos.

This example uses the QIF mean-field model with bi-exponential spike-frequency adaptation (qif_biexp_sfa.yaml next to this script). The bi-exponential kernel is a strict generalisation of the alpha-kernel QIF-SFA from Hopf Bifurcation and Limit Cycle Continuation: both kernels satisfy \(\tau_a^2 A'' + 2 \tau_a A' + A = \alpha r \tau_a\) when tau_r == tau_d == tau_a, so the codim-2 structure in the \((\bar\eta,\, \Delta)\) plane (generalised-Hopf, Bogdanov-Takens, cusps) is identical to the alpha-kernel case at that parameter point — verified by side-by-side bifurcation analysis. Picking the bi-exponential model now means we can also explore the period-doubling regime that opens up when tau_r is taken much smaller than tau_d (Section 2.1 below).

References

[1]

R. Gast, H. Schmidt, T.R. Knösche (2020) A Mean-Field Description of Bursting Dynamics in Spiking Neural Networks with Short-Term Adaptation. Neural Computation 32 (9): 1615-1634.

Step 1: Load the model

Bi-exponential QIF-SFA from the co-located YAML, with rise and decay adaptation time constants both set to 10 — the alpha-kernel-equivalent parameter point. Coupling J = 15 sqrt(2), alpha = 1.0, Delta = 2.0 to match the QIF-SFA example.

from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt

from pycobi import ODESystem
from pycobi.automated_continuation import (
    codim2_search,
    continue_period_doubling_bf,
)

here = Path(__file__).resolve().parent
yaml_path = str(here / 'qif_biexp_sfa' / 'qif_biexp_sfa')

ode = ODESystem.from_yaml(
    yaml_path,
    auto_dir="~/PycharmProjects/auto-07p",
    node_vars={
        'p/qif_biexp_sfa_op/Delta': 2.0,
        'p/qif_biexp_sfa_op/alpha': 0.8,
        # Start at tau_r = 11 (slightly above tau_d=10) and continue down in
        # tau_r before any eta work — see Step 2 below.
        'p/qif_biexp_sfa_op/tau_r': 11.0,
        'p/qif_biexp_sfa_op/tau_d': 10.0,
        'p/qif_biexp_sfa_op/eta': -8.0,
    },
    edge_vars=[('p/qif_biexp_sfa_op/r', 'p/qif_biexp_sfa_op/r_in',
                {'weight': 15.0 * np.sqrt(2.0)})],
    init_cont=True, NPR=100, NMX=30000,
)

# Auto-07p doesn't ship a built-in style for the cusp (CP) marker; add one
# so it shows up in the codim-2 diagram below alongside BT and GH.
ode.update_bifurcation_style('CP', marker='d', color='#7F4FBF')

Step 2: Pre-scan in \(\\tau_r\) to anchor the two regimes

We’ll cover two parameter regimes of this model in one example:

1. \(\\tau_r = \\tau_d = 10\) (alpha-kernel-equivalent): clean codim-2 structure in \((\\bar\\eta,\\, \\Delta)\) — BT / GH / CP, no PD.

2. \(\\tau_r = 0.1 \\ll \\tau_d\): same codim-2 backbone plus a period-doubling cascade on the LC born at the upper Hopf.

Rather than re-instantiate the model twice, we continue the steady state in \(\\tau_r\) from the initial \(\\tau_r = 11\) down to small values, planting user points at \(\\tau_r = \\tau_d = 10\) (UZ1) and at \(\\tau_r = 0.1\) (UZ2). Each UZ then serves as the starting point for an \(\\bar\\eta\) continuation in its own regime.

tau_sols, tau_cont = ode.run(
    starting_point='EP2', name='tau_r_branch',
    ICP='p/qif_biexp_sfa_op/tau_r',
    IPS=1, ILP=1, ISP=2, ISW=1, NTST=400, NCOL=4,
    NMX=5000, NPR=100, DS=-1e-3, DSMIN=1e-9, DSMAX=5e-2,
    UZR={'p/qif_biexp_sfa_op/tau_r': [10.0, 0.1]},
    UZSTOP={'p/qif_biexp_sfa_op/tau_r': [0.01, 10.0]},
)
print(f"tau_r pre-scan bifurcations: {dict(tau_sols['bifurcation'].value_counts())}")

Shared bookkeeping

We run the same analysis pipeline at two values of \(\\tau_r\) (Steps 3-4 and 5-6 below), so it pays to factor the codim-1 / codim-2 bookkeeping into a small helper. codim1_points lists the four starting points we expect on every \(\\bar\\eta\) scan (the model has two folds and two Hopfs in the parameter regime we’re exploring); CURVE_COLORS keys colours by bifurcation type so the codim-2 diagrams collapse their legends to two entries even though we draw eight curves.

codim1_points = [
    ('LP1', 'fold'),
    ('LP2', 'fold'),
    ('HB1', 'Hopf'),
    ('HB2', 'Hopf'),
]
CURVE_COLORS = {'fold': '#1F77B4', 'Hopf': '#FF7F0E'}
PD_COLOR = '#D62728'


def run_codim2_in_delta(origin_name, name_prefix):
    """Run the 8-call codim-2 search (LP1/LP2/HB1/HB2 × DS=±1e-3) in
    :math:`(\\bar\\eta,\\, \\Delta)` and return ``[(key, bif_type)]``.
    Forward + reverse continuations are kept separate (rather than one
    bidirectional call) because bidirectional gets stuck bouncing at
    the cusp adjacent to LP1 / LP2.
    """
    shared_kwargs = dict(
        pyauto_instance=ode,
        params=['p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/Delta'],
        origin=origin_name,
        max_recursion_depth=0,
        NMX=1500, NPR=10,
        DSMIN=1e-9, DSMAX=5e-2,
        RL0=-15.0, RL1=5.0,
        bidirectional=False,
        UZSTOP={'p/qif_biexp_sfa_op/Delta': [0.0, 4.0]},
    )
    curves = []
    for sp, bif_type in codim1_points:
        for ds in (1e-3, -1e-3):
            try:
                result = codim2_search(
                    starting_points=[sp], DS=ds,
                    name=f'{name_prefix}_{sp}_{"pos" if ds > 0 else "neg"}',
                    **shared_kwargs,
                )
                curves.append((list(result.keys())[0], bif_type))
            except Exception as exc:
                print(f"{name_prefix} {sp} DS={ds:+g}: skipped "
                      f"({type(exc).__name__}: {exc})")
    return curves


def plot_2d_diagram(ax, codim2_curves, eta_origin, *, pd_names=(),
                     title, xlim=(-6.0, 2.0), ylim=(0.0, 2.5)):
    """Plot a 2D codim-2 diagram in :math:`(\\bar\\eta,\\, \\Delta)` from
    the curves produced by :func:`run_codim2_in_delta`. Optionally overlay
    one or more PD continuations (rendered in red). Fold curves all share
    the blue colour, Hopf curves all share orange — so the legend has at
    most three entries (fold curve, Hopf curve, PD curve)."""
    labels_used = set()
    for key, bif_type in codim2_curves:
        color = CURVE_COLORS[bif_type]
        label = f'{bif_type} curve' if bif_type not in labels_used else None
        labels_used.add(bif_type)
        ode.plot_continuation(
            'p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/Delta',
            cont=key, ax=ax,
            line_color_stable=color, line_color_unstable=color,
            line_style_stable='solid', line_style_unstable='solid',
            bifurcation_legend=False, get_stability=False,
            ignore=['LP', 'HB', 'UZ'], label=label,
        )
    for pd_name in pd_names:
        label = 'PD curve' if 'PD' not in labels_used else None
        labels_used.add('PD')
        ode.plot_continuation(
            'p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/Delta',
            cont=pd_name, ax=ax,
            line_color_stable=PD_COLOR, line_color_unstable=PD_COLOR,
            line_style_stable='solid', line_style_unstable='solid',
            bifurcation_legend=False, get_stability=False,
            ignore=['PD', 'UZ'], label=label,
        )
    for sp, bif_type in codim1_points:
        try:
            sol, _, _ = ode.get_solution(point=sp, cont=eta_origin)
            eta_p = float(sol['eta'])
            delta_p = float(sol['Delta'])
            ax.scatter(eta_p, delta_p, marker='*', s=120,
                        c=CURVE_COLORS[bif_type],
                        edgecolor='k', linewidth=0.5, zorder=10)
            ax.annotate(sp, (eta_p, delta_p), xytext=(5, 5),
                         textcoords='offset points', fontsize=9)
        except Exception:
            pass
    ax.set_xlim(*xlim)
    ax.set_ylim(*ylim)
    ax.set_xlabel(r'$\bar\eta$')
    ax.set_ylabel(r'$\Delta$')
    ax.set_title(title)
    ax.legend(loc='best')


def plot_1d_diagram(ax, eta_cont_name, lc_cont_name, *, title,
                     xlim=(-8.0, 2.0), cascade_names=()):
    """Plot a 1D bifurcation diagram in :math:`(\\bar\\eta,\\, r)` with the
    equilibrium branch (blue), the limit cycle branch (orange), and
    optionally a list of cascade LC continuations (red — each doubled
    cycle from a PD branch switch) overlaid. Codim-1 markers (HB / LP /
    PD) appear automatically; UZ / BP / EP are filtered out as cosmetic
    noise."""
    ode.plot_continuation(
        'p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/r',
        cont=eta_cont_name, ax=ax,
        line_color_stable='#1F77B4', line_color_unstable='#1F77B4',
        bifurcation_legend=False, label='equilibrium branch',
        ignore=['UZ', 'BP', 'EP'],
    )
    ode.plot_continuation(
        'p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/r',
        cont=lc_cont_name, ax=ax,
        line_color_stable='#FF7F0E', line_color_unstable='#FF7F0E',
        bifurcation_legend=False, label='limit cycle branch',
        ignore=['UZ', 'BP', 'EP'],
    )
    cascade_label_done = False
    for cascade_name in cascade_names:
        label = None if cascade_label_done else 'PD cascade LCs'
        cascade_label_done = True
        try:
            ode.plot_continuation(
                'p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/r',
                cont=cascade_name, ax=ax,
                line_color_stable=PD_COLOR, line_color_unstable=PD_COLOR,
                bifurcation_legend=False, label=label,
                ignore=['UZ', 'BP', 'EP'],
            )
        except Exception as exc:
            print(f"  warning: could not plot cascade LC {cascade_name!r}: {exc}")
    ax.set_xlim(*xlim)
    ax.set_xlabel(r'$\bar\eta$')
    ax.set_ylabel(r'$r$')
    ax.set_title(title)
    ax.legend(loc='best')

Step 3: Continuation + codim-2 search at \(\\tau_r = \\tau_d = 10\)

This is the alpha-kernel-equivalent regime — clean codim-2 structure (BT / GH / CP), no period doublings. A bidirectional \(\\bar\\eta\) scan from UZ1 exposes the codim-1 starting points; we then branch-switch to the limit cycle born at HB2 and run codim2_search() from each codim-1 point.

eta_sols, eta_cont = ode.run(
    origin='tau_r_branch', starting_point='UZ1', name='eta_branch',
    ICP='p/qif_biexp_sfa_op/eta', bidirectional=True,
    RL0=-8.0, RL1=2.0,
    IPS=1, ILP=1, ISP=2, ISW=1, NTST=400, NCOL=4,
    NMX=2000, NPR=10, DS=1e-4, DSMIN=1e-8, DSMAX=5e-2,
    ITMX=40, ITNW=40, NWTN=12,
)
print("\neta scan (tau_r=10):", dict(eta_sols['bifurcation'].value_counts()))

lc_sols, lc_cont = ode.run(
    origin='eta_branch', starting_point='HB2', name='lc_branch',
    IPS=2, ISP=2, ISW=-1,
    ICP=['p/qif_biexp_sfa_op/eta', 11],
    NMX=2000, NPR=10, DS=1e-3, DSMIN=1e-9, DSMAX=5e-2,
    bidirectional=True, get_period=True, STOP=["BP3", "LP5"]
)
print("LC (tau_r=10):", dict(lc_sols['bifurcation'].value_counts()))

codim2_curves = run_codim2_in_delta('eta_branch', name_prefix='delta10')
print("codim-2 curves at tau_r=10:")
for key, bif_type in codim2_curves:
    print(f"  {key} ({bif_type}): "
          f"{dict(ode.get_summary(key)['bifurcation'].value_counts())}")

Step 4: Figures at \(\\tau_r = 10\)

Two figures: a 1D bifurcation diagram in \((\\bar\\eta,\\, r)\) with the equilibrium and limit-cycle branches overlaid, and a 2D codim-2 diagram in \((\\bar\\eta,\\, \\Delta)\) showing the BT / CP / GH structure.

fig, ax = plt.subplots(figsize=(7, 4.5))
plot_1d_diagram(ax, 'eta_branch', 'lc_branch',
                 title=r'1D bifurcation diagram, $\tau_r = 10$')
plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(7, 5))
plot_2d_diagram(ax, codim2_curves, 'eta_branch',
                 title=r'codim-2 bifurcation diagram, $\tau_r = 10$')
plt.tight_layout()
plt.show()

Reading the 2D diagram:

• Blue fold curves traced from LP1 / LP2, both directions each. The cusp CP where these collide marks the closing of the bistable wedge.

• Orange Hopf curves traced from HB1 / HB2. The BT (Bogdanov-Takens) point is where the Hopf curve from HB1 meets the fold manifold — two fold branches pass through BT and one Hopf branch terminates at it (a homoclinic curve also emerges from BT, but tracking it requires auto-07p’s HomCont package IPS=9 and is left as a manual follow-up).

• The GH (generalised-Hopf / Bautin) point on the Hopf curve from HB2 is where the first Lyapunov coefficient changes sign — supercritical Hopfs on one side, subcritical on the other.

Step 5: Continuation + codim-2 + PD cascade at \(\\tau_r = 0.1\)

Same pipeline at the second user point UZ2, with the addition of a continue_period_doubling_bf() call after the LC continuation — at \(\\tau_r = 0.1 \\ll \\tau_d\) the LC born at HB2 carries at least one period-doubling point.

eta_lo_sols, eta_lo_cont = ode.run(
    origin='tau_r_branch', starting_point='UZ2', name='eta_branch_lo',
    ICP='p/qif_biexp_sfa_op/eta', bidirectional=True,
    RL0=-8.0, RL1=2.0,
    IPS=1, ILP=1, ISP=2, ISW=1, NTST=400, NCOL=4,
    NMX=2000, NPR=10, DS=1e-4, DSMIN=1e-8, DSMAX=5e-2,
    ITMX=40, ITNW=40, NWTN=12,
)
print("\neta scan (tau_r=0.1):", dict(eta_lo_sols['bifurcation'].value_counts()))

lc_lo_sols, lc_lo_cont = ode.run(
    origin='eta_branch_lo', starting_point='HB2', name='lc_branch_lo',
    IPS=2, ISP=2, ISW=-1,
    ICP=['p/qif_biexp_sfa_op/eta', 11],
    NMX=2000, NPR=20, DS=1e-3, DSMIN=1e-9, DSMAX=5e-2,
    bidirectional=True, get_period=True, STOP=["BP3", "LP5"]
)
print("LC (tau_r=0.1):", dict(lc_lo_sols['bifurcation'].value_counts()))

codim2_curves_lo = run_codim2_in_delta('eta_branch_lo', name_prefix='delta01')
print("codim-2 curves at tau_r=0.1:")
for key, bif_type in codim2_curves_lo:
    print(f"  {key} ({bif_type}): "
          f"{dict(ode.get_summary(key)['bifurcation'].value_counts())}")

# Trace the locus of the *first PD point* on the LC in (eta, Delta), the
# same way we'd trace any other codim-1 manifold — `codim2_search` treats
# PD as a codim-1 starting point when `periodic=True` (so ICP gets PAR(11)
# appended for the period). This gives the red PD curve in the 2D diagram.
pd_locus = codim2_search(
    pyauto_instance=ode,
    params=['p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/Delta'],
    starting_points=['PD1'], origin='lc_branch_lo',
    max_recursion_depth=0,
    periodic=True,
    NMX=1500, NPR=20,
    DSMIN=1e-9, DSMAX=5e-2,
    bidirectional=True,
    UZSTOP={'p/qif_biexp_sfa_op/Delta': [0.0, 4.0]},
    name='pd_locus',
)
print("PD-locus continuations:", list(pd_locus.keys()))

# Chase the period-doubling cascade itself with the dedicated helper.
# Unlike `codim2_search` (which traces a locus in 2 params), this one
# does the auto-07p "branch-switch onto each new doubled LC and repeat"
# workflow — same recipe as the c.lor.2 / c.lor.3 manual cascade in
# auto-07p's lor demo. Each cascade step is a separate 1-parameter LC
# continuation in `eta` (with PAR(11) for the doubled period).
cascade_names, _ = continue_period_doubling_bf(
    solution=ode.results[ode.get_continuation('lc_branch_lo').key],
    continuation=lc_lo_cont,
    pyauto_instance=ode,
    icp='p/qif_biexp_sfa_op/eta',
    max_iter=4,
    NMX=2000, NPR=20, DS=1e-3, DSMIN=1e-9, DSMAX=5e-2,
    STOP=["BP3", "LP5"],
)
print("PD cascade LCs:", cascade_names)

Step 6: Figures at \(\\tau_r = 0.1\)

Same shape as Step 4. The 1D diagram now overlays each LC born along the PD cascade (each doubled cycle is its own continuation in eta). The 2D diagram overlays the PD locus from codim2_search() on the usual codim-2 backbone.

fig, ax = plt.subplots(figsize=(7, 4.5))
plot_1d_diagram(ax, 'eta_branch_lo', 'lc_branch_lo',
                 title=r'1D bifurcation diagram, $\tau_r = 0.1$',
                 cascade_names=cascade_names)
plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(7, 5))
plot_2d_diagram(ax, codim2_curves_lo, 'eta_branch_lo',
                 pd_names=list(pd_locus.keys()),
                 title=r'codim-2 + PD-locus bifurcation diagram, '
                       r'$\tau_r = 0.1$')
plt.tight_layout()
plt.show()

Step 7: Failure modes and what to expect on unfamiliar models

Both helpers wrap a series of nested ODESystem.run calls. If any individual sub-run raises inside auto-07p (common on low-quality parameter regimes or aggressive step sizes), the failure surfaces as a UserWarning and the search continues with the remaining starting points rather than aborting. Read the warnings carefully: they cite the auto-07p exception type and message, the sub-run label that failed, and the kwargs hook (kwargs_1D_lc_cont, kwargs_2D_cont, …) you can use to override the default constants for that path.

Step 8: Clean up

ode.close_session(clear_files=True)