Bifurcation Diagrams

Bifurcation diagrams can be produced from Catalyst generated models through the use of the BifurcationKit.jl package. This tutorial gives a simple example of how to create such a bifurcation diagram.

First, we declare our reaction model. For this example we will use a bistable switch, but one which also contains a Hopf bifurcation.

using Catalyst
rn = @reaction_network begin
    (v0 + v*(S * X)^n / ((S*X)^n + (D*A)^n + K^n), d), ∅ ↔ X
    (X/τ, 1/τ), ∅ ↔ A
end S D τ v0 v K n d

\[ \begin{align*} \require{mhchem} \ce{ \varnothing &<=>[$v0 + \frac{v \left( S X \right)^{n}}{K^{n} + \left( A D \right)^{n} + \left( S X \right)^{n}}$][$d$] X}\\ \ce{ \varnothing &<=>[$\frac{X}{\tau}$][$\frac{1}{\tau}$] A} \end{align*} \]

Next, we specify the system parameters for which we wish to plot the bifurcation diagram. We also set the parameter we wish to vary in our bifurcation diagram, as well as the interval to vary it over. Finally, we set which variable we wish to plot the steady state values of in the bifurcation plot.

p = Dict(:S => 1., :D => 9., :τ => 1000., :v0 => 0.01,
         :v => 2., :K => 20., :n => 3, :d => 0.05)
bif_par = :S           # bifurcation parameter
p_span = (0.1, 20.)    # interval to vary S over
plot_var = :X          # we will plot X vs S
:X

When creating a bifurcation diagram, we typically start at some point in parameter phase-space. We will simply select the beginning of the interval over which we wish to compute the bifurcation diagram, p_span[1]. We thus create a modified parameter set where S = .1. For this parameter set, we also make a guess for the steady-state of the system. While a good estimate could be provided through an ODE simulation, BifurcationKit does not require the guess to be very accurate.

p_bstart = copy(p)
p_bstart[bif_par] = p_span[1]
u0 = [:X => 1.0, :A => 1.0]
2-element Vector{Pair{Symbol, Float64}}:
 :X => 1.0
 :A => 1.0

Finally, we extract the ODE derivative function and its jacobian in a form that BifurcationKit can use:

oprob = ODEProblem(rn, u0, (0.0,0.0), p_bstart; jac = true)
F = (u,p) -> oprob.f(u, p, 0)
J = (u,p) -> oprob.f.jac(u, p, 0)
#3 (generic function with 1 method)

In creating an ODEProblem an ordering is chosen for the initial condition and parameters, and regular Float64 vectors of their numerical values are created as oprob.u0 and oprob.p respectively. BifurcationKit needs to know the index in oprob.p of our bifurcation parameter, :S, and the index in oprob.u0 of the variable we wish to plot, :X. We calculate these as

# get S and X as a symbolic variables
@unpack S, X = rn

# find their indices in oprob.p and oprob.u0 respectively
bif_idx  = findfirst(isequal(S), parameters(rn))
plot_idx = findfirst(isequal(X), species(rn))
1

Now, we load the required packages to create and plot the bifurcation diagram. We also bundle the information we have compiled so far into a BifurcationProblem.

using BifurcationKit, Plots, LinearAlgebra, Setfield

bprob = BifurcationProblem(F, oprob.u0, oprob.p, (@lens _[bif_idx]);
                           recordFromSolution = (x, p) -> x[plot_idx], J = J)
┌─ Bifurcation Problem with uType Vector{Float64}
├─ inplace:  false
├─ Symmetric: false
└─ Parameter: p

Next, we need to specify the input options for the pseudo-arclength continuation method (PACM) which produces the diagram.

bopts = ContinuationPar(dsmax = 0.05,          # Max arclength in PACM.
                        dsmin = 1e-4,          # Min arclength in PACM.
                        ds=0.001,              # Initial (positive) arclength in PACM.
                        maxSteps = 100000,     # Max number of steps.
                        pMin = p_span[1],      # Min p-val (if hit, the method stops).
                        pMax = p_span[2],      # Max p-val (if hit, the method stops).
                        detectBifurcation = 3) # Value in {0,1,2,3}
BifurcationKit.ContinuationPar{Float64, BifurcationKit.DefaultLS, BifurcationKit.DefaultEig{typeof(real)}}
  dsmin: Float64 0.0001
  dsmax: Float64 0.05
  ds: Float64 0.001
  θ: Float64 0.5
  a: Float64 0.5
  pMin: Float64 0.1
  pMax: Float64 20.0
  maxSteps: Int64 100000
  finDiffEps: Float64 1.0e-9
  newtonOptions: BifurcationKit.NewtonPar{Float64, BifurcationKit.DefaultLS, BifurcationKit.DefaultEig{typeof(real)}}
  η: Float64 150.0
  saveToFile: Bool false
  saveSolEveryStep: Int64 0
  nev: Int64 3
  saveEigEveryStep: Int64 1
  saveEigenvectors: Bool true
  plotEveryStep: Int64 10
  tolStability: Float64 1.0e-10
  detectFold: Bool true
  detectBifurcation: Int64 3
  dsminBisection: Float64 1.0e-16
  nInversion: Int64 2
  maxBisectionSteps: Int64 15
  tolBisectionEigenvalue: Float64 1.0e-16
  detectEvent: Int64 0
  tolParamBisectionEvent: Float64 1.0e-16
  detectLoop: Bool false

Here detectBifurcation determines to what extent bifurcation points are detected and how accurately their values are determined. Three indicates to use the most accurate method for calculating their values.

We are now ready to compute the bifurcation diagram:

bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> bopts)
[Bifurcation diagram]
 ┌─ From 0-th bifurcation point.
 ├─ Children number: 0
 └─ Root (recursion level 1)
      ┌─ Number of points: 778
      ├─ Curve of EquilibriumCont
      ├─ Type of vectors: Vector{Float64}
      ├─ Parameter p starts at 0.1, ends at 20.0
      ├─ Algo: PALC
      └─ Special points:

If `br` is the name of the branch,
ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`

- #  1,       bp at p ≈ +9.06308828 ∈ (+9.06308828, +9.06308828), |δp|=2e-09, [converged], δ = ( 1,  0), step = 138, eigenelements in eig[139], ind_ev =   1
- #  2,       bp at p ≈ +4.22810281 ∈ (+4.22810281, +4.22810821), |δp|=5e-06, [converged], δ = ( 1,  0), step = 226, eigenelements in eig[227], ind_ev =   2
- #  3,     hopf at p ≈ +11.18866014 ∈ (+11.18528123, +11.18866014), |δp|=3e-03, [converged], δ = (-2, -2), step = 579, eigenelements in eig[580], ind_ev =   2
- #  4, endpoint at p ≈ +20.00000000,                                                                     step = 777

Finally, we can plot it:

plot(bf, xlabel = string(bif_par), ylabel = string(plot_var))

Here, the Hopf bifurcation is marked with a red dot and the fold bifurcations with blue dots. The region with a thinner line width corresponds to unstable steady states.

This tutorial demonstrated how to make a simple bifurcation diagram where all branches are connected. However, BifurcationKit.jl is a very powerful package capable of a lot more. For more details, please see that package's documentation: BifurcationKit.jl.