Higher-order Lie symmetries in identifiability and predictability analysis of dynamic models

Parameter estimation in ordinary differential equations (ODEs) has manifold applications not only in physics but also in the life sciences. When estimating the ODE parameters from experimentally observed data, the modeler is frequently concerned with the question of parameter identifiability. The source of parameter nonidentifiability is tightly related to Lie group symmetries. In the present work, we establish a direct search algorithm for the determination of admitted Lie group symmetries. We clarify the relationship between admitted symmetries and parameter nonidentifiability. The proposed algorithm is applied to illustrative toy models as well as a data-based ODE model of the NFκB signaling pathway. We find that besides translations and scaling transformations also higher-order transformations play a role. Enabled by the knowledge about the explicit underlying symmetry transformations, we show how models with nonidentifiable parameters can still be employed to make reliable predictions.


I. INTRODUCTION
Modeling of dynamic systems by ordinary differential equations (ODEs) has always been concerned with the problem of unknown parameters in the model equations.As opposed to other fields, in the life sciences a priori unknown parameters like reaction velocities, association rates, dissociation rates, etc., can often not be measured directly.The system under investigation is only observable as a whole and cannot be dissected into isolated reactions.Therefore, the parameters need to be determined by parameter estimation, which relies on a certain amount of time-resolved measurements for at least a subset of system components.
In this procedure, the modeler is frequently confronted with the problem of nonidentifiable parameters, i.e., given the observed data there exists a submanifold in parameter space describing, in terms of an objective value, the data equally or almost equally well: these cases are denoted as structural or practical nonidentifiability [1].
In Ref. [2], a data-based approach is utilized to identify submanifolds of equal objective values, thus being able to detect nonidentifiabilities.A complementary data-based approach employing the profile likelihood is presented in Ref. [1], giving a precise definition of structural and practical identifiability in terms of parameter profiles.On the other hand, a priori data-independent analytical methods can be applied to directly test the model equations for structural nonidentifiability.The authors of [3,4] introduced two methods which can be applied to study the identifiability of ODE models.However, these do not yield algorithms to automatically detect structural nonidentifiabilities.The approach presented in Refs.[5][6][7] relies on differential algebra methods and can be employed to study the input-output relation between the parameters and the predicted observation.If a unique solution exists, the model is structurally identifiable.This approach is purely algebraic.However, as the method relies on Buchberger's algorithm to solve the equations, its complexity grows rapidly with the size of the model and quickly becomes infeasible.Additionally, the model is expected to be in minimal form, which is often not fulfilled.In such cases, the model has to be manually transformed before the analysis can be applied.
A drawback of all the above techniques is that they do not elucidate the fundamental source of the structural nonidentifiabilities. in Ref. [8] it was shown that the admittance of Lie transformations is equivalent to the existence of similarity transformations [4] and thus equivalent to structural nonidentifiability of the model.To concretize this statement, let ẋ(t) = f (x(t),θ dyn ,u(t)), x(0) = θ ini , y(t) = g(x(t),θ obs ) (1) be a dynamic system with states x ∈ R m , dynamic parameters θ dyn ∈ R p , inputs u ∈ R q , and initial conditions θ ini ∈ R m .In general, the observation y ∈ R n * is an arbitrary function of the states x and observational parameters θ obs ∈ R p .Then, an admitted Lie symmetry is a continuous group of transformations of the states x and parameters = (θ dyn ,θ ini ,θ obs ) ∈ R p ,p = m + p + p : : R m+p+q × G −→ R m+p+q (x, ,u; ) −→ (x * , * ,u * ), such that the observation y is unchanged for each element in the Lie group G. I.e., it must hold that where x(t) and x * (t) implicitly depend on the parameters (θ ini ,θ dyn ) and (θ * ini ,θ * dyn ).This illustrates the relation between Lie groups of transformations and structural nonidentifiability: Although the parameters are changed by the transformation , the observation stays invariant.Equation (2) can also be expressed using the infinitesimal generators X of , forming the Lie algebra g: The advantage of the latter formulation is that it yields a condition for admitted symmetries which can be exploited to explicitly calculate the transformations.
Generally, the challenge of finding symmetries is approached by two methods.The first one, introduced by [9], is based on differential algebra combined with probabilistic sampling.The sampling is employed to investigate the rank of the variational system of the ODEs in order to identify the set of nonidentifiable parameters.This identification is computationally efficient as it is polynomial in time.Furthermore, having identified the subset of nonidentifiable parameters, it is possible to determine the associated Lie algebra generators from the maximal singular minor of the Jacobian.A Mathematica implementation of the algorithm was presented in Ref. [10].
The second approach is a more direct method employing an explicit polynomial Ansatz for the generator.This leads to a determining system for the coefficients which specify the Ansatz.The authors of [11,12] applied this approach to compute the reduced system of ODEs, especially for scaling and Möbius transformations, and to define coordinates leaving the steady states of a dynamic system invariant.Translational and scaling symmetries were also investigated by [13], pointing out the close relation between nondimensionalization and Lie group symmetries.Both groups work with Maple implementations.Translation, scaling, and affine transformations were also studied by [14] and used to compute a minimal set of observations which make the model identifiable.In Ref. [8] the determining system is solved by methods available in Mathematica with applications to observation functions with scalings.
The approach presented in this work has the same starting point as the methods used in Refs.[8,[11][12][13][14]; however, it differs in three major aspects.First, instead of using random specializations, we develop a method to solve the determining equations by directly transforming the rational expressions to a linear system.Second, by applying the proposed algorithm [15] to three examples, we demonstrate the significance of higher-order generators which to our knowledge were not discussed before.Third, we show that the knowledge of the explicit underlying symmetries that we gain by the proposed algorithm enables us to determine reliable model predictions although model parameters are nonidentifiable.
The work is structured as follows: Sec.II recapitulates Lie group theory and illustrates the algorithm as being implemented in Python with the libraries SymPy [16] for symbolic mathematics and SciPy and NumPy for numerical calculations.Section III introduces several examples demonstrating the relevance of higher-order transformations.In this context, we also study a data-based ODE model of the NFκB signaling pathway.Besides determining and discussing the intrinsic symmetries of the dynamic system, we infer quantities which can be reliably predicted despite the structural nonidentifiabilities in the model.

II. METHODS
For a simplified notation, Eq. ( 1) is rewritten in a form where the input u(t) and parameters θ are absorbed in the state vector x(t): In the following, we call n = m + p + q the total number of states.Hereby, the ODE model becomes If the biological system is in steady state at the start of the experiment, this has to be accounted for by the additional condition where x 0 = x(t = 0), which introduces an additional constraint for the parameters.

A. Theory of Lie groups of transformations
We adopt the definition of [17] for a Lie group of transformations.
Definition 1. (i) For all ∈ S, the mapping (•, ) is one to one.
(ii) S together with φ forms a group with identity element zero.
(v) S is an interval in R and φ is analytic.(vi) (x, ) is smooth in x and analytic in .A Lie group of transformations is said to be admitted by ODE model (4) if for every solution x(t) and every ∈ S the transformed function x * (t) = (x(t), ) is also a solution of the system.
For every Lie group of transformations, the infinitesimal generator is a differential operator defined as . The functions η i (x) are called infinitesimals of the transformation.
Theorem 1 [17].By utilizing the infinitesimal generator, the corresponding transformation can be written as Hence, the infinitesimal generator is an equivalent representation of a Lie group of transformations.Each component of the state vector x(t) in system (4) is a function of time.Therefore, a transformation of x(t) also induces a transformation of the time derivatives ẋ(t).In order to compute the explicit transformation of ẋ(t), the infinitesimal generator is extended such that the transformation (6) also acts on the derivatives: The transformation of the derivatives is now represented by the infinitesimals η i (x) which are calculated by ẋj ∂η i ∂x j (x) [18].This prolonged infinitesimal generator enables the formulation of an explicit criterion for admitted Lie groups of transformations.
Theorem 2 [17,18].The system of ordinary differential equations (4) admits a Lie group of transformations defined by the infinitesimal generator whenever ẋ = f (x) and g(x) = 0.The steady state equation ( 5) imposes additional constraints for the parameter values.Therefore, it could, in principle, also restrict the admitted symmetries.This has to be accounted for by the third symmetry condition: However, it can be shown (see Appendix A, Theorem A1) that this is fulfilled whenever Eq. ( 8) holds and, consequently, the steady state equation ( 5) does not require additional considerations.If the initial concentrations are given by general algebraic expressions, these can be treated similarly to the observation functions.For details see [14].
If X i and X j are two infinitesimal generators satisfying Eq. ( 8), then also every linear combination μ i X i + μ j X j is admitted by the system of ODEs.Therefore, its solution set defines a linear vector space of admitted generators.Furthermore, it can be shown [18] that the Lie bracket also yields an infinitesimal generator which is admitted.The vector space of admitted generators together with the Lie bracket forms a Lie algebra, denoted by g.

Rational differential equations
In the following, a method for solving the determining equations ( 8) is derived.Theorem 2 states that a system of the form of Eq. ( 4) admits a Lie group of transformations x = (x, ) if and only if holds for all solutions x(t) of the system where η i (x) are the infinitesimals of the transformations.Equation (9) defines a system of coupled partial differential equations for these infinitesimals.If a solution to this system is found, the admitted transformations can be computed by Eq. ( 6).However, Eq. ( 9) cannot be solved analytically for general differential equations.Therefore, in the following, we consider rational differential equations and observation functions of the form where , and S k (x) are multivariate polynomial expressions.Hereby, a large class of ODE models is covered.In particular, mass action, Michaelis-Menten, and Hill kinetics are included.The derivatives of the input u(t) are generally unknown.Therefore, the terms in Eq. ( 9) containing u have to vanish and, consequently, it must hold that ∂η k ∂x j (x) = 0 for k = 1, . . .,m and j = m + p + 1, . . .,n.Taken together, if Eq. ( 10) is plugged into Eq.( 9), the conditions become where the subscripts x i denote partial derivatives with respect to that variable.Note that, due to the structure of these conditions, whenever a set of infinitesimals η i ,i = 1, . . .,n is admitted, also the infinitesimals η i = θ j η i ,i = 1, . . .,n are admitted.The next step is to choose an Ansatz for the infinitesimals η i (x),i = 1, . . .,n.This corresponds to a restriction to a subclass of symmetries the admitted transformations are sought in.

Polynomial generators
Many elementary transformations are generated by infinitesimals which are simple powers of the variable.Examples are Although the structures of these common transformations are rather different, their infinitesimals are given by similar expressions of the form x d i ,d ∈ N.This motivates the polynomial Ansatz The goal is then to determine the coefficients r i,d for which the corresponding transformation is admitted.If plugged into Eq.(11a), the first sum reduces to a single term because ∂η k ∂x j = 0 for j = k.Consequently, the conditions can be multiplied by the denominators Q k 2 and S k 2 , respectively, to obtain With Ansatz (12), these equations define polynomial expressions which can always be rearranged into the form where r is a vector containing all parameters r i,d and the coefficients c i 1 ,...,i n are linear in r.Since the polynomial vanishes if and only if all its coefficients vanish, condition ( 9) is equivalent to the linear system where the entries of the matrix C are computed by the differential equations and observation functions via Eq.( 13).
Consequently, the problem of finding structural nonidentifiabilities in an ODE model was turned into the problem of solving a linear system of size N r = n(d max + 1).However, there can be admitted transformations which are not covered by Ansatz (12).In order to also detect these symmetries using the presented approach, a more general Ansatz is required.The intuitive generalization of the infinitesimals is to allow multivariate polynomials where the infinitesimal of state x i can also depend on other variables x j =i .Hereby, the most general polynomial Ansatz is given by The necessary calculations for solving Eq. ( 9) are analogous to the case of univariate polynomial generators.For details, the reader is referred to Appendix B.
In summary, we have derived an algorithm to automatically and efficiently search for admitted symmetries in rational ODE models.Unlike the methods used in Refs.[11][12][13][14], no random specializations of the determining equations ( 8) and ( 11) are utilized but they are converted to a linear system which can be solved directly.Once the form of the infinitesimal generator is chosen and the expressions are combined into the form of Eq. ( 13), the algorithm performs two major steps.First, the terms in the polynomials are rearranged and combined such that the coefficients which are required to vanish can be extracted.Although these calculations are a priori symbolic, they can be performed in a linear vector space because each term in the polynomials is fully defined by its coefficients and exponents.Because the coefficients are always linear in the parameters, both can be represented by a real vector.This step can be performed independently for each differential equation and observation function and, therefore, allows for parallel evaluation.In the second step, the system of linear equations as being defined by the coefficients is solved.This yields the set of infinitesimals which generate admitted transformations.Note that these infinitesimals only generate a basis and every linear combination is also admitted.The solution for the parameters is then plugged into the infinitesimal Ansatz which, in turn, fully defines the admitted transformations.The computational complexity of the two algorithm steps is bounded by O(N 2 r ) and O(N 3 r ), respectively, where N r is the total number of constants in the Ansatz.Further details on the computational complexity with univariate and multivariate generators are presented in Appendix C.

III. APPLICATION
In the following, the algorithm introduced above is applied to example ODE models highlighting different aspects of the general theory.In particular, the significance of higher-order generators is demonstrated.For details on how these examples are analyzed with our Python implementation see [15].

As a simple example, consider the chemical reaction
We assume that A is measurable on a relative scale and the observation saturates for high values: If the univariate polynomial Ansatz (12) is used to solve the symmetry condition (13) for this system, two distinct solutions are found.

Solution 1
The first solution is given by the infinitesimals which results in the infinitesimal generator It corresponds to a scaling of species A by together with a reversed scaling of the parameters The scaled species A * satisfies the same ODE because all the scaling factors e cancel out in the differential equation (14).Furthermore, in the observation function (15), the scaling is compensated by the transformation of the parameters s 1 and s 2 leaving A obs invariant.The admittance of this transformation expresses the fact that the unit in which A is measured is not determined due to the relative measurement.

Solution 2
The second solution of the determining equations ( 13) is given by the infinitesimals which generate the Möbius transformation In this case, it is not immediately obvious that the transformed concentration A * solves the same differential equation.However, if also the time derivative of the transformation is plugged into Eq.( 14) it is indeed fulfilled.Note that by Theorems 1 and 2 this is equivalent to calculating Ȧ * = exp[ X ] Ȧ, where X is the extended infinitesimal generator (7).On the other hand, the observation A obs is not altered because the transformation of A and the second factor in Eq. ( 15) have the same structure and if the transformations (16) are plugged in the transformation parameter cancels out.In fact, the second factor in Eq. ( 15) itself can be viewed as an element of the transformation group with transformation parameter = −s 2 .This example pinpoints one source of second order symmetries.The quadratic rate expression −2kA 2 of Eq. ( 14), generated by the bimolecular reaction, is preserved by the Möbius transformation because of its quadratic generator.On the other hand, the observation function allows the interpretation as a transformation, corresponding to a measurement technique that is affected by saturation effects.
In general, observation functions can also compensate transformations which do not have the same structure, as demonstrated in the following example taken from [9]: Here, the input function u(t) is assumed to be known and x 2 and x 3 are measurable on a relative scale.This results in the following linear ODE model: Because u(t) is known, it cannot be transformed and η u = 0 must hold.The system does not admit symmetries which are generated by univariate polynomial infinitesimals.However, the application of one of the more general multivariate Ansatz functions reveals an admitted transformation which is determined by the infinitesimal generator The necessary computations take about 15 s on a standard PC, which is due to the high number of constants required in the multivariate polynomial Ansatz.In order to construct the explicit transformation x * i = exp[ X]x i for this symmetry, it has to be slightly modified and integrated, which can be found in Ref. [9].However, even without the explicit transformation it can be seen from the infinitesimals η x 2 = k 2 x 2 and η x 3 = −k 1 x 3 that the transformations of x 2 and x 3 are not of the scaling type but also depend on the parameters k 1 and k 2 .

B. Symmetries of a complex signaling pathway
In a next step, the presented method is applied to an ODE model of the NFκB signaling pathway.The parameters occurring in the ODE model, i.e., initial conditions, rate constants, and observation parameters, are a priori unknown.Therefore, in a modeling application, they would need to be determined from experimental data.However, every parameter involved in a symmetry transformation of the ODE cannot be uniquely determined from this data because the observation is left invariant.Note that the scope of this paper is not to make conclusions about the function of the NFκB pathway.Here, it is solely used as a complex, interesting example to demonstrate the relationship between nonidentifiability and the various underlying symmetry transformations.Figure 1 shows a graphical representation of the analyzed NFκB model.NFκB is a transcription factor that is prevented from exerting its activity by its inhibitor IκBα, with which it is bound in a complex in the cell cytoplasm.Upon stimulation, the complex is phosphorylated by pIKK and dissociates into phosphorylated NFκB and IκBα.Whereas IκBα is degraded, NFκB shuttles into the nucleus, where it induces gene expression.Besides NFκB target genes associated with apoptosis, NFκB induces also gene expression of its own inhibitor, IκBα, leading to the production of IκBα mRNA and translation into the IκBα protein.These IκBα proteins shuttle into the nucleus, forming complexes with NFκB which are transported back into the cytoplasm.Depending on the duration of the activating pIKK signal, this cycle begins again, resulting in characteristic NFκB oscillations.In the model, each of the reactions is associated with a mathematical expression for the reaction flow, indicated as arrow labels in Fig. 1.See Ref. [15] for the comprehensive list of differential equations.
Protein levels of NFκB are measured independently in the cytoplasm and the nucleus but the measurement technique does not distinguish between single proteins and protein complexes.Phosphorylation of NFκB and IκBα is measured in the cytoplasm only.Again, the measurement technique does not distinguish complexes from single molecules.This results in the following observation functions: NFκB obs cyt : y 1 = s 1 (x 1 + x 2 + x 3 ) + I 0 cyt , NFκB obs nuc : y 2 = s 2 (x 10 + x 5 + x 6 ) + I 0 nuc , pNFκB obs cyt : y 3 = s 3 (x 2 + x 3 ), pIκBα obs cyt : Here, the species x i are numbered according to Fig. 1.The observation parameters s 1 ,s 2 ,s 3 ,s 4 ,I 0 cyt , and I 0 nuc in Eq. ( 19) are unknown.The analysis with Ansatz (12) and p max = 2 yields five admitted transformations while the computation time is on the order of 1 s.The first four are regular scaling symmetries.
(i) Because there are only relative observations in Eq. ( 19), it is possible to scale all the internal variables as long as this is compensated by an inverse transformation of the observation parameters: x * i = e x i , i = 1, . . .,10, s * j = e − s j , j = 1, . . .,4.Here, x i stands for every concentration in the model except pIKK.In addition, because complex formation in the nucleus described by the reaction flow v = k 10 x 9 x 6 involves a product of species, one scaling factor has to be compensated by k * 10 = e − k 10 .(ii) The second scaling symmetry is introduced by the factor ρ vol .Because the only measurement in the nucleus has a separate scaling factor, the species in the nucleus can be scaled independently from those in the cytoplasm, i.e., the transformation (iii) Another symmetry is It is admitted because neither is mIκBα observed nor is its scale coupled to any other concentrations.
(iv) Finally, a fourth admitted scaling transformation is k 1 because these variables only occur as products with each other.
In addition to the scalings, the model admits the higherorder transformation This solution pinpoints another source for second order transformations.The Möbius transformation generated by φ(x) = x 2 and the Michaelis-Menten-like rate expression coincide structurally and, thus, are able to compensate each other.Because the transformation is admitted, a change of the input by transformation (20) results in the same model response and does not change the observation.To illustrate this, Fig. 2(a) shows the time course of pIKK employed as input function in the ODE model.The activation starts with a strong peak in the first 30 min, which is followed by a slower decay starting at about 25% of the peak activation and requiring approximately 6 h for complete deactivation.In Fig. 2(b), transformation (20) is applied to the input for different values of .This illustrates that the observation remains invariant for inputs which are qualitatively very different.For negative values of , the plateau is almost removed and the activation is dominated by a sharp peak.On the other hand, for positive values of , the input becomes more sustained.

C. Reliable model predictions despite structural nonidentifiability
As described in Ref. [14], in case of a model with nonidentifiable variables, Lie group analysis can be utilized to find additional measurements which would make these variables identifiable.However, even if the nonidentifiabilities are not removed by additional experiments, the model can still be employed to make accurate predictions because for every structural nonidentifiability the presented method also yields the underlying transformation.To this end, it is tested if a candidate prediction y pred (t) = h(x(t)) admits all transformations which are admitted by the model, i.e., it must hold that X • [y pred − h(x)] = 0 for all infinitesimal generators X.If this if fulfilled, the predicted quantity is not affected by the admitted transformations, similar to the case of observation functions.Furthermore, if we allow transformations of y pred by assuming η y pred = 0, it can be detected if a prediction is possible only modulo some specific transformation.
As an example, consider again the NFκB model of Fig. 1 together with the quantity N pred rel : z 1 = ρ vol (x 1 + x 2 + x 3 ) x 5 + x 6 + x 10 , which gives the ratio between the absolute numbers of NFκB molecules in the cytoplasm and nucleus.If this is analyzed with the presented method, it is revealed that all transformations are admitted by the expression with η N pred rel = 0. Therefore, the quantity can be predicted by the model although all the occurring variables are structurally nonidentifiable.Another example is that of the quantities pIκBα pred : z 2 = x 2 + x 4 , mIκBα pred : z 3 = x 7 , (21) which describe the total concentration of phosphorylated IκBα molecules and of IκBα mRNA.These predictions are not possible on an absolute scale.However, by introducing the additional infinitesimals η pIκBα pred = z 2 , η mIκBα pred = z 3 , the symmetries are preserved.These infinitesimals generate scaling transformations, which means that it is possible to predict the relative time courses of the quantities in Eq. ( 21) despite the structural nonidentifiabilities in the model.

IV. DISCUSSION
Structural nonidentifiabilities in ODE models are generated by admitted transformations, meaning that parameters can be simultaneously transformed in a way that the observed quantities are not altered.Therefore, in order to detect structural nonidentifiability, the search for admitted transformations in the model equations is a promising approach.For this purpose, Lie group theory is a suitable framework as it is capable of handling a large variety of classes of transformations in a uniform way.This is achieved by reformulating every transformation in terms of its infinitesimal generator.Hereby, the necessary calculations can be performed in a linear space, the Lie algebra, which enables the formulation of a compact criterion to guarantee invariance of the observed quantities.The set of coupled partial differential equations arising from this criterion cannot be solved in general.However, in many practical applications, the ODE model consists only of rational equations.In particular, this holds for chemical reaction networks based on mass action, Michaelis-Menten, or Hill kinetics.Furthermore, all of the most frequent transformations are generated by polynomial infinitesimal generators.Based on these assumptions, the criterion for admitted transformations can be reformulated as a set of linear equations.The algorithm we developed to obtain the admitted infinitesimal generators exploits these assumptions to automatically and efficiently search for structural nonidentifiabilities in rational ODE models.To this end it does not rely on random specializations as previous methods but directly converts the conditions to a linear system which can then be solved.As the presented method is independent of any measured data, it is not capable of finding practical nonidentifiabilities which arise due to uncertainties in a given set of measurements.In order to detect also this type of nonidentifiability, a data-based method like the one presented in Ref. [1] has to be utilized.
In the present work, we employed the proposed algorithm in two ways as being illustrated on an ODE model of the NFκB signal transduction pathway with ten state variables, four observables, one input, and 20 parameters.First, given the observation functions, the algorithm is applied to determine admitted symmetry transformations.The parameters involved in these transformations are structurally nonidentifiable.Most of the found symmetries are of scaling type.Since the determination of protein levels occurs mostly on a relative scale, these observables are compatible with the intrinsic symmetries, yielding nonidentifiable parameters.However, we demonstrated that also Möbius transformations play an important role.This is due to its possible interpretation as saturating observation function or Michaelis-Menten-like rate expression.While saturation of measurement signals frequently occurs due to technical limitations, Michaelis-Menten rates are employed for modeling enzymatic reactions.Both scenarios are frequently found in various applications, underlining the practical relevance of the second order symmetries.
We also tested the performance of the differential algebra for identifiability of systems [7] and exact arithmetic rank [10] methods on our examples.Although they consistently found the same parameters as being nonidentifiable, they give no hint to the underlying transformation.However, this information is crucial for the second application of the proposed method, where we detect quantities which can be predicted by the model even if structural nonidentifiabilities are present.Such a prediction is possible if the corresponding expression is invariant under all admitted transformations or if it is varied with some specific, known transformation.
The limitation of the presented method in terms of feasibility depends on two factors, which are the size of the analyzed model and the size of the infinitesimal generator.While the number of necessary constants quickly rises for multivariate generators, it only grows linearly for the most common transformations like translations, scalings, and Möbius transformations.Therefore, in the latter case, the algorithm can be applied to large models with as much as hundreds of variables.

ACKNOWLEDGMENT
This work was supported by the Mechanism Based Integrated Systems for the Prediction of Drug Induced Liver Injury project, Innovative Medicines Initiative Joint Undertaking under Grant Agreement No. 115336.

APPENDIX A: SYMMETRIES AND STEADY STATE CONDITIONS
Theorem A1 (symmetries of steady state conditions).Let η i (x), i = 1, . . .,n be a set of infinitesimals and let the transformation . .,n be admitted by the differential equations ẋ = f (x(t)).Then the transformation is admitted by the steady state condition f (x 0 ) = 0.
Proof.Because the transformation is admitted by the differential equations, Eq. ( 8) of the main text holds, which means that Now let x 0 ∈ R n be a fixed point with f (x 0 ) = 0 and hence ẋ0 = 0.If this is plugged into Eq.(A1), we see that

APPENDIX B: MULTIVARIATE POLYNOMIAL GENERATORS
In the most general polynomial Ansatz, each infinitesimal η i (x) is given by a multivariate polynomial in all variables.However, as mentioned in Sec.II, the infinitesimals of the concentrations x i , i = 1, . . .,m cannot depend on the inputs u(t).In addition, parameters are time-independent states.In order to conserve this independence, the infinitesimals of the transformations can only depend on other parameters.Consequently, the most general polynomial Ansatz is With these infinitesimals, the first sum in the symmetry condition (11a) does not reduce to a single term.Therefore, the equation has to be multiplied by Q k n l=1 Q l to obtain a polynomial condition: Analogously to the case of univariate polynomial generators, this can be translated into a linear system for the constants r by rearranging the terms of the polynomial.The number of constants r i,p in the Ansatz is given by where ( • • ) denotes binomial coefficients.This shows that for this general Ansatz the size of the linear system C grows rapidly with the number of variables and eventually becomes infeasible for large systems.
A set of infinitesimals which is smaller than Eq.(B1) but still covers a large class of transformations is given by i.e., each infinitesimal η i (x) is given by a multivariate polynomial in the corresponding variable x i and all parameters x m+1 , . . .,x m+p .Similarly to the case of univariate polynomial generators, condition (9) can be reduced to Eq. ( 13) because either ∂η i ∂x j = 0 or f j (x) = 0 for j = k in the first sum.Therefore, the condition for admitted transformations can also be reduced to a linear system by multiplication of Q k 2 .In this case, the number of constants r i,d in the Ansatz is and therefore only grows with the number of parameters p and not with the total number of variables n = m + p + q.
R with 0 ∈ S and let φ : S × S → S be a law of composition.A mapping x * = (x, ) ∈ D with x ∈ D and ∈ S is called a Lie group of transformations if it satisfies the following.

FIG. 1 .
FIG. 1. Reaction scheme of the NFκB signaling pathway.NFκB is bound in a complex with IκBα in the cytoplasm.Complex phosphorylation by pIKK leads to complex dissociation, IκBα degradation, and NFκB shuttling to the nucleus.The associated reactions are indicated by arrows labeled by reaction flow expressions.

)FIG. 2 .
FIG. 2. Different time courses of pIKK.(a) pIKK activation as expected from the literature.(b) Transformation (20) applied to the expected activation for different values of (renormalized to a peak height of 1).
f (x) = 0 ∀k is fulfilled and the transformation is admitted by the steady state condition.