Structural Identifiability for Interpolation

[MaAl13]

Hello, i have the problem that i have a interpolation function on the righthand-side of one of my ODE terms. In order to have an analytical expression for it i used Lagrange polynomials. However, it seem that the StructuralIdentifiability doesn't really work in this case. Here is some example code:

using DifferentialEquations
using StructuralIdentifiability
using ModelingToolkit

days = [0, 2, 4, 6, 8, 10, 12, 14]
surface = [0, 2, 4, 6, 8, 10, 12, 14] 
diam = [0, 2, 4, 6, 8, 10, 12, 14].*2  

@variables t IB(t)  # Define the variables
@parameters k gamma thres pow  # Define the parameters
D = Differential(t)

ode = @ODEmodel(
    IB'(t) = k * sum(
        (
            begin
                (
                    prod((t .- days[1:end .!= i])) / 
                    prod((days[i] .- days[1:end .!= i]))
                ) * diam[i]
            end
        ) ^ pow / (
            begin
                (
                    prod((t .- days[1:end .!= i])) / 
                    prod((days[i] .- days[1:end .!= i]))
                ) * diam[i]
            end
        ) ^ pow + thres ^ pow
        * (
            begin
                (
                    prod((t .- days[1:end .!= i])) / 
                    prod((days[i] .- days[1:end .!= i]))
                ) * surface[i]
            end
        )
        for i in 1:length(days)
    ) - gamma * IB,
    y(t) = IB(t)
)

measured = [IB]

# Assess structural identifiability
identifiability_result = assess_identifiability(ode, measured_quantities=measured)

[ChrisRackauckas]

Yes structural identifiability can only be performed on a reduced set of models (rational functions).

[pogudingleb]

The problem is that ODEmodel macros does not really evaluate your Lagrange polynomials. I would suggest the following workarounds:

• compute the polynomials separately and put them in explicit form in the macros
• (perhaps less tedious) use ModelingToolkit interface (as described here), then you should be able to write such expressions.

Let me know if neither of these helps.

[MaAl13]

Hello, i tried both approaches, but neither of them seems to work unfortunately. Do you maybe have another idea ho to fix it? This is my attempt at the moment:

using DifferentialEquations
using StructuralIdentifiability
using ModelingToolkit
using Symbolics

# Data points
days = [0, 2, 4, 6, 8, 10, 12, 14]
diam = [0, 2, 4, 6, 8, 10, 12, 14].*2
surface = [0, 2, 4, 6, 8, 10, 12, 14] 

# Compute the Lagrange polynomial and return as a symbolic expression
function lagrange_poly(t, x_points, y_points)
    n = length(x_points)
    sum_expr = 0
    for i in 1:n
        term = y_points[i]
        for j in 1:n
            if j != i
                term *= (t - x_points[j]) / (x_points[i] - x_points[j])
            end
        end
        sum_expr += term
    end
    return sum_expr
end

# Define the ODE using the Lagrange polynomial
@variables t IB(t)
@parameters k gamma thres pow
D = Differential(t)

interp_expr = lagrange_poly(t, days, diam)  # Use the computed polynomial here

# Ensure the polynomial expression is compatible with ModelingToolkit
interp_expr_sym = Symbolics.build_function(interp_expr, t)
interp_expr_fn = eval(interp_expr_sym[1])

ode = @ODEmodel(
    IB'(t) = k * (interp_expr_fn(t)^pow / (interp_expr_fn(t)^pow + thres^pow)) - gamma * IB,
    y(t) = IB(t)
)

# Rest of the code for problem setup and analysis
measured = [IB]
identifiability_result = assess_identifiability(ode, measured_quantities=measured)

[pogudingleb]

Hi @MaAl13 ,

I would change several things in your script. First of all, you are using Symbolics and then pass the expression to the ODEmodel macros while these are in fact two different ways of defining the model, see the docs.

Second, some rewriting is needed:

• it is not possible to have t directly in the equations. You can though introduce a function T(t) with T' = 1.
• your equations are non-rational due to pow in the exponents. I would suggest the following variable transformation to rationalise the model: we introduce new state F(t) := 1 / (1 + interp_expr_fn(t)^pow / thres^pow). This new state satisfies F'(t) = pow * F(t) * (1 - F(t)) / interp_expr_fn(t) * interp_expr_fn'(t) and the equation for IB becomes IB'(t) = k * (1 - F(t)) - gamma * IB(t).

After these transformations (if I did not make any mistakes), the script becomes

using DifferentialEquations
using StructuralIdentifiability
using ModelingToolkit
using Symbolics

# Data points
days = [0, 2, 4, 6, 8, 10, 12, 14]
diam = [0, 2, 4, 6, 8, 10, 12, 14].*2
surface = [0, 2, 4, 6, 8, 10, 12, 14] 

# Compute the Lagrange polynomial and return as a symbolic expression
function lagrange_poly(t, x_points, y_points)
    n = length(x_points)
    sum_expr = 0
    for i in 1:n
        term = y_points[i]
        for j in 1:n
            if j != i
                term *= (t - x_points[j]) / (x_points[i] - x_points[j])
            end
        end
        sum_expr += term
    end
    return sum_expr
end

# Define the ODE using the Lagrange polynomial
@variables t IB(t) T(t) F(t)
@parameters k gamma thres pow
D = Differential(t)

interp_expr = lagrange_poly(t, days, diam)  # Use the computed polynomial here

tmp = substitute(expand_derivatives(D(interp_expr)) / interp_expr, t => T)

eqs = [D(IB) ~ k * (1 - F) - gamma * IB, D(F) ~ pow * F * (1 - F) * tmp, D(T) ~ 1]
ode = ODESystem(eqs, t, name = :name)

assess_identifiability(ode, measured_quantities = [IB, T])

The results is that all the parameters and states are globally identifiable. Since thres can be expressed via F(t) and pow, it is identifiable as well. Let me know if you have any further questions.

[MaAl13]

Hey @pogudingleb,

thank you so much for your incredibly helpful response! Absolutely awesome that you made it work here. I am just trying to run your script with the following package specifications:

[0c46a032] DifferentialEquations v7.11.0 [961ee093] ModelingToolkit v8.73.1 [91a5bcdd] Plots v1.39.0 [220ca800] StructuralIdentifiability v0.5.1 [0c5d862f] Symbolics v5.10.0

But weirdly i get the following error message:

ERROR: LoadError: MethodError: no method matching getindex(::Expr, ::Int64) Stacktrace: [1] top-level scope @ ~/Documents/Identifiability.jl:36

This is weird, because it seems that it is running through in your case.

[pogudingleb]

@MaAl13 This isstrange, I have almost the same versions, it works fine for me. Your stacktrace seem to refer to getindex at the line eqs = ... which I do not see there. Could it be a typo while copying? Like some extra square brackets somewhere instead of parenthesis?

[pogudingleb]

I am closing the issue. @MaAl13 feel free to reopen if you will have any further questions.

[MaAl13]

It was indeed a typo i think. It now runs for me :) This is absolutely genius, thank you very much. Also for your detailed explanation in the answers!

[MaAl13]

@pogudingleb Do you also know a clever way to write dF/dt = k/(1+(G(t)/a)^b+(H(t)/c)^d) in rational form? I so far only managed to do it by including 3 ODEs instead of the original one: dF/dt = -(bA(t)G'(t)+dB(t)H'(t))F^2(t) dA/dt = (b-1)G'(t)A(t)G(t)/a dB/dt = (d-1)H'(t)B(t)*H(t)/c

Also if i have instead of the Lagrange polynomial only a fraction of two polynomials that is given and only depends on x, not on x_points and y_points, would you include it directly in the equation, or do the same as with the Lagrange polynomial?

[pogudingleb]

@pogudingleb Do you also know a clever way to write dF/dt = k/(1+(G(t)/a)^b+(H(t)/c)^d) in rational form? I so far only managed to do it by including 3 ODEs instead of the original one: dF/dt = -(bA(t)_G'(t)+d_B(t)H'(t))F^2(t) dA/dt = (b-1)G'(t)A(t)G(t)/a dB/dt = (d-1)H'(t)B(t)*H(t)/c

3 ODEs sounds a fair number in this case but I am not sure I understand the particular transformation you have performed. What are your A and B? I would take A(t) = (G(t)/a)^b and B(t) = (H(t) / c)^d. Note that in this case a and b should disappear from the equations. And this is how it should work: you "trade" two parameters for two new initial conditions A(0) and B(0), so the total dimension of the model stays the same, and you will be able to transfer identifiability results back and forth.

Also if i have instead of the Lagrange polynomial only a fraction of two polynomials that is given and only depends on x, not on x_points and y_points, would you include it directly in the equation, or do the same as with the Lagrange polynomial?

I would just include it explicitly.

[MaAl13]

Alright, thanks. So if i take A(t) = (G(t)/a)^b and B(t) = (H(t) / c)^d i get for the original equation k/(1+A(t)+B(t)) which is then rational you mean? Also just to be sure, since i am loosing the two parameters in my new system i get two new initial conditions for the two new equations. Do they have then a one to one correspondance to the two parameters and if so, which initial condition the corressponds to which parameter, is there a general way to find this out?

Also two additional question:

• If i have a summation where the number of summands are changing during the course of time how would i include that in the identifiability. For example if i have dA/dt = ksum_i=1^N Polynom(t)/Polynom(t)Vol_i(t)-gamma*A. So the N changes over time aswell as the volume

• If i know the initial conditions from some components how can i adress this? Especially in the case where i know the Ics from the original system but not from the additional system i get from transforming the original system to have rational functions on the right hand side. Then i would know some Ics and some i want to assess identifiability because this translates to the parameters

[pogudingleb]

There is indeed a some sort of correspondence here: for example, a can be found if b and A(0) are known from the formula A(0) = (G(0) / a)^b, so identifiability of b and A would imply identifiability of a.

For the summation, I am not sure this can be incorporated in such generality. You can do computation for some fixed N or, if you can get a closed formula for the sum, use the closed formula.

For known initial conditions, I have an experimental branch in which you can specify some initial conditions to be assumed to be know under the assumption that their values are generic. Here is the relevant function.

[MaAl13]

Okay, i see. So i need to solve the equations with the initial conditions for the unknown parameters.

Thank you for providing me with the experimental branch for the initial conditions! I would the put just and array like [A, B, C] in the known_ic argument if A, B and C were states? Also can is still provide then the measured quantities?

Also i get from the identifiability the identifiability of the states back, can you maybe tell me how i can interpret this? Does this mean that i will be able to recover the true trajectory given that we don't have any noise in the data? But that would also depend on the number of observations right? Or does this only tell me that the initial conditions are identifiable?

[MaAl13]

Also i have massive problems to run the identifiability for the following code

using DifferentialEquations
using StructuralIdentifiability
using ModelingToolkit
using Symbolics

@variables t X1(t) X2(t) X3(t) X4(t) X5(t) X6(t) X7(t) X8(t) X9(t) X10(t) X11(t) X12(t) X13(t) X14(t) 
@parameters P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13 P14 P15 P16 P17 P18 P19 P20 
D = Differential(t)

eqs = [D(X1) ~ P1 * X14 - P2 * X1,
       D(X2) ~ P3 + P4 * X3 - P5 * X2,
       D(X4) ~ P6 + P7 * X5 * π *  X10 * X10 - P8 * X4,
       D(X6) ~ P9 * 1/(1+X7 + X8) + P10 * X9 - P11 * X6,
       D(X10) ~ P12 * X11 + P13 * X12 + P14 * X13,
       D(X3) ~ P6 * (P6 + P7 * X5 * π *  X10 * X10 - P8 * X4)/X4 * (1 - X3) * X3,
       D(X5) ~ P16 * (P12 * X11 + P13 * X12 + P14 * X13)/X10 * (1 - X5) * X5,
       D(X7) ~ P17 * X7 * (P1 * X14 - P2 * X1)/X1,
       D(X8) ~ P18 * X8 * (P6 + P7 * X5 * π *  X10 * X10 - P8 * X4)/X4,
       D(X9) ~ P19 * X9 * (1 - X9) * (P6 + P7 * X5 * π *  X10 * X10 - P8 * X4)/X4,
       D(X11) ~ P20 * X11 * (1- X11) * (P9 * 1/(1+X7 + X8) + P10 * X9 - P11 * X6)/X6,
       D(X12) ~ P15 * X12 * (1- X12) * (P3 + P4 * X3 - P5 * X2)/X2,
       D(X13) ~ P16 * X13 * (1- X13) * (P12 * X11 + P13 * X12 + P14 * X13)/X10]

ode = ODESystem(eqs, t, name = :ANONYMOUS)

identifiability_result = assess_identifiability(ode, measured_quantities = [X4, X6, X2, X10, X14], p=0.999)

# Print out the results
println("Identifiability analysis results:")
println(identifiability_result)

The local identifiability works, but the global one seems to blow up like crzy the memory and doesn't produce a result. Is this even possible in this case?