Example: Parameter continuation and bifurcation detection

In this tutorial, you will learn how to perform a 1D numerical parameter continuation in PyCoBi with automatic fold bifurcation detection. Furthermore, you will learn how to plot a simple bifurcation diagram. Throughout this example, we will use the quadratic integrate-and-fire population model [1], a detailed introduction of which is given in the model introductions example gallery. The dynamic equations of this model read the following:

\[\begin{split}\tau \dot r &= \frac{\Delta}{\pi\tau} + 2 r v, \\ \tau \dot v &= v^2 +\bar \eta + I(t) + J r \tau - (\pi r \tau)^2,\end{split}\]

where \(r\) is the average firing rate and \(v\) is the average membrane potential of the QIF population. It is governed by 4 parameters:

• \(\tau\) –> the population time constant

• \(\bar \eta\) –> the mean of a Lorentzian distribution over the neural excitability in the population

• \(\Delta\) –> the half-width at half maximum of the Lorentzian distribution over the neural excitability

• \(J\) –> the strength of the recurrent coupling inside the population

In this tutorial, we will demonstrate how to (1) perform a simple 1D parameter continuation in \(\bar \eta\), and (2) plot the corresponding bifurcation diagram. The latter has also been done in [1], so you can compare the resulting plot with the results reported by Montbrió et al.

References

[1] (1,2)

E. Montbrió, D. Pazó, A. Roxin (2015) Macroscopic description for networks of spiking neurons. Physical Review X, 5:021028, https://doi.org/10.1103/PhysRevX.5.021028.

[2]

E.J. Doedel, T.F. Fairgrieve, B. Sandstede, A.R. Champneys, Y.A. Kuznetsov and W. Xianjun (2007) Auto-07p: Continuation and bifurcation software for ordinary differential equations. Technical report, Department of Computer Science, Concordia University, Montreal, Quebec.

Part 1: Creating an ODESystem Instance

In this first part, we will be concerned with how to create a model representation that is compatible with auto-07p, which is the software that is used for parameter continuations and bifurcation analysis in PyCoBi [2]. To this end, we make use of the code generation software PyRates, which allows us to generate the Fortran files required to run auto-07p from YAML model definition files. The QIF mean-field model comes pre-implemented with PyRates, so we can simply point to the ready-to-use YAML definition file.

from pycobi import ODESystem

# path to YAML model definition
model = "model_templates.neural_mass_models.qif.qif"

# installation directory of auto-07p
auto_dir = "~/projects/auto-07p"

# ODESystem initialization: pass `init_cont=True` to run an initial
# time-integration ("IVP") at instantiation, which is what most
# equilibrium-continuation workflows want — see note below.
ode = ODESystem.from_yaml(model, auto_dir=auto_dir, init_cont=True)

If you haven’t manually configured your system environment variables such that the Python installation of auto-07p is recognized by your Python interpreter, you can simply add the installation directory of auto-07p to the ODESystem instantiation (as above) and it will take care of these environment variables for you. With init_cont=True, an integration of the QIF model over time is automatically performed during initialization to ensure that the system converged to a steady-state. This is required for equilibrium-style parameter continuations: the starting point must be an equilibrium solution, i.e. \(\dot r = 0\) and \(\dot v = 0\) in our example. Since PyCoBi >= 0.10.0 the IVP is opt-in — if you already have a converged state (e.g. from a hand-written c.eq file or a previous run loaded via ODESystem.from_file), leave init_cont=False (the default) and skip the IVP.

You can check whether the system indeed converged to a steady-state solution by plotting the results of the time integration:

ode.plot_continuation("t", "p/qif_op/r", cont=0)
plt.show()

The above code plots the first state variable \(r\) against time. The column name "p/qif_op/r" is the namespaced PyRates identifier — PyRates >= 1.1 also writes auto-07p unames / parnames into the generated c.* files so the column is additionally available under its bare local name "r". The auto-07p-native "U(1)" form continues to work too. Alternatively, you can look at the auto-07p output in the terminal, where the values of U(1) and U(2) converge to specific numbers — these are the steady-state values of \(r\) and \(v\) the QIF system converged to.

Part 2: Performing Parameter Continuations

In this part, we will demonstrate how to perform simple 1D parameter continuations via the ODESystem.run() method. Since our model converged to an equilibrium and we are now save to perform the continuation in our parameter of interest: \(\bar \eta\). We can do this via the ODESystem.run method:

eta_sols, eta_cont = ode.run(
    origin=0, starting_point='EP2', name='eta', bidirectional=True,
    ICP=4, RL0=-20.0, RL1=20.0, IPS=1, ILP=1, ISP=2, 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, EPSL=1e-06, EPSU=1e-06, EPSS=1e-04,
    DS=1e-4, DSMIN=1e-8, DSMAX=5e-2, IADS=1, THL={}, THU={}, UZR={}, STOP={}
)

In this call, we specified the full set of auto-07p constants. Don’t worry, usually, you do not have to bother with most of them. It is common practice to specify most of them in constants files that you would pass to ODESystem.run via the keyword argument c=name. In such a case, you would specify all auto-07p constants that do not change between calls to the .run() method in a file with the name c.name and only provide the constants that need to be altered between .run() calls directly to the .run() method.

Note that we use starting_point='EP2' rather than 'EP1' — the IVP that ran when we instantiated the ODESystem produced two labeled solutions: 'EP1' is the \(t = 0\) initial condition, while 'EP2' is the converged steady state. Equilibrium continuations want the latter.

We did not pass a JAC=... value here because from_yaml defaults to analytical_jacobian=True, which instructs PyRates to symbolically differentiate the vector field and write a DFDU/DFDP block into the generated Fortran. The generated c.* file then sets JAC=1 so auto-07p picks up the analytical Jacobian — usually faster and more accurate than the default finite-difference fallback. Pass analytical_jacobian=False at instantiation, or JAC=0 on the per-call run(), to opt out.

Checking the terminal output of auto-07p, you will realize that the output in column TY shows LP for two of the solutions we computed along our branch in \(\bar\eta\). These indicate the detection of limit point or fold bifurcations. We can visualize the full bifurcation diagram via the following call:

ode.plot_continuation('PAR(4)', 'U(1)', cont='eta')
plt.show()

The curve in this plot represents the value of \(r\) (y-axis) at the equilibrium solutions that exist for each value of \(\bar \eta\) (x-axis). A solid line indicates that the equilibrium is stable, whereas a dotted line indicated that the equilibrium is unstable. The triangles mark the points at which auto-07p detected fold bifurcations. At a fold bifurcation, the critical eigenvalue of the vector field defined by the right-hand sides of our model’s ODEs crosses the imaginary axis (i.e. its real part changes the sign). This indicates a change of stability of an equilibrium solution, which happens because an unstable and a stable equilibrium approach and annihilate each other. This behavior can be read from the plot as well. The solid and the dotted line approach each other towards the fold bifurcation marks. After they meet at the fold bifurcation, both cease to exist.