Including initial conditions in equations

[TorkelE]

Is it possible to include initial conditions in the equations (in some sense, these are essentially parameters, right?). E.g something like:

si_ode = @ODEmodel(
    X1'(t) = p - k1*X1(t) + X1(0),
    X2'(t) = X1(t) -  X2(t),
    y1(t) = X2(t)
)

The reason I specifically am wondering about this is that, when reaction system have conservation laws (such as X1 + X2 = C), this can be used to eliminate a variable, e.g. X2. However, this introduces a new parameter C, which can be expressed X1(0) + X2(0) = C. Also, if X2 (which is eliminated) happens to be a measured quantity, this would need to be replaced with something like X1(0) + X2(0) - X1(t).

[pogudingleb]

But you can just express C as C = X1(t) + X2(t) and express everywhere X2 with C - X1. Am I missing something?

[TorkelE]

Yes, ish. The problem is that then in the output you get identifiability for C, and not for X1(0) and X2(0). Although I guess I can do a transformation on the output as well to put in identifiability for X2(0) given what I know of C and `X1(0). However, StructuralIdentfiabiltiy does not given the output in a pure symbol => symbol form, but the keys in the disc is a bit weird, so I am not sure how to generate these.

[pogudingleb]

I think that one can use the funcs_to_check argument which allows to provide custom functions to evaluate identifiability for. Say, you had a system

ode = @ODEmodel(
    x1'(t) = x1(t)^2 + x2(t)^2,
    x2'(t) = -x1(t)^2 - x2(t)^2,
    y(t) = x1(t) - x2(t)
)

Here you can set C = x1(t) + x2(t) and transform

ode2 = @ODEmodel(
    x1'(t) = x1(t)^2 + (c - x1(t))^2,
    y(t) = 2 * x1(t) - c
)

If you wanted to assess identifiability of x1, x2, you now call assess_identifiabilit(ode2, funcs_to_check = [x1, c - x1]).

[TorkelE]

That's a good point, would it be possible to name the output of these things (so that it looks like c - x1 is x2, here my c is a parameter that the user has never seen and would not necessarily know what it is about)? Unfortunatley, and this is a side thing, I was not able to get funcs_to_check to work for my ModelingToolkit type input.

[pogudingleb]

Indeed, there should be a code for renaming. Does the function you have for reducing wrt first integrals return you the correspondence between old and new variables?

Sorry to hear about the issue with funds_to_check! Could you tell me what is the problem you have? The following code works fine on my computer

using StructuralIdentifiability, ModelingToolkit

@parameters r1, r2, c1, c2, beta1, beta2, chi1, chi2
@variables t, x1(t), x2(t), y(t), u(t)

D = Differential(t)
eqs = [
    D(x1) ~ r1 * x1 * (1 - c1 * x1) + beta1 * x1 * x2 / (chi1 + x2) + u,
    D(x2) ~ r2 * x2 * (1 - c2 * x2) + beta2 * x1 * x2 / (chi2 + x1),
]

measured_quantities = [y ~ x1]

ode_mtk = ODESystem(eqs, t, name = :mutualist)
assess_identifiability(ode_mtk, measured_quantities = measured_quantities, funcs_to_check = [c1 + c2])

[TorkelE]

I think it might also be related to going through Catalyst. I've forgotten the exact detailed but remembered that I looked at options and decided supporting it was not going to be with the coding mess it'd produce. I will have a look again and message you if I can narrow it down properly.

[TorkelE]

So, I think I remember now. Basically, internally we run a StructuralIdentifiability.preprocess_ode(ode_mtk, measured_quantities = measured_quantities), and then use that as input to the assess_identifiability function. Here, I can only give funcs_to_check information in MTK form, but that doesn't match when the input ODE is in SI form.

The obvious answer is to not do the preprocess_ode conversion but instead use the MTK system. I think I had some good reason for this as well (which currently eludes me, but when tell you when I recall).

[isaacsas]

What is the reason to not just use the MTK system?

[TorkelE]

See above. I meant to implement it that was, however, there were some issues and I decided it would be a really big mess. That was last week, forgotten what it was. Might make a try to move back to MTK (which would be preferably, but also depends on MTK keeping up to date with SI).

[isaacsas]

I don't think you ever stated how funcs_to_check was failing going through Catalyst.

Can you point me to the comment that states why you can't go via MTK?

[TorkelE]

Sorry, not sure I fully understand.

The input to funcs_to_check has to be of the same type as the model. Given that we use SI-type odes, funcs_to_check must have SI-type expressions (which we cannot generate, hence we cannot support it).

If we used MTK-type systems internally, funcs_to_check could take MTK-style inputs, which we could generate, and we could support it. Ideally I'd want to support funcs_to_check, but at the time I decided it was no worth it/possible since we then would have to use MTK internally. If we can use MTK internally, I would want to do this.

[TorkelE]

Earlier in this conversation (just above your reply):

[isaacsas]

You said:

If we used MTK-type systems internally, funcs_to_check could take MTK-style inputs, which we could generate, and we could support it. Ideally I'd want to support funcs_to_check, but at the time I decided it was no worth it/possible since we then would have to use MTK internally. If we can use MTK internally, I would want to do this.

and

The obvious answer is to not do the preprocess_ode conversion but instead use the MTK system. I think I had some good reason for this as well (which currently eludes me, but when tell you when I recall).

I'm still missing why we can't use the MTK internally.

[TorkelE]

"The obvious answer is to not do the preprocess_ode conversion but instead use the MTK system. I think I had some good reason for this as well (which currently eludes me, but when tell you when I recall)."

As in:

• Using MTK internally seems like a good idea.
• I looked at implementing it, but it had implications which made it not work out well, so I didn't.
• I don't remember what these reasons were, but when I do, I will elaborate.

Sorry, I realise this was not directly obvious from my previous comment.

I am working on again trying to use MTK internally, hoping that I will remember the problems when/if I encounter them. This PR depends on https://github.com/SciML/Catalyst.jl/pull/713 in any case, and we won't have to worry about merging this one until that one gets merged.

[isaacsas]

No worries, I was just trying to understand what the limitation in using MTK is.

As for the other PR, we need to settle the issue I raised about where to do (optional) expansion. I still prefer doing it as a system transformation at the ReactionSystem level rather than trying to support doing it within every system conversion function.

[TorkelE]

Yeah, me and Gleb had discussions at like 3 different places, didn't realise first that it was pretty impossible to follow that from an outsider.

For the other point, my general feeling is that I fully agree. I think that PR is better than what we currently got, but if you'd rather go directly to something else, I am happy to skip it. Since that issue is not a big hurry, I think next time we meet we just quickly go through so that we are on the same page, I am certain I know what you mean, and I know how we want it implemented. Then I can make a new PR to implement that.