Notebook

A differential example from Hong et al 2020¶

Let us consider the following differential system with rational dynamics (see example 6.4 in Hong et al):

$$\left\{\begin{array}{ll} \dot x_1 & = a_1(x_2-x_1) - \frac{k_aV_mx_1}{k_ck_a + k_cx_3 + k_ax_1}\\ \dot x_2 & = a_2(x_1-x_2)\\ \dot x_3 & = b_1(x_4 - x_3) - \frac{k_cV_mx_3}{k_ck_a + k_cx_3 + k_ax_1}\\ \dot x_4 & = b_2(x_3-x_4) \end{array}\right.$$

This is a model arising from pharmacokinetics (see Demignot et al) and the only output in the system is the fucntion $y = x_1$. Let us transform this problem so we can use dalgebra to solve the identifiability problem.

In [1]:

%display latex
from dalgebra import *
from timeoutcontext import timeout

The system has 4 state variables $x_1(t), x_2(t), x_3(t)$ and $x_4(t)$ and the differential system can be written using the parameters $k_a, k_c, V_m, a_1, a_2, b_1$ and $b_2$. We need to create these variables and set the differential ring such that everything is a constant:

In [2]:

R.<a_1,a_2,b_1,b_2,k_a,k_c,V_m> = QQ[]
DR = DifferentialRing(R, lambda p : 0)
## We update the variables
a_1,a_2,b_1,b_2,k_a,k_c,V_m = DR.gens()
## We check that all these are constants
print("All are constants -> ", all(el.derivative() == 0 for el in DR.gens()))

DR

All are constants ->  True

Out[2]:

$\displaystyle \newcommand{\Bold}[1]{\mathbf{#1}}\left(\Bold{Q}[a_{1}, a_{2}, b_{1}, b_{2}, k_{a}, k_{c}, V_{m}], 0\right)$

In order to build the system, we need to create the four differential variables $x_1(t), x_2(t), x_3(t)$ and $x_4(t)$. In the code, these will be represented by $\texttt{x1}$, $\texttt{x2}$, $\texttt{x3}$ and $\texttt{x4}$ such that $\texttt{xi[k]} = x_i^{(k)}(t)$. We do that with the class DifferentialPolynomialRing:

In [3]:

DPR.<x1,x2,x3,x4> = DifferentialPolynomialRing(DR)
DPR

Out[3]:

$\displaystyle \newcommand{\Bold}[1]{\mathbf{#1}}\left(\Bold{Q}[a_{1}, a_{2}, b_{1}, b_{2}, k_{a}, k_{c}, V_{m}], 0\right) \{ x1, x2, x3, x4 \}$

In [4]:

x1,x2,x3,x4

Out[4]:

$\displaystyle \left(\mathit{x1}_{\ast}, \mathit{x2}_{\ast}, \mathit{x3}_{\ast}, \mathit{x4}_{\ast}\right)$

At this point we can create a differential system with the appropiate equations. Since $x_1(t)$ is our output variable that is the variable we do not need to eliminate, so the variables for teh system are the remaining $x_2(t), x_3(t)$ and $x_4(t)$:

In [5]:

denom = k_c*k_a + k_c*x3[0] + k_a*x1[0]
system = DifferentialSystem([
    x1[1]*denom - a_1*(x2[0]-x1[0])*denom + k_a*V_m*x1[0],
    x2[1] - a_2*(x1[0] - x2[0]),
    x3[1]*denom - b_1*(x4[0]-x3[0])*denom + k_c*V_m*x3[0],
    x4[1] - b_2*(x3[0]-x4[0])
], DPR, [x2,x3,x4])
system

Out[5]:

\(\displaystyle \newcommand{\Bold}[1]{\mathbf{#1}}\text{System over } \left(\Bold{Q}[a_{1}, a_{2}, b_{1}, b_{2}, k_{a}, k_{c}, V_{m}], 0\right) \{ x1, x2, x3, x4 \} \text{ with variables } \mathit{x2}_{\ast}, \mathit{x3}_{\ast}, \mathit{x4}_{\ast} : \left\{\begin{array}{ll} k_{a} \mathit{x1}_{1} \mathit{x1}_{0} + k_{c} \mathit{x1}_{1} \mathit{x3}_{0} + \left(a_{1} k_{a}\right) \mathit{x1}_{0}^{2} + \left(-a_{1} k_{a}\right) \mathit{x1}_{0} \mathit{x2}_{0} + \left(a_{1} k_{c}\right) \mathit{x1}_{0} \mathit{x3}_{0} + \left(-a_{1} k_{c}\right) \mathit{x2}_{0} \mathit{x3}_{0} + \left(k_{a} k_{c}\right) \mathit{x1}_{1} + \left(a_{1} k_{a} k_{c} + k_{a} V_{m}\right) \mathit{x1}_{0} + \left(-a_{1} k_{a} k_{c}\right) \mathit{x2}_{0} & = 0 \\ \left(-a_{2}\right) \mathit{x1}_{0} + 1 \mathit{x2}_{1} + a_{2} \mathit{x2}_{0} & = 0 \\ k_{a} \mathit{x1}_{0} \mathit{x3}_{1} + \left(b_{1} k_{a}\right) \mathit{x1}_{0} \mathit{x3}_{0} + \left(-b_{1} k_{a}\right) \mathit{x1}_{0} \mathit{x4}_{0} + k_{c} \mathit{x3}_{1} \mathit{x3}_{0} + \left(b_{1} k_{c}\right) \mathit{x3}_{0}^{2} + \left(-b_{1} k_{c}\right) \mathit{x3}_{0} \mathit{x4}_{0} + \left(k_{a} k_{c}\right) \mathit{x3}_{1} + \left(b_{1} k_{a} k_{c} + k_{c} V_{m}\right) \mathit{x3}_{0} + \left(-b_{1} k_{a} k_{c}\right) \mathit{x4}_{0} & = 0 \\ \left(-b_{2}\right) \mathit{x3}_{0} + 1 \mathit{x4}_{1} + b_{2} \mathit{x4}_{0} & = 0 \\ \end{array}\right.\)

Let us try to execute blindly the differential resultant algorithm:

In [8]:

try:
    with timeout(3):
        res = system.diff_resultant(alg_res="iterative", verbose=True)
except TimeoutError:
    print("Timeout reached")

Checking if there is any linear variable...
Found linear variables [x4_*, x2_*]: we remove them first...
##########################################################
The system is linear: we use Macaulay resultant
Getting the appropriate extension for having a SP2 system...
Trying the extension (0, 0, 0)
Trying the extension (1, 0, 0)
Trying the extension (0, 1, 0)
Trying the extension (0, 0, 1)
Trying the extension (2, 0, 0)
Trying the extension (1, 1, 0)
Trying the extension (1, 0, 1)
Trying the extension (0, 2, 0)
Trying the extension (0, 1, 1)
Trying the extension (0, 0, 2)
Trying the extension (3, 0, 0)
Trying the extension (2, 1, 0)
Trying the extension (2, 0, 1)
Trying the extension (1, 2, 0)
Trying the extension (1, 1, 1)
Trying the extension (1, 0, 2)
Found the valid extension (1, 0, 2)
Getting the homogenize version of the equations...
Computing the Macaulay resultant...
Timeout reached

In [ ]: