Equilibrate

[TorkelE]

This contains various changes for finding steady states and plotting bifurcation diagrams. I will try to list everything that is going on:

General Changes

• I have created a file "test_networks.jl" which creates 3 sets of 10 different reaction networks to be used in tests (currently only used in the steady states tests, but could potentially be used for other stuff). The first set is just standard reaction networks (of these there exists a minimal reaction network version), the second set contains reaction networks with hill functions (so were some parameters is required to be integer). The final set contains reaction networks were some reactants are neither created nor destroyed (like X1 + X2 ↔ X3). This will have an infinite amount of potential steady states, which can become finite by introducing additional equations or initial conditions.
• The macro, instead of returning the reaction network as last argument returns the reaction network as input to the manage_reaction_network!() function (which then outputs the reaction network. This have no real changes, but allows to make some post-creation manipulations of the reaction network without dealing with meta programming. This possibility is currently only used by the equilibration functionality. The function is currently not called by minimum reaction network (because they currently does not need any such processing). However, since the addodes!, addsdes!... functions uses eval these can not use the manage_reaction_network!() function, which have to be used separately. As a result of this, the way to get equilibrium features to you minimal reaction network is:

rn = @min_reaction_network begin
....
end p1 p2 p3
p = [1.,2.5,4.5]
addequi!(rn)
manage_reaction_network!(rn)
steady_states(rn,p)

• I have also added a separate file, equilibrate_utils.jl, containing all the steady state stuff.

New required packages
DiffEqBiological now requires HomotopyContinuation (does the actual solving), Plots (just by the functions which plots bifurcation diagrams), and DynamicPolynomials (since HC deals in polynomials). In addition DynamicPolynomials is now @reexport'ed by DiffEqBiological (I tried to solve stuff without this, but failed. Might still be possible).

Homotopy Continuation
Currently the only solver for finding steady states and doing bifurcations I have implemented is homotopy continuations. Basically it works very well for solving multivariate polynomial systems of n variables and n and n equations. So as long as the reaction rates only depends on other reactants as michaelis-menten functions and hill functions it should be fine (or some other rational polynomial expression). I believe this covers a very large portion of biochemical reaction network models.

While homotopy continuation is not really designed for the purpose it can be used to track the steady states solutions from one system to another, essentially creating bifurcation diagrams. Some times it can have problems tracking solutions through bifurcations (I had some problem with some saddle-node bifurcations, but sometimes these works as well). Right now I have a simple heuristics for this (checking if any solutions are lost while tracking from p1 to p2. If that is the case I track from p2 to p1 as well to get the missing paths). Right now this has worked for the system I have tried with. If we get problem with new systems it should be easy to introduce further heuristics (It is quite easy to realise when something funny is going on, probably a bifurcation. Just move the tracker forward a little bit, past the bifurcation, and continue).

I also looked at other methods for tracking bifurcations in Julia, but the docs were kinda confusing and I knew how to get it done with homotopy continuation so went with that. In the future we might want to switch this though.

Awesome work by the guys at https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl for actually implementing this algorithm. We#ll have to add their paper to the reference list.

How the update work

• To find steady states using HC we need a polynomial. Ideally we would want to create this polynomial when the reaction network is first created. However, parameters are represented as variables in the polynomial until they are set. This means that if one have a reaction rate with a reactant concentration to the power of a parameter (typically a hill equation) a polynomial cannot be created initially since it would have a variable to a power of a variable. Instead the macro creates a function, make_polynomial which takes parameters as input and creates the polynomial (this function is primarily accessed by another function). However, if no such parameters are present in the network a polynomial will be created and stored in the equilibratium_polynomial field.
• If some parameters are exponents which needs to be set afterwards this is done with the fix_parameter function, e.g.

rn = @reaction_network begin
(Y^n1,X^n1), X ↔ Y
end n1 n2
fix_parameter(rn,n1=2,n2=3)
steady_states(rn,p)

It can also be done by inserting a whole parameter vector (but then the whole vector needs to be inserted), e.g.

rn = @reaction_network begin
(k1*Y^n1,k2*X^n1), X ↔ Y
end k1 k2 n1 n2
p= [2.,3.,2,2]
fix_parameter(rn,p)
steady_states(rn)

• Certain systems, like (k1,k2), X ↔ Y, have an infinite amount of steady states (since neither reactant is created nor degraded). Here the equation system to solve is just the same equation two times over. To solve this additional information is needed. I have added the @fixed_concentration macro for this. The code to solve the system becomes:

rn = @reaction_network begin
(k1,k2), X ↔ Y
end k1 k2 
p = [1.0, 2.0]
@fixed_concentration rn X+Y=5.
steady_states(rn,p)

The fixed concentrations can also include parameters (which should be declared as a part of the network):

rn = @reaction_network begin
(k1,k2), X ↔ Y
end k1 k2 xyTot
p = [1.0, 2.0, 5.]
@fixed_concentration rn X+Y=xyTot
steady_states(rn,p)

The fixed concentrations are stored in a dictionary in the field fixed_concentrations. When the fixed concentrations are set the corresponding equations are replaced in the equilibrium polynomial. If the equilibrium polynomial is created after the fixed concentrations the fixed concentrations equations will be fetched from the fixed_concentrations fields and inserted automatically.

It should be possible to compute these relations from the reactions itself, and then the steady states should be computable by just giving some initial conditions of the system, instead of giving the algorithms this extra information. I would be surprised if someone already have looked into how this is done (or if anyone of you knows about it). In the future we might want to add that in, but this works ok for now.

• HC works much faster if it has a solution to a similar system to work with. So if you solve it for some parameter values it then goes very fast to solve it for another set. For some reason this works best if you just select random (complex) parameters instead of just taking the first set of parameters solved by the system. This HC template is stored in the homotopy_continuation_template field and is initially nothing. The first time you solve a system a template is created and stored here (a tuple of a random parameter set and the corresponding solutions). Alternatively the template can be initialised using a macro or function (@make_hc_template rn, make_hc_template(rn)).
• Once everything is prepared the function steady_states(reaction_network::DiffEqBase.AbstractReactionNetwork, params=Vector{Float64}()::Vector{Float64}, solver=HcSteadyStateSolver::Function). For the future one can define other solvers, but right now there is only the homotopy continuation one.
• Finally there is also a function stability(solutions,params,reaction_network) which returns true if the inserted fix point is stable).

Bifurcations

• If everything is ready to calculate steady states everything is also there for bifurcations. The function is simply: bifurcations(reaction_network,parameters,param,range). param is the bifurcation parameter and range is a tuple with the range were you want to do the diagram.
  (there is also an option for selecting an algorithm, with only one available). The output is a bifurcation_diagram structure (containing the parameter in question, the range, and a set of bifurcation_paths). The bifurcation_paths are a set of parameter value and corresponding fixed points, as well as the value of the corresponding eigenvalues.
• I have also added the option of making grid like bifurcation diagrams (bifurcations_grid and bifurcations_grid_s2). These does not actually track paths but simply solves the system at various grid points.

Plotting

• I initially wanted to simply overload the plot function, but was a bit unsure how to do it properly (my initial attempts failed). Now I have simple added two functions, bif_plot and bif_plot! which take bifurcation diagrams as input and plots them.
• Wanted to push this out now, but figure proper functions can be set of soon-ish.
• Rights now the color in the bifurcation diagrams depends on the eigenvalues of the corresponding fixed point (red if unstable with real eigenvalues, orange if unstable with imaginary ones, cyan if stable with imaginary ones, blue if stable with real ones).
• Some work might be needed to improve the plot recipes (e.g. better labelling), but that is for alter.

Finally
Still think the tests return some errors (which I think is due to a minor bug, will fix when I get time). The test got a bit long, would rewrite it a little bit more but got to go now, and figured it might be nice to have this up finally. Will update with more details.

[ChrisRackauckas]

why?

[TorkelE]

This is possibly me working with things I am not expert on, but as we require Plots in the tests, but no where else, I figure it belonged here.
(then whenever we should test plots is another question raised)

[ChrisRackauckas]

why reexport?

[TorkelE]

It has to do with the polyvars being created in the scope in which the macro is called, and some irritating issues that appeared were polyvars could have same names, but be different.

However, after your comment I had a look at this again. Due to the revised approach were I do a lot of stuff in the EquilibrateContent constructor I was able to move were the polyvars were created and could remove this reexport.

[ChrisRackauckas]

Please do "clean up" PRs separate from big changes...

[TorkelE]

You are ofcourse right. Will avoid this in the future.

[ChrisRackauckas]

do we need to plot in the test?

[TorkelE]

Basically these were to try the plot recipes, ensuring these generated no errors. Could remove them if they are unnecessary though.

[ChrisRackauckas]

why so many API functions? Can this not just work by dispatch?

[TorkelE]

@korsbo suggested a method were one gave a network topology type to designate which one to use, in the same spirit the solvers are given. I thought this would generate some messy since each method have requires different types of input

We have:
bifurcations requires AbstractReactionNetwork, Vector{Float64}, Symbol, Tuple{Number,Number} (technically the last one can be anything, but something like this which contains two numbers)
bifurcation_grid requires AbstractReactionNetwork, Vector{Float64}, Symbol, AbstractRange
bifurcation_grid_2d requires AbstractReactionNetwork, Vector{Float64}, Symbol, AbstractRange, Symbol, AbstractRange
bifurcation_diagram_grid requires AbstractReactionNetwork, Vector{Float64}, Symbol, AbstractRange, Symbol, Tuple{Number,Number}

I thought it would be clearer to have a separate function name for each, but it would be possible, and maybe a bit simpler, to have a single function , bifurcations, and depending on what form the input have it can have a dispatch for the corresponding type of diagram.

[ChrisRackauckas]

this should be with the plot recipes

[TorkelE]

Good catch, moved.

[ChrisRackauckas]

x?

[TorkelE]

I have removed all mentions of x and replaced them so that only u should be used.

[ChrisRackauckas]

Is x this the parameter being changed?

[TorkelE]

internal___polyvar___x[1:n] are the polyvars corresponding to the n reactants of the system- I've changed it so u is used instead of x to designate polyvars.

[ChrisRackauckas]

If we add constraints we should get a DAE, but I don't see a mass matrix being added to the ODEProblem?

[isaacsas]

That is exactly what I had in mind when I suggested using add_constraints instead of fix_concentrations for the macro name. Right now constraints are only added into the equilibration data structures, but ultimately we should extend this to allow the generation of DAE models (or at least allow them when using min reaction networks where ODE construction is deferred).

[TorkelE]

I have never really been involved with DAEs, but do we really get one? The simplest examplee would be X <--> Y, with the equations
dX/dt = Y - X
dY/dt = X - Y
Adding the information X+Y=1 doesn't really make the system different, since it is already defined by one differential equation for each variable, and an initial condition.

However, when we want to find steady states we get the system
Y - X = 0
X -Y = 0
which have an infinite amount of solutions, we need some additional information. Here a human "knows" the relation X+Y=1 (or = some parameter).@add_constraint adds this information, making the system
Y - X = 0
X -Y = 0
X+Y-1 =0
which we can solve.

It is likely that one can, given some initial condition, compute this constraint for a reaction network. Adding such an algorithm would save one from having from defining the constraints (although one would still have to supply some initial condition, or name for a parameter designating the combined concentration).

All that said, adding a mass matrix to the ODE problem would still make sense (I presume this could increase efficiency?). Or would creating a DAE problem from it also be a way to solve the system, given this additional information, in a more efficient way?

[ChrisRackauckas]

Suggested change

	Makes a template for solving the riven reaction networks steady states using homotopy continuation. Generally not needed as this is genrated automatically when netowkr is first solved *if one does not already exists). However, it can be sued to replace the current template, for whatever reason.
	Makes a template for solving for the given reaction networks steady states using homotopy continuation. Generally not needed as this is generated automatically when network is first solved *if one does not already exists*. However, it can be used to replace the current template, for whatever reason.

[TorkelE]

Fixed.

[ChrisRackauckas]

network then solver to match conventions?

[TorkelE]

Yes, I have shifted this.

[ChrisRackauckas]

don't need a second function. This is just a different dispatch of bifurcations_grid

[TorkelE]

Currently there only exists a single solver to solve steady states. It was written like this to prepare for if/when more was added. But yes, in the meantime we could remove this and have only a single function.

[ChrisRackauckas]

if you mean u, then use u. We shouldn't have one file in DiffEq use a different naming convention from all of the others.

[TorkelE]

Yes, I've went over the files to ensure there is not more x used.

[ChrisRackauckas]

Does it all need to be Float64s?

[TorkelE]

Probably me being overzealous defining types more narrowly then they should. used ::Vector would be better?

[korsbo]

By not restricting it at all you increase the users ability to draw on the rest of Julias ecosystem. You don't have to sacrifice performance if you do something like:

struct BifurcationPath{T} where T
    p_vals::T
    ...

However, just using ::Vector would be bad for performance since the compiler cannot figure out what you are giving it.

[ChrisRackauckas]

The issue is Int64. You should almost never use Int64, always Int.

[ChrisRackauckas]

Change that last thing to Int and it's good to go. Docs should follow up.

[ChrisRackauckas]

Let's continue to improve, but well done.

[isaacsas]

@TorkelE @korsbo Awesome job on this one! Thanks!