Implement equations in DSL planning

[TorkelE]

Ok, I hope to start on this soon, just wanted to go through a few things before I begin implementation.

Currently, I plan to implement two different types of equations (over two separate PRs):

• Observables.
• Differential equations. Would any other cases make sense?

Observables

I think these are the easiest. We would add an option @observable to the DSL, and the syntax be something like:

@reaction_network begin
    @observable Xtot = X1 + X2
   ...
end

Then, X (calculated according to the rule) would be made an observable of the system. Several observables could be computed in a block:

@reaction_network begin
    @observables begin
        Xtot = X1 + X2
        Ytot = Y1 + Y2 + Y3
    end
   ...
end

Considerations

• The current plan is to store these in the observed field. This has some other uses in MTK when it performs structural_simplify. I would need to check that this is indeed ok.
• The observables should not create equations in between each other (and each should only depend on the equations). It might be possible though to allow some simple inter-dependency like

    @observables begin
        Ptot = P1 + P2
        prot_tot = Ptot + Q
    end

• Observables could depend on species, parameters, and variables. What if an unknown symbol is encountered in the expression? We could either default to creating it as a parameter, or throwing an error. I think I would prefer the later, but for differential equations I would prefer the latter, and consistency is good. Generally I feel stronger about wanting to default to parameters for differential equations (than wanting to throw errors for observables).

Differential equations

The idea is to permit something like this:

@reaction_network begin
    @equation D(V) ~ r*G # Here r could be a parameter and G a species.
   ...
end

Again, blocks of several equations could be created using

    @equations begin
        D(V) ~ r*G # 
        D(N) ~ r*e*G 
    end

Considerations

• The left-hand side would be forced to be a single variable and its differential.
• The left-hand side would be automatically be detected and created as a variable and added to the system (no @variables V(t) needed).
• It would automatically be assumed to be time-dependant (we could force something like D(V(t)) ~ r*G to make this explicit).
• For the right-hand side expression, I would want to default to making undeclared component parameters (rather than requiring the user to declare these).
• Currently, no @equation or @equations macros actually exist. So in actual code we would only run D(V) ~ r*G internally. Technically this mean we could use the term @diff_equation or @diff_eq or something. This would make it explicit that we are expecting a differential equation.
• Here, the D variable is reserved for a differential. It is possible to (internally) replace it with the full word Differential in the code (mutating D(V) ~ r*G to Differential(V) ~ r*G before it is run. That way, we could still permit D to be a species/parameter/variable.
• Technically, we could use syntax like V' ~ r*G instead (or even V' = r*G). While I would have preferred that, keeping as similar to modeling toolkit as possible does make sense.

[ChrisRackauckas]

Is there a reason to exclude algebraic equations?

[TorkelE]

Not really. Just that I have little experience with them, and afraid of what would happen if you defined something like

@reaction_network begin
    @eqaution X = Y + 12
     (p,d), 0 <--> X
     (p,d), 0 <--> Y
end

which is non-solvable.

As long as stuff like this don't happen I guess it would be fine? However, I figured most (but not all) algebraic equations people would want to use would be covered by the observables system.

[sebapersson]

Observables could depend on species, parameters, and variables. What if an unknown symbol is encountered in the expression? We could either default to creating it as a parameter, or throwing an error. I think I would prefer the later, but for differential equations I would prefer the latter, and consistency is good. Generally I feel stronger about wanting to default to parameters for differential equations (than wanting to throw errors for observables).

Some input on this, I think throwing an error might be nice. If the observable is used in parameter estimation, being able to detect non-dynamic parameters that only appear in the observable is beneficial (which I do not know is possible if you add them as a parameter), because for non-dynamic parameters we can efficiently compute the gradient given an ODE solution.

[TorkelE]

In this case, it might be possible to add this kind of information as metadata to these parameters?

[sebapersson]

As metadata in catalyst - or downstream packages?

[TorkelE]

ModelingToolkit (and in extent, Catalyst) variables and parameters have a metadata field where additional information on them can be stored (accessed through e.g. my_param.val.metadata, although each metadata typically has its own getters). We could store stuff like this (about a parameter) here.

Another possibility is that stuff like priors on a parameter's value (currently declared when a PEtabParameter is created) could be declared as parameter metadata (if desired) already when the Catalyst ReactionSystem is created. PEtab could then check this, and if it exist, use it.

[isaacsas]

Generally our goal should be to maintain consistency with the conventions of MTK as much as possible. This makes it easier for Catalyst users to use MTK and vice-versa.

Along those lines,

1. Observables should just be equations of the form A ~ B + C and declared as such. I have to double check, but I believe this is how they are stored (i.e. they are assumed to define the value of something on the left side and are stored as MTK equations), and will presumably be what MTK would use if/when they add observables into their DSL.
2. Outside of reactions, for which we've already established our conventions, for all other new blocks we should require pre-declaring any parameters, species, or variables that are used in an expression. So this would apply to observables, equations, events, compounds, etc. If the block is defining something, such as an observable A via A ~ B + C then we could allow A to be defined via the block, however, B and C should need an explicit pre-declaration. Note that some care will then have to be taken to figure out if A is a species or variable. (edit: Maybe we should just make all observables a variable, they wouldn't be included in the species list anyways.)
3. By default, the derivative operator D(A) should be with respect to DEFAULT_IV, however, if a user manually specifies another default iv as part of the DSL it should then be with respect to this iv. Likewise, we should ultimately give the ability for a user to define and use different derivative operators to support PDEs (but this can be a followup). We do currently support multiple ivs in Catalyst to enable conversion to PDESystems for MethodOfLines, so we don't want to introduce too many non-interoperable features in adding these blocks.
4. @equations should allow algebraic equations too. I believe ReactionSystem supports this, but if it doesn't we should fix that up (it should be straightforward since we just forward them to the MTK system type). The notation for defining them should just be what MTK would use. D(A) ~ -k*A for example or 0 ~ -k*A.

[isaacsas]

5. We should allow higher derivatives and any other equation types MTK allows, so really we just want to assume users have typed a valid MTK equation, and then pass it into the ReactionSystem constructor as part of the equation and reaction list.

[isaacsas]

6. We should reuse as much as possible the @equations and other parsing code that has been put into MTK. This will make it easier for Catalyst to pickup any new functionality added to the MTK DSL and simplify what we have to maintain.

[TorkelE]

now works on master