example in docs using @reaction_network

[gszep]

Would be good to show users that its possible to use this library with DiffEqBiological

using DiffEqBiological,BifurcationKit
using Setfield,LinearAlgebra,Plots

system = @reaction_network begin
     k₁,        Y → 2X
     2k₂,       2X → X + Y
     k₃,        X + Y → Y
     k₄,        X → 0
end k₁ k₂ k₃ k₄

# get rate function and jacobian, for use in continuation
F! = oderhsfun(system)
function F(u::Vector{T},p::NamedTuple) where T<:Number
    f = similar(u)
    F!(f,u,[ p[key] for (key, value) ∈ paramsmap(system) ],Inf)
    return f
end

∂Fu! = jacfun(system)
function ∂Fu(u::Vector{T},p::NamedTuple) where T<:Number
    J = zeros(T,length(u),length(u))
    ∂Fu!(J,u,[ p[key] for (key, value) ∈ paramsmap(system) ],Inf)
    return J
end

# define the parameter direction that you want to explore
pRange = range(0,12,length=200)
parameter_name = (@lens _.k₁)
p = (k₁=0.0, k₂=1.0, k₃=1.0, k₄=1.5)

# options for branch continuation along parameter direction
options = ContinuationPar( maxSteps = 5000,
    dsmax=10*step(pRange), dsmin=step(pRange), ds=step(pRange),
    pMin = minimum(pRange), pMax=maximum(pRange),
    newtonOptions = NewtonPar(tol = 1e-8),
    saveEigenvectors = false
)

# options for detecting multiple branches
deflation = DeflationOperator(2.0, dot, 1.0, [ zeros(length(species(system))) ])
termination_condition(x, f, J, residual, iteration, niter, options; kwargs...) = residual <1e3
root_search(x,p,id) = x .+ rand(length(x))

##################################### main continuation method
@assert( det(∂Fu(randn(length(species(system))),p)) ≠ 0, "species must be linearly independent")
branches, = continuation( F, ∂Fu, p, parameter_name, options, deflation; showplot=false,
    callbackN = termination_condition, perturbSolution = root_search )

#################################### display results
plot(size=(500,500),ylabel="steady states", xlabel="parameter $(first(typeof(parameter_name).parameters))")
plot!(branches,lw=3,label="")

[rveltz]

Nice! I can add bifurcation detection to this.

@ChrisRackauckas would you prefer to see this example?

[gszep]

note that if the determinant of the jacobian is zero det(∂Fu(u,p)) = 0 the continuation method encounters a singularity exception. This becomes a problem any time the user wants to add conservation laws such as

@add_constraints system begin
           X + Y = 1
end

or in cases where the reactive species are not linearly independent, or more generally where there are local conservation laws at play, i.e. when the vector field F(u,p) is solonoidal - has zero divergence - in the vicinity of F(u,p)=0

It seems that jacfun doesn't return the correct jacobian that takes into account the use of @add_constraints. We would expect in the example above that one of the columns of the jacobian is zero due to the linear dependence between X and Y

[rveltz]

Do you have such example where the constraint induces singularity? For deflated continuation, this should not a problem (too much)

[rveltz]

This

# options for branch continuation along parameter direction
options = ContinuationPar( maxSteps = 5000,
    dsmax=10*step(pRange), dsmin=step(pRange), ds=step(pRange),
    pMin = minimum(pRange), pMax=maximum(pRange),
    newtonOptions = NewtonPar(tol = 1e-8),
    saveEigenvectors = false, detectBifurcation = 3
)

##################################### main continuation method
@assert( det(∂Fu(randn(length(species(system))),p)) ≠ 0, "species must be linearly independent")
branches, = continuation( F, ∂Fu, p, parameter_name, options, deflation; showplot=false,
    callbackN = termination_condition, perturbSolution = root_search )

#################################### display results
plot(size=(500,500),ylabel="steady states", xlabel="parameter $(first(typeof(parameter_name).parameters))")
plot!(branches...,label="")

gives the stability:

[isaacsas]

This looks great! It would probably be better to build the example using Catalyst instead (which replaced DiffEqBiological). It generates a ModelingToolkit.ODESystem from which one can still get Jacobian and rhs functions (and probably directly get out of place versions if that is needed).

[isaacsas]

If it would be helpful I can try to adapt this to Catalyst and add here.

[rveltz]

Yes, please do! Also, I was hoping to get an example where the jacobian is singular but maybe that's not possible

[gszep]

A simple example of a singular jacobian would be adding any reaction that includes "auxiliary" species P that does not change the dynamics, but is merely proportional to it.

system = @reaction_network begin
     k₁,        Y → 2X
     2k₂,       2X → X + Y
     k₃,        X + Y → Y
     k₄,        X → 0
     1,         X → P
end k₁ k₂ k₃ k₄

[rveltz]

👍

[TorkelE]

Starting this up again.

While the new version of DiffEqBiological.jl looks similar to the previous one on the surface, the internals has been redone more or less from scratch (we also took the time to rename it Catalyst.jl). While it previously was more or less self-contained, it is now based on ModelingToolkit.jl (MTK), https://github.com/SciML/ModelingToolkit.jl. The @reaction_network macro, while taking the same kind of input as before, now generates a new type of structure, a ReactionSystem (this structure is implemented in MTK). MTK also have many similar structures implemented, like ODESystem, SDESystem, JumpSystem, NonlinearSystem, etc. It also contains the means to transform a ReactionSystem to some desired system (this is done automatically when e.g. using a ReactionSystem as input to solve an ODE).

The advantage of this is that, while previously, we had to construct methods to act specifically on the quite specific systems generated by DiffEqBiological.jl, now we can make those act on the generalised systems produced by MTK. Now several different modelling packages can be implemented via MTK, and they will all automatically get access to loads of features. At the same time, lots of the stuff we created for DiffEqBiological.jl (in addition to some new stuff) can now be used by other modelling packages. This means that it will be a lot easier to implement new packages like Catalyst.jl. The people at SciML are also preparing lots of cool stuff MTK.

So what does this mean? For a starter, some things look a little bit different now. The system function and Jacobian are no longer stored in the ReactionSystem. Instead, we must generate an ODESystem, from this we can generate an ODEFunction which contains these (also, this ODEFunction, as opposed to your example, takes a third argument, t, the time, which must be accounted for). The example thus becomes

using Catalyst,BifurcationKit
using Setfield,LinearAlgebra,Plots

system = @reaction_network begin
     k₁,        Y → 2X
     2k₂,       2X → X + Y
     k₃,        X + Y → Y
     k₄,        X → 0
end k₁ k₂ k₃ k₄

# get rate function and jacobian, for use in continuation
odefun = ODEFunction(convert(ODESystem,system),jac=true)
function F(u::Vector{T},p::NamedTuple) where T<:Number
    odefun(u,p,0.)
end
function ∂Fu(u::Vector{T},p::NamedTuple) where T<:Number
    odefun.jac(u,p,0.)
end

# define the parameter direction that you want to explore
pRange = range(0,12,length=200)
parameter_name = (@lens _.k₁)
p = (k₁=0.0, k₂=1.0, k₃=1.0, k₄=1.5)

# options for branch continuation along parameter direction
options = ContinuationPar( maxSteps = 5000,
    dsmax=10*step(pRange), dsmin=step(pRange), ds=step(pRange),
    pMin = minimum(pRange), pMax=maximum(pRange),
    newtonOptions = NewtonPar(tol = 1e-8),
    saveEigenvectors = false
)

# options for detecting multiple branches
deflation = DeflationOperator(2.0, dot, 1.0, [ zeros(length(species(system))) ])
termination_condition(x, f, J, residual, iteration, niter, options; kwargs...) = residual <1e3
root_search(x,p,id) = x .+ rand(length(x))

##################################### main continuation method
@assert( det(∂Fu(randn(length(species(system))),p)) ≠ 0, "species must be linearly independent")
branches, = continuation( F, ∂Fu, p, parameter_name, options, deflation; showplot=false,
    callbackN = termination_condition, perturbSolution = root_search )

#################################### display results
plot(size=(500,500),ylabel="steady states", xlabel="parameter $(first(typeof(parameter_name).parameters))")
plot!(branches,lw=3,label="")

Secondly, it would be cool to implement so that BifurcationKit.jl could work directly on some structure generated by MTK. That way, anything which uses MTK would get a shortcut into using BifurcationKit.jl. I am not sure what would be most convenient, either the ODESystem which I have shown here, or maybe the NonlinearSystem, it would depend on what was most convenient. Possibly, one could also create a BifrucationSystem, and then some way of generating those from the above-mentioned system. I'm happy to help, as well as do some implementation in MTK.

[rveltz]

Thank you a lot for the detailed explanations, it is very useful!!

I already have DiffEqBase used in BK which export ODEFunction and NonlinearSystem. So, if I use your example, it would be very easy to generate the functions from a system. But, adding Catalyst as a dependency for BK would probably be too big.

[isaacsas]

@TorkelE Thanks for putting the example together! Unfortunately I just didn't have the time this month so far...

[isaacsas]

@rveltz There's no reason to add Catalyst as dependency here. I think if we just add a tutorial in Catalyst showing the code above that should be good. That would give users the idea of how to interface with BifurcationKit.

[TorkelE]

Yeah, makes little sense to add Catalyst as a dependency (also an advantage of it being built of MTK and general structures). It is possible to make a NonlinearSystem from a ReactionSystem, might redo the example to fo down that path instead, might be the most natural opne.

[isaacsas]

Yeah, though we need some tests for NonlinearSystems in MT before I’d rely on them; we haven’t verified the conversion is working well since the MT updates.

[TorkelE]

Ahh, I see. Maybe keep it to ODEFunction then. But aiming to have NonlinearSystem as the long term medium for these kinds of things might make sense?

[rveltz]

yes, why not

[TorkelE]

Hi again.

I was playing around with this for one of my systems, and was comparing your code here with a similar script I threw together for myself during the summer. I was trying to interpret the results, but they look... odd.

First I declare my model:

using Catalyst,BifurcationKit
using Setfield,LinearAlgebra,Plots
using OrdinaryDiffEq

system = @reaction_network begin
    v0+v*(S*σ)^n/((S*σ)^n+(D*A)^n+K^n), ∅ → σ
    d, σ → ∅
    (τ*σ,τ), ∅ ↔ A
end S D τ v v0 K n d;
parameters = [10.,9.,0.005,2.,0.01,20.,3,0.05]
pRange = (0.1,20.)
p = :S;

Next, I use my previous script to generate a bifurcation diagram:

p_init = setindex!(copy(parameters),pRange[1],findfirst(Symbol.(system.ps) .== p))
fun = ODEFunction(convert(ODESystem,system),jac=true)
x0 = solve(ODEProblem(system,[1.,1.],(0.,2000),p_init),Rosenbrock23()).u[end];
options = ContinuationPar( maxSteps=50000, 
    dsmax=(pRange[2]-pRange[1])*1e-3, dsmin=(pRange[2]-pRange[1])*1e-9, ds=(pRange[2]-pRange[1])*1e-9, 
    pMin=pRange[1], pMax=pRange[2], detectBifurcation=1, 
    newtonOptions = NewtonPar(tol = 1e-8, verbose = false));
path = continuation((u,p)->fun(u,p,0), (u,p)->fun.jac(u,p,0), x0, p_init, (@lens _[findfirst(Symbol.(system.ps) .== p)]), options; tangentAlgo=BorderedPred());

plot(path[1],grid=false,xguide="$p",plotfold=false,label="",framestyle=:box)

which generates: which is as it should be (a bistable switch where the unstable region continuous a bit).

Now I try the code used here:

p_init = (S = 0.001, D = 9., τ = 0.005, v = 2., v0 = 0.01, K = 20., n = 3, d = 0.05);
fun = ODEFunction(convert(ODESystem,system),jac=true)
pRange = range(0,20,length=200)
parameter_name = (@lens _.S);

options = ContinuationPar( maxSteps=100000, 
    dsmax=10*step(pRange), dsmin=step(pRange), ds=step(pRange),
    pMin=pRange[1], pMax=pRange[end], 
    newtonOptions = NewtonPar(tol = 1e-8, verbose = false));
deflation = DeflationOperator(2.0, dot, 1.0, [ zeros(length(species(system))) ])
termination_condition(x, f, J, residual, iteration, niter, options; kwargs...) = residual <1e3
root_search(x,p,id) = x .+ rand(length(x));
@assert(det(fun.jac(randn(length(species(system))),p_init,0.)) ≠ 0, "species must be linearly independent")
branches, = continuation((u,p)->fun(u,p,0), (u,p)->fun.jac(u,p,0), p_init, parameter_name, options, deflation; showplot=false,callbackN = termination_condition, verbosity = 0, perturbSolution = root_search);

plot(branches,lw=3,label="")

which yields this (a bit confusing) diagram: If I reduce the dsmin it gets a bit better, but there's still another weird, extra, line which I do not understand:

options = ContinuationPar( maxSteps=150000, 
    dsmax=10*step(pRange), dsmin=step(pRange)/1000., ds=step(pRange)/1000.,
    pMin=pRange[1], pMax=pRange[end], 
    newtonOptions = NewtonPar(tol = 1e-8, verbose = false));

Also, the different colours, are they different paths, or do they signify some other property?

I have a feeling there's something here which I am missing, but not really sure what.

[rveltz]

Hi,

I think there is nothing wrong. The lines blue-brown-magenta correspond to the diagram you first computed (with 3 folds) with continuation. Now, deflated continuation found another branch which is not connected to the main branch (the blue one you first computed), this is cool!!

Note that there is an elaboration of the deflated continuation algorithm that I did not push. In a second pass, it connects together the branches that should be. So it would return two branches and not 4 as now.

[rveltz]

By the way, using the options:

options = ContinuationPar( maxSteps=50000, 
           dsmax=(pRange[2]-pRange[1])*1e-3, dsmin=(pRange[2]-pRange[1])*1e-9, ds=(pRange[2]-pRange[1])*1e-9, 
           pMin=pRange[1], pMax=pRange[2], detectBifurcation=3,nInversion=4, 
           newtonOptions = NewtonPar(tol = 1e-8, verbose = false));

it gives bifurcation points, in particular a Hopf one:

julia> path
(Branch number of points: 1941
Branch of Equilibrium
Parameters from 0.1 to 20.0
Bifurcation points:
 (ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`)
- #  1,    bp at p ≈ +9.06308826 ∈ (+9.06308826, +9.06308848), |δp|=2e-07, [converged], δ = ( 1,  0), step = 358, eigenelements in eig[359], ind_ev =   1
- #  2,    bp at p ≈ +4.22810078 ∈ (+4.22810078, +4.22810078), |δp|=4e-10, [converged], δ = ( 1,  0), step = 566, eigenelements in eig[567], ind_ev =   2
- #  3,  hopf at p ≈ +10.75597251 ∈ (+10.75597184, +10.75597251), |δp|=7e-07, [converged], δ = (-2, -2), step = 1408, eigenelements in eig[1409], ind_ev =   2
Fold points:
- #  1, fold at p ≈ 9.06308826, [    guess], step = 355, eigenelements in eig[355], ind_ev =   0
- #  2, fold at p ≈ 4.22810078, [    guess], step = 563, eigenelements in eig[563], ind_ev =   0
, BorderedArray{Array{Float64,1},Float64}([36.858740221512086, 36.858740221512086], 20.0), BorderedArray{Array{Float64,1},Float64}([0.6414478443904162, 0.6414478443904162], 1.396422020836672))

[TorkelE]

I think there is nothing wrong. The lines blue-brown-magenta correspond to the diagram you first computed (with 3 folds) with continuation. Now, deflated continuation found another branch which is not connected to the main branch (the blue one you first computed), this is cool!!

This might be the moment when I reveal my ignorance, but what exactly is the green branch? Unless I am very mistaken, the blue-brown-magenta branch(es) is the only equilibrium solutions to the system. Unless the green is (a possible projected) branch with negative/imaginary values, I am not really sure what it is?

I try extracting the values by setting saveSolEveryStep=true in ContinuationPar(), and then plotting the green branch, as well as the extracted values (or their negation, as they are negative).

plot(branches[2],lw=3,label="That branch",color=:green)
y = getproperty.(branches[2].sol,:x);
x = getproperty.(branches[2].sol,:p);
plot!(x,-first.(y),lw=3,color=:red,label="Absolute value of the first sol value")
plot!(x,-last.(y),lw=3,color=:blue,linestyle=:dash,label="Absolute value of the first sol value")

This seems related to some negative solution, but not exactly?

(The blue/red line seems to be vaguely a solution. If I select some value of my parameter and insert the corresponding value of (sigma,A), as picked from the blue/red line, I get values like [ 4.222571785623963e-6, 0.0], which is not zero, but at least small).

[rveltz]

Should we close this?

[TorkelE]

Yes, let's do.