Parameter Continuation and Bifurcation Detection with TVBO

In this tutorial, you will learn how to perform a 1D numerical parameter continuation with automatic fold bifurcation detection using TVBO for model specification and PyRates/PyCoBi for the continuation analysis.

import os

import matplotlib.pyplot as plt
from pycobi import ODESystem
from pyrates import CircuitTemplate
from IPython.display import Markdown
from tvbo import Dynamics

qif_model = Dynamics.from_string(
    """
name: QIF
description: Quadratic Integrate-and-Fire (QIF) population mean-field model.
parameters:
    tau:
        value: 1.0
        description: Population time constant
    eta:
        value: -5.0
        description: mean of a Lorentzian distribution over the neural excitability
    Delta:
        value: 1.0
        description: half-width of a Lorentzian distribution over the neural excitability
    J:
        value: 15.0
        description: strength of the recurrent coupling inside the population

state_variables:
    r:
        description:  Average firing rate
        equation:
            rhs: Delta/(tau*pi) + 2*r*v
        initial_value: 0.1
    v:
        description: Average membrane potential
        equation:
            rhs: v**2 + eta + J*r*tau - (pi*r*tau)**2
        initial_value: -2.0

references:
    - Montbriot2015
"""
)

Markdown(qif_model.generate_report())

QIF

Quadratic Integrate-and-Fire (QIF) population mean-field model.

autonomous: True; modes: 1; state variables: 2; parameters: 4.

State Equations

\[ \dot{r} = 2*r*v + \frac{\Delta}{\pi*\tau} \] \[ \dot{v} = \eta + v^{2} + J*r*\tau - \pi^{2}*r^{2}*\tau^{2} \]

State Variables

Variable Initial Value Equation Flags Description
\(r\) 0.1 differential (order 1) recorded Average firing rate
\(v\) -2 differential (order 1) recorded Average membrane potential

Parameters

Parameter Value Description
\(\tau\) 1 Population time constant
\(\eta\) -5 mean of a Lorentzian distribution over the neural excitability
\(\Delta\) 1 half-width of a Lorentzian distribution over the neural excitability
\(J\) 15 strength of the recurrent coupling inside the population

References

Citation key ‘Montbriot2015’ not found.

# Export to PyRates YAML format
os.makedirs('_qif_model', exist_ok=True)
with open('_qif_model/__init__.py', 'w') as f:
    f.write('')

qif_model.to_yaml(format="pyrates", filepath='_qif_model/qif.yaml')
print("Model exported to PyRates format")
Model exported to PyRates format

Step 2: Generate Fortran Routines for AUTO-07p

Next, we translate the model into a Fortran file containing all subroutines required by AUTO-07p. This requires using the Fortran backend of PyRates with vectorize=False:

# Load the TVBO-exported model into PyRates
qif = CircuitTemplate.from_yaml("_qif_model.qif.QIF_circuit")

# Generate Fortran file with AUTO-07p subroutines
# Setting auto=True creates two files: .f90 (Fortran code) and c.ivp (AUTO constants)
qif.get_run_func(
    func_name='qif_rhs',
    file_name='qif',
    step_size=1e-4,
    auto=True,
    backend='fortran',
    solver='scipy',
    vectorize=False,
    float_precision='float64'
)
Compilation Progress
--------------------
    (1) Translating the circuit template into a networkx graph representation...
        ...finished.
    (2) Preprocessing edge transmission operations...
        ...finished.
    (3) Parsing the model equations into a compute graph...
        ...finished.
    Model compilation was finished.
(<fortran function qif_rhs>,
 (array(0.),
  array([ 0.1, -2. ]),
  array([0., 0.]),
  array(1.),
  array(-5.),
  array(1.),
  array(15.)),
 ('t',
  'y',
  'dy',
  'p/QIF_op/tau',
  'p/QIF_op/eta',
  'p/QIF_op/Delta',
  'p/QIF_op/J'),
 {'p/QIF_op/r': 0, 'p/QIF_op/v': 1})

Calling get_run_func with auto=True creates two files. The first is a .f90 file containing the Fortran subroutines:

f = open('qif.f90', 'r')
print(f.read())
f.close()
module qif

double precision :: PI = 4.0*atan(1.0)
complex :: I = (0.0, 1.0)

contains


subroutine qif_rhs(t,y,dy,tau,eta,Delta,J)

implicit none

double precision :: r
double precision :: v
double precision, intent(in) :: tau
double precision, intent(in) :: eta
double precision, intent(in) :: Delta
double precision, intent(in) :: J
double precision, intent(in) :: t
double precision, intent(in) :: y(2)
double precision, intent(inout) :: dy(2)


r = y(1)
v = y(2)

dy(1) = Delta/(pi*tau) + 2*r*v
dy(2) = J*r*tau + eta - pi**2*r**2*tau**2 + v**2

end subroutine


end module


subroutine func(ndim,y,icp,args,ijac,dy,dfdu,dfdp)

use qif
implicit none
integer, intent(in) :: ndim, icp(*), ijac
double precision, intent(in) :: y(ndim), args(*)
double precision, intent(out) :: dy(ndim)
double precision, intent(inout) :: dfdu(ndim,ndim), dfdp(ndim,*)

call qif_rhs(args(14), y, dy, args(1), args(2), args(3), args(4))

if (ijac .eq. 0) return

dfdu(1,1) = 2*y(2)
dfdu(1,2) = 2*y(1)
dfdu(2,1) = args(1)*args(4) - 2*y(1)*args(1)**2*pi**2
dfdu(2,2) = 2*y(2)

if (ijac .eq. 1) return

dfdp(1,1) = -args(3)/(args(1)**2*pi)
dfdp(1,3) = 1/(args(1)*pi)
dfdp(2,1) = args(4)*y(1) - 2*args(1)*y(1)**2*pi**2
dfdp(2,2) = 1
dfdp(2,4) = args(1)*y(1)

end subroutine func


subroutine stpnt(ndim, y, args, t)

implicit None
integer, intent(in) :: ndim
double precision, intent(inout) :: y(ndim), args(*)
double precision, intent(in) :: t

args(1) = 1.0  ! tau
args(2) = -5.0  ! eta
args(3) = 1.0  ! Delta
args(4) = 15.0  ! J
y(1) = 0.1  ! r
y(2) = -2.0  ! v

end subroutine stpnt



subroutine bcnd
end subroutine bcnd


subroutine icnd
end subroutine icnd


subroutine fopt
end subroutine fopt


subroutine pvls
end subroutine pvls

The second file contains the AUTO-07p parameters that determine how it performs parameter continuations:

f = open('c.ivp', 'r')
print(f.read())
f.close()
parnames = {1: 'tau', 2: 'eta', 3: 'Delta', 4: 'J'}
unames = {1: 'r', 2: 'v'}
NDIM = 2
NPAR = 4
IPS = -2
ILP = 0
ICP = [14]
NTST = 1
NCOL = 3
IAD = 0
ISP = 0
ISW = 1
IPLT = 0
NBC = 0
NINT = 0
NMX = 9000
NPR = 20
MXBF = 10
IID = 2
ITMX = 2
ITNW = 5
NWTN = 2
JAC = 1
EPSL = 1e-06
EPSU = 1e-06
EPSS = 0.0001
IRS = 0
DS = 0.0001
DSMIN = 1e-08
DSMAX = 0.01
IADS = 1
THL = {}
THU = {}
UZR = {}
UZSTOP = {}

The default parameters allow solving the initial value problem (performing simulations). For detailed explanation of these parameters, see the AUTO-07p documentation.

Step 3: Generate a PyCoBi Instance

Now that the model equations are compiled, we generate an instance of pycobi.ODESystem, which provides an interface to AUTO-07p:

qif_auto = ODESystem(eq_file="qif", working_dir=None, init_cont=False)

Now we can use all tools provided by AUTO-07p to investigate how the model reacts to changes in its parameterization.

Part 2: Performing Parameter Continuations

Step 1: Time Continuation

In parameter continuations, we must start from an equilibrium or periodic orbit. We run a time continuation to let the system converge to an equilibrium:

# Time continuation: integrate until system converges to equilibrium
t_sols, t_cont = qif_auto.run(
    c='ivp', name='time',
    DS=1e-4, DSMIN=1e-10, DSMAX=1e-2,
    EPSL=1e-08, EPSU=1e-08, EPSS=1e-06,
    NMX=1000,
    UZR={14: 4.0},   # Create marker when time (PAR 14) reaches 4.0
    STOP={'UZ1'}     # Stop at first user marker
)

# Plot time evolution
qif_auto.plot_continuation('PAR(14)', 'r', cont='time')
plt.xlabel('Time')
plt.ylabel('r (firing rate)')
plt.title('Time Evolution to Equilibrium')
plt.show()
/opt/homebrew/bin/gfortran -arch arm64 -O -fopenmp -O -c qif.f90 -o qif.o
/opt/homebrew/bin/gfortran -arch arm64 -O -fopenmp -O qif.o -o qif.exe /Applications/auto-07p/lib/*.o
Starting qif ...
ld: warning: duplicate -rpath '/Library/Developer/CommandLineTools/SDKs/MacOSX26.4.sdk/usr/lib' ignored

  BR    PT  TY  LAB    PAR(14)       L2-NORM          r             v       
   1     1  EP    1   0.00000E+00   2.00250E+00   1.00000E-01  -2.00000E+00
   1    20        2   2.97235E-02   1.99165E+00   9.77522E-02  -1.98925E+00
   1    40        3   7.34782E-02   1.97917E+00   9.49818E-02  -1.97689E+00
   1    60        4   1.25882E-01   1.96849E+00   9.23394E-02  -1.96633E+00
   1    80        5   1.93602E-01   1.95964E+00   8.97448E-02  -1.95759E+00
   1   100        6   2.90740E-01   1.95324E+00   8.71506E-02  -1.95129E+00
   1   120        7   4.58082E-01   1.95095E+00   8.45402E-02  -1.94912E+00
   1   140        8   6.58061E-01   1.95319E+00   8.29707E-02  -1.95143E+00
   1   160        9   8.58037E-01   1.95619E+00   8.21742E-02  -1.95446E+00
   1   180       10   1.05802E+00   1.95860E+00   8.17431E-02  -1.95690E+00
   1   200       11   1.25801E+00   1.96029E+00   8.14983E-02  -1.95860E+00
   1   220       12   1.45801E+00   1.96141E+00   8.13547E-02  -1.95972E+00
   1   240       13   1.65801E+00   1.96212E+00   8.12688E-02  -1.96043E+00
   1   260       14   1.85801E+00   1.96257E+00   8.12168E-02  -1.96088E+00
   1   280       15   2.05801E+00   1.96284E+00   8.11850E-02  -1.96116E+00
   1   300       16   2.25801E+00   1.96302E+00   8.11656E-02  -1.96134E+00
   1   320       17   2.45801E+00   1.96312E+00   8.11536E-02  -1.96145E+00
   1   340       18   2.65801E+00   1.96319E+00   8.11462E-02  -1.96151E+00
   1   360       19   2.85801E+00   1.96323E+00   8.11417E-02  -1.96155E+00
   1   380       20   3.05801E+00   1.96326E+00   8.11389E-02  -1.96158E+00
   1   400       21   3.25801E+00   1.96327E+00   8.11372E-02  -1.96159E+00
   1   420       22   3.45801E+00   1.96328E+00   8.11361E-02  -1.96160E+00
   1   440       23   3.65801E+00   1.96329E+00   8.11355E-02  -1.96161E+00
   1   460       24   3.85801E+00   1.96329E+00   8.11351E-02  -1.96161E+00
   1   475  UZ   25   4.00000E+00   1.96329E+00   8.11349E-02  -1.96162E+00

 Total Time    0.111E-01
qif ... done

Key parameters explained:

  • DS=1e-4: Initial step-size (in time units)
  • DSMIN=1e-10: Minimum step-size
  • DSMAX=1e-2: Maximum step-size
  • NMX=1000: Maximum number of continuation steps
  • UZR={14: 4.0}: Create user marker when parameter 14 (time) reaches 4.0
  • STOP={'UZ1'}: Stop at first user marker

The values for r and v converged to certain values, representing the firing rate \(r\) and membrane potential \(v\) at equilibrium.

Step 2: Continuation of \(\bar\eta\)

Now we perform the parameter continuation in \(\bar\eta\):

# Parameter continuation in eta (parameter 4)
eta_sols, eta_cont = qif_auto.run(
    origin=t_cont, starting_point='UZ1', name='eta',
    bidirectional=True,
    ICP=4,              # Continue parameter 4 (eta)
    RL0=-20.0,          # Lower bound for eta
    RL1=20.0,           # Upper bound for eta
    IPS=1,              # Equilibrium continuation
    ILP=1,              # Detect fold bifurcations
    ISP=2,              # Full automatic bifurcation detection
    ISW=1, NTST=400, NCOL=4, IAD=3, IPLT=0, NBC=0, NINT=0,
    NMX=2000, NPR=10, MXBF=5, IID=2, ITMX=40, ITNW=40, NWTN=12, JAC=0,
    EPSL=1e-06, EPSU=1e-06, EPSS=1e-04,
    DS=1e-4, DSMIN=1e-8, DSMAX=5e-2, IADS=1,
    THL={}, THU={}, UZR={}, STOP={}
)
Starting qif ...

  BR    PT  TY  LAB       J          L2-NORM          r             v       
   1    10       26   1.50384E+01   1.96237E+00   8.11728E-02  -1.96069E+00
   1    20       27   1.55267E+01   1.95050E+00   8.16686E-02  -1.94879E+00
   1    30       28   1.60265E+01   1.93808E+00   8.21939E-02  -1.93634E+00
   1    40       29   1.65264E+01   1.92537E+00   8.27385E-02  -1.92359E+00
   1    50       30   1.70262E+01   1.91235E+00   8.33039E-02  -1.91053E+00
   1    60       31   1.75260E+01   1.89900E+00   8.38917E-02  -1.89715E+00
   1    70       32   1.80258E+01   1.88530E+00   8.45037E-02  -1.88341E+00
   1    80       33   1.85256E+01   1.87123E+00   8.51420E-02  -1.86929E+00
   1    90       34   1.90254E+01   1.85674E+00   8.58091E-02  -1.85476E+00
   1   100       35   1.95252E+01   1.84181E+00   8.65077E-02  -1.83978E+00
   1   110  EP   36   2.00250E+01   1.82640E+00   8.72410E-02  -1.82431E+00

 Total Time    0.427E-02
qif ... done
Starting qif ...

  BR    PT  TY  LAB       J          L2-NORM          r             v       
   1    10       26   1.49616E+01   1.96422E+00   8.10962E-02  -1.96254E+00
   1    20       27   1.44733E+01   1.97581E+00   8.06189E-02  -1.97416E+00
   1    30       28   1.39734E+01   1.98744E+00   8.01457E-02  -1.98582E+00
   1    40       29   1.34736E+01   1.99884E+00   7.96871E-02  -1.99725E+00
   1    50       30   1.29737E+01   2.01002E+00   7.92423E-02  -2.00846E+00
   1    60       31   1.24738E+01   2.02100E+00   7.88105E-02  -2.01946E+00
   1    70       32   1.19739E+01   2.03178E+00   7.83910E-02  -2.03027E+00
   1    80       33   1.14740E+01   2.04238E+00   7.79830E-02  -2.04089E+00
   1    90       34   1.09742E+01   2.05280E+00   7.75861E-02  -2.05133E+00
   1   100       35   1.04743E+01   2.06305E+00   7.71996E-02  -2.06160E+00
   1   110       36   9.97436E+00   2.07314E+00   7.68229E-02  -2.07171E+00
   1   120       37   9.47446E+00   2.08307E+00   7.64557E-02  -2.08166E+00
   1   130       38   8.97456E+00   2.09285E+00   7.60975E-02  -2.09146E+00
   1   140       39   8.47465E+00   2.10248E+00   7.57478E-02  -2.10112E+00
   1   150       40   7.97474E+00   2.11198E+00   7.54062E-02  -2.11063E+00
   1   160       41   7.47483E+00   2.12135E+00   7.50725E-02  -2.12002E+00
   1   170       42   6.97491E+00   2.13058E+00   7.47462E-02  -2.12927E+00
   1   180       43   6.47500E+00   2.13970E+00   7.44271E-02  -2.13840E+00
   1   190       44   5.97508E+00   2.14869E+00   7.41148E-02  -2.14741E+00
   1   200       45   5.47516E+00   2.15757E+00   7.38091E-02  -2.15630E+00
   1   210       46   4.97523E+00   2.16634E+00   7.35097E-02  -2.16509E+00
   1   220       47   4.47531E+00   2.17499E+00   7.32164E-02  -2.17376E+00
   1   230       48   3.97538E+00   2.18355E+00   7.29288E-02  -2.18233E+00
   1   240       49   3.47546E+00   2.19200E+00   7.26469E-02  -2.19080E+00
   1   250       50   2.97553E+00   2.20036E+00   7.23704E-02  -2.19917E+00
   1   260       51   2.47559E+00   2.20862E+00   7.20991E-02  -2.20745E+00
   1   270       52   1.97566E+00   2.21679E+00   7.18329E-02  -2.21563E+00
   1   280       53   1.47573E+00   2.22487E+00   7.15714E-02  -2.22372E+00
   1   290       54   9.75790E-01   2.23287E+00   7.13147E-02  -2.23173E+00
   1   300       55   4.75853E-01   2.24078E+00   7.10625E-02  -2.23965E+00
   1   310       56  -2.40855E-02   2.24860E+00   7.08146E-02  -2.24749E+00
   1   320       57  -5.24025E-01   2.25635E+00   7.05710E-02  -2.25525E+00
   1   330       58  -1.02397E+00   2.26402E+00   7.03315E-02  -2.26293E+00
   1   340       59  -1.52391E+00   2.27161E+00   7.00959E-02  -2.27053E+00
   1   350       60  -2.02385E+00   2.27913E+00   6.98642E-02  -2.27806E+00
   1   360       61  -2.52380E+00   2.28658E+00   6.96362E-02  -2.28552E+00
   1   370       62  -3.02374E+00   2.29396E+00   6.94118E-02  -2.29291E+00
   1   380       63  -3.52369E+00   2.30127E+00   6.91909E-02  -2.30023E+00
   1   390       64  -4.02364E+00   2.30851E+00   6.89734E-02  -2.30748E+00
   1   400       65  -4.52358E+00   2.31569E+00   6.87592E-02  -2.31467E+00
   1   410       66  -5.02353E+00   2.32281E+00   6.85482E-02  -2.32179E+00
   1   420       67  -5.52348E+00   2.32986E+00   6.83404E-02  -2.32886E+00
   1   430       68  -6.02343E+00   2.33685E+00   6.81356E-02  -2.33586E+00
   1   440       69  -6.52339E+00   2.34378E+00   6.79337E-02  -2.34280E+00
   1   450       70  -7.02334E+00   2.35066E+00   6.77347E-02  -2.34968E+00
   1   460       71  -7.52329E+00   2.35748E+00   6.75385E-02  -2.35651E+00
   1   470       72  -8.02325E+00   2.36424E+00   6.73450E-02  -2.36328E+00
   1   480       73  -8.52320E+00   2.37095E+00   6.71541E-02  -2.37000E+00
   1   490       74  -9.02316E+00   2.37760E+00   6.69659E-02  -2.37666E+00
   1   500       75  -9.52311E+00   2.38420E+00   6.67801E-02  -2.38327E+00
   1   510       76  -1.00231E+01   2.39076E+00   6.65968E-02  -2.38983E+00
   1   520       77  -1.05230E+01   2.39726E+00   6.64159E-02  -2.39634E+00
   1   530       78  -1.10230E+01   2.40371E+00   6.62374E-02  -2.40280E+00
   1   540       79  -1.15229E+01   2.41011E+00   6.60611E-02  -2.40921E+00
   1   550       80  -1.20229E+01   2.41647E+00   6.58870E-02  -2.41557E+00
   1   560       81  -1.25229E+01   2.42278E+00   6.57151E-02  -2.42189E+00
   1   570       82  -1.30228E+01   2.42905E+00   6.55453E-02  -2.42817E+00
   1   580       83  -1.35228E+01   2.43527E+00   6.53776E-02  -2.43439E+00
   1   590       84  -1.40227E+01   2.44145E+00   6.52120E-02  -2.44058E+00
   1   600       85  -1.45227E+01   2.44758E+00   6.50483E-02  -2.44672E+00
   1   610       86  -1.50227E+01   2.45368E+00   6.48865E-02  -2.45282E+00
   1   620       87  -1.55226E+01   2.45973E+00   6.47267E-02  -2.45888E+00
   1   630       88  -1.60226E+01   2.46574E+00   6.45687E-02  -2.46489E+00
   1   640       89  -1.65226E+01   2.47171E+00   6.44125E-02  -2.47087E+00
   1   650       90  -1.70225E+01   2.47764E+00   6.42580E-02  -2.47681E+00
   1   660       91  -1.75225E+01   2.48354E+00   6.41053E-02  -2.48271E+00
   1   670       92  -1.80225E+01   2.48939E+00   6.39544E-02  -2.48857E+00
   1   680       93  -1.85224E+01   2.49521E+00   6.38050E-02  -2.49439E+00
   1   690       94  -1.90224E+01   2.50099E+00   6.36573E-02  -2.50018E+00
   1   700       95  -1.95224E+01   2.50674E+00   6.35113E-02  -2.50593E+00
   1   710  EP   96  -2.00223E+01   2.51245E+00   6.33667E-02  -2.51165E+00

 Total Time    0.207E-01
qif ... done
Merge done

Key parameters explained:

  • origin=t_cont: Start from time continuation branch
  • starting_point='UZ1': Start from first user marker
  • bidirectional=True: Continue in both positive and negative directions
  • ICP=4: Parameter index for \(\bar\eta\) (4th parameter in Fortran file)
  • RL0=-20.0, RL1=20.0: Parameter bounds
  • IPS=1: Indicates equilibrium continuation of an ODE system
  • ILP=1: Enable fold bifurcation detection
  • ISP=2: Full automatic bifurcation detection

The output shows LP in column TY for some solutions, indicating limit point (fold) bifurcations.

Part 3: Bifurcation Diagram

We can visualize the full bifurcation diagram:

# Plot bifurcation diagram
fig, axes = plt.subplots(1, 2, layout="compressed", figsize=(10, 5))

# Firing rate vs eta
ax = axes[0]
qif_auto.plot_continuation("eta", "r", cont="eta", ax=ax)
ax.set_xlabel("η (Background drive)")
ax.set_ylabel("Firing rate r")
ax.set_title("QIF Bifurcation Diagram: Firing Rate")

# Membrane potential vs eta
ax = axes[1]
qif_auto.plot_continuation("eta", "v", cont="eta", ax=ax)
ax.set_xlabel("η (Background drive)")
ax.set_ylabel("Membrane potential v")
ax.set_title("QIF Bifurcation Diagram: Voltage")
plt.show()

Interpreting the Diagram

The curve represents the value of \(r\) (y-axis) at equilibrium solutions for each value of \(\bar\eta\) (x-axis):

  • Solid line: Stable equilibrium
  • Dotted line: Unstable equilibrium
  • Triangles (LP): Fold (limit point) bifurcations

At a fold bifurcation, the critical eigenvalue of the vector field crosses the imaginary axis (its real part changes sign), indicating a change of stability. A stable and unstable equilibrium approach and annihilate each other.

import shutil

# Clean up temporary files
qif.clear()
shutil.rmtree('_continuation', ignore_errors=True)

# Remove generated files
for f in ['qif.f90', 'c.ivp']:
    if os.path.exists(f):
        os.remove(f)

print("Cleaned up temporary files")
Cleaned up temporary files

References