{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# A differential example from Hong et al 2020\n", "\n", "Let us consider the following differential system with rational dynamics (see example 6.4 in [Hong et al](https://arxiv.org/abs/1801.08112)):\n", "\n", "$$\\left\\{\\begin{array}{ll}\n", " \\dot x_1 & = a_1(x_2-x_1) - \\frac{k_aV_mx_1}{k_ck_a + k_cx_3 + k_ax_1}\\\\\n", " \\dot x_2 & = a_2(x_1-x_2)\\\\\n", " \\dot x_3 & = b_1(x_4 - x_3) - \\frac{k_cV_mx_3}{k_ck_a + k_cx_3 + k_ax_1}\\\\\n", " \\dot x_4 & = b_2(x_3-x_4)\n", "\\end{array}\\right.$$\n", "\n", "This is a model arising from pharmacokinetics (see [Demignot et al](https://pubmed.ncbi.nlm.nih.gov/2855567/)) 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." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%display latex\n", "from dalgebra import *\n", "from timeoutcontext import timeout" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "All are constants -> True\n" ] }, { "data": { "text/html": [ "\$\\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)\$" ], "text/latex": [ "$\\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)$" ], "text/plain": [ "Differential Ring [[Multivariate Polynomial Ring in a_1, a_2, b_1, b_2, k_a, k_c, V_m over Rational Field], (0,)]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "R. = QQ[]\n", "DR = DifferentialRing(R, lambda p : 0)\n", "## We update the variables\n", "a_1,a_2,b_1,b_2,k_a,k_c,V_m = DR.gens()\n", "## We check that all these are constants\n", "print(\"All are constants -> \", all(el.derivative() == 0 for el in DR.gens()))\n", "\n", "DR" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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`:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\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 \\}\$" ], "text/latex": [ "$\\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 \\}$" ], "text/plain": [ "Ring of operator polynomials in (x1, x2, x3, x4) over Differential Ring [[Multivariate Polynomial Ring in a_1, a_2, b_1, b_2, k_a, k_c, V_m over Rational Field], (0,)]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "DPR. = DifferentialPolynomialRing(DR)\n", "DPR" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\left(\\mathit{x1}_{\\ast}, \\mathit{x2}_{\\ast}, \\mathit{x3}_{\\ast}, \\mathit{x4}_{\\ast}\\right)\$" ], "text/latex": [ "$\\displaystyle \\left(\\mathit{x1}_{\\ast}, \\mathit{x2}_{\\ast}, \\mathit{x3}_{\\ast}, \\mathit{x4}_{\\ast}\\right)$" ], "text/plain": [ "(x1_*, x2_*, x3_*, x4_*)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x1,x2,x3,x4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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)$:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\\(\\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} :\n", "\n", " \\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 \\\\\n", "\\left(-a_{2}\\right) \\mathit{x1}_{0} + 1 \\mathit{x2}_{1} + a_{2} \\mathit{x2}_{0} & = 0 \\\\\n", "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 \\\\\n", "\\left(-b_{2}\\right) \\mathit{x3}_{0} + 1 \\mathit{x4}_{1} + b_{2} \\mathit{x4}_{0} & = 0 \\\\ \n", "\\end{array}\\right.\\)" ], "text/latex": [ "$\\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} :\n", "\n", " \\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 \\\\\n", "\\left(-a_{2}\\right) \\mathit{x1}_{0} + 1 \\mathit{x2}_{1} + a_{2} \\mathit{x2}_{0} & = 0 \\\\\n", "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 \\\\\n", "\\left(-b_{2}\\right) \\mathit{x3}_{0} + 1 \\mathit{x4}_{1} + b_{2} \\mathit{x4}_{0} & = 0 \\\\ \n", "\\end{array}\\right.$" ], "text/plain": [ "System over [Ring of operator polynomials in (x1, x2, x3, x4) over Differential Ring [[Multivariate Polynomial Ring in a_1, a_2, b_1, b_2, k_a, k_c, V_m over Rational Field], (0,)]] with variables [(x2_*, x3_*, x4_*)]:\n", "{\n", "\tk_a*x1_1*x1_0 + k_c*x1_1*x3_0 + a_1*k_a*x1_0^2 + (-a_1*k_a)*x1_0*x2_0 + a_1*k_c*x1_0*x3_0 + (-a_1*k_c)*x2_0*x3_0 + k_a*k_c*x1_1 + (a_1*k_a*k_c + k_a*V_m)*x1_0 + (-a_1*k_a*k_c)*x2_0 == 0\n", "\t(-a_2)*x1_0 + x2_1 + a_2*x2_0 == 0\n", "\tk_a*x1_0*x3_1 + b_1*k_a*x1_0*x3_0 + (-b_1*k_a)*x1_0*x4_0 + k_c*x3_1*x3_0 + b_1*k_c*x3_0^2 + (-b_1*k_c)*x3_0*x4_0 + k_a*k_c*x3_1 + (b_1*k_a*k_c + k_c*V_m)*x3_0 + (-b_1*k_a*k_c)*x4_0 == 0\n", "\t(-b_2)*x3_0 + x4_1 + b_2*x4_0 == 0\n", "}" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "denom = k_c*k_a + k_c*x3[0] + k_a*x1[0]\n", "system = DifferentialSystem([\n", " x1[1]*denom - a_1*(x2[0]-x1[0])*denom + k_a*V_m*x1[0],\n", " x2[1] - a_2*(x1[0] - x2[0]),\n", " x3[1]*denom - b_1*(x4[0]-x3[0])*denom + k_c*V_m*x3[0],\n", " x4[1] - b_2*(x3[0]-x4[0])\n", "], DPR, [x2,x3,x4])\n", "system" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us try to execute blindly the differential resultant algorithm:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Checking if there is any linear variable...\n", "Found linear variables [x4_*, x2_*]: we remove them first...\n", "##########################################################\n", "The system is linear: we use Macaulay resultant\n", "Getting the appropriate extension for having a SP2 system...\n", "Trying the extension (0, 0, 0)\n", "Trying the extension (1, 0, 0)\n", "Trying the extension (0, 1, 0)\n", "Trying the extension (0, 0, 1)\n", "Trying the extension (2, 0, 0)\n", "Trying the extension (1, 1, 0)\n", "Trying the extension (1, 0, 1)\n", "Trying the extension (0, 2, 0)\n", "Trying the extension (0, 1, 1)\n", "Trying the extension (0, 0, 2)\n", "Trying the extension (3, 0, 0)\n", "Trying the extension (2, 1, 0)\n", "Trying the extension (2, 0, 1)\n", "Trying the extension (1, 2, 0)\n", "Trying the extension (1, 1, 1)\n", "Trying the extension (1, 0, 2)\n", "Found the valid extension (1, 0, 2)\n", "Getting the homogenize version of the equations...\n", "Computing the Macaulay resultant...\n", "Timeout reached\n" ] } ], "source": [ "try:\n", " with timeout(3):\n", " res = system.diff_resultant(alg_res=\"iterative\", verbose=True)\n", "except TimeoutError:\n", " print(\"Timeout reached\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "SageMath 9.7", "language": "sage", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.5" } }, "nbformat": 4, "nbformat_minor": 4 }