Symmetries in Dynamic Models of Biological Systems: Mathematical Foundations and Implications

Abstract

Symmetries are ubiquitous in nature. Almost all organisms have some kind of “symmetry”, meaning that their shape does not change under some geometric transformation. This geometrical concept of symmetry is intuitive and easy to recognize. On the other hand, the behavior of many biological systems over time can be described with ordinary differential equations. These dynamic models may also possess “symmetries”, meaning that the time courses of some variables remain invariant under certain transformations. Unlike the previously mentioned symmetries, the ones present in dynamic models are not geometric, but infinitesimal transformations. These mathematical symmetries can be related to certain features of the system’s dynamic behavior, such as robustness or adaptation capabilities. However, they can also arise from questionable modeling choices, which may lead to non-identifiability and non-observability. This paper provides an overview of the types of symmetries that appear in dynamic models, the mathematical tools available for their analyses, the ways in which they are related to system properties, and the implications for biological modeling.

1. Introduction

Symmetrical features are abundant in the biological world. Practically all organisms have some kind of symmetry in three-dimensional space. Due to developmental noise, these symmetries are approximate; the deviation from perfect symmetry is called fluctuating asymmetry [1]. Such geometrical symmetries can often be easily noticed.

On the other hand, symmetries may also appear in equations that do not represent a geometric shape. For example, an equation may contain two parameters that have exactly the same influence in its outcome, e.g., if they always multiply each other. More generally, these symmetries refer to an invariance of some variables under certain transformations of other variables.

The dynamic behavior of many biological processes can be accurately described by systems of differential equations, which are often ordinary differential equations (ODEs). Such systems of ODEs, along with the definition of the inputs and outputs, constitutes a dynamic model [2]. It is often the case that said dynamic models of biosystems have symmetries, whose significance may vary depending on the context. For example, for the purpose of system identification, if two parameters can be interchanged without affecting the model output, this means that they cannot be distinguished by observing (measuring) the output—and thus, they are unidentifiable. However, from a different point of view, this feature may be interpreted as biological robustness—variations in the value of one parameter, due for example to environmental fluctuations, can be compensated by variations in other parameter, leaving the output invariant. The former view is relevant for model building and the latter for the system behavior, but both of them are connected.

The symmetries in dynamic models are less obviously noticed than symmetric shapes. They can be of different types, and some of them defy intuition—but they can nevertheless play an important role. Fortunately, there exists a well established theory for its study. The symmetry analysis of differential equations was founded in the nineteenth century by the Norwegian mathematician Sophus Lie, and there are many review articles [3] and books [4,5] on the topic of Lie groups of transformations.

The purpose of the present review is threefold. First, it presents the fundamental concepts of Lie symmetries and related mathematical tools (Section 2). Second, it describes a number of features found in dynamic models of biological systems, and their connections to symmetries, as well as the biological and modeling implications (Section 3). Third, it discusses some subtleties and relationships between properties, providing some remarks and lines for future research (Section 4).

2. Mathematical Concepts and Analysis Tools

We consider ODE models of the general form:

$$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)

where $x\left(t\right)\in {\mathbb{R}}^{n_x}$ is the state vector, $\theta \in {\mathbb{R}}^{n_{\theta }}$ is the parameter vector, $u\left(t\right)\in {\mathbb{R}}^{n_u}$ is the known input vector, $w\left(t\right)\in {\mathbb{R}}^{n_w}$ is the unknown input vector, $y\left(t\right)\in {\mathbb{R}}^{n_y}$ is the output vector, and f and g are vectors of (possibly nonlinear) analytical functions.

To simplify the notation, we will usually omit the dependency on time. Furthermore, despite the fact that the distinction made in (1) between states, parameters, and inputs is the usual practice in dynamic modeling, we will group them together in an augmented state vector $x\in {\mathbb{R}}^n$ with dimension $n=n_x+n_u+n_w+n_{\theta }$. The reason for this is that all those variables can be involved in the symmetries that we consider. Thus, the Equation (1) may be rewritten as:

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

(2)

It should be noted that the simplicity of this notation conceals certain subtleties related to the inputs, both known (u) and unknown (w). Since the aim of this Section is not to describe those aspects in depth, but instead to provide an overview of the fundamental math, we will proceed with this abuse of notation. The reader interested in a more rigorous treatment is referred to, e.g., [6,7,8,9].

The remainder of this section discusses concepts related to symmetries in models such as (2). The definitions in Section 2.1 are taken from [3,4,5], and the procedure in Section 2.2 mostly follows [6].

2.1. Lie Symmetries: Definitions

Definition 1.

A transformation of variables $x^{*}=X\left(x\right)$ is called a “symmetry transformation”, a “Lie symmetry”, or a “point symmetry” of a system of differential equations (2) if

$$\begin{array}{cc}\hfill \dot{x^{*}}& =f\left(x^{*}\right),\hfill \\ \hfill y& =g\left(x^{*}\right),\hfill \\ \hfill x^{*}{|}_{t=t_0}& =x_0^{*}.\hfill \end{array}$$

A symmetry transformation maps the set of all solutions of (2) into itself. If the equations admit such a symmetry, the model variables can be transformed by it without altering the dynamics nor the output.

A closely related concept is equivariance, which can be defined as follows [10]:

Definition 2.

Given a group of operators Γ acting on $\mathbb{R}$, a vector field $f:{\mathbb{R}}^n\times \mathbb{R}\to {\mathbb{R}}^n$ is said to be “Γ invariant” if $f(\gamma x,t)=\gamma f(x,t)$ for all $\gamma \in \Gamma$ and $x\in {\mathbb{R}}^n$.

Note that $\gamma$ can be a Lie group, and it does not need to be a linear operator. If f from (2) is $\Gamma$ invariant, then $\gamma$ is a symmetry of (2).

Definition 3.

A “one-parameter Lie group of transformations” is a mapping of the form

$$x^{*}=X(x,\epsilon ),$$

(3)

where ε is the parameter, if it constitutes a group — i.e., it satisfies the four group axioms: (1) closure, (2) associativity, and existence of (3) an identity and (4) an inverse—and, additionally, the transformation is (5) smooth in x and (6) analytic in ε, and (7) the identity element can be chosen as$\epsilon =0$.

Note that, expanding (3) in a neighborhood of $\epsilon =0$, we obtain the following expression for a Lie group:

$$x^{*}=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)+...=x+\epsilon \left({\displaystyle \frac{\partial X(x,\epsilon )}{\partial \epsilon }}{|}_{\epsilon =0}\right)+O\left({\epsilon }^2\right),$$

(4)

where the following term appears:

Definition 4.

The “infinitesimal” of a Lie group (3) is:

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

(5)

Definition 5.

The “infinitesimal transformation” of a Lie group (3) is:

$$x^{*}=x+\epsilon ·\eta \left(x\right)+O\left({\epsilon }^2\right)$$

(6)

Definition 6.

The “infinitesimal generator” of a Lie group (3) is:

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

(7)

The infinitesimal generator X is an equivalent representation of a Lie group. Using X, the infinitesimal transformation (6) can be written as

$$x^{*}=X(x,\epsilon )=exp\left[\epsilon X\right]x$$

(8)

In order to calculate the transformation of $\dot{x}$, we can extend the definition of X so that (8) acts on the derivatives:

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

(9)

where

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

(10)

Types of Lie symmetries. We can classify symmetries according to their expressions. Below we list some of the most common ones, along with their infinitesimal generators.

Translation:

$$x_i^{*}=x_i+\epsilon ,X=\frac{\partial }{\partial x_i}$$

(11)

Scaling:

$$x_i^{*}=e^{\epsilon }{x}_i,X=x_i\frac{\partial }{\partial x_i}$$

(12)

Moebius:

$$x_i^{*}=\frac{x_i}{1-\epsilon x_i},X=x_i^2\frac{\partial }{\partial x_i}$$

(13)

Higher order:

$$x_i^{*}=\frac{x_i}{{[1-(p-1)\epsilon x_i^{p-1}]}^{\frac{1}{p-1}}},X=x_i^p\frac{\partial }{\partial x_i}$$

(14)

The symmetries most commonly found in dynamic models of biological systems are traslation and scaling.

2.2. Finding Lie Symmetries

Lie symmetries can be found by solving a system of differential equations, whose solutions are the expressions of the transformations admitted by the model. To formulate the problem, we use the following result:

Theorem 1.

An ODE model of the form (2) admits the one-parameter Lie group of transformations defined by the infinitesimal generator (7) if and only if [3,4]

$$\begin{array}{cc}\hfill X^{\prime }·({\dot{x}}_k-f_k\left(x\right))=0,& k=1,...,n\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill X·(y_l-g_l\left(x\right))=0,& l=1,...,n_y\hfill \end{array}$$

(16)

where $X^{\prime }$ was defined in (9) and (10).

Application of Theorem 1 to the equations of (2) yields the following set of partial differential equations (PDEs):

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

(17)

By solving these equations we obtain the infinitesimals, ${\eta }_i\left(x\right)$, and by replacing them in (6), we calculate the transformations that we are looking for.

Unfortunately, in the general case it is not possible to solve equations (17) analytically. A way of simplifying the problem is to assume that $f_k\left(x\right)$ and $g_l\left(x\right)$ are rational functions, that is:

$$\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,...,n,\hfill \\ \hfill y_l& =g_l\left(x\right)={\displaystyle \frac{R^l\left(x\right)}{S^l\left(x\right)}},l=1,...,n_y,\hfill \end{array}$$

(18)

where $P^k\left(x\right)$, $Q^k\left(x\right)$, $R^l\left(x\right)$, and $S^l\left(x\right)$ are polynomial expressions. Most dynamic models encountered in biology are of this form. If we replace (18) in (17), the PDEs are transformed into ODEs. Next, we introduce an assumption about the shape of the solutions (or ansatz), which amounts to restricting the search to certain classes of symmetries. To choose a suitable ansatz, we may note that the infinitesimals of the generators (11)–(14) are polynomial expressions. Thus, a possible ansatz could be:

$${\eta }_i\left(x\right)=\sum _{d_i=0}^{d_{max}}{r}_{i,d}{x}_1^{d_1}{x}_2^{d_2}···x_n^{d_n},i=1,...,n.$$

(19)

Note, however, that in this expression every variable is allowed to depend on any other variable. In reality, parameters (which are constant) cannot have symmetries depending on states or inputs (which can be time-varying), so the corresponding coefficients must be zero. Using the ansatz above, after some modifications Equation (17) become:

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

(20)

When the infinitesimals are polynomials as in (19), the equations (20) are of the form:

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

(21)

where r is the vector of all $r_{i,d}$, and the coefficients $c_{i_1,...,i_n}$ are linear in r. Equation (21) holds if the following condition holds:

$$C·r=0$$

(22)

where all the coefficients have been grouped in a matrix C. Thus, if a model has symmetries defined by (19), said symmetries can be found by solving a system of linear equations, which amounts to computing the kernel of C.

2.3. Structural Identifiability and Observability

Identifiability is a mathematical property that describes the ability to infer the parameters in a model from observations of its output vector over time. It is common to distinguish between structural and practical identifiability; the former is determined exclusively by the model ODEs and output equations, while the latter is also affected by the quantity and quality of the available data. Thus, structural identifiability is necessary but not sufficient for practical identifiability. A recent review of both concepts is [11].

Structural identifiability is related to another property, observability, which refers to states instead of parameters [12]. These concepts are defined in the following paragraphs.

Definition 7.

A parameter ${\theta }_i$ (that is, the ith element of the parameter vector) is structurally locally identifiable (SLI) if, for almost any parameter vector ${\theta }^{*}\in {\mathbb{R}}^{n_{\theta }}$, there is a neighbourhood $\mathcal{N}\left({\theta }^{*}\right)$ where the following condition holds [2]:

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

(23)

If (23) is not true in any neighborhood of ${\theta }^{*}$, the parameter ${\theta }_i$ is structurally unidentifiable (SU). If (23) is true for all model parameters, the model is said to be SLI as well, and SU otherwise.

Definition 8.

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

A model is observable if all its states are observable.

We will use the acronym SIO to refer to Structural Identifiability and Observability jointly.

3. Connections: Symmetries and Other Properties

3.1. Symmetries and SIO

The presence of symmetries in the equations of a model is a possible source of structural non-identifiability and non-observability [13]. It has been shown that symmetries allow for similarity transformations, i.e., transformations of parameters and states that do not change the model output [14,15,16]. If that is the case, it is not possible to infer the values of the parameters and states involved in such symmetries by observing the output. Therefore, those parameters are unidentifiable and those states are unobservable. Since structural identifiability and observability are desirable properties for the purpose of system identification, it may be useful to break the symmetries that prevent them. Symmetries can be removed by reparameterizing the model equations, or by measuring additional outputs when possible. Several efforts have been made in this regard.

Meshkat and Sullivant addressed this problem for a class of linear compartment models, which are used in systems biology and pharmacokinetics [17]. They provided partial results for determining which models allow for scaling reparameterizations that make them identifiable, and an algorithm for finding them. A different solution was sought in [18], where an algorithm for finding the minimal output sets for obtaining structural identifiability (i.e., which subsets of state variables must be measured) was presented. The method relies on finding translation, scaling, and affine symmetries, and then breaking these symmetries by adding states to the output vector. A Mathematica implementation of the identifiability analysis method is available. Exact simplification methods based on Lie symmetries, including scaling and affine transformations, were presented in [19,20,21] and implemented in the Maple package MABSys.

An algorithm capable of detecting higher-order Lie symmetries, along with its Python implementation, was provided in [6], where the implications for predictability of the dynamic model were also discussed. A related algorithm, which finds symmetries and their symmetry-breaking transformations, was presented in [9] and included in the Matlab toolbox STRIKE-GOLDD [22]. Shi and Chatzis studied this problem in the context of monitoring the state of mechanical structures, for which it is important to consider unknown inputs. Initially, they assumed that the symmetries are translations or scalings [7]; their results were later extended so as to relax those assumptions [8]. A computational method called AutoRepar was presented in [23] as a follow-up to [9]. AutoRepar extends the preceding algorithm by calculating the reparameterizations associated to the symmetry breaking transformations, thus rendering a model structurally identifiable and observable.

There are many software packages available for finding symmetries in dynamic models. A few of them are focused on structural identifiability, some of which are listed in

Table 1. Some methods for finding symmetries in dynamic models and analysing their SIO.

Name	Language	Reference	Types of Detected Symmetries
MABSys	Maple	[26]	scaling
MinimalOutputSets	Mathematica	[18]	translation, scaling, affine
symmetryDetection	Python	[6]	up to higher order
AutoRepar	Matlab	[23]	up to higher order
SIM	–	[24]	scaling

. The first four rows in Table 1 list algorithms that have computational implementations in several symbolic environments. Software implementations are convenient due to the complexity of these methods, which makes computations by hand infeasible for all but the simplest models. In contrast, the last entry on the table, the Scaling Invariance Method (SIM), is a simple technique that can be applied manually [24]. It detects scaling symmetries, which are arguably the most frequently found in biological models. However, it must be taken into account that, if a model contains other types of symmetries, the results of SIM do not inform about structural identifiability [25].

3.2. Symmetries and Biological Robustness

The term robustness refers to the ability of a biological system to maintain its behavior in varying environmental conditions. It is a desirable property, since organisms must be able to tolerate uncertainty—both in external inputs and in internal components—in order to survive. From a modeling point of view, one way of achieving robustness is insensitivity. A model may be largely insensitive to changes in an external input or parameter [27]. However, biological systems also need to be sensitive to certain signals or variations, so as to sense the environment and adapt to changes in it. Thus, there must be a balance between sensitivity and robustness. Organisms must change their behavior in response to certain (types of) changes, and keep it unaltered in response to others. A number of mechanisms exist which can achieve this goal.

The concept of fold-change detection (FCD) describes the behavior of sensory systems whose response depends only on relative changes to the input, instead of on its absolute values [28] (that is, the output is the same when the input changes from 1 to 2 and when it changes from 3 to 6). FCD entails two properties: exact adaptation (i.e., the steady state is the same) and the Weber-Fechner law (i.e., the maximal initial response is the same), but it is more general than both of them: FCD concerns the whole dynamic response. Examples of biological systems that exhibit FCD include bacterial chemotaxis, vision, and cellular signaling. A few mechanisms have been shown to produce FCD, including incoherent feed-forward loops (IFFL) and nonlinear integral feedback (NLIFL), which are shown in Figure 1. FCD is a scale invariance: it corresponds to the existence of a scaling symmetry involving input and state [29]. FCD is a particular case of symmetry invariance; Shoval et al. [30] characterized systems whose transient behaviors are invariant under a set of symmetries.

More recently, a related concept, dynamical compensation (DC), was introduced [31]. It was presented as a design principle that provides robustness to physiological circuits, in the face of uncertainties in certain parameters. It was shown by Sontag [32] that DC and FCD are manifestations of the same phenomenon.

The interplay between symmetries and stability of dynamical systems was discussed in [10]. Both concepts were linked by means of a third property, contraction; an approach that is able to explain FCD. It was also shown how the invariance under input scaling can be exploited for network synchronization, thus playing a key role in biological applications such as quorum-sensing networks in bacteria.

4. Discussion and Conclusions

In the preceding section two different viewpoints have been discussed. On the one hand, Section 3.1 has noted that the presence of symmetries in dynamic models can lead to non-identifiability and non-observability, which hampers system identification and is a possible source of wrong predictions (that is bad!). On the other hand, Section 3.2 has presented symmetries as providers of robust behavior, which allow for adaptation to changes in external and internal conditions (that is good!).

Albeit conflicting, these two viewpoints are not mutually exclusive and they can be subtly related. Here we illustrate them with a case study. Figure 2A shows a graphical representation of a plausible physiological mechanism for regulating plasma glucose levels. An external input representing ingestion (u) increases glucose concentration (G), and indirectly influences the $\beta$ cells’ mass ($\beta$). In turn, $\beta$ cells secrete insulin (I), and insulin decreases glucose concentration. In this way, plasma glucose concentration tends to be pushed back to a reference level after each meal, as can be seen in Figure 2B. The model ODEs are:

$$\begin{array}{cc}\hfill \dot{G}& =u-(c+s_{i}I)G,\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill \dot{\beta }& =\beta ({\lambda }^{+}-{\lambda }^{-}),\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill \dot{I}& =p\beta {\rho }_G-\gamma I,\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill y& =G,\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\rho }_G& ={\displaystyle \frac{G^2}{{\alpha }^2+G^2}},{\lambda }^{+}={\displaystyle \frac{1.4583·{10}^{-5}}{1+{\left({\displaystyle \frac{8.4}{G}}\right)}^{1.7}}},{\lambda }^{-}={\displaystyle \frac{1.7361·{10}^{-5}}{1+{\left({\displaystyle \frac{G}{4.8}}\right)}^{8.5}}}\hfill \end{array}$$

(28)

This model was presented in [31] to illustrate the DC property, which in this case can be summarized as follows: if the rate at which insulin removes glucose from the system (which is given by parameter $s_i$) changes, the system goes through an adaptation process which eventually leads to the glucose response being identical to the pre-change one. This is shown in Figure 2B, where it can be seen that the two time courses of glucose are indistinguishable. That is, the model exhibits the DC property, which, as explained in Section 3.2, is related to the existence of a scaling symmetry. Indeed, it can be seen [23] that the model admits the following transformations:

$$\begin{array}{cc}\hfill \tilde{\beta }& =\beta ·s_i·p\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill \tilde{I}& =I·s_i\hfill \end{array}$$

(30)

which leave the model output, $y=G$, invariant. The resulting equations are:

$$\begin{array}{cc}\hfill \dot{G}& =u-(c+\tilde{I})G,\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill \dot{\tilde{\beta }}& =\tilde{\beta }({\lambda }^{+}-{\lambda }^{-}),\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill \dot{\tilde{I}}& =\tilde{\beta }{\rho }_G-\gamma \tilde{I}\hfill \end{array}$$

(33)

The mechanism that gives rise to these symmetries is a feedforward loop, shown in green in Figure 2A, that resembles the IFFL depicted in Figure 1A.

Let us now adopt the other viewpoint. It was noted in [32,33,34] that, at least according to its initial definition, DC amounts to non-identifiability. Indeed, the two parameters involved in the symmetries, $s_i$ and p, are structurally unidentifiable. This means that their values cannot be determined by measuring the system output (glucose). In system identification terms, this model is over-parameterized, since it contains more unknowns that can be inferred from data. In contrast, in the reparameterized equations given by (31)–(33) $s_i$ and p have been removed, and the model is identifiable.

Importantly, the non-identifiability of the model given by (24)–(28) also entails non-observability: insulin is an non-observable state variable, and hence it cannot be estimated from output measurements. This is shown in Figure 2C, where it can be seen that the two values of $s_i$ yield two different time courses of I; since it is not possible to know $s_i$ from data, it is also impossible to determine I. This limitation may be problematic, e.g., if we want to implement an artificial pancreas controller that requires knowledge of I. In this regard, it should be noted that the reparameterization given by (29)–(33) transforms I, so even though the transformed variable becomes observable, it no longer corresponds to insulin concentration. Fortunately, the transformation given by (24)–(28) is not the only admissible one. As was discussed in [23], if the input u is considered unknown it is possible to lump other terms in it, obtaining a model reformulation that leaves not only G, but also I invariant. To this end, we perform the following transformations:

$$\begin{array}{cc}\hfill u^{*}& =u+GI(-s_i+1),\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill \tilde{\beta }& =\beta p,\hfill \end{array}$$

(35)

obtaining the following reparameterized model:

$$\begin{array}{cc}\hfill \dot{G}& =u^{*}-(c+I)G,\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill \dot{\tilde{\beta }}& =\tilde{\beta }({\lambda }^{+}-{\lambda }^{-}),\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill \dot{I}& =\tilde{\beta }{\rho }_G-\gamma I.\hfill \end{array}$$

(38)

The model above would be suitable for implementing the previously mentioned type of controller.

It should be stressed that structural non-identifiability and/or non-observability should not be invariably interpreted as signs of biological robustness. While these features may sometimes be related to adaptation mechanisms, they always reveal deficiencies in the model formulation for the purpose of system identification and monitoring. As the example in the present Section has shown, the absence of SIO hampers the information that a model can provide, and, if unnoticed, it can lead to wrong predictions and conclusions. Hence, it is essential to perform a proper SIO analysis and to address any issues found by it. In this regard, we remark that a biosystem’s adaptation capabilities (such as, e.g., the DC property) do not necessarily entail lack of SIO; indeed, it is possible to use alternative definitions of DC that preserve its intended mechanistic meaning while avoiding the implication of non-identifiability, as was proposed in [34].

The example discussed above contains scaling symmetries. While they are arguably the most commonly found type of symmetry, dynamic models of biological systems may admit symmetry transformations of other types, instead or in addition to scaling. Since the absence of scaling symmetries does not exclude the existence of other types, caution must be exercised when interpreting the results of analyses based exclusively on scalings [25].

The study of symmetries in dynamic models of biosystems presents multiple opportunities for future research. The extent to which symmetry invariances such as FCD are present in different biological contexts is still an open question. So is the determination of all the mechanisms capable of producing said behavior, as well as the comparison of their relative advantages. The implementation of these mechanisms in synthetic biology circuits is a challenging problem, which is likely to inspire many developments. Other avenues for future research are discussed in [28].

From the mathematical point of view, the study of Lie symmetries has a long and rich history. However, there is still a need for computationally efficient algorithms for finding symmetries, especially of higher order and non-conventional types. Likewise, their systematic exploitation for the analysis of identifiability, observability, and other control-theoretic properties such as stability and controllability still poses many challenges.

In order to succeed, these research efforts will necessitate a truly interdisciplinary collaboration between biologists, mathematicians, and engineers.