Checklist for MT conversion

[isaacsas]

Thought it might be good to get a TODO list for converting to MT. @TorkelE, @ChrisRackauckas please feel free to add.

• [x] ODE support

• [x] SDE support

• [x] Jump support

• [x] Jump dep graphs

• [ ] Bifurcation support

• [ ] Check for performance regressions (e.g. on BCR example for ODEs and jumps)

• [x] Update API as appropriate.

• [x] Remake tests.

• [x] Nonlinnear system support

• [x] Support reaction network expansion

• [x] Add function for getting field values

• [x] Decide which fields to keep in reaction_network

• [x] Update API Docs

• [x] Update README

• [x] Update reaction model introduction in SciML docs.

• [ ] Update tutorials as appropriate.

• [x] Bring back Latexify support

[TorkelE]

Good idea, I ticked off ODE and SDE Support, and added a few things myself:

• Remake tests: We will need to remake the tests to some extent anyway. Previously I created (and just now updated) a set of test reaction networks for the steady state functionality. It would be possible to automatically use these in more cases. First, we should ensure that all current tests still applicable work, before scrambling things up.
• Nonlinear system support: It seems relatively easy to add in, and it might come in handy in some cases. One advantage of using Modelling toolkit IR should be that we can get these kinds of features vaguely for free.
• Support reaction network expansion: Adding a macro which adds reactions to a previously made reaction network (and updates the internal reaction_system function).
• Add function for getting field values: There seems to functions like states() for fetching attributes of MTK structures. We should add similar for reaction_networks.
• Decide which fields to keep in reaction_network: The structure have previously contained fields like params and p_matrix. Since we will be introducing breaking changes we might as well sort out exactly what to keep here (and if we wish to change some fields names). It seems like much of it will exist in reaction_system (or some other _system), so we might remove these and defer fetching these values to functions like states().

[ChrisRackauckas]

I'm not sure what to do about bifurcation support, but I don't know if we've gone the right route with it. I wonder if https://rveltz.github.io/PseudoArcLengthContinuation.jl/dev/ is an easier target to maintain. The current system is a bit wonky at times.

(@rveltz)

[rveltz]

My package is a bit wonky, I agree, especially for computing periodic orbits. The methods are here though but I am improving the interface a lot these days. I will publish a major update soon. The part concerning equilibria is polished in my opinion though.

I also have not had time to polish the interface to DiffEq

[ChrisRackauckas]

Oh no, the current package we're using is HomotopyContinuation.jl with DynamicPolynomials.jl, so it only supports polynomial rate functions (mass action kinetics). There's a lot of good there but I'm not entirely convinced it's the right route given a lot of the nonlinear stuff we do. PseudoArcLengthContinuation.jl is probably more appropriate.

Also I think they lead to upper bounding issues. https://travis-ci.org/github/SciML/ODE.jl/jobs/682337046#L274 depending on packages that aren't well-maintained is just a PITA, especially when DiffEqBio shows up in the testing suites.

[isaacsas]

Or bifurcation functionality could be split to a separate package that depends on DiffEqBio (or perhaps just a ReactionSystem). We can still advertise it here in the readme and tutorials.

[isaacsas]

@TorkelE we have functions for getting all the various field values in the current API, they just need to be updated to extract that info from the ReactionSystem as appropriate.

[TorkelE]

I think the current bifurcation support is good for making some "good-enough" plots, and quite useful in practical work (at least for my own). But ultimately, I don't think HomotopyContinuation as an algorithm is very suited for bifurcations (although I think it will remain useful for steady-state finding at specific parameter values).

PseudoArcLengthContinuation.jl is more probably more in the direction we would want. When I looked through alternatives while working at this a while ago that and some other wasn't very well documented enough, and also seemed like there was still some development to do.

While Julia is up to quality in many aspects with more established languages, I believe that we still lack a bifurcation tool at an even similar level to what you can find elsewhere. Hopefully, we will get that in the not too distant future. In the meantime, we probably should remain flexible, and no consider the functionality "done". If PseudoArcLengthContinuation.jl is having a major update soon though, maybe that will be good enough to settle the issue?

[rveltz]

You should also consider Bifurcations.jl for small systems. It seems very well written and in the spirit of Auto07p. The branch switching is far less complete though but given the base already coded, it is "just" a matter of time. It boils down to coding switching to limit cycles from codim 2 points.

[ChrisRackauckas]

I think our design here is a little bit off. I agree with @isaacsas here that what we probably need to do is just make sure MTK can generate the various functions these packages need and leave it at that.

[TorkelE]

Or bifurcation functionality could be split to a separate package that depends on DiffEqBio (or perhaps just a ReactionSystem). We can still advertise it here in the readme and tutorials.

Eventually, it would be nice to try and put bifurcation things in other places. Ideally, I think it would be nice to have a polynomial_system in MTK. There are quite a lot of powerful things one can do if one knows that a system is polynomial, and I think there are many cases not just in Systems Biology where this is the case. If this is possible, then one main part of the equilibrated functionality (generating something which is recognized as a polynomial) could be put elsewhere.

Then when we have a reliable bifurcation package (preferably with mutual support for and from MTK) the second main part of the equilibrate functionality (the heuristics for generating bifurcation diagrams using HomotopyContinuation.jl) could be removed. Potentially from there, we would have everything elsewhere, or just be left with a minimal set of convenience functions (without any actual heavy lifting in DiffEqBio).

You should also consider Bifurcations.jl for small systems. It seems very well written and in the spirit of Auto07p. The branch switching is far less complete though but given the base already coded, it is "just" a matter of time. It boils down to coding switching to limit cycles from codim 2 points.

At the time I skipped by this one for pretty much the same reasons as PseudoArcLengthContinuation.jl, I am still not sure to what extent these two actually overlap in functionality. Hopefully, it will clear with tim.

[rveltz]

What do you want to see ? What does the current bifurcation packages lack?

[TorkelE]

I was going to say "extensive documentation", since what mostly kept me away before was that I didn't actually know what the packages would do, and how. Giving it a second look now it seems to be much better for PseudoArcLengthContinuation.jl, and there is also some of it Bifurcations.jl.

I probably should read through what it says now before I say much more.

[TorkelE]

So I think one thing which is good with the current approach is that it is quick to use, basically the syntax is:

# This first part generates the equations for the system, it could be used to generate whenever have to be solved by the bifurcation algorithm.
rn = @reaction_network begin 
  (p,d), 0 ↔ X               
end p d 
params = [1.,2.] # Some parameter values.
bif = bifurcations(rn, params, :p, (0.1,5.)) # bifurcation(system,parameter values, target parameter, interval to vary it over)
plot(bif) #Plots the bifurcation diagram

Which gives you a diagram of the steady-state values as a function of p (in this case p/d = p/2.). Similarly, you can do for the Brusselator:

brusselator_network = @reaction_network begin
    A, ∅ → X
    1, 2X + Y → 3X
    B, X → Y
    1, X → ∅
end A B;
brusselator_p = [1.,4.]
bif = bifurcationsbrusselator_network, brusselator_p, :A, (0.1,8.)) 
plot(bif) 

It is not able to actually pinpoint the bifurcation types and points themself, nor the amplitude of periodic orbits (or existence, but it is sometimes implied from stability information). It will give you a plot of the steady-state values and their parameter values. Because you don't really need much fine-tuning from the user, you can spit out the diagram from your system without further considerations.

Is it actually possible to have something similar, but with the proper methods for bifurcation diagrams, but where the user is not required to "work" for it? Basically a somewhat reliable pipeline which takes you from the system to the diagram.

[rveltz]

For the equilibria and for the branching from the Hopf points, I would say yes. However, if the branch of periodic solutions has a period which increases without bounds (like close to homoclinic points) solutions, it requires a shooting method (PALC has it) or adaptive collocations (not in Bifurcations.jl).

I would say that both PALC and bifurcations should be able to handle most of the (simple) cases though.

For PALC, there are many options to disable the automation you describe because for the targeted systems (PDE,..._), already finding the equilibria or periodic orbits can be really tough., see this simple example to see that periodic orbits can be tough to compute

Your examples are nice! The analogous of https://tutorials.sciml.ai/html/models/04b-diffeqbio_III_steadystates.html is brusselator 1d. There, you will see the two methods for computing the periodic orbits (shooting and finite differences). This example (which is one the most complete in PALC) will soon be simplified a lot with automatic switching to periodic orbits.

Note that brusselator 1d is a simple example for the PDE targetted in PALC.

[TorkelE]

Sounds promising. I had a look at the Brusselator example, and at seems like PALC has a lot of other functionality baked in as well.

[TorkelE]

@isaacsas The JumpSystem Seems to work now. Have you started on a function Base.convert(::Type{<:ODESystem},rs::ReactionSystem) in the "reactionsystem.jl" file? I was thinking of attempting to do that today, so that the JUMP things could then be added to the updated DiffEqBio. I wanted to check with you first, however. To ensure that you had not already started, or preferred to do it yourself.

[isaacsas]

I just finished up JumpSystems yesterday, so haven't had a chance to add conversion from a ReactionSystem yet. If you want to take a first pass feel free, but it might be a bit of work since we'll want to analyze which reactions can be preferentially represented as MassActionJumps then ConstantRateJump and then VariableRateJumps (in decreasing priority). JumpSystems also don't have dependency graph generation, which means we can really only use Direct right now too. So there is some work to do. I'm at a (virtual) conference and grading the next few days, but hope to return to this over the weekend.

[TorkelE]

Sounds good. I will have a look at things, and while you are busy, and if things end up not working for me you can take over when available.

[TorkelE]

We have finally sent in that paper, and now I am fixing for the bifurcation support.

I have started looking at PALC, and think I should, relatively simply, be able to assemble something to make bifurcation diagram for a parameter p, over an interval (p1,p2). I have asked @rveltz for some help explaining certain features (and will probably ask some more, before this is done). The idea would be something like:

• Find all steady states for values p1 and p2.
• Loop over all those values, tracking the solution until it reaches the boundary (p1 or p2).
• Remove the end value from the steady state list (thus avoiding tracking the same path twice). That should cover those bifurcation diagrams we could expect to see from sensible reaction networks? (I have not started looking at periodic orbits and such yet, just reproducing what we already have)

What is needed is still a way to find steady states for a certain parameter value. I was thinking we should create a covert(NonlinearSystem,rs::ReactionSystem). That should enable some cases (but I am unsure how got it is at finding all the available states?)/

Regarding polynomials, it doesn't seem to be supported to create polynomials with parameters (in e.g. exponents), see https://github.com/JuliaAlgebra/MultivariatePolynomials.jl/issues/141. What would work is to create a function, which for a specified set of parameter values, creates a polynomial from a NonlinearSystem (or, slightly messier, an ODESystem). This should be good enough for all situation when one does not try to do it for a very large number of parameter values.

[TorkelE]

So it seems that there is a possibility of HC starting to support MTK stuff in the future, in which case we could hold off thinking about Polynomials for a while. Maybe settle for generating NonlinearSystem for now, and wait and see where that develops?

[isaacsas]

That sounds good to me!

[TorkelE]

Ok, so I should have something for making naive bifurcation diagrams right now. However, does anyone have a good heuristic for finding steady states of an ODESystem/NonlinearSystem at some parameter value? I would use it at the edge values, and then we just track those solutions to find the bifurcation paths. Currently using

function defaultSteadyStates(os::ODESystem,p_vals)
    ss_prob = SteadyStateProblem(ODEProblem(os,Pair.(os.states,ones(length(os.states))),(0.,1000.),Pair.(os.ps,p_vals)))
    return [solve(ss_prob).u]
end

but there is probably something better. It would be especially useful with something that can find several steady states, possibly simply by scanning phase space. If/when HC make it to MTK we should have something good for a large class of systems. But it would be good to have something until then, as well as for non-polynomial systems.

[ChrisRackauckas]

That is fine, but it should take in args...;kwargs... that go to solve so the method can be chosen. Also, using ones for the initial condition doesn't really make sense, that probably needs user input.

[TorkelE]

I have written some basic code for doing bifurcation analysis. I have only checked it on the reaction networks we used for testing in DiffEqBio, so will need some more checking. It should contain most of the features we already have (finding the initial steady states could be done both better and faster though).

In principle, it should take any ODESystem, and analyze it, but I have only really kept biochemical reaction networks in mind while doing this. Periodic orbits currently not supported.

I have also made a new plot receipt for the output. This one could use some improvement, especially by building on the receipts already present in PALC (I ran into some troubles when doing so initially, and just wanted to get something out. Will try to make an issue on it later).

I am a little bit unsure though exactly where this code belongs. If we aim to remove most things no directly related to the DSL from DiffEqBio, then maybe not here?

While it is nice to have something now, I am again unsure how much effort we want to put into it. This considering that both HC and PALC have updated lined up which might make things a lot easier.

I am attaching the code. The basic syntax should be:

using ModelingToolkit, DiffEqBiological, Plots

v0 = 0.005; v = 0.1; K = 2.8; n = 4;
kD = 10; kB = 100; kC = 0.1;
deg = 0.01; L = 0.;
#Normal σV model, without inducable production.
σV_p = [v0, v, K, n, kD, kB, kC, deg, L];
σV_network = @MT_reaction_network begin
    v0 + hill(σV,v,K,n), ∅ → (σV+A)
    deg, (σV,A,AσV) → ∅
    (kB,kD), A + σV ↔ AσV
    L*kC, AσV → σV
end v0 v K n kD kB kC deg L;
ode_sys = convert(ODESystem,σV_network.reaction_system)
bif = MT_bifurcations(ode_sys,σV_p,:L,(0.,1.5))
plot(bif)

Cannot attach files in GitHub, so here is the code:

using DifferentialEquations, ModelingToolkit, DiffEqBiological, PseudoArcLengthContinuation, Setfield, LinearAlgebra, Plots, RecipesBase

function MT_bifurcations(os::ODESystem,p_vals,p,pspan,args...;
        steadystateAlg=defaultSteadyStatesAuto,
        steaty_state_args=(), steaty_state_kwargs=(),
        init_guess_min = [ones(length(os.states))],
        init_guess_max = [ones(length(os.states))],
                verbosity = 1, tangentAlgo = BorderedPred(),
                theta=0.5,
                contParams=ContinuationPar(
                        dsmax = (pspan[2]-pspan[1])/500.,
                        dsmin = min((pspan[2]-pspan[1])/10000.,1e-5),
                        ds = min((pspan[2]-pspan[1])/500., 0.001),
                        maxSteps = 10000,
                        pMin = pspan[1],
                        pMax = pspan[2],
                        computeEigenValues = true,
            newtonOptions = NewtonPar(tol = 1e-8, verbose = true, maxIter = 15),
                        theta = theta,
            saveSolEveryNsteps = 1
                        ),
                kwargs...)
        contParamsBackwards = ContinuationPar(contParams,ds=-contParams.ds)
                p_idx = findfirst(p .== Symbol.(os.ps))
                ss_pmin = steadystateAlg(os,setindex!(copy(p_vals),pspan[1],p_idx),init_guess_min,steaty_state_args...;steaty_state_kwargs...)
                ss_pmax = steadystateAlg(os,setindex!(copy(p_vals),pspan[2],p_idx),init_guess_max,steaty_state_args...;steaty_state_kwargs...)
        ode_function = ODEFunction(os,jac=true)
        F = (x,p) -> ode_function.f(x,p,0.)
        J = (x,p) -> ode_function.jac(x,p,0.)
        fps = [ss_pmin...,ss_pmax...]; ignore = fill(false,length(fps));
        bifurcation_paths = Vector{ContResult}()
        for (i,ss) in enumerate(fps)
            ignore[i] && continue
            p_vals_here = setindex!(copy(p_vals),(i<=length(ss_pmin)) ? pspan[1] : pspan[2],p_idx)
            new_path = continuation(F, J, ss, p_vals_here, (@lens _[p_idx]), (i<=length(ss_pmin)) ? contParams  : contParamsBackwards ; tangentAlgo=tangentAlgo)
            push!(bifurcation_paths,new_path[1])
            ignore[argmin(map(s -> norm(s .- new_path[2].u),fps))] = true
        end
        return bifurcation_paths
end

function defaultSteadyStatesAuto(os::ODESystem,p_vals,init_guesses,args...;kwargs...)
    output = []
    for init_guess in init_guesses
        for factor in [1e-3,1e-2,1e-1,1,1e1,1e2,1e3]
            ss_prob = SteadyStateProblem(ODEProblem(os,Pair.(os.states,factor*init_guess),(0.,100.),Pair.(os.ps,p_vals)))
            sol = solve(ss_prob,args...;kwargs...).u
            (minimum([1.,map(o -> norm(o .- sol)/norm(o), output)...]) > 0.001) && push!(output,sol)
        end
    end
    return output
end

@recipe function f(contreses::Vector{ContResult}; var = 1,linewidth=3)
    label --> ""
    linewidth --> linewidth
    for (i, contres) in enumerate(contreses)
        @series begin
            primary --> i==1
            stabilities = maximum.(real.(getproperty.(contres.eig,:eigenvals))).<0
            color --> map(b -> b ? :blue : :red, stabilities)
            linestyle --> map(b -> b ? :solid : :dot ,stabilities)
            getproperty.(contres.sol,:p),getindex.(getproperty.(contres.sol,:x),var)
        end
    end
end

[rveltz]

Nice!

[TorkelE]

If anyone wants to actually test it (and not have to get this branch of DiffEqBio), here is a small example)

using ModelingToolkit

@parameters t µ
@variables x(t)
@derivatives D'~t

eqs = [D(x) ~ x^2 - μ^2 + 0.5*μ^4]
de = ODESystem(eqs,t,[x],[μ])

bif = MT_bifurcations(de,[1.],:μ,(-1.,2.))
plot(bif)

also, swapping to

bif = MT_bifurcations(de,[1.],:μ,(-2.,2.))

shows one weakness of the methodology (since there are no steady states for either µ=-2 or µ=2 the method cannot initiate the tracking of anything).

[rveltz]

What was the previous way doing in this case?

Also, I would like to mention the following which are two byproducts of a project I am working on:

1. I will soon push a method to compute a complete bifurcation diagram recursively (only for equilibria in a first case). This will require an initial guess though
2. I will push a deflated continuation method which does basically like 1. but dont need an initial guess

[TorkelE]

That sounds really cool!

Yes, this is very much a temporary solution. Basically, a part of making this update to the package would include scrapping the previous methods for bifurcation analysis (using HC for tracking as well as steady-state solving). It would be lots of work to enable that with MTK base internals, and it wouldn't make sense to put that effort in when we can do much better using PALC.

The ideal approach would be to wait a couple of months, until the updates of PALC and HC, which should be able to do some really cool stuff. However, that would mean also pushing this change forward substantially. If we wish to update DiffEqBio before then, I think it would be nice to have some kind of "temporary" solution for people to use in the meantime.

Alternatively, we could remove the functionality while waiting, and just tell people to use the previous version for stability analysis

[isaacsas]

I’ll start on API and readme updates tomorrow. After that I can update the tutorials, and/or temporarily take them down if we want to release.

[TorkelE]

Could you also have a look at the inactive files, especially under tests? To see which can be thrown away, and if there are any that we should keep.

[rveltz]

1. I will soon push a method to compute a complete bifurcation diagram recursively (only for equilibria in a first case). This will require an initial guess though

Did you see this update and the examples on the README?

[TorkelE]

Missed it, but looks good! Will have a closer look over the weekend