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.