TODO per discussion of 11/25

[isaacsas]

• [x] Eliminate @reaction_func and how Hill functions work. (Niklas)
• [x] Fix ReactionNetworkImporters. (Sam)
• [x] Fix BCR loading and update BCR benchmarks. (Sam)
• [x] Work towards previous BCR performance. (Sam)

Tutorials:

• [x] Bifurcations (adapting https://github.com/rveltz/BifurcationKit.jl/issues/27) (Torkel and Alex)
• [x] Adding forcing to a network. (Sam)
• [x] Parameter estimation. (Chris and Alex)
• [ ] Discovery of model components. (Chris and Alex)
• [ ] Discovery and elimination of conservation laws and such.
• [ ] Diagram of the IR (Chris)
• [ ] Diagram of the user features (Sam and Torkel)
• [ ] Structural transformations?

[isaacsas]

@ChrisRackauckas @TorkelE @korsbo The non-tutorial issues that we had listed here are now basically done aside from getting BCR back to where we were at before moving to MT. (@shashi and I have been working more on that, and will hopefully meet up to make some more progress this month.)

I think with those tutorials we listed done we could use them in building the paper.

[yewalenikhil65]

Discovery and elimination of conservation laws and such. @ChrisRackauckas @isaacsas , I am thinking along following lines for the point in to-do list

The dynamics of most CRNs is written as dx/dt= S. v(x(t)) , where S is a net-stoichiometric matrix, and v(x(t)) is vector of oderatelaws.(rn.eqs) If κ ∈ ker(Sᵀ) ,We can see that κᵀ. dx/dt == 0 i.e. κᵀ.x = constant is satisfied. So κ could be easily found from nullspace of Sᵀ (scaled-to-integer version of columns)

Now to find out what is the conserved quantity, we can easily write this as, κᵀ.x = κᵀ.x(t=0) (conservation law, where κᵀ.x(t=0) is the conserved quantity ) So, the non-zero quantities in κ vector are the species involved in conservation laws. ?

example

using Catalyst, ModelingToolkit

rn = @reaction_network begin
           k₁, A --> B
           k₂, B + B --> B + C
           k₃, B + C --> A + C
end k₁ k₂ k₃

u0 = Float32[1; 0; 0]; # initial condition

# stoichiometric matrix
S = -substoichmat(rn)' + prodstoichmat(rn)';

using LinearAlgebra
κ = nullspace(S')

# conservation law
lhs_law = κ'*states(rn);
rhs_law = κ'*u0;
conservation_law = [lhs_law rhs_law]

julia> conservation_law
1×2 Array{Any,2}:
 0.57735A(t) + 0.57735B(t) + 0.57735C(t)  0.57735

[isaacsas]

If you add the conservation law to the generated ODEsystem does structural simplification then eliminate one of the ODEs?

[ChrisRackauckas]

If it's a linear conservation law it should. Nonlinear is much harder, won't right now, and only cases would in the future (since of course symbolic nonlinear solving is very nontrivial).

[yewalenikhil65]

If you add the conservation law to the generated ODEsystem does structural simplification then eliminate one of the ODEs?

I tried replacing simply the last ODE in above mentioned code..as follows, and but couldn't create an ODEProblem out of it.

osys = convert(ODESystem,rn)
osys.eqs[3] = lhs_law[1] - rhs_law[1] ~ 0 
oprob = ODEProblem(osys,[1.0,0.0,0.0],(0.0,1e5),(0.04,3e7,1e4))
ERROR: Only semi-explicit constant mass matrices are currently supported. Faulty equation: 0.5773502691896258A(t) + 0.5773502691896258B(t) + 0.5773502691896256C(t) - 0.5773502691896258 ~ 0.

[isaacsas]

Don't try to change an equation, just add the constraint as a new equation (i.e. push it onto osys.eqs) and then run structural_simplify to get a new system. Take look at the latest ModelingToolkit docs to see how to mix ODEs and constraints in component-based systems, and then use structural_simplify.

[yewalenikhil65]

Don't try to change an equation, just add the constraint as a new equation (i.e. push it onto osys.eqs) and then run structural_simplify to get a new system. Take look at the latest ModelingToolkit docs to see how to mix ODEs and constraints in component-based systems, and then use structural_simplify.

julia> push!(osys.eqs, (lhs_law[1] ~ rhs_law[1]))
4-element Array{Equation,1}:
 Differential(t)(A(t)) ~ k₃*B(t)*C(t) - (k₁*A(t))
 Differential(t)(B(t)) ~ k₁*A(t) - (0.5k₂*(B(t)^2)) - (k₃*B(t)*C(t))
 Differential(t)(C(t)) ~ 0.5k₂*(B(t)^2)
 0.5773502691896258A(t) + 0.5773502691896258B(t) + 0.5773502691896256C(t) ~ 0.5773502691896258

julia> new_sys = structural_simplify(osys)
ERROR: 

... Internal error in Tearing.jl: vs[1] = 0.

[isaacsas]

Here is a minimal example that doesn't work:

using Catalyst, ModelingToolkit
rn = @reaction_network begin
       k1, A+B --> C
       k2, C --> A + B
       end k1 k2
osys = convert(ODESystem, rn)
@parameters t; @variables  A(t) B(t) C(t)
eq = 0 ~ 1 - (A + C)
push!(osys.eqs, eq)
sys = structural_simplify(osys)

This gives the normal tearing error:

ERROR: 

... Internal error in Tearing.jl: vs[1] = 0.

@ChrisRackauckas should this work?

[ChrisRackauckas]

@YingboMa you got another one today 😆

[YingboMa]

Woah, Jung's theory of collective unconscious confirmed.

[YingboMa]

But, what is the expected behavior when adding another algebraic equation 0 ~ 1 - (A + C)? Do we assume that's true and eliminate differential equations? I am not sure what to do when the system itself is not consistent.

[isaacsas]

@YingboMa See above, where I asked "If you add the conservation law to the generated ODEsystem does structural simplification then eliminate one of the ODEs?" and Chris indicated that it would...

So currently linear conservation laws are not detected / reduced? i.e. without the algebraic equation structural_simplify leaves the system unchanged.

[YingboMa]

We reduce them if the system is actually consistent.

[isaacsas]

What is your definition of a consistent ODE system? If I leave out the algebraic equation and just include the three ODEs is this system "not consistent" in your definition?

[YingboMa]

 Differential(t)(A(t)) ~ k2*C(t) - (k1*A(t)*B(t))
 Differential(t)(B(t)) ~ k2*C(t) - (k1*A(t)*B(t))
 Differential(t)(C(t)) ~ k1*A(t)*B(t) - (k2*C(t))

already uniquely define A, B, and C. If 0 ~ 1 - A(t) - C(t) is added, then C = 1 - A, so it conflicts with

  Differential(t)(C(t)) ~ k1*A(t)*B(t) - (k2*C(t))

So, we have a choice here. We can eliminate the differential equation or the added algebraic equation. It's not clear which one should "win".

[yewalenikhil65]

I believe the issue was that of eliminating a differential equation corresponding to some species and replacing it by a (linear) conservation law if there exists any such in the system https://github.com/SciML/Catalyst.jl/issues/280#issuecomment-804021518 found in similar manner.

It's like a automatic conversion from a ODEProblem to DAEProblem? Or replacing the C variable with 1 - A in first two differential equations so that the first two ODEs have only A and B as variables. ?

P. S :- I see application of this for finding automatically a proper constraint while fitting NeuralODEs to chemical reactions as demonstrated in some of the example in DiffEqFlux

[isaacsas]

Ignore the algebraic equation, that was from my misunderstanding of Chris’s reply. I thought the algebraic equation would have precedence and lead to the elimination of one ODE and conversion of one variable to an observable...

The independent system of three ODEs has an implicit conservation relation that should allow us to eliminate one of the three variables (knowing the initial condition). That is what we were discussing. We can easily solve for this relation ourselves, but I was hoping that the structural simplification could somehow handle that reduction already.

[YingboMa]

I misunderstood the question as well. So you want to find conservation laws from ODEs and convert them into DAEs, so in this case, naturally, the algebraic equations take precedence. Unfortunately, structural_simplify won't understand that for now as it views all equations equally.

[YingboMa]

I think by using the matching algorithm, one can find what differential equations to replace with added algebraic equations. I am not sure if it is fully correct. I have to think more about this.

[isaacsas]

@YingboMa Once we get @kaandocal's PR merged (i.e. adding tests for accuracy/correctness, cleaning up and such), it would be great to work on making the conservation law detection and reduction code more efficient. This might be a relevant reference:

http://web.math.ku.dk/~wiuf/journalWiuf/SIAMApplMath72.pdf

More generally though, I have to think there are standard methods for large systems for detecting and eliminating conservation laws in the chemical engineering literature... Maybe next week I'll spend a little time searching around for more on what is used there.

[isaacsas]

Another reference discussing methods for calculating conservation laws:

https://www.sciencedirect.com/science/article/pii/S0301462203002503?casa_token=lwJ3ukNRX6cAAAAA:oA1T3kKDFf9ebVYhPkO-j3v7XneNP2YhnfEYBmrOhAinjkUpmOEDGVqmp03--EnyjMJndAMg#BIB8

[isaacsas]

A library with mostly BSD code for doing structural analysis (including detecting conservation laws):

Paper: https://doi.org/10.1016/j.biosystems.2018.05.008 Library: https://github.com/sys-bio/Libstructural

[yewalenikhil65]

@ChrisRackauckas @TorkelE @korsbo The non-tutorial issues that we had listed here are now basically done aside from getting BCR back to where we were at before moving to MT. (@shashi and I have been working more on that, and will hopefully meet up to make some more progress this month.)

I think with those tutorials we listed done we could use them in building the paper.

@isaacsas , I hope writing a paper on catalyst is an ongoing issue. I would also like to be a contributor if that's okay with sciml organization (still need to add a tutorial based on newly added API functions) . ReactionNetworkImporters will be treated as part of this , right?

[ChrisRackauckas]

It's happening. In fact, the best thing would be to get someone to do some benchmarks of gillespy2 and such so we can have some comprehensive benchmarks for the paper. If you're willing to do those, you're in. Let's try to get all of that onto SciMLBenchmarks.

[yewalenikhil65]

It's happening. In fact, the best thing would be to get someone to do some benchmarks of gillespy2 and such so we can have some comprehensive benchmarks for the paper. If you're willing to do those, you're in. Let's try to get all of that onto SciMLBenchmarks.

@ChrisRackauckas I was considering on networkAPI functions that I added in past some PRs, but I would take a look at gillespy2 too. Is this gillespy2 you were talking about? https://stochss.github.io/GillesPy2/ There also seems a julia wrapper for it available

[ChrisRackauckas]

Yes. The steps to benchmark would be to first create a benchmark of the Julia wrapper vs Python directly, and understand the overhead. If the Julia wrapper is then good, we can create a SciMLBenchmark weave file benchmarking against Catalyst. Maybe throw in some of the Gillespie.jl and R GillespieSSA results too.

[isaacsas]

I think the network API stuff needs more testing before we start advertising it. I'd be hesistant right now to say more about it in a paper than that we are currently adding such functionality with maybe a couple sentences.

What Chris is mentioning would be good to have in the paper though as a subsection. If not for Catalyst then for a DiffEqJump paper. So it would be used either way.

[yewalenikhil65]

Yes. The steps to benchmark would be to first create a benchmark of the Julia wrapper vs Python directly, and understand the overhead. If the Julia wrapper is then good, we can create a SciMLBenchmark weave file benchmarking against Catalyst. Maybe throw in some of the Gillespie.jl and R GillespieSSA results too.

I can try this

[yewalenikhil65]

Yes. The steps to benchmark would be to first create a benchmark of the Julia wrapper vs Python directly, and understand the overhead. If the Julia wrapper is then good, we can create a SciMLBenchmark weave file benchmarking against Catalyst. Maybe throw in some of the Gillespie.jl and R GillespieSSA results too.

I can try this and build upon the lines of https://benchmarks.sciml.ai/html/Jumps/Diffusion_CTRW.html

[isaacsas]

This is pretty dated at this point. If someone wants to add the missing tutorials it would be appreciated but I don't think the issue is warranted further at this point.