BifKit extension, conservation laws, and changing initial conditions

[isaacsas]

We need to check if it is actually working when one is using an initial condition as the bifurcation parameter. I suspect that it has no way to know it needs to update the conserved constant in this case (i.e. it is just changing that parameter in the underlying parameter vector), and hence will give incorrect results.

[isaacsas]

This doesn't work, but at least it errors:

rn = @reaction_network begin
              @parameters M
              @species A(t) = M
              (1.0,1.0), A <--> B
              end

osys = complete(convert(ODESystem, rn; remove_conserved = true))
u_guess = [osys.B => 2.0]
p_start = [osys.M => 3.0]
bprob = BifurcationProblem(osys, u_guess, p_start, osys.M; plot_var = osys.A)

giving

ERROR: ArgumentError: SymbolicUtils.BasicSymbolic{Real}[M] are either missing from the variable map or missing from the system's unknowns/parameters list.
Stacktrace:
 [1] throw_missingvars_in_sys(vars::Vector{SymbolicUtils.BasicSymbolic{Real}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/XHRvt/src/utils.jl:765
 [2] promote_to_concrete(vs::Vector{SymbolicUtils.BasicSymbolic{Real}}; tofloat::Bool, use_union::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/XHRvt/src/utils.jl:784
 [3] varmap_to_vars(varmap::Vector{…}, varlist::Vector{…}; defaults::Dict{…}, check::Bool, toterm::Function, promotetoconcrete::Nothing, tofloat::Bool, use_union::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/XHRvt/src/variables.jl:171
 [4] BifurcationProblem(::NonlinearSystem, ::Vector{…}, ::Vector{…}, ::Num; plot_var::Num, record_from_solution::Function, jac::Bool, kwargs::@Kwargs{})
   @ MTKBifurcationKitExt ~/.julia/packages/ModelingToolkit/XHRvt/ext/MTKBifurcationKitExt.jl:104
 [5] BifurcationProblem(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any}; kwargs::@Kwargs{plot_var::Num})
   @ MTKBifurcationKitExt ~/.julia/packages/ModelingToolkit/XHRvt/ext/MTKBifurcationKitExt.jl:152
 [6] top-level scope
   @ REPL[35]:1
Some type information was truncated. Use `show(err)` to see complete types.

[isaacsas]

Continuing with respect to Γ[1], the conservation constant, seems to (silently) give wrong answers:

using Catalyst, BifurcationKit, Plots
rn = @reaction_network begin
              (1.0,1.0), A <--> B
              end
osys = structural_simplify(convert(ODESystem, rn; remove_conserved = true))
u_guess = [osys.A => 1.0, osys.B => 2.0]
p_start = [osys.Γ[1] => 3.0]  # A = 1.0 then
bprob = BifurcationProblem(osys, u_guess, p_start, osys.Γ[1]; plot_var = osys.A)
p_span = (3.0, 20.) # A goes from 1.0 to 18.0
opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps = 1000)
bif_dia = bifurcationdiagram(bprob, PALC(), 2, opts_br; bothside = true)
plot(bif_dia; xlabel = "A(0) + B(0)", label = "A(t)")

which gives:

Actually this is right and I'm just being dumb. So that is good.

[isaacsas]

Seems like we need to disallow changing initial conditions and/or the conservation constant in the current extension setup.

[isaacsas]

@rveltz is there any way you could use remake for updating parameters and initial conditions in the bifurcation solvers (at least when using problems coming via MTK)? I guess right now you just modify the bifurcation parameter, but in cases like this one needs to hook into the SciML remake functionality for everything to update correctly. If you have an initial condition given by a symbolic parameter you want to make the bifurcation parameter, or a conservation law over which you want to vary the conserved constant, you need to have any other dependent parameters/ICs also update accordingly.

[rveltz]

Do you want the same name because there is BK.re_make

I was worried about type piracy so that s why the new name. I can name it as in SciMLBase if needed

[isaacsas]

Does that hook into the symbolic SciMLBase.remake?

[rveltz]

no, it does not extend it

[isaacsas]

Yeah, so I suspect that if one has parameters that are functions of initial conditions, initial conditions that are functions of parameters, or parameters that are functions of other parameters, one may be likely to get incorrect results (because these may need to get updated as the bifurcation parameter changes).

[rveltz]

If this https://github.com/SciML/SciMLBase.jl/blob/0f8ec1f501e916e5f8c39d94ddc02757d574cc9c/src/remake.jl#L27 is enough, it wont take me long.

But if you want this one

https://github.com/SciML/SciMLBase.jl/blob/0f8ec1f501e916e5f8c39d94ddc02757d574cc9c/src/remake.jl#L80

then I have to think

[isaacsas]

Maybe this would be easier by having a reeval type function that a user can pass in (in this case Catalyst/MTK)? That could handle updating all other parameters when one changes during a solve? Then we could perhaps pass a wrapped version of SciMLBase.remake? @AayushSabharwal is the expert on this functionality in SciMLBase, but it would be nice to get this working and updating correctly.

[isaacsas]

I guess trying to worry about changing a parameter in an initial condition is a bit much, but I'll try to work up a simple example where there are dependent parameters but things aren't updating correctly when the bifurcation parameter changes.

[rveltz]

OK!

[vyudu]

    rn = @reaction_network begin
        @parameters k ksq = k^2
        (k, ksq), A <--> B
    end

    rn = complete(rn)
    u0_guess = [:A => 1., :B => 1.] 
    p_start = [:k => 2.]

    bprob = BifurcationProblem(rn, u0_guess, p_start, :k; plot_var = :A, u0 = [:A => 5., :B => 3.])
    p_span = (0.1, 6.0)
    opts_br = ContinuationPar(dsmin = 0.0001, dsmax = 0.001, ds = 0.0001, max_steps = 10000, p_min = p_span[1], p_max = p_span[2], n_inversion = 4)
    bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true)
    plot(bif_dia, xlabel = "k", ylabel = "A", xlims = (0, 6), ylims=(0,8))

will give the following plot. It looks like ksq is fixed at the original value of k^2 = 2^2 = 4, and does not change when k is varied. (Notice that the equilibrium where [A] = [B] occurs at k = 4).

[isaacsas]

Yeah, this is the kind of issue I was worried about. Thanks for the simple example.

[rveltz]

Hi,

Can you tell me what you would like to see implemented? Did you settle on this?