Homoclinic Continuation in a QIF Mean-Field Model

This example demonstrates PyCoBi’s support for auto-07p’s HomCont extension (\(IPS=9\), AUTO manual Ch. 20) via three new entry points:

• The 'hom' scenario in auto_constants=(...) generates a c.hom file pre-configured with HomCont’s eight specific constants (NUNSTAB, NSTAB, IEQUIB, ITWIST, ISTART, IREV, IFIXED, IPSI). These flow through from_template(...) / run(...) as ordinary kwargs.

• ODESystem.extract_orbit_to_dat() writes a labelled limit-cycle profile from auto-07p’s solution object straight into the .dat format HomCont reads when ISTART=1.

• ODESystem.continue_homoclinic() orchestrates the full pipeline: extract the near-homoclinic orbit, seed a HomCont continuation in two parameters, append the requested PSI test-function PARs (PAR(20 + IPSI[j])) to ICP so they land in the summary, then post-process the result to flag every PSI zero-crossing as a 'SNIC' bifurcation in the bifurcation column.

  Note

  The 'SNIC' label is the historical PyCoBi name and not the standard codim-1 SNIC bifurcation (saddle-node on invariant cycle). With the default IPSI=(15, 16) we are detecting the codim-2 non-central SNIC of Nechyporenko, Ashwin & Tsaneva-Atanasova 2026 (arXiv:2412.12298): a homoclinic orbit to a saddle-node equilibrium whose return trajectory enters the saddle-node along its stable manifold rather than along the central (zero-eigenvalue) direction. This is a vertex terminating a curve of codim-1 SNIC bifurcations, not the codim-1 SNIC itself.

The model is the bi-exponential QIF mean field with spike-frequency adaptation from the Automated Codim-2 Search example. At \(\tau_r = \tau_d = 10\) (the alpha-kernel limit, Gast 2020 [1]) the upper Hopf (HB2) gives birth to a stable limit cycle whose period diverges as \(\bar\eta\) is pushed toward a saddle-loop homoclinic. We pick the near-homoclinic LC up at its endpoint and continue the homoclinic curve in \((\bar\eta,\, \Delta)\).

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.

[2]

K. Nechyporenko, P. Ashwin, K. Tsaneva-Atanasova (2026) A novel route to oscillations via non-central SNICeroclinic bifurcation: unfolding the separatrix loop between a saddle-node and a saddle. SIAM J. Appl. Dyn. Syst., to appear. arXiv:2412.12298.

Step 1: load the QIF-SFA model + run the equilibrium continuation

Same setup as automated_continuation.py: alpha-kernel-equivalent regime, four-dimensional state (r, v, A, B). auto_constants now also requests 'hom' so the c.hom file is ready for use.

from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt

from pycobi import ODESystem

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.5,
        'p/qif_biexp_sfa_op/tau_r': 10.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_constants=('ivp', 'eq', 'lc', 'hom'),
)

# Equilibrium continuation in eta — locates HB1 / HB2 / LP1 / LP2.
# UZR pins :math:`\bar\eta = -5.394` (a value slightly different from the
# LC's homoclinic terminus at :math:`\bar\eta \approx -5.39358`, see the
# note in Step 4) so the saddle equilibrium's state values are available
# as labelled UZ solutions later when seeding HomCont.
eta_sols, _ = ode.run(
    starting_point='EP2', name='eta_branch',
    c='eq', ICP='p/qif_biexp_sfa_op/eta', bidirectional=True,
    RL0=-8.0, RL1=2.0,
    NMX=2000, NPR=10, DS=1e-4, DSMIN=1e-8, DSMAX=5e-2,
    ITMX=40, ITNW=40, NWTN=12, NTST=400, NCOL=4,
    UZR={'p/qif_biexp_sfa_op/eta': [-5.394]},
)
print("eta scan:", dict(eta_sols['bifurcation'].value_counts()))

Step 2: continue the limit cycle from HB2 toward the homoclinic

The Hopf at the upper edge of the spiking regime gives birth to an unstable limit cycle whose period grows as \(\bar\eta\) walks the branch toward a saddle-loop homoclinic. At \(\alpha=0.5,\ \Delta=2.0, \ J=15\sqrt{\Delta},\ \tau_r=\tau_d=10\) the LC reaches the homoclinic at \(\bar\eta \approx -5.394\) with a period of ~388.

Note

For high-period limit-cycle continuations approaching a homoclinic, auto-07p’s Newton tolerances and bifurcation-detection tolerance need to be tighter than the global defaults — otherwise the test functions used to flag LP / PD / BP along the LC become dominated by numerical noise and auto-07p reports spurious bifurcations. PyCoBi’s 'lc' scenario ships with EPSL = EPSU = 1e-7 and EPSS = 1e-5. For this LC the progression is:

• Defaults (NTST=50, EPSL=1e-7): ~24 LPs, ~8 BPs, ~4 PDs.

• NTST=400, EPSL=1e-9: ~24 LPs drop to <10, but 3 PDs remain.

• NTST=800, EPSL=1e-9 (the values below): all PDs disappear.

The remaining ~800 BP labels concentrate at the homoclinic \(\bar\eta\) and are not numerical noise — they’re auto-07p detecting the LC branch nearly intersecting the saddle equilibrium’s branch, which is exactly what a saddle-loop homoclinic is. We keep them and ignore them in the 1D plot.

lc_sols, _ = ode.run(
    starting_point='HB2', name='lc_branch',
    c='lc', origin='eta_branch',
    ICP=['p/qif_biexp_sfa_op/eta', 11],
    NMX=8000, NPR=20, DS=1e-3, DSMIN=1e-12, DSMAX=5e-2,
    EPSL=1e-9, EPSU=1e-9, EPSS=1e-7,
    NTST=800, NCOL=4,
    bidirectional=False, get_period=True,
)
periods = np.asarray(lc_sols[('PAR(11)', '')], dtype=float)
print(f"LC bifurcations (raw): {dict(lc_sols['bifurcation'].value_counts())}")
print(f"LC period grows from {periods.min():.1f} to {periods.max():.1f}")

# Heuristically tag the homoclinic terminus as 'HC' on the LC summary:
# the max-to-min period ratio is ~22 here, well above the default 5x
# threshold for the heuristic to fire.
hc_idx = ode.label_homoclinic_terminus('lc_branch')
if hc_idx is not None:
    _eta_col = ('eta', '') if ('eta', '') in lc_sols.columns \
        else ('p/qif_biexp_sfa_op/eta', '')
    print(f"HC marker placed at row {hc_idx}: "
          f"eta = {float(lc_sols[_eta_col].iloc[hc_idx]):+.5f}, "
          f"period = {periods[hc_idx]:.1f}")

# 1D bifurcation diagram of the eta scan + LC envelope + HC marker.
fig, ax = plt.subplots(figsize=(9, 4))
ode.plot_continuation('p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/r',
                      cont='eta_branch', ax=ax,
                      line_color_stable='#1F77B4', line_color_unstable='#1F77B4',
                      bifurcation_legend=False,
                      ignore=['UZ', 'BP', 'EP'],
                      label='equilibrium')
ode.plot_continuation('p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/r',
                      cont='lc_branch', ax=ax,
                      line_color_stable='#FF7F0E', line_color_unstable='#FF7F0E',
                      bifurcation_legend=False,
                      ignore=['UZ', 'BP', 'EP', 'RG'],
                      label='limit cycle')
# Plot the HC marker directly — plot_continuation's bifurcation table
# doesn't know about our custom 'HC' label yet, so we drop a marker
# manually at the labelled row.
if hc_idx is not None:
    eta_col = ('eta', '') if ('eta', '') in lc_sols.columns \
        else ('p/qif_biexp_sfa_op/eta', '')
    r_col = [c for c in lc_sols.columns
             if isinstance(c, tuple) and c[0] == 'r'][0]
    eta_hc = float(lc_sols[eta_col].iloc[hc_idx])
    r_hc = lc_sols[r_col].iloc[hc_idx]
    if hasattr(r_hc, '__len__'):
        r_hc = float(np.max(r_hc))
    else:
        r_hc = float(r_hc)
    ax.scatter([eta_hc], [r_hc], marker='X', s=120,
               c='#D62728', edgecolor='k', linewidth=0.8, zorder=10,
               label=f'HC (period ≈ {periods[hc_idx]:.0f})')
ax.set_xlabel(r'$\bar\eta$')
ax.set_ylabel(r'$r$')
ax.set_title('1D bifurcation diagram — LC born at HB2 terminates at HC')
ax.legend(loc='best')
plt.tight_layout()
plt.show()

Step 3: extract the near-homoclinic orbit and write it as a .dat file

Auto-07p’s HomCont uses an existing orbit as the initial guess; the orbit lives on disk as a whitespace-separated .dat file of (time, state, state, …) columns. extract_orbit_to_dat() reads the per-mesh state values straight off auto-07p’s solution object and writes the .dat for us — no need to have run the LC continuation with get_timeseries=True.

Inspecting the seed profile confirms it has the characteristic “slow near saddle, fast excursion” shape of a near-homoclinic orbit.

seed_path = ode.extract_orbit_to_dat(
    cont='lc_branch', point='EP1', path='lc_seed', n_points=201,
)
data = np.loadtxt(seed_path)
print(f"seed written to {seed_path.name}: shape {data.shape}")
print(f"r profile: min={data[:,1].min():.3f}, max={data[:,1].max():.3f}, "
      f"median={np.median(data[:,1]):.3f}")

fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(data[:, 0], data[:, 1], color='#FF7F0E', lw=1.5)
ax.set_xlabel(r'$t / T$')
ax.set_ylabel(r'$r(t)$')
ax.set_title(f'Near-homoclinic orbit profile (period $\\approx${periods.max():.0f})')
plt.tight_layout()
plt.show()

Step 4: find the saddle equilibrium (the homoclinic’s target)

A homoclinic orbit by definition leaves and returns to a saddle equilibrium. At \(\bar\eta \approx -5.394\) the QIF-SFA has three equilibria: a stable down-state, an unstable middle saddle, and a stable up-state. We pull the equilibrium values straight off the UZ labels we planted on the eta_branch continuation in Step 1 — the saddle is identified as the unstable middle equilibrium whose r coordinate is between the LC orbit’s min and median r (i.e. the “slow-point” the orbit hovers near).

Note

The UZR \(\bar\eta = -5.394\) doesn’t exactly match the LC’s EP1 \(\bar\eta = -5.39358\) — and that small mismatch is intentional. Reading the saddle off the LC’s orbit slow-point puts us AT the codim-2 non-central SNIC vertex (PSI(15) ≈ 0), where HomCont has no direction to advance from. A saddle taken from the eta-branch at a nominally-nearby \(\bar\eta\) gives a slightly off-vertex starting orbit (PSI(15) of order 0.3), and HomCont can walk the codim-1 homoclinic curve forward.

saddle_state = None
for uz_label in ('UZ1', 'UZ2', 'UZ3'):
    try:
        s, _, _ = ode.get_solution(cont='eta_branch', point=uz_label)
        if hasattr(s, 'b') and isinstance(getattr(s, 'b', None), dict):
            s = s.b['solution']
        coords = {c: float(np.asarray(s[c]).ravel()[0]) for c in s.coordnames}
        # the saddle is the unstable middle equilibrium; r matches LC's min r
        if abs(coords['r'] - data[:, 1].min()) < 0.05:
            saddle_state = coords
            print(f"saddle equilibrium at {uz_label}: r={coords['r']:.4f}, "
                  f"v={coords['v']:.4f}, A={coords['A']:.4f}, B={coords['B']:.4f}")
            break
    except (KeyError, AttributeError):
        pass

Step 5: continue the homoclinic curve in \((\bar\eta,\, \Delta)\)

continue_homoclinic() now wraps the full LC-to-HomCont seed preparation pipeline:

1. extract_orbit_to_dat reads the labelled LC profile straight off auto-07p’s solution object.

2. phase_shift_to=saddle_state re-rolls the orbit so the point closest to the saddle sits at \(t = 0\) — without this, the LC’s arbitrary phase choice leaves HomCont’s stable / unstable manifold projections misaligned and Newton MXs at step 2.

3. saddle_state also populates PAR(11 + i) with the saddle’s coordinates (one per state variable, in uname order). Combined with IEQUIB = 0, this pins HomCont’s equilibrium at the precomputed saddle rather than asking Newton to find it.

4. The seed’s parameter values (from the auto solution’s PAR vector) are forwarded via PAR={...} so the homoclinic starts at the correct \((\bar\eta, \Delta, \alpha, \dots)\).

5. After the run, _flag_psi_zero_crossings scans PAR(35) = PSI(15) and PAR(36) = PSI(16) for sign changes and marks each as 'SNIC' in the bifurcation column. Per the docstring note above: these are codim-2 non-central SNIC bifurcations (Nechyporenko et al. 2026, arXiv:2412.12298) — endpoints of a curve of standard SNIC bifurcations — not the codim-1 SNIC itself.

Note

Eigenvalue-split kwargs (NUNSTAB, NSTAB) still need to match the model’s saddle structure; for this 4-D QIF-SFA the saddle has one unstable + three stable directions, hence NUNSTAB=1, NSTAB=3. Adapting the recipe to a different model means picking these from the local eigenvalue spectrum of the saddle.

snic_sols, _ = ode.continue_homoclinic(
    origin='lc_branch', starting_point='EP1',
    ICP=['p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/Delta'],
    NUNSTAB=1, NSTAB=3,
    IEQUIB=0, ITWIST=0,
    IPSI=(15, 16),
    saddle_state=saddle_state,
    n_points=201,
    NTST=100, NCOL=4, IAD=1, ISP=0, ILP=0,
    NMX=1000, NPR=5, IID=2, ITMX=10, ITNW=7, NWTN=3,
    DS=0.001, DSMIN=1e-5, DSMAX=0.05,
    UZSTOP={'p/qif_biexp_sfa_op/Delta': [0.01, 4.0]},
    name='homoclinic',
)
bifs = dict(snic_sols['bifurcation'].value_counts())
print(f"\nHomCont curve in (eta, Delta): {bifs}, {len(snic_sols)} points")
n_snic = bifs.get('SNIC', 0)
if n_snic:
    snic_col = ('bifurcation', '')
    print(f"\n{n_snic} non-central SNIC bifurcation(s) detected on the "
          f"homoclinic curve (Nechyporenko et al. 2026, arXiv:2412.12298):")
    snic_rows = snic_sols[snic_sols[snic_col] == 'SNIC']
    eta_col = [c for c in snic_sols.columns
               if isinstance(c, tuple) and 'eta' in c[0]][0]
    delta_col = [c for c in snic_sols.columns
                 if isinstance(c, tuple) and 'Delta' in c[0]][0]
    for _, row in snic_rows.iterrows():
        print(f"  eta = {float(row[eta_col]):+.5f}, Delta = {float(row[delta_col]):.4f}")
hom_curve_available = len(snic_sols) > 2

Step 6: codim-2 continuation of HB / LP curves in \((\bar\eta,\, \Delta)\)

PyCoBi’s codim2_search() wraps the standard auto-07p 2-parameter continuation of codim-1 starting points. Tracking the four equilibrium bifurcations we located on eta_branch (HB1, HB2, LP1, LP2) in \((\bar\eta,\, \Delta)\) gives the codim-1 backbone we need to plot alongside the homoclinic locus from Step 5.

Note

NPR=5 is deliberately tighter than the typical NPR=20 you’d use for a 1-parameter scan. NPR is auto-07p’s stored-point spacing: with the default 20, the first stored row on a codim-2 branch lies ~20 continuation steps (up to ~1 unit of arclength at DSMAX=5e-2) past the LP/HB starting point on eta_branch, leaving a visible gap in the \((\bar\eta,\, \Delta)\) plot between each codim-1 source on eta_branch and the start of its own codim-2 curve. Lowering to NPR=5 shrinks the first-step offset to ~0.01 units — visually contiguous on the plot.

from pycobi.automated_continuation import codim2_search

codim2_curves = []
for sp in ('LP1', 'LP2', 'HB1', 'HB2'):
    bif_type = 'fold' if sp.startswith('LP') else 'Hopf'
    for ds in (1e-3, -1e-3):
        try:
            res = codim2_search(
                pyauto_instance=ode,
                params=['p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/Delta'],
                starting_points=[sp],
                origin='eta_branch',
                name=f'codim2_{sp}_{"pos" if ds > 0 else "neg"}',
                max_recursion_depth=0,
                NMX=3000, NPR=5, DSMIN=1e-9, DSMAX=5e-2, DS=ds,
                RL0=-15.0, RL1=5.0,
                bidirectional=False,
                UZSTOP={'p/qif_biexp_sfa_op/Delta': [0.0, 4.0]},
            )
            curve_name = next(iter(res))
            codim2_curves.append((curve_name, bif_type))
        except Exception as exc:
            print(f"  codim-2 {sp} DS={ds:+g} skipped: "
                  f"{type(exc).__name__}: {str(exc)[:80]}")

Step 7: 2D bifurcation diagram with all the codim-1 curves + homoclinic

Stack everything in the \((\bar\eta,\, \Delta)\) plane: fold curves (blue), Hopf curves (orange), homoclinic curve (red), and any detected non-central SNIC points (red stars). The result is the complete codim-1 bifurcation portrait of the QIF-SFA model in this slice.

FOLD_COLOR = '#1F77B4'
HOPF_COLOR = '#FF7F0E'
HOM_COLOR  = '#D62728'

fig, ax = plt.subplots(figsize=(9, 6))
fold_labelled, hopf_labelled = False, False
for curve_name, bif_type in codim2_curves:
    color = FOLD_COLOR if bif_type == 'fold' else HOPF_COLOR
    if bif_type == 'fold':
        label = None if fold_labelled else 'fold of equilibria'
        fold_labelled = True
    else:
        label = None if hopf_labelled else 'Hopf'
        hopf_labelled = True
    try:
        ode.plot_continuation(
            'p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/Delta',
            cont=curve_name, 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', 'EP'],
            label=label,
        )
    except (KeyError, ValueError):
        pass

if hom_curve_available:
    ode.plot_continuation(
        'p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/Delta',
        cont='homoclinic', ax=ax,
        line_color_stable=HOM_COLOR, line_color_unstable=HOM_COLOR,
        line_style_stable='solid', line_style_unstable='dashed',
        bifurcation_legend=False, get_stability=False,
        ignore=['UZ', 'EP', 'RG'],
        label='homoclinic',
    )
    # Plant explicit markers at the codim-2 non-central SNIC points
    # (the PSI(15)/(16) zero-crossings flagged as 'SNIC' in Step 5).
    snic_col = ('bifurcation', '')
    eta_col = [c for c in snic_sols.columns
               if isinstance(c, tuple) and 'eta' in c[0]][0]
    delta_col = [c for c in snic_sols.columns
                 if isinstance(c, tuple) and 'Delta' in c[0]][0]
    snic_rows = snic_sols[snic_sols[snic_col] == 'SNIC']
    if len(snic_rows):
        ax.scatter(snic_rows[eta_col].to_numpy(dtype=float),
                   snic_rows[delta_col].to_numpy(dtype=float),
                   marker='*', s=200, c=HOM_COLOR,
                   edgecolor='k', linewidth=0.8, zorder=10,
                   label='non-central SNIC')

ax.set_xlabel(r'$\bar\eta$')
ax.set_ylabel(r'$\Delta$')
ax.set_xlim(-12.0, 5.0)
ax.set_ylim(0.0, 4.0)
ax.set_title('Complete codim-1 bifurcation portrait in '
             r'$(\bar\eta,\, \Delta)$')
ax.legend(loc='best')
plt.tight_layout()
plt.show()

Reference table — HomCont PSI test functions

The IPSI list above selects which “PSIHO” test functions auto-07p computes along the homoclinic continuation. Their values land in PAR(20 + IPSI[j]) and zero crossings flag codim-2 phenomena on the homoclinic curve. PyCoBi keeps a reference table so users picking IPSI=[...] have the meanings close at hand:

for j, desc in ODESystem.HOMCONT_PSI_NAMES.items():
    marker = " ←" if j in (15, 16) else ""
    print(f"  PSI({j:>2}) → PAR({20+j}): {desc}{marker}")