steady-state of reaction system

[yewalenikhil65]

Here is something is on mind for steadystate of reactionSystem along with conservationlaws

rs::ReactionSystem  # some reaction system

N = convert(Array{Int},conservationlaws(netstoichmat(rs)))
expr₁ = build_function(ratelaw,[species(rs);parameters(rs)])
f₁ = eval(expr₁[2]);
out₁ = zeros(Float64, nr);
S = netstoichmat(rs)'
function du_steady!(u)
    f₁(out₁, [u ;ps])
    du = S*out₁      # equivalent of ODESystem of ReactionSystem.. du/dt = S*ratevectors
    Conservation_Laws = N*u - N*x0;  # should follow conservation laws..here x0 is initial-condition 
    x✳ = [du; Conservation_Laws];     # x✳
end
using LsqFit     # steady state
result = LsqFit.lmfit(du_steady!, x0, Float64[])
x⃰ = result.param

Any ideas on this if anything is lacking ?

Specifically,

1. Is this enough to find steady-state(meaning anything more than conservationlaws is required) ?
2. Anything in SciML to use without LsqFit ?
3. Possibility of adding positivity constraint(concentration of species is positive) in du_steady! function ?

[ChrisRackauckas]

You want to use SteadyStateProblem instead. LsqFit isn't very stable.

[yewalenikhil65]

You want to use SteadyStateProblem instead. LsqFit isn't very stable.

Hi @ChrisRackauckas Is it good idea to construct separate SteadyState problem from reactionSystems in Catalyst which in addition to ODESystem also incorporates conservationlaws ?

[ChrisRackauckas]

I don't get what you mean.

[yewalenikhil65]

Here is something is on mind for steadystate of reactionSystem along with conservationlaws

function du_steady!(u)
    f₁(out₁, [u ;ps])
    du = S*out₁      # equivalent of ODESystem of ReactionSystem.. du/dt = S*ratevectors
    Conservation_Laws = N*u - N*x0;  # should follow conservation laws..here x0 is initial-condition 
    x✳ = [du; Conservation_Laws];     

This part I mean(lines with comments). Shouldn't this part be included by default in the steady-state calculation of reactionSystems always?

[ChrisRackauckas]

Yeah, the SteadyStateProblem infrastructure will do all of that in a better way. What's a case where it won't? Maybe we should just document more clearly NonlinearSystem stuff.

[isaacsas]

I don't think we have any SteadyStateProblem or NonlinearProblem examples currently. Are SteadyStateProblems recommended anymore, or should we just be pointing towards NonlinearProblem?

@yewalenikhil65 are you advocating for making the system over-determined by including the conservation laws in the set of equations? I don't think we'd want to do that -- instead it would be better to get the PR that adds system reduction via conservation laws resurrected and reduce the system size. I think it stalled looking for a more scalable linear algebra approach than the Smith normal form, and needed to get updated to use observables.

[yewalenikhil65]

Yeah, the SteadyStateProblem infrastructure will do all of that in a better way. What's a case where it won't? Maybe we should just document more clearly NonlinearSystem stuff.

This issue came from following thought.

## REGULATION OF NEURAL STEM CELL
using ModelingToolkit, Catalyst, DifferentialEquations
function funcsym(S::Symbol, args...)
    Num(Variable{ModelingToolkit.FnType{Tuple{Any},Real}}(S,args...))
end

@parameters t
X = [funcsym(:X,i)(t) for i = 1:21]  # species
k = [funcsym(:k,i)(t) for i = 1:20]  # parameters


rxs = [Reaction(k[1]/(X[1] + (1+X[4])*(1+X[5])), [X[1]], [X[4], X[5]],[1],[1,1]),
      Reaction(k[2]/(X[1] + (1+X[4])*(1+X[5])), [X[4], X[5]], [X[1]], [1,1], [1]),
      Reaction(k[3], [X[4], X[6]], [X[7]], [1, 1], [1]),
      Reaction(k[4], [X[7]], [X[4], X[6]], [1], [1, 1]),
      Reaction(k[5], [X[7]], [X[8]], [1], [1]),
      Reaction(k[6], [X[8]], [X[7]], [1], [1]),
      Reaction(k[7], [X[9] , X[10]], [X[11]], [1,1],[1]),
      Reaction(k[8], [X[11]], [X[9] , X[10]], [1],[1,1]),
      Reaction(k[9], [X[2], X[12]],[X[13]], [1,1], [1]),
      Reaction(k[10], [X[13]], [X[2], X[12]], [1], [1,1]),
      Reaction(k[11], [X[13]], [X[8]],[1], [1]),
      Reaction(k[12], [X[8]], [X[13]], [1], [1]),
      Reaction(k[13], [X[14], X[15]],[X[16]], [1,1],[1]),
      Reaction(k[14], [X[16]], [X[14], X[15]], [1],[1,1]),
      Reaction(k[15]*X[16]/(1 + X[17] + X[18] + X[19] + X[18]*X[19]), [X[17], X[18]],[X[19]], [1,1],[1]),
      Reaction(k[16]*X[16]/(1 + X[17] + X[18] + X[19] + X[18]*X[19]), [X[19]],[X[17],X[18]], [1],[1,1]),
      Reaction(k[17]/(1 + X[19]), [X[20], X[21]],[X[3]], [1,1], [1]),
      Reaction(k[18]/(1 + X[19]), [X[3]],[X[20], X[21]], [1], [1,1]),
      Reaction(k[19], [X[21]], [X[8]], [1], [1]),
      Reaction(k[20], [X[8]], [X[21]], [1], [1])
     ]

@named rs = ReactionSystem(rxs, t, X, k)

# Rate constants and initial condition
ps = [√2*5/101; √2*5/101; 1; 1; 1; 1; 1; 1;
    1/101; 1/101; 1; 1; 1; 1; 1; 1; 1/101; 1/101; 1;1];

x0 = [5.0; 5; 0; 0; 0; 5; 0; 0; 5; 5; 0;5; 0; 5; 5; 0; 5; 5; 0; 5; 5]

tspan = (0.0, 550.0)

First method...

odesys = convert(ODESystem, rs)
oprob = ODEProblem(odesys, x0, tspan, ps)
sprob = SteadyStateProblem(oprob)
ssol = solve(sprob,SSRootfind())
julia> ssol.u
21-element Vector{Float64}:
 -0.07162857178504112
  2.534820671446184
 -0.2873517061039479
 -0.020636518322810105
  3.470962042186102
  4.696581928732968
 -0.09692109910892555
 -0.09692109902499624
  5.20402236314834
  0.9566264965606641
  4.978305681754148
 -0.03823588006233081
 -0.09692109934408677
  1.6294904630882535e-20
 -0.971800461612748
 -9.505627145074941e-21
  4.999999999999997
  5.0
  1.9354238582636987
  2.964800316349802
 -0.09692109902235302

if we compare this with..second method

S = netstoichmat(rs)'
N = convert(Array{Int},conservationlaws(netstoichmat(rs)))
expr₁ = build_function(ratelaw,[species(rs);parameters(rs)])
f₁ = eval(expr₁[2]);
out₁ = zeros(Float64, nr);

function du_steady!(u)
    f₁(out₁, [u ;ps])
    du = S*out₁
    Conservation_Laws = N*u - N*x0;  # should follow conservation laws
    x✳ = [du; Conservation_Laws];     # x✳
end
using LsqFit     # steady state
result = LsqFit.lmfit(du_steady!, x0, Float64[])
x⃰ = result.param
julia> x⃰
21-element Vector{Float64}:
 2.0475481971503635
 1.378747087521816
 3.276422843405382
 0.6935077467301308
 2.952451802849641
 2.741055943880478
 1.9009435313008654
 1.900943531298466
 1.79128784747792
 1.79128784747792
 3.20871215252208
 1.3787470875218064
 1.9009435313004441
 1.79128784747792
 1.79128784747792
 3.20871215252208
 1.7912878474779201
 1.7912878474779201
 3.2087121525220805
 1.7235771565946263
 1.9009435312925895

The second solution, where i have explicitly added conservationlaws as constraint in the nonlinearproblem, seems more accurate steady-state solution.

[yewalenikhil65]

@yewalenikhil65 are you advocating for making the system over-determined by including the conservation laws in the set of equations? I don't think we'd want to do that -- instead it would be better to get the PR that adds system reduction via conservation laws resurrected and reduce the system size. I think it stalled looking for a more scalable linear algebra approach than the Smith normal form, and needed to get updated to use observables.

I wanted to know if there is a way to calculate, du/dt = f(u,p,t) = 0, with constraint that g(u,t) = 0 where du/dt = f(u,p,t) is formed from reactionSystem and g(u,t) is formed from conservationlaws Isn't this common in almost all reaction -networks ? I think adding conservationlaws as constraint gives us more unique and physically/chemically plausible steady-state(correct me please if this isn't correct)

SteadyStateProblem docs i think direct only for f(u,p,t)= 0 part.

Now i see NonlinearSolver.jl ,https://nonlinearsolve.sciml.ai/dev/ it seems little relatable

[ChrisRackauckas]

I don't think we have any SteadyStateProblem or NonlinearProblem examples currently. Are SteadyStateProblems recommended anymore, or should we just be pointing towards NonlinearProblem?

Well SteadyStateProblem is a wrapper over an ODEProblem to turn it into a NonlinearProblem with an interpretation on t. It may end up just being fully deprecated to NonlinearProblem though, that would be a good idea, but for now it's slightly different because of the function form still containing t. But we should probably have ReactionSystem have a conversion to a NonlinearSystem. structural_simplify on that will give a much more efficient solution due to the tearing process, so that will be a good route to have.

du/dt = f(u,p,t) = 0, with constraint that g(u,t) = 0

You don't want it overconstrained though. You want to eliminate some of the other equations first.

[ChrisRackauckas]

The second solution, where i have explicitly added conservationlaws as constraint in the nonlinearproblem, seems more accurate steady-state solution.

DynamicSS on a DAE solver should be a bit better.

[TorkelE]

There are methods (like in SteadyStateProblem?) that are dependent on an initial condition one give, right? In these cases conversation law constraints should not be necessary (since these are implied from the initial condition).

However, in other cases, it would make sense to have them, since they are required to find a steady state (unless the user gives additional information. But creating overdetermined systems by default seems dangerous (although I think some solvers should be able to do them, think e.g. HomotopycOntinuation.jl is able to). If we have some method for removing some unnecessary equations in the case of steady state problems things should be fine though.

[isaacsas]

We have an untested conversion to NonlinearProblem at

https://github.com/SciML/Catalyst.jl/blob/8e3e05c6d1d501c8994bcdee471cd2a19ea2e6b3/src/reactionsystem.jl#L476

but it doesn’t appear to drop any time dependence from states or equations.

[yewalenikhil65]

The second solution, where i have explicitly added conservationlaws as constraint in the nonlinearproblem, seems more accurate steady-state solution.

DynamicSS on a DAE solver should be a bit better.

This works ..Thanks @ChrisRackauckas . used DynamicSS(Tsit5()) and gives expected answer

[isaacsas]

I think this was settled, but @yewalenikhil65 feel free to reopen if there is something actionable in Catalyst still to do.