Rationalizing input model

[pogudingleb]

At the moment all the algorithms in the package require input models to have rational right-hand side (that is, the rhs must be a quotient of two polynomials). One way to lift this limitation is to transform the system into a rational or even polynomial one. It is know that this can be done by introducing new variables (detailed algorithm with implementation in Python).

However, there is an important subtlety. Consider an example:

x'(t) = a * e^x(t)
y(t) = x(t)

The standard idea would be to introduce a new variable z(t) := e^x(t), this will recast the system into

x'(t) = a * z(t)
z'(t) = a * z(t)^2
y(t) = x(t)

In the latter model, parameter a will be non identifiable. However, if we add a different new variable w'(t) := a * e^x(t), the system will become

x'(t) = w(t)
w'(t) = w(t)^2
y(t) = x(t)

Although a has disappeared from the system, knowing x(0) and w(0) is sufficient to find it. Since these initial conditions are identifiable (and they are), a is in fact identifiable as well.

The difference between the two transformations above is that the former increases the number of degrees of freedom in the system and the latter does not. Because of this, the former yield an incorrect result. Therefore, the overall goal would be to have a procedure for transforming the model into a rational one while keeping the number of degrees of freedom the same. For practical examples of "careful" transformation, see sections A.2 and A.3 in this paper.

One plausible approach is:

1. Implement any reasonable procedure for transforming a system into a rational one (e.g. adding all non rational terms as new variables until we are done). This may increase the number of degrees of freedom.
2. Reduce back the number of degrees of freedom by applying an algorithm for finding scaling transformation from this paper.

This approach will work nicely for the example above. Implementing any of the two steps above separately would be a nice and welcome contribution on its own.

UPD: A simpler procedure can be established for parameters only #190

[mirams]

Was hoping to use StructuralIdentifiability.jl to look at a simple ion channel model example, where some transition rates in a chemical reaction feature terms like exp(theta*V(t)) where V(t) is the time-dependent input/forcing function and theta is a model parameter. i.e.

ode = @ODEmodel(
    k1(t) = theta_1*exp(theta_2*V(t)),
    k2(t) = theta_3*exp(-theta_4*V(t)),

    a_inf(t) = k1(t) / (k1(t)+k2(t)),
    tau_a(t) = 1.0 / (k1(t)+k2(t)),
    
    a'(t) = (a_inf(t) - a(t)) / tau_a(t),

    Current(t) = Conductance * a(t) * (V(t)-NernstPotential),
)
results = assess_identifiability(ode)
display(results)

N.B. in this case sometimes V appears in an exponential, but on the last line it also appears linearly, meaning a simple substitution won't be sensible/possible?

At the moment I get exp listed as a parameter:

[ Info: Summary of the model:
[ Info: State variables: a
[ Info: Parameters: Conductance, NernstPotential, exp, theta_1, theta_2, theta_3, theta_4
[ Info: Inputs: V
[ Info: Outputs: k1, k2, a_inf, tau_a, Current

and the following error:

ERROR: Number of variables does not match number of values
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] (::Nemo.QQMPolyRingElem)(vals::Nemo.QQMPolyRingElem)
   @ Nemo ~/.julia/packages/Nemo/Dukvj/src/flint/fmpq_mpoly.jl:562
 [3] top-level scope
   @ ~/.julia/packages/StructuralIdentifiability/r7Tss/src/input_macro.jl:5

I guess even if this isn't supported at the moment it would be nice to pick up Julia functions like exp and throw a nice error rather than specify them as parameters.

Incidentally, it also leaves the julia REPL in a confusing state so that it can't evaluate exponentials any more

julia> exp(2)
ERROR: Number of variables does not match number of values
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] (::Nemo.QQMPolyRingElem)(vals::Int64)
   @ Nemo ~/.julia/packages/Nemo/Dukvj/src/flint/fmpq_mpoly.jl:557
 [3] top-level scope
   @ REPL[3]:1

(this is Julia 1.11.5 and StructuralIdentifiability.jl version 0.5.14)

[pogudingleb]

@mirams Thanks for an interesting model! I think one can rewrite the model in the format acceptable by StructuralIdentifiability.jl as follows:

1. As input we will take the derivative V'(t), denote it by Vp(t). Then V(t) can be regarded as a state satisfying V'(t) = Vp(t). In order to reflect the fact that V(t) is known, we add an extra output y(t) = V(t).
2. Then k1 satisfies k1'(t) = theta_2 * Vp(t) * k1(t). Note that this transformation does not change the number of degrees of freedom: we get an extra state variable k1 but get rid of a parameter theta_1 which can now be reconstructed from V(t), theta_2 and k1(t). Same story for k2.

Putting all these together, we get:

ode = @ODEmodel(
    k1'(t) = theta_2 * Vp(t) * k1(t),
    k2'(t) = theta_3 * Vp(t) * k2(t),
    V'(t) = Vp(t),
    a'(t) = (k1(t) / (k1(t)+k2(t)) - a(t)) * (k1(t) + k2(t)),
    y(t) = V(t),
    Current(t) = Conductance * a(t) * (V(t)-NernstPotential)
)

By running assess_identifiability(ode), we get that all the parameters and states are identifiable. Since V(t), theta_2 and k1(t) are identifiable, theta_1 also is. Same for theta_3.

Let me know if this makes sense.