Speeding things up via linear first integrals

[pogudingleb]

Main idea

Many models of interest (especially, in life sciences) have linear first integrals (for example, conservation of mass, etc.). If there is such a first integral a_1 * x_1 + ... + a_n * x_n, where x_1, ..., x_n are the states and a_1, ...., a_n are constants (maybe involving parameters), one can do the following:

1. introduce a new parameter C = a_1 * x_1 + ... + a_n * x_n
2. use x_1 = 1 / a_1 * (C - a_2 * x_2 - ... - a_n * x_n) to eliminate x_1 from the system

As a result, we reduce the number of states by one at the expense of having one more parameter. Since the typical bottleneck of the algorithm is differential elimination, this will very likely speed up the computation.

Motivating example

As a motivating/illustrating example, we consider this CRN model but with only one output:

ode = @ODEmodel(
    x1'(t) = -k1 * x1(t) * x2(t) + k2 * x4(t) + k4 * x6(t),
    x2'(t) = -k1 * x1(t) * x2(t) + k2 * x4(t) + k3 * x4(t),
    x3'(t) = k3 * x4(t) + k5 * x6(t) - k6 * x3(t) * x5(t),
    x4'(t) = k1 * x1(t) * x2(t) - k2 * x4(t) - k3 * x4(t),
    x5'(t) = k4 * x6(t) + k5 * x6(t) - k6 * x3(t) * x5(t),
    x6'(t) = -k4 * x6(t) - k5 * x6(t) + k6 * x3(t) * x5(t),
    y1(t) = x3(t)
)

This model is quite hard, it does not finish in an hour (and I remember trying it earlier, I think it just took all my memory in a couple of hours).

One can observe that x5(t) + x6(t) is constant. Thus, we can introduce C1 = x5(t) + x6(t) and eliminate x5(t):

ode = @ODEmodel(
    x1'(t) = -k1 * x1(t) * x2(t) + k2 * x4(t) + k4 * x6(t),
    x2'(t) = -k1 * x1(t) * x2(t) + k2 * x4(t) + k3 * x4(t),
    x3'(t) = k3 * x4(t) + k5 * x6(t) - k6 * x3(t) * (C1 - x6(t)),
    x4'(t) = k1 * x1(t) * x2(t) - k2 * x4(t) - k3 * x4(t),
    x6'(t) = -k4 * x6(t) - k5 * x6(t) + k6 * x3(t) * (C1 - x6(t)),
    y1(t) = x3(t)
)

Doing the same with first integrals C2 = x2(t) + x4(t) and C3 = x1(t) - x2(t) + x3(t) + x6(t), we arrive at the following model with only three states:

ode = @ODEmodel(
    x1'(t) = -k1 * x1(t) * x2(t) + k2 * (C2 - x2(t)) + k4 * x6(t),
    x2'(t) = -k1 * x1(t) * x2(t) + k2 * (C2 - x2(t)) + k3 * (C2 - x2(t)),
    x6'(t) = -k4 * x6(t) - k5 * x6(t) + k6 * (C3 - x1(t) + x2(t) - x6(t)) * (C1 - x6(t)),
    y1(t) = (C3 - x1(t) + x2(t) - x6(t))
)

Now global identifiability is successfully analyzed in 60 seconds (and computing the input-output equation takes only 15 seconds)

What to do

The substantial speedup above can be achieved automatically if we can

1. find all linear first integrals
2. perform a reduction using this integrals
3. add this as a preprocessing in assess_global_identifiability

[iliailmer]

If I understand correctly, would this require solving a system Sum(A[i] Diff(x[i], t^j))=0 to find the coefficients A[i]?

[pogudingleb]

What is Diff(x[i], t^j) in your notation?

[iliailmer]

What is Diff(x[i], t^j) in your notation?

j-th Lie derivative of x[i], the i-th state variable.

[pogudingleb]

I would not think that the Lie derivatives are necessary. Basically, we are talking about "one" equation sum( A[i] * Diff(x[i], t)) = 0 but since A[i] must depend only on parameters but not states or inputs, this equation will yield several equations (for the polynomial case via collecting with respect to the degree in x). Does this make sense?

[iliailmer]

Not entirely, sorry. I understand the starting idea: we have one linear combination for which we seek coefficients A[i]. We would need to differentiate w.r.t. time several times to obtain a square linear system for A[i], right? I guess I'm confused at this part:

or the polynomial case via collecting with respect to the degree in x

[pogudingleb]

Well, I think one could approach this via differentiation but this is very subtle (not to mention computationally demanding): the coefficients of the system will involve parameters and states while the solution will be sought in the parameters only.

Maybe discussing this on an example would be a better way:

x1'(t) = x1(t) + a * x2(t)
x2'(t) = x2(t) + x3(t)
x3'(t) = b * x3(t) + b/a * x1(t)

We would like to find A1, A2, A3 such that A1 * x1'(t) + A2 * x2'(t) + A3 * x3'(t) = 0, we can write this as follows:

A1 * (x1(t) + a * x2(t)) + A2 * (x2(t) + x3(t)) + A3 * (b * x3(t) + b/a * x1(t)) = 
= x1(t) * (A1 + A3 * b/a) + x2(t) * (A1 * a + A2) + x3(t) * (A2 + A3 * b) = 0

Since A1, A2, A3 are rational functions in the parameters a, b, the single equality above yields three equations by equating to zero each coefficient of the above expression considered as a polynomial in x1, x2, x3, that is:

A1 + A3 * b/a = 0
A1 * a + A2 = 0
A2 + A3 * b = 0

This is what I meant by "collecting with respect to the degree in x". Does this make more sense now?

[pogudingleb]

In progress with @StefanVaylBX2023