StructuralIdentifiability.jl icon indicating copy to clipboard operation
StructuralIdentifiability.jl copied to clipboard

*very* long identifiability runtime for relatively small model

Open TorkelE opened this issue 8 months ago • 5 comments

I have the following model:

ode = @ODEmodel(
    X'(t) = p - X(t),
    Y'(t) = -Y(t) + ((K1^3)*v1) / (K1^3 + X(t)^3),
    Z'(t) = -Z(t) + (v2*(Y(t)^3)*(X(t)^3)) / ((K2^3 + X(t)^3)*(K3^3 + Y(t)^3)),
    y(t) = Z(t)
)

for which the runtime of assess_identifiability is incredibly long. I.e.

assess_identifiability(ode)

takes at least 24 hours (I have not actually managed to complete it).

The model is essentially a incoherent feedforward loop (X deactives Y and activets Z, Y activates Z). In Catalyst it can be implemented like

rn = @reaction_network begin
    (p,1.0), 0 <--> X
    hillr(X, v1, K1, 3), 0 --> Y
    v2*hill(X, 1.0, K2, 3)*hill(Y, 1.0, K3, 3), 0 --> Z
    1.0, (Y,Z) --> 0
end

TorkelE avatar Jun 29 '25 13:06 TorkelE

Hi @TorkelE , thanks a lot for reporting! I have checked briefly, this seemingly innocent model indeed triggers complicated computations (because of high nonlinearity appearing in denominators), so this is in some sense "correct behavior". But, of course, this is annoying (if not embarrasing) to take so long on such a small model, I will think about this, there should be a way to get this work.

pogudingleb avatar Jun 29 '25 14:06 pogudingleb

Thanks for the reply! Just happy to know that the sluggishness is a property of the model and nothing else.

Generally, if I have something which takes like a day, is there hope in putting it on an HPC and let run for a week, or at that point it is probably better to give up?

TorkelE avatar Jun 29 '25 14:06 TorkelE

Yes, in some cases passing to a machine with more memory and waiting longer may allow to get ultimately result. Looking at the internal debugging info in this example, I am skeptical about this particular case.

However! I have found a workaround (but not sure how to automatize it). From the first equation, we know that p - X(t) is proportional to $e^{-t}$. So we can take any known output $y_2(t) = b e^{-t}$. Then, for some (unknown) constant a, we will have y2(t) = a(p - X(t)). This will give us the following system (equivalent identifiability-wise to the original one):

ode = @ODEmodel(
    X'(t) = p - X(t),
    Y'(t) = -Y(t) + ((K1^3)*v1) / (K1^3 + X(t)^3),
    Z'(t) = -Z(t) + (v2*(Y(t)^3)*(X(t)^3)) / ((K2^3 + X(t)^3)*(K3^3 + Y(t)^3)),
    y(t) = Z(t),
    y2(t) = a * (p - X(t)),
)

After this rewriting, the identifiability assessment is very quick.

PS. This seem to be also related to first integrals as we essentially use that $(p - X(t)) e^t$ is a first integral of the model.

pogudingleb avatar Jun 29 '25 16:06 pogudingleb

Thanks a lot for that extra!

Works flawlessly, could compute identifiability for my model, and a couple of variants of it, without any problems.

TorkelE avatar Jun 30 '25 15:06 TorkelE

A bit longer explanation/justification for myself, in order not to forget.

  1. Augment the system with an observed state satisfying W'(t) = -W(t):
ode = @ODEmodel(
    X'(t) = p - X(t),
    Y'(t) = -Y(t) + ((K1^3)*v1) / (K1^3 + X(t)^3),
    Z'(t) = -Z(t) + (v2*(Y(t)^3)*(X(t)^3)) / ((K2^3 + X(t)^3)*(K3^3 + Y(t)^3)),
    y(t) = Z(t),
    W'(t) = -W(t),
    y2(t) = W(t),
)

The subsystem consisting of the last two equations is not coupled with the original, so this transformation does not change identifiability.

  1. Remark that the new system has a first integral C = W(t) / (p - X(t)). Use it to eliminate W from the original model:
ode = @ODEmodel(
    X'(t) = p - X(t),
    Y'(t) = -Y(t) + ((K1^3)*v1) / (K1^3 + X(t)^3),
    Z'(t) = -Z(t) + (v2*(Y(t)^3)*(X(t)^3)) / ((K2^3 + X(t)^3)*(K3^3 + Y(t)^3)),
    y(t) = Z(t),
    y2(t) = C * (p - X(t)),
)

So we get the extended model as described in the post above.

pogudingleb avatar Jul 08 '25 19:07 pogudingleb