Issue with scoping rules of parameter (and variable) declaration in components (Catalyst vs ModelingToolkit)

Specific Domains Modelling & Simulations
,
1

Working with Catalyst v15.0.8 and ModelingToolkit v9.80.5,
I found a strange and misleading behavior w.r.t. parameter declarations in reaction components: suppose I declare a parameter (k) in two different components (n1 and n2) with two different default values (say 1.0 in n1 and 2.0 in n2), and in both components I also declare a variable x with k as initial condition. I would thus expect that the initial condition of n1.x is 1.0 and n2.x is 2.0. This is not the case, both variable have the same initial value! More concretely:

n1 = @network_component n1 begin
    @parameters k = 1.0
    @variables x(t) = k
end

n2 = @network_component n2 begin
    @parameters k = 2.0
    @variables x(t) = k
end

@named sys1  = System([D(n2.x) ~ 0, D(n1.x) ~ 0], t)
prob1 = ODEProblem(complete(sys1))

I get:

ODEProblem with uType Vector{Float64} and tType Nothing. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: false
timespan: (nothing, nothing)
u0: 2-element Vector{Float64}:
 1.0
 1.0

So both n1.x and n2.x have 1.0 as initial condition. If I swap both equations:

@named sys2  = System([D(n1.x) ~ 0, D(n2.x) ~ 0], t)
prob2 = ODEProblem(complete(sys2))

I get:

ODEProblem with uType Vector{Float64} and tType Nothing. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: false
timespan: (nothing, nothing)
u0: 2-element Vector{Float64}:
 2.0
 2.0

both n1.x and n2.x have now 1.0 as initial condition.

If I translate this code in ModelingToolkit as follows:

@mtkmodel Decl1 begin
    @parameters begin
        k = 1.0
    end
    @variables begin
        x(t) = k
    end
end

@mtkmodel Decl2 begin
    @parameters begin
        k = 2.0
    end
    @variables begin
        x(t) = k
    end
end

@mtkmodel Eq begin
    @components begin
        n1 = Decl1()
        n2 = Decl2()
    end
    @equations begin
        D(n2.x) ~ 0
        D(n1.x) ~ 0
    end
end

@mtkbuild sys3 = Eq()
prob3 = ODEProblem(sys3)

then I get the expected initial conditions:

ODEProblem with uType Vector{Float64} and tType Nothing. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: false
timespan: (nothing, nothing)
u0: 2-element Vector{Float64}:
 1.0
 2.0

Is it a bug or a feature of @reaction_component ?

2

I wouldn’t recommend mixing Catalyst ReactionSystems and ModelingToolkit.System as you are doing. These were not designed to work together and have not been tested to work together in any way. What is it you actually want to do? If you make your model purely in Catalyst is there some problem you are encountering?

3

Note you can make a model that includes both network components in Catalyst, and then everything seems fine:

@named rn = ReactionSystem([], default_t(); systems = [n1, n2])
oprob = ODEProblem(complete(rn), [], (0.0, 10.0))

which gives

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: false
timespan: (0.0, 10.0)
u0: 2-element Vector{Float64}:
 1.0
 2.0
4

Thanks a lot for your answer!
Following your suggestion, if add a reaction to the reaction system, the problem arises:

n1 = @network_component n1 begin
    @parameters k = 1.0
    @species x(t) = k
end

n2 = @network_component n2 begin
    @parameters k = 2.0
    @species x(t) = k
end

@named rn = ReactionSystem([Reaction(1.0,[n1.x],[n2.x])], default_t(); systems = [n1, n2])

which gives

Model rn:
Subsystems (2): see hierarchy(rn)
  n1
  n2
Unknowns (2): see unknowns(rn)
  n1₊x(t) [defaults to k]
  n2₊x(t) [defaults to k]
Parameters (3): see parameters(rn)
  k [defaults to 2.0]
  n1₊k [defaults to 1.0]
  n2₊k [defaults to 2.0]

I don’t know why their are 3 parameters here… And

oprob = ODEProblem(complete(rn), [], (0.0, 10.0))

gives:

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: false
timespan: (0.0, 10.0)
u0: 2-element Vector{Float64}:
 2.0
 2.0

Maybe I did something wrong. If not, the problem may come from MTK because I could actually also reproduce the same thing:

@mtkmodel Decl1 begin
    @parameters begin
        k = 1.0
    end
    @variables begin
        x(t) = k
    end
end

@mtkmodel Decl2 begin
    @parameters begin
        k = 2.0
    end
    @variables begin
        x(t) = k
    end
end

@named m1 = Decl1()
@named m2 = Decl2()
@named sys1  = ODESystem([D(m2.x) ~ 0, D(m1.x) ~ 0], t)
prob1 = ODEProblem(complete(sys1))

gives:

ODEProblem with uType Vector{Float64} and tType Nothing. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: false
timespan: (nothing, nothing)
u0: 2-element Vector{Float64}:
 1.0
 1.0
5

Catalyst uses the same code to discover the unknowns and parameters as MTK (or at least it did prior to recent MTK changes). If you update your MTK example to

@named sys1  = ODESystem([D(m2.x) ~ 0, D(m1.x) ~ 0], t; systems = [m1, m2])

(which I think you need to do to tell it the subsystems), then you will again see that it adds an extra parameter

julia> @named sys1  = ODESystem([D(m2.x) ~ 0, D(m1.x) ~ 0], t; systems = [m1, m2])
Model sys1:
Subsystems (2): see hierarchy(sys1)
  m1
  m2
Equations (2):
  2 standard: see equations(sys1)
Unknowns (2): see unknowns(sys1)
  m2₊x(t) [defaults to k]
  m1₊x(t) [defaults to k]
Parameters (3): see parameters(sys1)
  k [defaults to 1.0]
  m1₊k [defaults to 1.0]
  m2₊k [defaults to 2.0]

I’d suggest you open an issue on the ModelingToolkit Github with your MTK example from above, but using

@named sys1  = ODESystem([D(m2.x) ~ 0, D(m1.x) ~ 0], t; systems = [m1, m2])
6

As a temporary workaround to the bug you can use the four argument constructor like:

@named rn = ReactionSystem([Reaction(1.0, [n1.x], [n2.x])], default_t(), [], []; systems=[n1, n2])

giving

Model rn:
Subsystems (2): see hierarchy(rn)
  n1
  n2
Unknowns (2): see unknowns(rn)
  n1₊x(t) [defaults to n1₊k]
  n2₊x(t) [defaults to n2₊k]
Parameters (2): see parameters(rn)
  n1₊k [defaults to 1.0]
  n2₊k [defaults to 2.0]

Then

julia> oprob = ODEProblem(complete(rn), [], (0.0, 10.0))
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: false
timespan: (0.0, 10.0)
u0: 2-element Vector{Float64}:
 1.0
 2.0
7

I’ve opened an issue on MTK here.

8

Thanks again! In the constructor

@named rn = ReactionSystem([Reaction(1.0, [n1.x], [n2.x])], default_t(), [], []; systems=[n1, n2])

according to the source code of Catalyst the two extra parameters ([] and []) stand for the lists of unknowns et parameters. For my understanding, how does it come that this solves the issue ?
Also, do you know any workaround that would work in “pure” MTK ?

9

If you use the four argument constructor, you disable automatic discovery of unknowns and parameters from the reactions/equations. Instead Catalyst (and I presume MTK) assume that you will tell it via the two arrays which parameters and unknowns should be owned by that system. Hence in this case, when passing [], the top level system assumes that it has no unknowns/parameters that it owns.