Finding and Breaking Lie Symmetries: Implications for Structural Identifiability and Observability in Biological Modelling

Abstract

A dynamic model is structurally identifiable (respectively, observable) if it is theoretically possible to infer its unknown parameters (respectively, states) by observing its output over time. The two properties, structural identifiability and observability, are completely determined by the model equations. Their analysis is of interest for modellers because it informs about the possibility of gaining insight into a model’s unmeasured variables. Here we cast the problem of analysing structural identifiability and observability as that of finding Lie symmetries. We build on previous results that showed that structural unidentifiability amounts to the existence of Lie symmetries. We consider nonlinear models described by ordinary differential equations and restrict ourselves to rational functions. We revisit a method for finding symmetries by transforming rational expressions into linear systems. We extend the method by enabling it to provide symmetry-breaking transformations, which allows for a semi-automatic model reformulation that renders a non-observable model observable. We provide a MATLAB implementation of the methodology as part of the STRIKE-GOLDD toolbox for observability and identifiability analysis. We illustrate the use of the methodology in the context of biological modelling by applying it to a set of problems taken from the literature.

1. Introduction

The present work is motivated mainly by problems arising in identification and modelling of biological systems, although its results are applicable in other fields. We consider nonlinear dynamic models defined by ordinary differential equations. This framework is sufficiently powerful to model a wide range of biological processes, from intracellular networks to whole ecosystems, with the appropriate level of detail.

When compared to other applications, biological models often pose specific challenges due to the combination of nonlinearity, uncertainty about the underlying system, and experimental limitations regarding the possibility of perturbations and of measurements [1]. These features make the identification of biological models particularly challenging, and call for new methodological developments and computational tools. Indeed, many theoretical advances in nonlinear systems identification have been motivated by biological problems, even though the type of problems being considered are often common to other scientific areas, which makes the resulting methodologies generally applicable [2,3]. An example is the study of structural identifiability, a property that was introduced in the context of physiological modelling [4] and has since then been adopted in many scientific and technological areas; since its study is of particular interest for biological models, many related theoretical developments have been motivated by biological problems [5,6,7]. Structural identifiability refers to the possibility of inferring the values of the unknown (constant) parameters in the model equations from observations (measurements) of the model output. It is closely related to—and in fact can be considered a particular case of—another property, observability, which describes the possibility of determining the (time-varying) state variables of the model [8]. By considering the parameters as constant state variables, structural (local) identifiability amounts to observability (In accordance with the usual terminology, we speak of structural identifiability to distinguish it from practical identifiability. The latter quantifies the uncertainty in parameter values taking into account the information content of the available data, which may be limited due to experimental errors or low sampling rates. We also specify that it is structural local identifiability to distinguish it from structural global identifiability. Although observability is also a structural local property, for historical reasons we do not add those adjectives to its name, and refer to it simply as observability.).

The cause of lack of structural identifiability or observability can be traced back to the existence of symmetries in the model equations. In the present work we use Lie’s theory on symmetry analysis of differential equations [9,10,11]. It has been shown [12] that the existence of Lie symmetries is equivalent to the existence of similarity transformations, i.e., transformations of parameters and variables that leave the output(s) unchanged [13,14]. This means that the existence of symmetries amounts to lack of structural identifiability and/or observability.

Here we revisit a method for finding symmetries by transforming rational expressions into linear systems [15]. Determining the existence of symmetries in a model is a way of analysing its structural identifiability and observability. Furthermore, if symmetries exist, their mathematical expressions provide information about the relationships between model variables that cause loss of identifiability and/or observability. One way of exploiting these insights is by fixing one or more parameters involved in a symmetry, in order to render the remaining ones identifiable. Another way is by using the symmetry-breaking transformations to reformulate the model, applying the transformations so that the symmetries disappear and the new model is identifiable and observable. In this paper, we extend the method by enabling it to provide symmetry-breaking transformations, which allows for a semi-automatic model reformulation that renders a non-observable model observable. Thus, the approach can be used not only for characterizing the identifiability and observability of a model, but also for suggesting alternative reformulations if the original model does not possess those properties. We illustrate the usefulness and applicability of this approach in biological applications with four models of biochemical processes. Furthermore, we provide an implementation of the methodology as part of a new version of the STRIKE-GOLDD software https://github.com/afvillaverde/strike-goldd_2.1. STRIKE-GOLDD is a MATLAB toolbox that analyses structural identifiability and observability using a differential geometry approach [16,17,18].

The organization of this paper is as follows: we begin by providing an overview of the methodological aspects in Section 2. Then we illustrate its application to a number of modelling problems in Section 3. Finally, we discuss the implications and provide some conclusions in Section 4.

2. Methods

2.1. Structural Identifiability and Observability

For the study of observability and identifiability we consider the following type of models:

$$M_{NL}:=\left\{\begin{array}{cc}\dot{x}\left(t\right)\hfill & =f\left(x\right(t),\theta ,u(t),w(t\left)\right),\hfill \\ y\left(t\right)\hfill & =g\left(x\right(t),\theta ,u(t),w(t\left)\right),\hfill \\ x\left(t_0\right)\hfill & =x^0\left(\theta \right)\hfill \end{array}\right.$$

(1)

with a parameter vector $\theta \in {\mathbb{R}}^q$, known input functions $u\left(t\right)\in {\mathbb{R}}^{m_{\mathrm{u}}}$, unknown input functions $w\left(t\right)\in {\mathbb{R}}^{m_{\mathrm{w}}}$, state vector $x\left(t\right)\in {\mathbb{R}}^m$, output vector $y\left(t\right)\in {\mathbb{R}}^n$ and, f and g vectors constituted by analytical functions. From this point on, the time dependency will be omitted to simplify the notation.

Definition 1.

A parameter ${\theta }_i$, denoting the ith component, of the $M_{NL}$ model is structurally locally identifiable (SLI) if for almost any parameter vector ${\theta }^{\ast }\in {\mathbb{R}}^q$ there is a neighbourhood $\mathcal{N}\left({\theta }^{\ast }\right)$ in which the following relation holds:

$$\widehat{\theta }\in \mathcal{N}\left({\theta }^{\ast }\right)andy(t,\widehat{\theta })=y(t,{\theta }^{\ast })\Rightarrow {\widehat{\theta }}_i={\theta }_i^{\ast },i=1,\dots ,q.$$

(2)

If the above definition does not hold for any neighbourhood of ${\theta }^{\ast }$, ${\theta }_i$ is structurally unidentifiable (SU). If all model parameters are SLI, the model is SLI too. Otherwise, the model is SU.

Similarly, the observability of a state can be defined as:

Definition 2.

Given an output vector $y\left(t\right)$ and a known input vector $u\left(t\right)$, the state $x_i\left(\tau \right)$ is observable if it can be determined from $y\left(t\right)$ and $u\left(t\right)$ in the interval $t_0\le \tau \le t\le t_f$, where $t_f$ is a finite time. Otherwise, it is unobservable.

A model is observable if all states are observable. The same definition applies to the unknown entries vector but instead of observable, it is called reconstructible.

The term Full Input-State-Parameter Observability (FISPO) has been recently proposed [8] to refer to a model that is observable, identifiable and reconstructible. Its formal definition is given by:

Definition 3.

Given the model (1), the unknown quantities vector is considered $z\left(t\right)=[x\left(t\right),\theta ,w\left(t\right)]\in {\mathbb{R}}^{m+q+m_w}$. Denoting each component of z in time τ as $z_i\left(\tau \right)$ with $t_0\le \tau \le t\le t_f$, for a finite $t_f$. The model $M_{NL}$ is FISPO if the following condition holds:

$$\widehat{z}\left(\tau \right)\in \mathcal{N}\left(z^{\ast }\left(\tau \right)\right)andy(t,\widehat{z}\left(\tau \right))=y(t,z^{\ast }\left(\tau \right))\Rightarrow {\widehat{z}}_i\left(\tau \right)=z_i^{\ast }\left(\tau \right).$$

(3)

Observability and identifiability can be studied jointly if the parameters are considered as constant states. Thus, the observability of these states is equivalent to the local structural identifiability. Therefore, the augmented state vector is [18]:

$$\tilde{x}\left(t\right)=\left[\begin{array}{c}x\left(t\right)\\ \theta \\ w\left(t\right)\end{array}\right],\dot{\tilde{x}}\left(t\right)=\left[\begin{array}{c}f(\tilde{x}\left(t\right),u\left(t\right))\\ 0\\ \dot{w}\left(t\right)\end{array}\right].$$

(4)

In this way, the model (1) is transformed into:

$$M_{NL}:=\left\{\begin{array}{cc}\dot{\tilde{x}}\left(t\right)\hfill & =f(\tilde{x}\left(t\right),u\left(t\right)),\hfill \\ y\left(t\right)\hfill & =g(\tilde{x}\left(t\right),u\left(t\right)).\hfill \end{array}\right.$$

(5)

To determine the observability of a model it is necessary to calculate the nonlinear observability matrix ${\mathcal{O}}^{NL}$ using Lie derivatives. In the case of time-dependent entries, i.e., $u\left(t\right)$, Lie derivatives are computed as follows:

Definition 4.

The first order Lie derivative of $g\left(\tilde{x}\right)$ with respect to $f\left(\tilde{x}\right)$ is:

$$L_{f}g\left(\tilde{x}\right)={\displaystyle \frac{\partial g\left(\tilde{x}\right)}{\partial \tilde{x}}}f(\tilde{x},u)+\sum _{j=0}^{\infty }{\displaystyle \frac{\partial g\left(\tilde{x}\right)}{\partial u^{\left(j\right)}}}{u}^{(j+1)},$$

(6)

where $u^{\left(j\right)}$ is the jth derivative of u.

For higher orders, $i\ge 2$, the calculation is done through a recursive procedure:

$$L_f^{i}g\left(\tilde{x}\right)={\displaystyle \frac{\partial L_f^{i-1}g\left(\tilde{x}\right)}{\partial \tilde{x}}}f(\tilde{x},u)+\sum _{j=0}^{\infty }{\displaystyle \frac{\partial L_f^{i-1}g\left(\tilde{x}\right)}{\partial u^{\left(j\right)}}}{u}^{(j+1)}.$$

(7)

The observability-identifiability matrix of the previous model is defined as:

$${\mathcal{O}}_I(\tilde{x},u)=\left(\begin{array}{c}{\displaystyle \frac{\partial }{\partial \tilde{x}}}g(\tilde{x},u)\\ \frac{\partial }{\partial \tilde{x}}\left(L_{f}g(\tilde{x},u)\right)\\ \frac{\partial }{\partial \tilde{x}}\left(L_f^{2}g(\tilde{x},u)\right)\\ ⋮\\ {\displaystyle \frac{\partial }{\partial \tilde{x}}}\left(L_f^{n_{\tilde{x}-1}}g(\tilde{x},u)\right)\end{array}\right).$$

(8)

Theorem 1.

Nonlinear Observability-Identifiability Condition (OIC). If the model (5) satisfies $\mathrm{rank}\left({\mathcal{O}}_I({\tilde{x}}_0,u)\right)=m+q+m_w$, with ${\mathcal{O}}_I({\tilde{x}}_0,u)$ given by (8), where ${\tilde{x}}_0$ is a point in the augmented state space, then the model is locally observable and locally structurally identifiable in a neighbourhood $\mathcal{N}\left({\tilde{x}}_0\right)$ of ${\tilde{x}}_0$ [18].

If the OIC is not fulfilled, there is at least one unobservable (respectively, unidentifiable) state (respectively, parameter). In that case, full observability may sometimes be achieved by measuring more variables or functions. However, sometimes it is not possible to perform more measurements due to the characteristics of the experiments.

2.2. Lie Symmetries and Structural Unidentifiability

Yates et al. [12] showed that the existence of Lie symmetries entails the existence of similarity transformations, and therefore denotes lack of structural identifiability [13,14]. Similarity transformations allow the existence of parameters and variables transformations that leave the output(s) unchanged.

In Lie algebra, similarity transformations are one-parameter Lie group morphisms that map a solution of the differential equation onto themselves in terms of state variables. There are an infinite number of ways to represent this morphism, however the representation is unique when independent variables are fixed. The uniqueness of representation is a key property for the purpose of implementing the algorithm.

The expression of the one-parameter Lie group of transformations is:

$$x^{\ast }=X(x;\epsilon ).$$

(9)

Expanding the above expression in some neighborhood of $\epsilon =0$:

$$x^{\ast }=x+\epsilon \left({\displaystyle \frac{\partial X(x;\epsilon )}{\partial \epsilon }}{|}_{\epsilon =0}\right)+{\displaystyle \frac{1}{2}}{\epsilon }^2\left({\displaystyle \frac{{\partial }^{2}X(x;\epsilon )}{\partial {\epsilon }^2}}{|}_{\epsilon =0}\right)+\dots =x+\epsilon \left({\displaystyle \frac{\partial X(x;\epsilon )}{\partial \epsilon }}{|}_{\epsilon =0}\right)+O\left({\epsilon }^2\right),$$

(10)

where

$$\eta \left(x\right)={\displaystyle \frac{\partial X(x;\epsilon )}{\partial \epsilon }}{|}_{\epsilon =0}$$

(11)

is the infinitesimal of (9) and $x+\epsilon \eta \left(x\right)$ is the infinitesimal transformation of the Lie group of transformations.

Definition 5

([9]). The infinitesimal generator of the one-parameter Lie group of transformations is the differential operator:

$$X=X\left(x\right)=\eta \left(x\right)·▽=\sum _{i=1}^n{\eta }_i\left(x\right){\displaystyle \frac{\partial }{\partial x_i}},$$

(12)

where ▽ defines the gradient

$$▽=\left({\displaystyle \frac{\partial }{\partial x_1}},{\displaystyle \frac{\partial }{\partial x_2}},\dots ,{\displaystyle \frac{\partial }{\partial x_n}}\right).$$

(13)

Theorem 2.

First Fundamental Theorem of Lie [9]. Given an initial value problem (IVP) for a system of first-order ODEs:

$${\displaystyle \frac{d{x}^{\ast }}{d\tau }}=\eta \left(x^{\ast }\right),$$

(14)

$$x^{\ast }=x\mathit{when}\tau =0.$$

(15)

The Lie group of transformations (9) is equivalent to the above IVP with the parametrization $\tau \left(\epsilon \right)$:

$$\tau \left(\epsilon \right)={\int }_0^{\epsilon }\Gamma \left({\epsilon }^{\prime }\right)d{\epsilon }^{\prime },$$

where

$$\Gamma \left(\epsilon \right)={\displaystyle \frac{\partial \varphi (a,b)}{\partial b}}{|}_{(a,b)=({\epsilon }^{-1},\epsilon )},$$

$$\Gamma \left(0\right)=1.$$

$\varphi (a,b)$ is the law of composition and ${\epsilon }^{-1}$ denotes the inverse element to ε.

The above theorem shows how infinitesimal transformations contain the essential information that determines the uniparametric Lie group of transformations. From Theorem (2) and without loss of generality, it is assumed that the Lie group of transformations has as its law of composition $\varphi (a,b)=a+b$ with ${\epsilon }^{-1}=-\epsilon$ and $\Gamma \left(\epsilon \right)=1$ [9]. Therefore, the uniparametric Lie group of transformations (9) is rewritten, in terms of its infinitesimals $\eta \left(x^{\ast }\right)$, as:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{x}^{\ast }}{d\epsilon }}& =\eta \left(x^{\ast }\right),\hfill \\ \hfill x^{\ast }& =x\mathrm{at}\epsilon =0.\hfill \end{array}$$

(16)

The above expression defines an Initial Value Problem (IVP) for the Lie group of transformations in terms of its infinitesimals generators.

The exponential parametrization of the Lie group around the identity is:

Theorem 3

([9]). The one-parameter Lie group of transformations (9) is equivalent to:

$$x^{\ast }=exp\left[\epsilon X\right]x=x+\epsilon Xx+{\displaystyle \frac{1}{2}}{\epsilon }^2{X}^{2}x+\dots =$$

(17)

$$=\left(1+\epsilon X+{\displaystyle \frac{1}{2}}{\epsilon }^2{X}^2+\dots \right)x=\sum _{k=0}^{\infty }{\displaystyle \frac{{\epsilon }^k}{k!}}{X}^{k}x,$$

where X is given by (12) and $X^k=X{X}^{k-1},k=1,2\dots$ with $X^{0}x=x$.

For a more detailed background, we refer the reader to [9,11].

2.2.1. Computing Symmetries

Let us consider the same ODE system as in (1). The state vector will be augmented with parameters and unknown inputs, as mentioned earlier:

$$\begin{array}{cc}\hfill {\dot{x}}_i\left(t\right)& =f_i(x\left(t\right),u\left(t\right)),i=1,\dots ,m\hfill \\ \hfill x_i\left(t\right)& =\theta ,i=m+1,\dots ,m+q\hfill \\ \hfill x_i\left(t\right)& =w_i\left(t\right),i=m+q+1,\dots ,n^{\ast }\hfill \end{array}$$

where $n^{\ast }=m+q+m_w$.

We will study three different infinitesimals. Considering $d_{max}$, $d\in \mathbb{N}$ as the maximum degree for polynomials and $r_{i,d}$ unknown constants to determine, the expressions of infinitesimals are:

• Univariate:

  $${\eta }_i\left(x\right)=\sum _{d=0}^{d_{max}}{r}_{i,d}{x}_i^d,i=1,\dots ,n^{\ast }.$$

  (18)
• Partially variate:

  $$\begin{array}{cc}\hfill {\eta }_i\left(x\right)& =\sum _{d_i,d_{m+1},\dots ,d_{m+q}=0}^{|d|=d_{max}}{r}_{i,d}{x}_i^{d_i}{x}_{m+1}^{d_{m+1}}···x_{m+q}^{d_{m+q}},i=1,\dots m,\hfill \\ \hfill {\eta }_i\left(x\right)& =\sum _{d_{m+1},\dots ,d_{m+q}=0}^{|d|=d_{max}}{r}_{i,d}{x}_{m+1}^{d_{m+1}}···x_{m+q}^{d_{m+q}},i=m+1,\dots m+q,\hfill \\ \hfill {\eta }_i\left(x\right)& =\sum _{d_i,d_{m+1},\dots ,d_{m+q}=0}^{|d|=d_{max}}{r}_{i,d}{x}_i^{d_i}{x}_{m+1}^{d_{m+1}}···x_{m+q}^{d_{m+q}},i=m+q+1,\dots n^{\ast }.\hfill \end{array}$$

  (19)
• Multivariate:

  $$\begin{array}{cc}\hfill {\eta }_i\left(x\right)& =\sum _{d_1,\dots ,d_{m+q}=0}^{|d|=d_{max}}{r}_{i,d}{x}_1^{d_1}···x_{m+q}^{d_{m+q}},i=1,\dots m,\hfill \\ \hfill {\eta }_i\left(x\right)& =\sum _{d_{m+1},\dots ,d_{m+q}=0}^{|d|=d_{max}}{r}_{i,d}{x}_{m+1}^{d_{m+1}}···x_{m+q}^{d_{m+q}},i=m+1,\dots m+q,\hfill \\ \hfill {\eta }_i\left(x\right)& =\sum _{d_1,\dots ,d_{n^{\ast }}=0}^{|d|=d_{max}}{r}_{i,d}{x}_1^{d_1}···x_{n^{\ast }}^{d_{n^{\ast }}},i=1,\dots n^{\ast }.\hfill \end{array}$$

  (20)

The derivative of infinitesimal generators is also defined, so that it can act on the $\dot{x}\left(t\right)$:

$$X^{\prime }=\sum _{i=1}^{n^{\ast }}{\eta }_i\left(x\right){\displaystyle \frac{\partial }{\partial x_i}}+\sum _{i=1}^{n^{\ast }}{\eta }_i^{\prime }\left(x\right){\displaystyle \frac{\partial }{\partial \dot{x_i}}},$$

(21)

where

$${\eta }_i^{\prime }\left(x\right)=\sum _{j=1}^{n^{\ast }}{\dot{x}}_j{\displaystyle \frac{\partial {\eta }_i}{\partial x_j}}.$$

Using the above formulation for infinitesimals generators, the following explicit criterion for admittance of a Lie group of transformations is obtained:

Theorem 4

([9,10]). The system of ordinary differential equations admits a one-parameter Lie group of transformations defined by the infinitesimal generator (12) if and only if:

$$X^{\prime }·({\dot{x}}_k-f_k\left(x\right))=0,k=1,\dots ,m$$

(22)

$$X·(y_l-g_l\left(x\right))=0,l=1,\dots ,n$$

(23)

Applying the previous theorem to the initial system (1), we obtain an explicit criterion:

$$\begin{array}{cc}\hfill \sum _{j=1}^{n^{\ast }}{\dot{x}}_j{\displaystyle \frac{\partial {\eta }_k}{\partial x_j}}\left(x\right)-\sum _{i=1}^{n^{\ast }}{\eta }_i\left(x\right){\displaystyle \frac{\partial f_k}{\partial x_i}}\left(x\right)=& 0,k=1,\dots ,m\hfill \\ \hfill \sum _{i=1}^{n^{\ast }}{\eta }_i\left(x\right){\displaystyle \frac{\partial g_l}{\partial x_i}}\left(x\right)=& 0,l=1,\dots ,n\hfill \end{array}$$

(24)

An admitted Lie symmetry is a continuous group of transformations X such that the observation (observed datum) is unchanged:

$$g(x^{\ast }\left(t\right),u^{\ast }\left(t\right))=g(x\left(t\right),u\left(t\right))$$

The output map should not be modified. The above system defines a system of partial differential equations. If we consider the rational form of $f_k\left(x\right)$ and $g_l\left(x\right)$:

$$\begin{array}{cc}\hfill {\dot{x}}_k& =f_k\left(x\right)={\displaystyle \frac{P^k\left(x\right)}{Q^k\left(x\right)}},k=1,\dots ,m,\hfill \\ \hfill y_l& =g_l\left(x\right)={\displaystyle \frac{R^l\left(x\right)}{S^l\left(x\right)}},l=1,\dots ,n,\hfill \end{array}$$

(25)

then, the system of PDEs (24) is converted into a system of ODEs.

There are different expressions of (24) depending on the infinitesimal generator:

• Univariate and partial:

  $$\begin{array}{cc}\hfill P^k{Q}^k{\displaystyle \frac{\partial {\eta }_k}{\partial x_k}}-\sum _{i=1}^{n^{\ast }}{\eta }_i[P_{x_i}^k{Q}^k-P^k{Q}_{x_i}^k]=& 0,k=1,\dots m,\hfill \\ \hfill \sum _{i=1}^{n^{\ast }}{\eta }_i[R_{x_i}^l{S}^l-R^l{S}_{x_i}^l]=& 0,l=1,\dots n.\hfill \end{array}$$

  (26)
• Multivariate:

  $$\begin{array}{cc}\hfill \sum _{j=1}^m{P}^j{Q}^k\left(\prod _{b\ne j}{Q}^b\right){\displaystyle \frac{\partial {\eta }_k}{\partial x_j}}-\sum _{i=1}^{n^{\ast }}{\eta }_i\left(\prod _{b\ne k}{Q}^b\right)[P_{x_i}^k{Q}^k-P^k{Q}_{x_i}^k]=& 0,k=1,\dots m,\hfill \\ \hfill \sum _{i=1}^{n^{\ast }}{\eta }_i[R_{x_i}^l{S}^l-R^l{S}_{x_i}^l]=& 0,l=1,\dots n.\hfill \end{array}$$

  (27)

Each of these equations can be reordered based on the combinations among the components of x. Let r be a vector containing all $r_{i,d}$:

$$\sum _{i_1,\dots ,i_n}{c}_{i_1,\dots ,i_n}\left(r\right)x_1^{i_1}···x_n^{i_n}=0.$$

(28)

The coefficients $c_{i_1,\dots ,i_n}$ are linear in r and its matrix form expression allows to reformulate condition (26) and (27) into:

$$C·r=0.$$

The problem of finding symmetries is equivalent to solving a linear system of equations with numeric entries. However, C is a non-square matrix and, in order to find all the solutions, it is necessary to compute its kernel. It is possible that the obtained solutions are not independent of each other, as a result of linear combinations or multiplication by $x_i$; for this reason, we will only consider solutions that are independent of each other and in their greatest degree of simplification.

The next step is to build the expression of $x^{\ast }$ with the infinitesimal generators as in Theorem (3). When the infinitesimal transformation is given by powers of the variable, the exact transformation is known and it is classified as “elementary”. Some examples are:

$$\begin{array}{cc}\hfill x_i^{\ast }& =x_i+\epsilon ,X={\displaystyle \frac{\partial }{\partial x_i}}\left(\mathrm{translation}\right)\hfill \\ \hfill x_i^{\ast }& =exp\left(\epsilon \right)x_i,X=x_i{\displaystyle \frac{\partial }{\partial x_i}}\left(\mathrm{scaling}\right)\hfill \\ \hfill x_i^{\ast }& ={\displaystyle \frac{x_i}{1-\epsilon x_i}},X=x_i^2{\displaystyle \frac{\partial }{\partial x_i}}\left(\mathrm{Mobius}\right)\hfill \\ \hfill x_i^{\ast }& ={\displaystyle \frac{x_i}{{[1-(p-1)\epsilon x_i^{p-1}]}^{\frac{1}{p-1}}}},X=x_i^p{\displaystyle \frac{\partial }{\partial x_i}}(\mathrm{higher}\mathrm{order})\hfill \end{array}$$

(29)

The most common types of symmetries in biochemical models are translation and scaling. However, those elementary transformations cover only a part of the possible symmetries that a model can contain. The others must be calculated using Lie series or solving the IVP (16).

It is possible to maximize the number of elementary transformations given by an infinitesimal generator before applying Lie series. This process starts by searching for all parameter combinations that forbid an elementary transformation, and dividing the initial infinitesimal generator by them. Each new infinitesimal generator provides a number of elementary transformations that can be greater than, lower than or equal to the initial one. Maximizing the total amount, we will obtain an infinitesimal generator that provides the largest possible number of elementary transformations.

2.2.2. Initial Conditions

The definition of a dynamic model may include specific initial conditions (ICS), which can be numeric, parametric (known or unknown) or a combination of both [19,20]. Perturbations in ICS must produce changes in the output $y\left(t\right)$ for the model to be observable.

If the symmetries are studied without taking into account the ICS, it may happen that the generators do not fulfill them. It is important to consider only the generators that satisfy the ICS [19]; by including them as output vectors, a symmetry is admitted by the system if and only if:

$$X·{\left(x_k\left(t_0\right)-x^0\left(\theta \right)\right)}_{{|}_{x=x^0\left(\theta \right)}}=0,k=1,\dots ,m,$$

(30)

where $x^0\left(\theta \right)$ are the ICS (parametric, numerical or a combination of both). Expanding the above equation:

$$\sum _{i=1}^{n^{\ast }}{\eta }_i{\displaystyle \frac{\partial }{\partial x_i}}{x_k\left(t_0\right)}_{{|}_{x=x^0\left(\theta \right)}}-\sum _{i=1}^{n^{\ast }}{\eta }_i{\displaystyle \frac{\partial }{\partial x_i}}{x^0\left(\theta \right)}_{{|}_{x=x^0\left(\theta \right)}}=0,k=1,\dots ,m.$$

(31)

Considering the rational form of $x^0\left(\theta \right)$:

$$x^0\left(\theta \right)={\displaystyle \frac{V^k\left(x\right)}{W^k\left(x\right)}},k=1,\dots ,m.$$

Expression (31) is reformulated as follow:

$$\sum _{i=1}^{n^{\ast }}{\eta }_i\left(x^0\left(\theta \right)\right)-\sum _{i=1}^{n^{\ast }}{\eta }_i{{\displaystyle \frac{V_{x_i}^k{W}^k-V^k{W}_{x_i}^k}{{\left(W^k\right)}^2}}}_{{|}_{x=x^0\left(\theta \right)}}=0,k=1,\dots ,m.$$

(32)

This new restriction allows us to consider only those symmetries that fulfill the initial conditions, and provides a tool for study the influence of them.

2.3. Implementation

An overview of the algorithm described in the preceding subsections is shown in Figure 1. We have made a MATLAB implementation available at https://github.com/GemmaMasFes6/Lie-Symmetries, and as part of a new version (v2.1.6) of the STRIKE-GOLDD toolbox (https://github.com/afvillaverde/strike-goldd_2.1); it will also be included in future releases of said toolbox. Our code represents an addition to the set of existing tools for studying differential equation symmetries, which include symmetryDetection [15] in Python, MinimalOutputSets [19] in Mathematica, and SADE [21] in Maple. Our software is open source and, to the best of our knowledge, it is the first tool of its characteristics available in MATLAB.

The code consists of a main MATLAB script, ‘Lie_Symmetry’, and ten auxiliary functions defined in separate files. Each of these functions performs one of the stages outlined above: calculation of infinitesimal polynomials (univariate, partially variate and multivariate); computation of polynomials for states (depending on the type of polynomial there are two possibilities), observations, and ICS (if they are specified); and obtaining the transformations. The last step in turn incorporates two other functions, corresponding to maximizing the number of elementary transformations and, if necessary, calculating the transformations using Lie series.

The algorithm has a number of options that can be specified by the user: the type of infinitesimal polynomial, its degree, and the number of terms in Lie series (in case it has to be used).

The input of the programme must include the following vectors of symbolic variables, declared through the MATLAB sym command: parameters, states, initial conditions, ODEs, observations, and inputs (known or unknown). For initial conditions, two vectors must be provided: a vector called known_ics with entries equal to 1 for a known initial condition and 0 otherwise, and a vector ics with the values of the known ICS, either numeric or (known or unknown) parametric value.

The output of the programme first reports whether there is any symmetry. If a symmetry exists, the programme prints the infinitesimal generator and the transformations on the screen.

3. Results

To illustrate the application of the method in biological and biomedical modelling, we use it to analyse a set of models taken from the literature. The models are listed in

Table 1. List of models analysed in this paper and summary of their features.

Model Name (and Acronym)	Reference	States	Parameters	Outputs
Simple chemical reaction (CR)	[15]	A	$k,s_1,s_2$	$A^{\mathrm{obs}}$
Pharmacokinetic model (PK)	[22]	$x_1,x_2,x_3,x_4$	$k_1,k_2,k_3,k_5,$	$x_2^{\mathrm{obs}},x_3^{\mathrm{obs}}$
			$k_6,k_7,s_2,s_3$
NF-$\kappa$B signalling pathway (NFKB)	[23]	$x_1,x_2,x_3,x_4,x_5,$	$k_0,k_1,k_{1p},k_2,k_3,k_4,$	$y_1,y_2,y_3,y_4$
		$x_6,x_7,x_8,x_9,x_{10}$	$k_5,k_6,k_7,k_8,k_9,k_{10},k_{11},$
			$s_1,s_2,s_3,s_4,{\rho }_{\mathrm{vol}},I_{0_{\mathrm{cyt}}},I_{0_{\mathrm{nuc}}}$
Glucose-insulin regulation (Bolie)	[24]	$q_1,q_2$	$p_1,p_2,p_3,p_4,V_p$	h

and their schematic representations are shown in Figure 2.

3.1. Simple Chemical Reaction

This model represents a bimolecular reaction described by one ODE and one observation [15]:

$$\begin{array}{cc}\hfill \dot{A}& =-2k{A}^2,\hfill \\ \hfill A^{\mathrm{obs}}& =s_1{\displaystyle \frac{A}{1+s_{2}A}}.\hfill \end{array}$$

(33)

It is used to provide a basic illustration of the methodology, due to its simplicity. Without considering initial conditions, and using an univariate polynomial of second order, the programme finds two infinitesimal generators:

$$\begin{array}{cc}\hfill X& =A{\displaystyle \frac{\partial }{\partial A}}-k{\displaystyle \frac{\partial }{\partial k}}-s_1{\displaystyle \frac{\partial }{\partial s_1}}-s_2{\displaystyle \frac{\partial }{\partial s_2}},\hfill \\ \hfill X& =A^2{\displaystyle \frac{\partial }{\partial A}}+{\displaystyle \frac{\partial }{\partial s_2}}.\hfill \end{array}$$

(34)

All the transformations are elementary:

$$\begin{array}{cc}\hfill & A^{\ast }=e^{\epsilon }A,k^{\ast }=e^{-\epsilon }k,s_1^{\ast }=e^{-\epsilon }{s}_1,s_2^{\ast }=e^{-\epsilon }{s}_2.\hfill \\ \hfill & A^{\ast }={\displaystyle \frac{A}{1-\epsilon A}},s_2^{\ast }=s_2+\epsilon .\hfill \end{array}$$

(35)

Our results coincide with those reported in [15]. It is possible to include ICS in order to study its influence in the model.

This example, because of its simplicity, allows us to check the results manually. Once the transformations are computed, it is easy to see that the second group of transformations solves the same ODE. The time derivative of $A^{\ast }$ is:

$${\dot{A}}^{\ast }={\displaystyle \frac{\dot{A}}{{(1-\epsilon A)}^2}}$$

Incorporating the above expression with $A^{\ast }$ in (33), the ODE is still fulfilled.

Below are two screenshots of the results of the observability and identifiability analysis obtained with STRIKE-GOLDD. In the first panel, Figure 3 (Page 1), corresponding to the initial model, all states and parameters are unobservable; in the second one, Figure 3 (Page 2), corresponding to the model with Lie transformations given by the second generator, states and parameters are observable.

3.2. Pharmacokinetic Model

A pharmacokinetic (PK) model describes the time course of the concentrations of a drug in different compartments, after entering an organism. This model includes one input, four ODEs, and two outputs [22,25]:

$$\begin{array}{cc}\hfill \dot{x_1}& =u-(k_1+k_2)x_1,\hfill \\ \hfill \dot{x_2}& =k_1{x}_1-(k_3+k_6+k_7)x_2+k_5{x}_4,\hfill \\ \hfill \dot{x_3}& =k_2{x}_1+k_3{x}_2-k_4{x}_3,\hfill \\ \hfill \dot{x_4}& =k_6{x}_2-k_5{x}_4,\hfill \\ \hfill x_2^{\mathrm{obs}}& =s_2{x}_2,\hfill \\ \hfill x_3^{\mathrm{obs}}& =s_3{x}_3.\hfill \end{array}$$

We use partial variate polynomials of second order, without ICS. Maximizing the number of elementary transformations leads to four of them, and the procedure yields the following infinitesimal generator:

$$X=k_1\left({\displaystyle \frac{\partial }{\partial k_1}}-{\displaystyle \frac{\partial }{\partial k_2}}\right)-{\displaystyle \frac{k_3(k_1+k_2)}{k_2}}\left({\displaystyle \frac{\partial }{\partial k_3}}-{\displaystyle \frac{\partial }{\partial k_7}}\right)-s_2{\displaystyle \frac{\partial }{\partial s_2}}+{\displaystyle \frac{k_1{s}_3}{k_2}}{\displaystyle \frac{\partial }{\partial s_3}}+x_2{\displaystyle \frac{\partial }{\partial x_2}}-{\displaystyle \frac{k_1{s}_3}{k_2}}{\displaystyle \frac{\partial }{\partial x_3}}+x_4{\displaystyle \frac{\partial }{\partial x_4}}.$$

The formulation of the IVP (16) for this infinitesimal generator, considering only the non-elementary transformations, is:

${\dot{k_2}}^{\ast }=-k_1^{\ast },$	$k_2^{\ast }\left(0\right)=k_2$ ,
${\dot{k_3}}^{\ast }=-k_3^{\ast }{\displaystyle \frac{(k_1^{\ast }+k_2^{\ast })}{k_2^{\ast }}},$	$k_3^{\ast }\left(0\right)=k_3$ ,
${\dot{k_7}}^{\ast }=k_3^{\ast }{\displaystyle \frac{(k_1^{\ast }+k_2^{\ast })}{k_2^{\ast }}},$	$k_7^{\ast }\left(0\right)=k_7$ ,
${\dot{s_3}}^{\ast }={\displaystyle \frac{k_1^{\ast }{s}_3^{\ast }}{k_2^{\ast }}},$	$s_3^{\ast }\left(0\right)=s_3$ ,
${\dot{x_3}}^{\ast }=-{\displaystyle \frac{k_1^{\ast }{x}_3^{\ast }}{k_2^{\ast }}},$	$x_3^{\ast }\left(0\right)=x_3$ .

The solution of the ODE system, replacing $k_1^{\ast }$ with its transformation, is:

$$\begin{array}{cc}\hfill k_2^{\ast }=& k_1+k_2-k_1{e}^{\epsilon },\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill k_3^{\ast }=& {\displaystyle \frac{k_3{e}^{-\epsilon }(k_1+k_2-k_1{e}^{\epsilon })}{k_2}},\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill k_7^{\ast }=& k_7+{\displaystyle \frac{k_3(k_1+k_2)}{k_2}}-{\displaystyle \frac{k_3{e}^{-\epsilon }(k_1+k_2)}{k_2}},\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill x_3^{\ast }=& {\displaystyle \frac{x_3(k_1+k_2-k_1{e}^{\epsilon })}{k_2}},\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill s_3^{\ast }=& {\displaystyle \frac{k_2{s}_3}{(k_1+k_2-k_1{e}^{\epsilon })}}.\hfill \end{array}$$

(40)

These new transformations coincide with those presented in [25]. Using Lie series in the infinitesimal generator, the new expressions are:

$$\begin{array}{cc}\hfill x_2^{\ast }& =x_2{e}^{\epsilon },x_4^{\ast }=x_4{e}^{\epsilon },k_1^{\ast }=k_1{e}^{\epsilon },s_2^{\ast }=s_2{e}^{-\epsilon },\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill x_3^{\ast }& =x_3-{\displaystyle \frac{\epsilon k_1{x}_3}{k_2}}-{\displaystyle \frac{{\epsilon }^2{k}_1{x}_3}{2k_2}}-{\displaystyle \frac{{\epsilon }^3{k}_1{x}_3}{6k_2}}-{\displaystyle \frac{{\epsilon }^4{k}_1{x}_3}{24k_2}},\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill k_2^{\ast }& =k_2-\epsilon k_1-{\displaystyle \frac{{\epsilon }^2{k}_1}{2}}-{\displaystyle \frac{{\epsilon }^3{k}_1}{6}}-{\displaystyle \frac{{\epsilon }^4{k}_1}{24}},\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill k_3^{\ast }& =k_3-{\displaystyle \frac{k_3(k_1+k_2)\epsilon }{k_2}}+{\displaystyle \frac{{\epsilon }^2{k}_3(k_1+k_2)}{2k_2}}-{\displaystyle \frac{{\epsilon }^3{k}_3(k_1+k_2)}{6k_2}}+{\displaystyle \frac{{\epsilon }^4{k}_3(k_1+k_2)}{24k_2}},\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill k_7^{\ast }& =k_7+{\displaystyle \frac{k_3(k_1+k_2)\epsilon }{k_2}}-{\displaystyle \frac{{\epsilon }^2{k}_3(k_1+k_2)}{2k_2}}+{\displaystyle \frac{{\epsilon }^3{k}_3(k_1+k_2)}{6k_2}}-{\displaystyle \frac{{\epsilon }^4{k}_3(k_1+k_2)}{24k_2}},\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill s_3^{\ast }& =s_3+{\displaystyle \frac{\epsilon k_1{s}_3}{k_2}}+{\displaystyle \frac{{\epsilon }^2{k}_1{s}_3(2k_1+k_2)}{2k_2^2}}+{\displaystyle \frac{{\epsilon }^3{k}_1{s}_3(6k_1^2+6k_1{k}_2+k_2^2)}{6k_2^3}}+\hfill \\ \hfill & +{\displaystyle \frac{{\epsilon }^4{k}_1{s}_3(24k_1^3+36k_1^2{k}_2+14k_1{k}_2^2+k_2^3)}{24k_2^4}}.\hfill \end{array}$$

(46)

These expressions seem to differ from the closed form, however, the result can be reformulated in order to obtain them. Taking as example $x_3^{\ast }$:

$$x_3^{\ast }=x_3-{\displaystyle \frac{k_1{x}_3}{k_2}}\left(\epsilon +{\displaystyle \frac{{\epsilon }^2}{2}}+{\displaystyle \frac{{\epsilon }^3}{6}}+{\displaystyle \frac{{\epsilon }^4}{24}}\right).$$

(47)

The last part of the previous expression consists of the first terms of the following series:

$$\sum _{n=0}^{\infty }{\displaystyle \frac{{\epsilon }^n}{n!}}-1=e^{\epsilon }-1.$$

(48)

Substituting (48) in (47), the transformation of $x_3^{\ast }$ is the same as the one obtained from IVP:

$$x_3^{\ast }=x_3-{\displaystyle \frac{k_1{x}_3}{k_2}}(e^{\epsilon }-1)={\displaystyle \frac{x_3{k}_2-k_1{x}_3(e^{\epsilon }-1)}{k_2}}={\displaystyle \frac{x_3(k_2-k_1(e^{\epsilon }-1))}{k_2}}.$$

The most complicated case appears to be $s_3^{\ast }$. The expression obtained from the IVP (40) can be rearranged as:

$$s_3^{\ast }={\displaystyle \frac{k_2{s}_3}{k_2+k_1(1-e^{\epsilon })}}={\displaystyle \frac{s_3}{1-\frac{k_1}{k_2}(e^{\epsilon }-1)}}=s_3\sum _{n=0}^{\infty }{\left(\frac{k_1}{k_2}(e^{\epsilon }-1)\right)}^n.$$

(49)

It is only necessary to prove that (46) is the same as (49) to show that both transformations are equal. Equation (46) can be reordered in terms of the powers of $k_1/k_2$:

$$s_3^{\ast }=s_3\left(1+{\displaystyle \frac{k_1}{k_2}}\left(\epsilon +{\displaystyle \frac{{\epsilon }^2}{2}}+{\displaystyle \frac{{\epsilon }^3}{6}}+{\displaystyle \frac{{\epsilon }^4}{24}}\right)+{\displaystyle \frac{k_1^2}{k_2^2}}\left({\epsilon }^2+{\epsilon }^3+{\displaystyle \frac{7{\epsilon }^4}{12}}\right)+{\displaystyle \frac{k_1^3}{k_2^3}}\left({\epsilon }^3+{\displaystyle \frac{3{\epsilon }^4}{2}}\right)+{\displaystyle \frac{k_1^4{\epsilon }^4}{k_2^4}}\right).$$

(50)

The coefficient of $k_1/k_2$ is $e^{\epsilon }-1$, as was proven before in (47) and (48).

The coefficient of the third term must be equal to

$${(e^{\epsilon }-1)}^2=e^{2\epsilon }+1-2e^{\epsilon }.$$

Considering the Taylor series of the previous expression until the fifth order,

$$1+2\epsilon +2{\epsilon }^2+{\displaystyle \frac{8{\epsilon }^3}{6}}+{\displaystyle \frac{16{\epsilon }^4}{24}}+1-2-2\epsilon -{\epsilon }^2-{\displaystyle \frac{{\epsilon }^3}{3}}-{\displaystyle \frac{2{\epsilon }^4}{24}}={\epsilon }^2+{\epsilon }^3+{\displaystyle \frac{7{\epsilon }^4}{12}},$$

(51)

we obtain that the final expression is the same as that presented in (50).

The coefficient of $k_1^3/k_2^3$ is the result of considering some terms of the Taylor series of ${(e^{\epsilon }-1)}^3$:

$$3\left(1+\epsilon +{\displaystyle \frac{{\epsilon }^2}{2}}+{\displaystyle \frac{{\epsilon }^3}{6}}+{\displaystyle \frac{{\epsilon }^4}{24}}\right)-3\left(1+2\epsilon +2{\epsilon }^2+{\displaystyle \frac{8{\epsilon }^3}{6}}+{\displaystyle \frac{16{\epsilon }^4}{24}}\right)+\left(1+3\epsilon +{\displaystyle \frac{9{\epsilon }^2}{2}}+{\displaystyle \frac{27{\epsilon }^3}{6}}+{\displaystyle \frac{81{\epsilon }^4}{24}}\right)-1=$$

$$={\epsilon }^3+{\displaystyle \frac{36{\epsilon }^4}{24}}={\epsilon }^3+{\displaystyle \frac{3{\epsilon }^4}{2}}.$$

It is necessary to consider more terms of the Lie series to prove that the coefficient of $k_1^4/k_2^4$ is ${(e^{\epsilon }-1)}^4$. In any case, regardless of the model under consideration, the transformations obtained with Lie series—that is, the output of the programme—are the expanding form of the solution of the IVP. In order to obtain a good approximation of the closed form (36)–(40) the programme needs to consider sufficient terms of the Lie series. For this case, two terms are sufficient to achieve full observability (i.e., FISPO).

This model had been previously analysed in [15,25]. In [15] a different infinitesimal generator was obtained; due to the lack of a method to compute non-elementary transformations, it was not possible to obtain all transformations. In [25] the same infinitesimal generator as reported here was obtained, and transformations were computed using Hermite-Padé polynomials. The resulting transformations in [25] are the same as those obtained from the IVP, as well as those obtained by our programme using Lie series.

3.3. NF-$\kappa$B Signalling Pathway

The model studied in this example was described in [23,26]. It represents a cellular signaling pathway found in most animal cells, corresponding to the NF-$\kappa$B transcription factor. It is depicted in Figure 2C, where black arrows indicate exit routes.

$$\begin{array}{cc}\hfill {\dot{x}}_1& =k_{11}{x}_{10}-\left({\displaystyle \frac{k_{1}u}{1+k_{0}u}}+k_{1p}\right)x_1,\hfill \\ \hfill {\dot{x}}_2& =\left({\displaystyle \frac{k_{1}u}{1+k_{0}u}}+k_{1p}\right)x_1-k_2{x}_2,\hfill \\ \hfill {\dot{x}}_3& =k_2{x}_2-k_3{x}_3,\hfill \\ \hfill {\dot{x}}_4& =k_2{x}_2-k_4{x}_4,\hfill \\ \hfill {\dot{x}}_5& =k_3{\rho }_{\mathrm{vol}}{x}_3-k_5{x}_5,\hfill \\ \hfill {\dot{x}}_6& =k_5{x}_5-k_{10}{x}_9{x}_6,\hfill \\ \hfill {\dot{x}}_7& =k_6{x}_6-k_7{x}_7,\hfill \\ \hfill {\dot{x}}_8& =k_8{x}_7-k_9{x}_8,\hfill \\ \hfill {\dot{x}}_9& =k_9{\rho }_{\mathrm{vol}}{x}_8-k_{10}{x}_9{x}_6,\hfill \\ \hfill {\dot{x}}_{10}& =k_{10}{x}_9{x}_6-k_{11}{\rho }_{\mathrm{vol}}{x}_{10},\hfill \\ \hfill y_1^{\mathrm{obs}}& =s_1(x_1+x_2+x_3)+I_{0_{\mathrm{cyt}}},\hfill \\ \hfill y_2^{\mathrm{obs}}& =s_2(x_{10}+x_5+x_6)+I_{0_{\mathrm{nuc}}},\hfill \\ \hfill y_3^{\mathrm{obs}}& =s_3(x_2+x_3),\hfill \\ \hfill y_4^{\mathrm{obs}}& =s_4(x_2+x_4).\hfill \end{array}$$

This model is the only one studied with ICS. In this case, they are parametric and include an additional parameter, $x_{1_0}$, that does not appear in the model equations:

$$\begin{array}{cc}\hfill x_1\left(0\right)& =x_{1_0},x_2\left(0\right)={\displaystyle \frac{k_{1p}{x}_1}{k_2}},x_3\left(0\right)={\displaystyle \frac{k_{1p}{x}_1}{k_3}},\hfill \\ \hfill x_4\left(0\right)& ={\displaystyle \frac{k_{1p}{x}_1}{k_4}},x_5\left(0\right)={\displaystyle \frac{k_3{\rho }_{\mathrm{vol}}{x}_3}{k_5}},x_6\left(0\right)={\displaystyle \frac{k_7{x}_7}{k_6}},\hfill \\ \hfill x_7\left(0\right)& ={\displaystyle \frac{k_9{x}_8}{k_8}},x_8\left(0\right)={\displaystyle \frac{k_3{x}_3}{k_9}},x_9\left(0\right)={\displaystyle \frac{k_5{x}_5}{x_6{k}_{10}}},x_{10}\left(0\right)={\displaystyle \frac{k_{1p}{x}_1}{k_{11}}}.\hfill \end{array}$$

Using a second order univariate polynomial, the following infinitesimal transformations were found:

• Scaling symmetry was found for the known input function and two parameters:

  $$u^{\ast }=u{e}^{-\epsilon };k_i^{\ast }=k_i{e}^{\epsilon }i=0,1.$$
• The second symmetry was also for the input function and one parameter, in this case, a Mobius and translation symmetry, respectively:

  $$u^{\ast }=-{\displaystyle \frac{u}{\epsilon u-1}};k_0^{\ast }=k_0+\epsilon .$$
• Another scaling symmetry involving one state and two parameters:

  $$x_7^{\ast }=x_7{e}^{-\epsilon };k_6^{\ast }=k_6{e}^{-\epsilon };k_8^{\ast }=k_8{e}^{\epsilon }.$$
• One scaling type symmetry is admitted using the parameter ${\rho }_{\mathrm{vol}}$. All the nucleus states, as well as four parameters, take part in the symmetry:

  $$x_i^{\ast }=x_i{e}^{\epsilon }i=5,6,9,10;s_2^{\ast }=s_2{e}^{-\epsilon };$$

  $$k_i^{\ast }=k_i{e}^{-\epsilon }i=6,10,11;{\rho }_{\mathrm{vol}}^{\ast }={\rho }_{\mathrm{vol}}{e}^{\epsilon }.$$
• The last symmetry is the only one that involves the initial condition parameter, $x_{1_0}$. All of the states have a scaling type symmetry, compensated by the scaling factor of $s_i$ and $k_{10}$:

  $$x_i^{\ast }=x_i{e}^{\epsilon }i=1,\dots ,10;s_i^{\ast }=s_i{e}^{-\epsilon }i=1,\dots ,4;$$

  $$k_{10}^{\ast }=k_{10}{e}^{-\epsilon };x_{1_0}^{\ast }=x_{1_0}{e}^{\epsilon }.$$

All of the symmetries are elementary transformations and it was not necessary to use Lie series. If ICS had not been considered, the symmetries would be elementary too.

This model was studied in [15]. The results of the first four transformations presented above coincide with those found in [15]; however, the last transformation includes the ICS parameter, unlike in the aforementioned article.

3.4. Glucose-Insulin Regulation

This model describes the regulation of blood glucose and insulin [24]. It has two states (glucose and insulin), one output (a glucose measurement), and a known input (the glucose entering from the digestive system):

$$\begin{array}{cc}\hfill \dot{q_1}& =u+p_1{q}_1-p_2{q}_2,\hfill \\ \hfill \dot{q_2}& =p_3{q}_2+p_4{q}_1,\hfill \\ \hfill y& ={\displaystyle \frac{q_1}{V_p}}.\hfill \end{array}$$

Here we analyse its symmetries without considering ICS. The model has two infinitesimal generators, which can be found using multivariate polynomials of second order:

$$\begin{array}{cc}\hfill X& =q_2{\displaystyle \frac{\partial }{\partial q_2}}-p_2{\displaystyle \frac{\partial }{\partial p_2}}+p_4{\displaystyle \frac{\partial }{\partial p_4}},\hfill \\ \hfill X& =-q_1{\displaystyle \frac{\partial }{\partial q_1}}+{\displaystyle \frac{u+p_3{q}_1}{p_2}}{\displaystyle \frac{\partial }{\partial q_2}}+p_3{\displaystyle \frac{\partial }{\partial p_1}}-p_2{\displaystyle \frac{\partial }{\partial p_2}}-p_3{\displaystyle \frac{\partial }{\partial p_3}}+\hfill \\ \hfill & +{\displaystyle \frac{p_1{p}_3-p_3^2+p_2{p}_4}{p_2}}{\displaystyle \frac{\partial }{\partial p_4}}-V_p{\displaystyle \frac{\partial }{\partial V_p}}.\hfill \end{array}$$

(52)

The first infinitesimal generator includes only elementary transformations:

$$\begin{array}{cc}\hfill q_2^{\ast }& =q_2{e}^{\epsilon },p_2^{\ast }=p_2{e}^{-\epsilon },p_4^{\ast }=p_4{e}^{\epsilon }.\hfill \end{array}$$

The second infinitesimal generator has seven transformations, four of which are elementary:

$$\begin{array}{cc}\hfill q_1^{\ast }& =q_1{e}^{-\epsilon },p_2^{\ast }=p_2{e}^{-\epsilon },p_3^{\ast }=p_3{e}^{-\epsilon },V_p^{\ast }=V_p{e}^{-\epsilon },\hfill \\ \hfill q_2^{\ast }& ={\displaystyle \frac{(u-p_3{q}_1){\epsilon }^4}{24p_2}}+{\displaystyle \frac{(u+p_3{q}_1){\epsilon }^3}{6p_2}}+{\displaystyle \frac{(u-p_3{q}_1){\epsilon }^2}{2p_2}}+{\displaystyle \frac{(u+p_3{q}_1)\epsilon }{p_2}}+q_2,\hfill \\ \hfill p_1^{\ast }& =-{\displaystyle \frac{p_3{\epsilon }^4}{24}}+{\displaystyle \frac{p_3{\epsilon }^3}{6}}-{\displaystyle \frac{p_3{\epsilon }^2}{2}}+p_3\epsilon +p_1,\hfill \\ \hfill p_4^{\ast }& ={\displaystyle \frac{(p_3^2+p_1{p}_3+p_2{p}_4){\epsilon }^4}{24p_2}}+{\displaystyle \frac{(-p_3^2+p_1{p}_3+p_2{p}_4){\epsilon }^3}{6p_2}}+\hfill \\ \hfill & +{\displaystyle \frac{(p_3^2+p_1{p}_3+p_2{p}_4){\epsilon }^2}{2p_2}}+{\displaystyle \frac{(-p_3^2+p_1{p}_3+p_2{p}_4)\epsilon }{p_2}}+p_4.\hfill \end{array}$$

Considering two terms of Lie series, the reparameterized model is FISPO, as classified by STRIKE-GOLDD.

The formulation of the IVP for the second infinitesimal generator is:

${\dot{q_1}}^{\ast }=-q_1^{\ast },$	$q_1^{\ast }\left(0\right)=q_1$ ,
${\dot{q_2}}^{\ast }={\displaystyle \frac{u+p_3^{\ast }{q}_1^{\ast }}{p_2^{\ast }}},$	$q_2^{\ast }\left(0\right)=q_2$ ,
${\dot{p_1}}^{\ast }=p_3^{\ast },$	$k_7^{\ast }\left(0\right)=k_7$ ,
${\dot{p_2}}^{\ast }=-p_2^{\ast },$	$p_2^{\ast }\left(0\right)=p_2$ ,
${\dot{p_3}}^{\ast }=-p_3^{\ast },$	$p_3^{\ast }\left(0\right)=p_3$ ,
${\dot{p_4}}^{\ast }={\displaystyle \frac{p_1^{\ast }{p}_3^{\ast }-{\left(p_3^{\ast }\right)}^2+p_2^{\ast }{p}_4^{\ast }}{p_2^{\ast }}},$	$p_4^{\ast }\left(0\right)=p_4$ .
${\dot{V_p}}^{\ast }=-V_p^{\ast },$	$V_p^{\ast }\left(0\right)=V_p$ .

The solution of the ODE system is:

$$\begin{array}{cc}\hfill q_1^{\ast }=& q_1{e}^{-\epsilon },\hfill \\ \hfill q_2^{\ast }=& q_2+{\displaystyle \frac{e^{-\epsilon }{p}_3{q}_1+e^{\epsilon }u}{p_2}},\hfill \\ \hfill p_1^{\ast }=& p_1+p_3-p_3{e}^{-\epsilon },\hfill \\ \hfill p_2^{\ast }=& p_2{e}^{-\epsilon },\hfill \\ \hfill p_3^{\ast }=& p_3{e}^{-\epsilon },\hfill \\ \hfill p_4^{\ast }=& {\displaystyle \frac{p_3^2{e}^{-\epsilon }}{p_2}}-{\displaystyle \frac{p_3^2+p_1{p}_3}{p_2}}+{\displaystyle \frac{e^{\epsilon }(p_1{p}_3+p_2{p}_4)}{p_2}},\hfill \\ \hfill V_p^{\ast }=& V_p{e}^{-\epsilon }.\hfill \end{array}$$

It is possible, using the procedure described for the pharmacokinetic model, to verify that the IVP solutions are the closed form of the solutions obtained through Lie series provided by the programme.

This model was proposed by Bolie [24] and its structural identifiability was analysed in [27] for the first time. However, its symmetries had not been studied until now.

4. Discussion

This article has addressed the relationship between non-observability and Lie algebra. Its main contribution is a computational method that searches for infinitesimal transformations in models composed of rational functions, in order to undo the symmetries that these may present. The procedure is based on expressing each transformation admitted by the ODE system according to its infinitesimal generator in polynomial form. In this way, the search for symmetries is equivalent to solving a system of linear equations, whose solution yields a transformation of the parameters that makes the model observable while leaving the observations invariant.

Our method builds on previous work [12,19,25], and especially on the procedure presented by Merkt et al. [15], with the addition of two features. The first one is the a priori maximization of the number of explicit transformations that can be obtained from the infinitesimal generator. The second one is the calculation of non-elementary transformations by means of Lie series. Increasing the number of explicit transformations is beneficial not only because it reduces the number of terms to consider from the Lie series, but also for calculating the solutions of the IVP. The complexity of the IVP solutions is inversely proportional to the number of explicit transformations. The pharmacokinetic model (PK) analysed here illustrates this point: without the use of the maximization of explicit transformations, MATLAB’s symbolic math toolbox did not manage to solve the IVP.

The algorithm allows to study the influence of the initial conditions in the model. The type of ICS (parametric, numeric or both) and the states that incorporate them may affect the number and type of symmetries of the models, varying from explicit to non-elementary transformations and reducing the number of infinitesimal generators.

We have implemented the method as a MATLAB programme that automates both the search for symmetries and the reconstruction of the model from the infinitesimal generators found. The programme has been integrated in the STRIKE-GOLDD toolbox for observability and identifiability analysis. The software has been tested with four previously published biomedical models, one of which—Bolie’s glucose-insulin regulation model—had not been tested for symmetries before. In the other cases our diagnoses mostly agree with those previously reported in the literature. An exception is the NF-$\kappa$B model, for which we found an infinitesimal generator that includes the parameter introduced by the initial conditions and that was not found in a previous analysis [15]. We observed another difference between our software and the one provided with [15]: when analysing the chemical reaction (CR) and pharmacokinetic (PK) models, the generators obtained with our code remained the same when varying the type of polynomial and degree; in contrast, the generators obtained with the programme of [15] changed when using partially varied and multivariate polynomials of order three or higher. These discrepancies may be due to implementation issues.

Our symmetry-detecting algorithm can be directly used to analyse structural identifiability and observability, providing an alternative to the OIC-checking algorithm already included in STRIKE-GOLDD for that purpose. More importantly however, this new code provides additional information about the relationships between model variables that cause loss of identifiability and/or observability. These insights can be exploited in two ways: (i) by fixing one or more parameters involved in a symmetry, in order to render the remaining ones identifiable, and (ii) by using the symmetry-breaking transformations to reformulate the model, yielding a modified model that is identifiable and observable. To facilitate the application of the latter procedure, we have implemented in our programme the semi-automatic transformation of a non-observable (respectively, non-identifiable) model into an observable (respectively, identifiable) model. It should be noted that, while said transformation may render a model fully observable, it also modifies the expression of the variables involved in its equations, which lose their original mechanistic meaning. Thus, while the results of the procedure can offer valuable insight about the model structure, they should be applied carefully for the purpose of model reformulation.

Our programme has some known limitations. First, while it considers high order generators and it can uncover a wide range of possible symmetries, it lacks procedures for determining a priori the type and total number of symmetries present in a model. Second, it does not provide a bound on the number of terms of the Lie series needed to obtain the infinitesimal transformations, when these are not given by explicit transformations. To the best of our knowledge, these limitations are shared with other existing methodologies. The possibility of overcoming them will be considered in future work.