FAQs

How to index solution objects using symbolic variables and observables?

One can directly use symbolic variables to index into SciML solution objects. Moreover, observables can also be evaluated in this way. For example, consider the system

using Catalyst, OrdinaryDiffEqTsit5, Plots
rn = @reaction_network ABtoC begin
  (k₊,k₋), A + B <--> C
end

Let's convert it to a system of ODEs, using the conservation laws of the system to eliminate two of the species:

osys = convert(ODESystem, rn; remove_conserved = true)
osys = complete(osys)

\[ \begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= \mathtt{k_-} \left( - A\left( t \right) + \Gamma_{2} \right) - \mathtt{k.} \left( A\left( t \right) + \Gamma_{1} \right) A\left( t \right) \end{align} \]

Notice the resulting ODE system has just one ODE, while algebraic observables have been added for the two removed species (in terms of the conservation law constants, Γ[1] and Γ[2])

observed(osys)

\[ \begin{align} B\left( t \right) &= A\left( t \right) + \Gamma_{1} \\ C\left( t \right) &= - A\left( t \right) + \Gamma_{2} \end{align} \]

Let's solve the system and see how to index the solution using our symbolic variables

u0 = [osys.A => 1.0, osys.B => 2.0, osys.C => 0.0]
ps = [osys.k₊ => 1.0, osys.k₋ => 1.0]
oprob = ODEProblem(osys, u0, (0.0, 10.0), ps)
sol = solve(oprob, Tsit5())

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 19-element Vector{Float64}:
  0.0
  0.06602162921791198
  0.16167383723177658
  0.27797398569864795
  0.42472768210212675
  0.5991542465684496
  0.8069250399875878
  1.0494129741596487
  1.3328285260259545
  1.6634553095820483
  2.0523543727756204
  2.514619744984259
  3.073985028187367
  3.7668739415643824
  4.65256513049854
  5.821585305489345
  7.328825483832256
  8.975283435170365
 10.0
u: 19-element Vector{Vector{Float64}}:
 [1.0]
 [0.8836569748526829]
 [0.7588242777314499]
 [0.6540331462534721]
 [0.5681302593629073]
 [0.5062418977317176]
 [0.46461886630918997]
 [0.43937914557022106]
 [0.42544919657534686]
 [0.4186148714289667]
 [0.4156790083953944]
 [0.4146117644861058]
 [0.4142972990536097]
 [0.4142270387749867]
 [0.4142160588486048]
 [0.4142152933849376]
 [0.4142196830866055]
 [0.41425212965225616]
 [0.4142262306892306]

Suppose we want to plot just species C, without having to know its integer index in the unknown vector. We can do this using the symbolic variable C, which we can get at in several ways

sol[osys.C]

19-element Vector{Float64}:
 0.0
 0.11634302514731709
 0.2411757222685501
 0.3459668537465279
 0.4318697406370927
 0.49375810226828243
 0.5353811336908101
 0.5606208544297789
 0.5745508034246531
 0.5813851285710333
 0.5843209916046056
 0.5853882355138942
 0.5857027009463903
 0.5857729612250133
 0.5857839411513952
 0.5857847066150623
 0.5857803169133945
 0.5857478703477439
 0.5857737693107694

or

@unpack C = osys
sol[C]

19-element Vector{Float64}:
 0.0
 0.11634302514731709
 0.2411757222685501
 0.3459668537465279
 0.4318697406370927
 0.49375810226828243
 0.5353811336908101
 0.5606208544297789
 0.5745508034246531
 0.5813851285710333
 0.5843209916046056
 0.5853882355138942
 0.5857027009463903
 0.5857729612250133
 0.5857839411513952
 0.5857847066150623
 0.5857803169133945
 0.5857478703477439
 0.5857737693107694

To evaluate C at specific times and plot it we can just do

t = range(0.0, 10.0, length = 101)
plot(sol(t, idxs = C), label = "C(t)", xlabel = "t")

If we want to get multiple variables we can just do

@unpack A, B = osys
sol(t, idxs = [A, B])

t: 0.0:0.1:10.0
u: 101-element Vector{Vector{Float64}}:
 [1.0, 2.0]
 [0.8342522596177404, 1.8342522596177404]
 [0.7195976232151855, 1.7195976232151855]
 [0.6383995987555826, 1.6383995987555826]
 [0.5799339465200811, 1.579933946520081]
 [0.5373276995576052, 1.537327699557605]
 [0.5060148770362881, 1.506014877036288]
 [0.4828472548133267, 1.4828472548133267]
 [0.46563399153226215, 1.465633991532262]
 [0.4527899376478868, 1.4527899376478868]
 ⋮
 [0.41422620725249526, 1.4142262072524954]
 [0.4142173055979239, 1.414217305597924]
 [0.41421185440153846, 1.4142118544015385]
 [0.41421025171366543, 1.4142102517136654]
 [0.4142122146278574, 1.4142122146278573]
 [0.41421677928089273, 1.4142167792808928]
 [0.4142223008527756, 1.4142223008527757]
 [0.4142264535667362, 1.4142264535667362]
 [0.4142262306892306, 1.4142262306892306]

Plotting multiple variables using the SciML plot recipe can be achieved like

plot(sol; idxs = [A, B])

How to disable rescaling of reaction rates in rate laws?

As explained in the Reaction rate laws used in simulations section, for a reaction such as k, 2X --> 0, the generated rate law will rescale the rate constant, giving k*X^2/2 instead of k*X^2 for ODEs and k*X*(X-1)/2 instead of k*X*(X-1) for jumps. This can be disabled when directly converting a ReactionSystem. If rn is a generated ReactionSystem, we can do

osys = convert(ODESystem, rn; combinatoric_ratelaws=false)

\[ \begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= \mathtt{k{_-}} C\left( t \right) - \mathtt{k.} B\left( t \right) A\left( t \right) \\ \frac{\mathrm{d} B\left( t \right)}{\mathrm{d}t} &= \mathtt{k{_-}} C\left( t \right) - \mathtt{k.} B\left( t \right) A\left( t \right) \\ \frac{\mathrm{d} C\left( t \right)}{\mathrm{d}t} &= - \mathtt{k{_-}} C\left( t \right) + \mathtt{k.} B\left( t \right) A\left( t \right) \end{align} \]

Disabling these rescalings should work for all conversions of ReactionSystems to other ModelingToolkit.AbstractSystems.

When creating a ReactionSystem using the DSL, combinatoric rate laws can be disabled (for the created system, and all systems derived from it) using the @combinatoric_ratelaws option (providing false as its only input):

rn = @reaction_network begin
    @combinatoric_ratelaws false
    k, 2X --> 0
end

How to use non-integer stoichiometric coefficients?

using Catalyst
rn = @reaction_network begin
  k, 2.5*A --> 3*B
end

\[ \begin{align*} 2.5 \mathrm{A} &\xrightarrow{k} 3.0 \mathrm{B} \end{align*} \]

or directly via

t = default_t()
@parameters k b
@species A(t) B(t) C(t) D(t)
rx1 = Reaction(k,[B,C],[B,D], [2.5,1],[3.5, 2.5])
rx2 = Reaction(2*k, [B], [D], [1], [2.5])
rx3 = Reaction(2*k, [B], [D], [2.5], [2])
@named mixedsys = ReactionSystem([rx1, rx2, rx3], t, [A, B, C, D], [k, b])
mixedsys = complete(mixedsys)
osys = convert(ODESystem, mixedsys; combinatoric_ratelaws = false)
osys = complete(osys)

\[ \begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= 0 \\ \frac{\mathrm{d} B\left( t \right)}{\mathrm{d}t} &= - 2 k B\left( t \right) - 5 \left( B\left( t \right) \right)^{2.5} k + \left( B\left( t \right) \right)^{2.5} k C\left( t \right) \\ \frac{\mathrm{d} C\left( t \right)}{\mathrm{d}t} &= - \left( B\left( t \right) \right)^{2.5} k C\left( t \right) \\ \frac{\mathrm{d} D\left( t \right)}{\mathrm{d}t} &= 5 k B\left( t \right) + 4 \left( B\left( t \right) \right)^{2.5} k + 2.5 \left( B\left( t \right) \right)^{2.5} k C\left( t \right) \end{align} \]

Note, when using convert(ODESystem, mixedsys; combinatoric_ratelaws=false) the combinatoric_ratelaws=false parameter must be passed. This is also true when calling ODEProblem(mixedsys,...; combinatoric_ratelaws=false). As described above, this disables Catalyst's standard rescaling of reaction rates when generating reaction rate laws, see also the Reaction rate laws used in simulations section. Leaving this keyword out for systems with floating point stoichiometry will give an error message.

For a more extensive documentation of using non-integer stoichiometric coefficients, please see the Symbolic Stochiometries section.

How to set default values for initial conditions and parameters?

How to set defaults when using the @reaction_network macro is described in more detail here. There are several ways to do this. Using the DSL, one can use the @species and @parameters options:

using Catalyst
sir = @reaction_network sir begin
    @species S(t)=999.0 I(t)=1.0 R(t)=0.0
    @parameters β=1e-4 ν=0.01
    β, S + I --> 2I
    ν, I --> R
end

Model sir:
Unknowns (3): see unknowns(sir)
  S(t) [defaults to 999.0]
  I(t) [defaults to 1.0]
  R(t) [defaults to 0.0]
Parameters (2): see parameters(sir)
  β [defaults to 0.0001]
  ν [defaults to 0.01]

When directly constructing a ReactionSystem, we can set the symbolic values to have the desired default values, and this will automatically be propagated through to the equation solvers:

using Catalyst, Plots, OrdinaryDiffEqTsit5
t = default_t()
@parameters β=1e-4 ν=.01
@species S(t)=999.0 I(t)=1.0 R(t)=0.0
rx1 = Reaction(β, [S, I], [I], [1,1], [2])
rx2 = Reaction(ν, [I], [R])
@named sir = ReactionSystem([rx1, rx2], t)
sir = complete(sir)
oprob = ODEProblem(sir, [], (0.0, 250.0))
sol = solve(oprob, Tsit5())
plot(sol)

One can also build a mapping from symbolic parameter/species to value/initial condition and pass these to the ReactionSystem via the defaults keyword argument:

@parameters β ν
@species S(t) I(t) R(t)
rx1 = Reaction(β, [S,I], [I], [1,1], [2])
rx2 = Reaction(ν, [I], [R])
defs = [β => 1e-4, ν => .01, S => 999.0, I => 1.0, R => 0.0]
@named sir = ReactionSystem([rx1, rx2], t; defaults = defs)

Finally, default values can also be added after creating the system via the setdefaults! command and passing a Symbol based mapping, like

sir = @reaction_network sir begin
    β, S + I --> 2I
    ν, I --> R
end
setdefaults!(sir, [:β => 1e-4, :ν => .01, :S => 999.0, :I => 1.0, :R => 0.0])

How to specify initial conditions and parameters values for ODEProblem and other problem types?

To explicitly pass initial conditions and parameters we can use mappings from Julia Symbols corresponding to each variable/parameter to their values, or from ModelingToolkit symbolic variables/parameters to their values. Using Symbols we have

using Catalyst, OrdinaryDiffEqTsit5
rn = @reaction_network begin
    α, S + I --> 2I
    β, I --> R
end
u0 = [:S => 999.0, :I => 1.0, :R => 0.0]
p  = (:α => 1e-4, :β => .01)
op1  = ODEProblem(rn, u0, (0.0, 250.0), p)

while using ModelingToolkit symbolic variables we have

t = default_t()
u0 = [rn.S => 999.0, rn.I => 1.0, rn.R => 0.0]
p  = (rn.α => 1e-4, rn.β => .01)
op2  = ODEProblem(rn, u0, (0.0, 250.0), p)

How to include non-reaction terms in equations for a chemical species?

One method to add non-reaction terms into an ODE or algebraic equation for a chemical species is to add a new (non-species) unknown variable that represents those terms, let it be the rate of zero order reaction, and add a constraint equation. I.e., to add a force of (1 + sin(t)) to $dA/dt$ in a system with the reaction k, A --> 0, we can do

using Catalyst
t = default_t()
@variables f(t)
rx1 = @reaction k, A --> 0
rx2 = @reaction $f, 0 --> A
eq = f ~ (1 + sin(t))
@named rs = ReactionSystem([rx1, rx2, eq], t)
rs = complete(rs)
osys = convert(ODESystem, rs)

\[ \begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= f\left( t \right) - k A\left( t \right) \\ f\left( t \right) &= 1 + \sin\left( t \right) \end{align} \]

In the final ODE model, f can be eliminated by using ModelingToolkit.structural_simplify

osyss = structural_simplify(osys)
full_equations(osyss)

\[ \begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= 1 + \sin\left( t \right) - k A\left( t \right) \end{align} \]

How to modify generated ODEs?

Conversion to other ModelingToolkit.AbstractSystems allows the possibility to modify the system with further terms that are difficult to encode as a chemical reaction or a constraint equation. For example, an alternative method to the previous question for adding a forcing function, $1 + \sin(t)$, to the ODE for dA/dt is

using Catalyst
rn = @reaction_network begin
    k, A --> 0
end
osys = convert(ODESystem, rn)
dAdteq = equations(osys)[1]
t      = ModelingToolkit.get_iv(osys)
dAdteq = Equation(dAdteq.lhs, dAdteq.rhs + 1 + sin(t))

# create a new ODESystem with the modified equation
@named osys2  = ODESystem([dAdteq], t)

\[ \begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= 1 + \sin\left( t \right) - k A\left( t \right) \end{align} \]

How to override mass action kinetics rate laws?

While generally one wants the reaction rate law to use the law of mass action, so the reaction

using Catalyst
rn = @reaction_network begin
    k, X --> ∅
end
convert(ODESystem, rn)

\[ \begin{align} \frac{\mathrm{d} X\left( t \right)}{\mathrm{d}t} &= - k X\left( t \right) \end{align} \]

occurs at the (ODE) rate $d[X]/dt = -k[X]$, it is possible to override this by using any of the following non-filled arrows when declaring the reaction: <=, ⇐, ⟽, =>, ⇒, ⟾, ⇔, ⟺. This means that the reaction

rn = @reaction_network begin
    k, X => ∅
end
convert(ODESystem, rn)

\[ \begin{align} \frac{\mathrm{d} X\left( t \right)}{\mathrm{d}t} &= - k \end{align} \]

will occur at rate $d[X]/dt = -k$ (which might become a problem since $[X]$ will be degraded at a constant rate even when very small or equal to 0).

Note, stoichiometric coefficients are still included, i.e. the reaction

rn = @reaction_network begin
    k, 2*X ⇒ ∅
end
convert(ODESystem, rn)

\[ \begin{align} \frac{\mathrm{d} X\left( t \right)}{\mathrm{d}t} &= - 2 k \end{align} \]

has rate $d[X]/dt = -2 k$.

How to specify user-defined functions as reaction rates?

The reaction network DSL can "see" user-defined functions that work with ModelingToolkit. e.g., this is should work

using Catalyst
myHill(x) = 2*x^3/(x^3+1.5^3)
rn = @reaction_network begin
    myHill(X), ∅ --> X
end

\[ \begin{align*} \varnothing &\xrightarrow{\frac{2 X^{3}}{3.375 + X^{3}}} \mathrm{X} \end{align*} \]

In some cases, it may be necessary or desirable to register functions with Symbolics.jl before their use in Catalyst, see the discussion here.

How does the Catalyst DSL (@reaction_network) infer what different symbols represent?

When declaring a model using the Catalyst DSL, e.g.

using Catalyst
rn = @reaction_network begin
 (p,d), 0 <--> X
end

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{p} \mathrm{X} \end{align*} \]

Catalyst can automatically infer that X is a species and p and d are parameters. In total, Catalyst can infer the following quantities:

• Species (from reaction reactants).
• Parameters (from reaction rates and stoichiometries).
• (non-species) Variables (from the @equations option).
• Differential functions (from the @equations option).
• Observables (from the @observables option).
• Compound species (from the @compounds option).

Inference of species, variables, and parameters follows the following steps:

1. Every symbol explicitly declared using the @species, @variables, and @parameters options are assigned to the corresponding category.
2. Every symbol not declared in (1) that occurs as a reaction reactant is inferred as a species.
3. Every symbol not declared in (1) or (2) that occurs in an expression provided after @equations is inferred as a variable.
4. Every symbol not declared in (1), (2), or (3) that occurs either as a reaction rate or stoichiometric coefficient is inferred to be a parameter.

Here, in

using Catalyst
rn = @reaction_network begin
    @parameters p1
    @equations V ~ X + p1
 X + V + p1 + p2, 0 --> X
end

\[ \begin{align*} \varnothing &\xrightarrow{X + \mathtt{p1} + \mathtt{p2} + V\left( t \right)} \mathrm{X} \\ V\left( t \right) &= \mathtt{p1} + X\left( t \right) \end{align*} \]

p is first set as a parameter (as it is explicitly declared as such). Next, X is inferred as a species. Next, V is inferred as a variable. Finally, p2 is inferred as a parameter.

Next, if any expression D(...) (where ... can be anything) is encountered within the @equations option, D is inferred to be the differential with respect to the default independent variable (typically t). Note that using D in this way, while also using it in another form (e.g. in a reaction rate) will produce an error.

Any symbol used as the left-hand side within the @observables option is inferred to be an observable. These are by default assumed to be variables. It is possible to simultaneously explicitly declare an observable using the @species or @variables options (in the former case, the observable will be treated as a species instead). Using observables within most other expressions (e.g. as a reactant) will produce an error.

Any symbol declared as a compound using the @compound option is automatically inferred to be a system species.

Symbols occurring within other expressions will not be inferred as anything. These must either occur in one of the forms described above (which enables Catalyst to infer what they are) or be explicitly declared. E.g. having a parameter which only occurs in an event:

using Catalyst
rn_error = @reaction_network begin
    @discrete_events 1.0 => [X ~ X + Xadd] 
 d, X --> 0
end

is not permitted. E.g. here Xadd must be explicitly declared as a parameter using @parameters:

using Catalyst
rn = @reaction_network begin
    @parameters Xadd
    @discrete_events 1.0 => [X ~ X + Xadd]
 d, X --> 0
end

\[ \begin{align*} \mathrm{X} &\xrightarrow{d} \varnothing \end{align*} \]

It is possible to turn off all inference (requiring all symbols to be declared using @parameters, @species, and @variables) through the @require_declaration option.

How can I turn off automatic inferring of species and parameters when using the DSL?

This option can be set using the @require_declaration option inside @reaction_network. In this case all the species, parameters, and variables in the system must be pre-declared using one of the @species, @parameters, or @variables macros. For more information about what is inferred automatically and not, please see the section on @require_declaration.

using Catalyst
rn = @reaction_network begin
    @require_declaration
    @species A(t) B(t)
    @parameters k1 k2
    (k1, k2), A <--> B
end

\[ \begin{align*} \mathrm{A} &\xrightleftharpoons[\mathtt{k2}]{\mathtt{k1}} \mathrm{B} \end{align*} \]

What to be aware of when using remake with conservation law elimination and NonlinearProblems?

When constructing NonlinearSystems or NonlinearProblems with remove_conserved = true, i.e.

# for rn a ReactionSystem
nsys = convert(NonlinearSystem, rn; remove_conserved = true)

# or 
nprob = NonlinearProblem(rn, u0, p; remove_conserved = true)

remake is currently unable to correctly update all u0 values when the conserved constant(s), Γ, are updated. As an example consider the following

using Catalyst, NonlinearSolve
rn = @reaction_network begin
    (k₁,k₂), X₁ <--> X₂
    (k₃,k₄), X₁ + X₂ <--> 2X₃
end
u0 = [:X₁ => 1.0, :X₂ => 2.0, :X₃ => 3.0]
ps = [:k₁ => 0.1, :k₂ => 0.2, :k₃ => 0.3, :k₄ => 0.4]
nlsys = convert(NonlinearSystem, rn; remove_conserved = true, conseqs_remake_warn = false)
nlsys = complete(nlsys)
equations(nlsys)

\[ \begin{align} 0 &= - \mathtt{k_1} \mathtt{X_1}\left( t \right) + \mathtt{k_2} \mathtt{X_2}\left( t \right) - \mathtt{k_3} \mathtt{X_1}\left( t \right) \mathtt{X_2}\left( t \right) + \frac{1}{2} \left( - \mathtt{X_1}\left( t \right) - \mathtt{X_2}\left( t \right) + \Gamma_{1} \right)^{2} \mathtt{k_4} \\ 0 &= \mathtt{k_1} \mathtt{X_1}\left( t \right) - \mathtt{k_2} \mathtt{X_2}\left( t \right) - \mathtt{k_3} \mathtt{X_1}\left( t \right) \mathtt{X_2}\left( t \right) + \frac{1}{2} \left( - \mathtt{X_1}\left( t \right) - \mathtt{X_2}\left( t \right) + \Gamma_{1} \right)^{2} \mathtt{k_4} \\ 0 &= - \mathtt{X_1}\left( t \right) - \mathtt{X_2}\left( t \right) - \mathtt{X_3}\left( t \right) + \Gamma_{1} \end{align} \]

If we generate a NonlinearProblem from this system the conservation constant, Γ[1], is automatically set to X₁ + X₂ + X₃ = 6 and the initial values are those in u0. i.e if

nlprob1 = NonlinearProblem(nlsys, u0, ps)

NonlinearProblem with uType Vector{Float64}. In-place: true
u0: 3-element Vector{Float64}:
 1.0
 2.0
 3.0

then

nlprob1[(:X₁, :X₂, :X₃)] == (1.0, 2.0, 3.0)

true

and

nlprob1.ps[:Γ][1] == 6.0

true

If we now try to change a value of X₁, X₂, or X₃ using remake, the conserved constant will be recalculated. i.e. if

nlprob2 = remake(nlprob1; u0 = [:X₂ => 3.0])

NonlinearProblem with uType Vector{Float64}. In-place: true
u0: 3-element Vector{Float64}:
 1.0
 3.0
 3.0

compare

println("Correct u0 is: ", (1.0, 3.0, 3.0), "\n", "remade value is: ", nlprob2[(:X₁, :X₂, :X₃)])

Correct u0 is: (1.0, 3.0, 3.0)
remade value is: (1.0, 3.0, 3.0)

and

println("Correct Γ is: ", 7.0, "\n", "remade value is: ", nlprob2.ps[:Γ][1])

Correct Γ is: 7.0
remade value is: 7.0

However, if we try to directly change the value of Γ it is not always the case that a u0 value will correctly update so that the conservation law is conserved. Consider

nlprob3 = remake(nlprob1; u0 = [:X₂ => nothing], p = [:Γ => [4.0]])

NonlinearProblem with uType Vector{Float64}. In-place: true
u0: 3-element Vector{Float64}:
 1.0
 2.0
 3.0

Setting [:X₂ => nothing] for other problem types communicates that the u0 value for X₂ should be solved for. However, if we examine the values we find

println("Correct u0 is: ", (1.0, 0.0, 3.0), "\n", "remade value is: ", nlprob3[(:X₁, :X₂, :X₃)])

Correct u0 is: (1.0, 0.0, 3.0)
remade value is: (1.0, 2.0, 3.0)

and

println("Correct Γ is: ", 4.0, "\n", "remade value is: ", nlprob3.ps[:Γ][1])

Correct Γ is: 4.0
remade value is: 4.0

As such, the u0 value for X₂ has not updated, and the conservation law is now violated by the u0 values, i.e,

(nlprob3[:X₁] + nlprob3[:X₂] + nlprob3[:X₃]) == nlprob3.ps[:Γ][1]

false

Currently, the only way to avoid this issue is to manually specify updated values for the u0 components, which will ensure that Γ updates appropriately as in the first example. i.e. we manually set X₂ to the value it should be and Γ will be updated accordingly:

nlprob4 = remake(nlprob1; u0 = [:X₂ => 0.0])

NonlinearProblem with uType Vector{Float64}. In-place: true
u0: 3-element Vector{Float64}:
 1.0
 0.0
 3.0

so that

println("Correct u0 is: ", (1.0, 0.0, 3.0), "\n", "remade value is: ", nlprob4[(:X₁, :X₂, :X₃)])

Correct u0 is: (1.0, 0.0, 3.0)
remade value is: (1.0, 0.0, 3.0)

and

println("Correct Γ is: ", 4.0, "\n", "remade value is: ", nlprob4.ps[:Γ][1])

Correct Γ is: 4.0
remade value is: 4.0

Finally, we note there is one extra consideration to take into account if using structural_simplify. In this case one of X₁, X₂, or X₃ will be moved to being an observed. It will then always correspond to the updated value if one tries to manually change Γ. Let's see what happens here directly

nlsys = convert(NonlinearSystem, rn; remove_conserved = true, conseqs_remake_warn = false)
nlsys = structural_simplify(nlsys)
nlprob1 = NonlinearProblem(nlsys, u0, ps)

NonlinearProblem with uType Vector{Float64}. In-place: true
u0: 2-element Vector{Float64}:
 2.0
 1.0

We can now try to change just Γ and implicitly the observed variable that was removed will be assumed to have changed its initial value to compensate for it. Let's confirm this. First we find the observed variable that was elminated.

obs_unknown = only(observed(nlsys)).lhs

\[ \begin{equation} \mathtt{X_3}\left( t \right) \end{equation} \]

We can figure out its index in u0 via

obs_symbol = ModelingToolkit.getname(obs_unknown)
obsidx = findfirst(p -> p[1] == obs_symbol, u0)

3

Let's now remake

nlprob2 = remake(nlprob1; u0 = [obs_unknown => nothing], p = [:Γ => [8.0]])

NonlinearProblem with uType Vector{Float64}. In-place: true
u0: 2-element Vector{Float64}:
 2.0
 1.0

Here we indicate that the observed variable should be treated as unspecified during initialization. Since the observed variable is not considered an unknown, everything now works, with the observed variable's assumed initial value adjusted to allow Γ = 8:

correct_u0 = last.(u0)
correct_u0[obsidx] = 8 - sum(correct_u0) + correct_u0[obsidx]
println("Correct u0 is: ", (1.0, 2.0, 5.0), "\n", "remade value is: ", nlprob2[(:X₁, :X₂, :X₃)])

Correct u0 is: (1.0, 2.0, 5.0)
remade value is: (1.0, 2.0, 5.0)

and Γ becomes

println("Correct Γ is: ", 8.0, "\n", "remade value is: ", nlprob2.ps[:Γ][1])

Correct Γ is: 8.0
remade value is: 8.0

Unfortunately, as with our first example, trying to enforce that a non-eliminated species should have its initial value updated instead of the observed species will not work.

Summary: it is not recommended to directly update Γ via remake, but to instead update values of the initial guesses in u0 to obtain a desired Γ. At this time the behavior when updating Γ can result in u0 values that do not satisfy the conservation law defined by Γ as illustrated above.