Bioinformatics Advance Access originally published online on July 28, 2007
Bioinformatics 2007 23(19):2612-2618; doi:10.1093/bioinformatics/btm382
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Data-based identifiability analysis of non-linear dynamical models
Physics Institute, University of Freiburg, Hermann Herder Strasse 3, 79104 Freiburg i.Br., Germany
*To whom correspondence should be addressed.
| ABSTRACT |
|---|
|
|
|---|
Motivation: Mathematical modelling of biological systems is becoming a standard approach to investigate complex dynamic, non-linear interaction mechanisms in cellular processes. However, models may comprise non-identifiable parameters which cannot be unambiguously determined. Non-identifiability manifests itself in functionally related parameters, which are difficult to detect.
Results: We present the method of mean optimal transformations, a non-parametric bootstrap-based algorithm for identifiability testing, capable of identifying linear and non-linear relations of arbitrarily many parameters, regardless of model size or complexity. This is performed with use of optimal transformations, estimated using the alternating conditional expectation algorithm (ACE). An initial guess or prior knowledge concerning the underlying relation of the parameters is not required. Independent, and hence identifiable parameters are determined as well. The quality of data at disposal is included in our approach, i.e. the non-linear model is fitted to data and estimated parameter values are investigated with respect to functional relations. We exemplify our approach on a realistic dynamical model and demonstrate that the variability of estimated parameter values decreases from 81 to 1% after detection and fixation of structural non-identifiabilities.
Availability: Our algorithm is written in Matlab and R. It is available from the authors on request. An implementation of ACE, written in Matlab as well as in C, is available online at www.stefanhengl.de
Contact: hengl{at}fdm.uni-freiburg.de
Supplementary information: Supplementary data are available at Bioinformatics online.
| 1 INTRODUCTION |
|---|
|
|
|---|
In the fast-growing field of Systems Biology, biological processes like signal transduction pathways and metabolic networks are often modelled mathematically based on systems of differential equations. They comprise parameters such as reaction rates, which have to be determined in accordance to measured data, e.g. by fitting to time course or dose response experiments. However, even for the simplest case of a reversible reaction of two species,
Two conceptually different sources for non-identifiability exist: first, the model structure itself may cause parameters to be functionally related. Second, since parameter values are estimated by fitting the model structure to experimental data even a structurally identifiable model may exhibit practical non-identifiabilities, because of an insufficient amount or quality of measurements. The noisier measurements are, and the lower the sampling frequency, the less information is contained in the measurement. Moreover, the dynamical response of the model depends on the input applied in the experiment, e.g. variable and constant stimuli may cause completely different dynamics. Therefore, the type of input may be decisive for parameters estimation (Faller et al., 2003).
Identifiability, however, is a necessary prerequisite for mathematical analysis of a model. Thus, the following question arises: how can non-identifiabilities be detected? Basically two approaches are used to handle non-identifiability: first, the model structure itself is investigated with respect to non-identifiabilities. If non-identifiabilities exist, they must be removed analytically by introduction of new parameters, representing, e.g. an identifiable combination of two non-identifiable parameters. This approach is referred to as a priori identifiability analysis, since the model is examined before simulating or fitting procedures. Second, a non-identifiable model structure manifests itself in functionally related parameters. Thus, non-identifiabilities may be detected by fitting a model repeatedly to data and investigating parameter estimates. Ideally, both methods are applied successively.
Numerous methods have been presented to deal with a priori identifiability analysis of linear and non-linear models. The Laplace transform or transfer function approach which may only be applied to linear models is thoroughly discussed in, Jacquez and Greif, 1985 and Godfrey and DiStefano 1987. However, when modelling biological systems, non-linear differential equations are ubiquitous, e.g. in Michaelis Menten kinetics and cooperative phenomena. The Similarity Transformation Approach (Vajda et al., 1989), the Power Series Expansion, (Pohjanpalo, 1978), the Volterra and Generating Power Series Approaches (Lecourtier et al., 1987) as well as identifiability tests derived from differential algebra (Ljung and Glad, 1994; Saccomani et al., 2003) are also applicable to non-linear models. Unfortunately, these methods become mathematically intractable with increasing model complexity. If we want to investigate which kind of non-identifiabilities of a model occur under realistic experimental conditions, a data based-method is to be applied. However, available data-based approaches for a posteriori identifiability analysis of non-linear models, like multivariate regression, require prior knowledge concerning explicit linear or non-linear functional relations (Quinn and Keough, 2002; Seber and Wild, 1988). In our approach, this is not necessary.
With our method, we provide a solution for following yet unsolved problems:
- structural identifiability analysis of large and complex non-linear dynamical models within reasonable time.
- automatic detection of non-linear dependencies between arbitrarily many parameters.
- practical identifiability analysis under consideration of the quality of data at disposal.
The article is organized as follows. First, we provide a definition of identifiability and introduce briefly the alternating conditional expectation (ACE) algorithm which estimates optimal transformations. Then, the test-function is introduced and its are discussed. Subsequently, we propose the algorithm which detects non-identifiabilities with help of the before defined test-function and apply it to a biological example of a non-linear dynamical model.
| 2 BACKGROUND |
|---|
|
|
|---|
2.1 Identifiability
Identifiability of a dynamical model depends on the model equations themselves, as well as on input and observation functions, initial conditions, constraints (Audoly et al., 2001; Godfrey and DiStefano, 1987) and the often unknown true parameter values (Vajda et al., 1989). A model together with all constraints is called a constrained structure.
We follow a definition given by Godfrey and DiStefano (Godfrey and DiStefano, 1987). Let
denote the state variables,
the externally given input signals,
the system parameters and
the observation function. The initial values
depend in general on the parameters
. Finally, let
denote all additional constraints mathematically formulated as explicit or implicit algebraic equations. A constrained structure is then given by
|
| (1) |
|
| (2) |
|
| (3) |
|
| (4) |
|
| (5) |
Biological data canonically comprises observational noise, generalizing Equation (2) to
|
|
Observational noise may render even structural identifiable parameters practically non-identifiable. Parameters which can be determined with a small enough SD are termed practical identifiable.
2.2 Alternating conditional expectation (ACE)
Initially, the ACE-algorithm has been proposed by Breiman and Friedman (Breiman and Friedman, 1985) for the purpose of regression analysis and has since been applied in various fields (Buja, 1990; Timer et al., 2000; Voss and Hurths, 1997; Wang and Murphy, 2005). In the bivariate case, ACE non-parametrically estimates optimal transformations
and
which maximize the linear correlation R between
and
,
|
|
Breiman and Friedman also provide weak conditions for convergence of the iterative algorithm. Since ACE itself is not the focus of this article, we just provide the basic knowledge needed to understand its further application.
The bivariate case can easily be extended to arbitrarily many parameters. Let
denote a (n x m) matrix of m parameters with n estimates for each parameter. Suppose the m parameters have an unknown linear or non-linear functional relation and let
and
j denote the true transformations that linearize the relation between the parameters,
|
|
is Gaussian noise. Then ACE estimates optimal transformations |
| (6) |
Note that ACE intrinsically distinguishes between left-and-right-hand side terms. It regards the left-hand side parameter to be the response and all others to be predictors.
The principle of the ACE algorithm for the multivariate case is as follows: It starts with an initial estimate for the optimal transformation of the response,
and of the predictors,
i = 0, i
{1.,n}/k. The transformation of the response,
, and the
i are estimated iteratively; new estimates of the
i serve as input for the estimation of
, and vice versa. For each
i it is calculated, how much variance of the response,
( pk), cannot yet be explained by the m – 2 other predictors,
. This unexplained variance is the next estimate for the predictor
i. In other words, the best estimate for a transformation
i minimizes the squared residuals of Equation (6). It can be shown, that this estimate is given by
|
|
The
i are updated successively until the fraction of unexplained variance of the predictor fails to decrease. The last updated estimates serve as input to estimate
( pk) like follows
|
|
(pk) does not cause further reduction of the fraction of unexplained variance. In practice, the conditional expectation is replaced by smoothing techniques.
2.2.1 Remarks
We use a modified version of ACE as implemented by H. Voss (Voss and Kurths, 1997), where data are ranked before optimal transformations are estimated. Nevertheless, our application does not depend on a special implementation of ACE. We also tested our algorithm with ACE as included in the acepack package of R.
ACE works non-parametrically, thus we do not have to make assumptions about underlying functional relations. However, it does not provide explicit functional expressions as output, but values of the estimated optimal transformations
and
at the parameter values which served as input.
| 3 APPROACH |
|---|
|
|
|---|
Our goal is to reveal, in a data-based way, non-identifiabilities of a non-linear dynamical model. Non-identifiability manifests itself in functionally related parameters, and we apply ACE to non-parametrically estimate these relations. However, to do this objectively, a test function quantifying this information is required. In order to motivate this test-function, we study the typical behaviour of the ACE algorithm by means of a simple example. Our findings can directly be applied to identifiability analysis, as will be outlined afterwards.
Let p2,p3 and p4 be three parameters uniformly distributed on the interval
. Suppose a fourth parameter p1 to be functionally related with p2,p3 by
. Real data is always corrupted with observational noise, i.e. random errors occurring during measurement, and we take this into account by adding Gaussian noise to the response p1
|
| (7) |
We independently draw 200 tuples (p2,p3,p4) and calculate p1 for each tuple, thus gaining a 200 x 4 matrix K, mimicking 200 parameter estimates. Only the first three columns of each row are functionally related, the fourth being independent. Then, we take K as input for ACE to estimate optimal transformations (see Fig. 1). Note that even p4 seems to have quite a smooth estimated optimal transformations significantly different from white noise. This may be mistaken for a real functional dependency. If, however, we draw a new 200 x 4 matrix the same way as outlined above, the transformations of the first three parameters will remain quite stable, while the transformation of the fourth parameter will in general look different, but still smooth. We can understand this volatility if we recall the iterative form of the ACE algorithm that estimates the optimal transformations. Since p2 and p3 suffice to explain all variance of p1 but Gaussian noise, the residuals will be Gaussian noise, smoothed by the filter employed by ACE. Noise will be distributed differently in each new sample we draw, therefore, we yield different estimates for the transformation of p4 every time.
|
This is exactly, what renders possible to distinguish between parameters with functional relations and those without; we calculate the average estimated transformation by repeatedly drawing new matrices K, estimating transformations each time. We expect optimal transformations of functionally related parameters to be invariant under averaging. In contrast, parameters without functional relations yield different estimates for each new drawn matrix K.
The connection to the identifiability analysis of a constraint structure follows immediately: a non-identifiable constraint structure causes parameters to be functionally related and we may employ the above findings to detect these relations. The process of drawing new matrices K is replaced by estimating parameters through repeated fitting with different initial guesses of the parameters of the constraint structure to experimental or simulated data. This means, to yield a single 200 x m matrix K, we fit the model 200 times to data; we term this a single fitting sequence. Thus, the problem of identifiability analysis is mapped onto the problem of detecting groups of functionally related parameters. In the following, the qualitative ideas of our approach are quantified within a test-function, which is used to decide for functional dependencies.
| 4 METHODS |
|---|
|
|
|---|
4.1 Construction of the test-function
A quantitative test-function based on the ACE algorithm is required to detect sets of functionally related parameters. In order to improve its robustness, all estimated optimal transformations for a given fitting sequence are ranked. Let
|
|
|
|

, i.e. zero variance. This motivates following definition of a test function
|
| (8) |
{2,...,m} used in the calculation of Hk, |
|
|
|
4.2 Properties
and Hk render it possible to distinguish between three different cases, see Figure 2.
: the response parameter pi1 on the left-hand side has no functional relation with any other parameter
: a given set of parameters contains not enough information, i.e. we need to add additional parameters to establish a strong functional relation.
: a given set of parameters contains enough information to establish a strong functional relation.
|
The three above stated cases correspond to three regions of variance marked off by two threshold values which are determined analytically (see supplementary Material). Let N denote the number of estimates, and NoB the number of bootstrap samples. The expectation value of the test function of a parameter which is not functionally related with a given set of other parameters is
|
| (9) |
The expectation value of the test-function of a parameter which has a strong monotone functional relation with a given set of parameters is
|
| (10) |
- Both results Equations (9) and (10) have to be equal in the limes of large N and NoB = 1 which is fulfilled (c.f. Fig. 2 inset).
- Equation (9) assumes independence of the drawn samples which is asymptotically fulfilled for infinitely many fits within one fitting sequence. Therefore, the obtained result is too small for finite N. If we replace resampling by new fitting sequences, Equation (9) holds for all N.
- Equations (9) and (10) assume ACE to estimate optimal transformations perfectly, which again is only asymptotically true for large N and no noise. Especially, noise results in a slight shift of the test-function for functionally related parameters to smaller variances (c.f. Fig. 2 inset).
4.3 The algorithm
The test-function Hk proposed in Equation (8) interprets the stability of estimated optimal transformations as a measure for functional relations of parameters. Our algorithm exploits the properties of Hk to identify those parameters pk of a given set
which have a functional relation and determines the relations non-parametrically.
ACE distinguishes between right-and left-hand side terms, which we want to make use of. We propose a step up algorithm, i.e. we start with parameter pi as left-hand side term and take a second parameter
, as right-hand side. Then, we calculate Hj(pi,pj) for each pj and choose the parameter pk with
which means we take the parameter which is likely to have the strongest functional relation with pi. If
, parameter pi is supposed to be independent of all others. Otherwise we calculate
. If
a strong functional relation is found. Else another parameter
, is added and
, calculated again. If T2 is never exceeded, even if we added all parameters, it is supposed to have no strong functional relation.
To increase power, we rerun the loop once again after having found a strong functional relation. If the
value increases with an additional parameter, we keep it, else we take the relation which has already been found before.
This is conducted successively for all parameters pi taken as left-hand side; Figure 8 in the Supplementary Material shows a flowchart of the proposed algorithm.
As a crosscheck, we determine the multiple r2, i.e. the fractional amount of variance of the response y explained by the predictor variables
.
4.4 Example
We work through an example to illustrate what kind of input is needed and which output is provided. Assume a non-identifiable constraint structure with seven parameters functionally related like follows:
|
| (11) |
are uniformley distributed on the interval [0 5]. The input used in our algorithm is |
|
4.5 Sensitivity and specificity
We determined sensitivity and specificity in dependence of the threshold values for Equation (11). In order to test robustness of our algorithm, noise was added to all left-hand side terms. Figure 3 confirms that defined threshold values yield high sensitivity as well as high specificity.
|
The inset of Figure 3 stresses the generality of the algorithm: sensitivity is largely independent of the actual functional relation. It only depends on the contribution strength of a predictor. The less a predictor pj on the right-hand side contributes to the response on the left-hand side, the noisier the estimated transformation
j( pj) gets, finally being indistinguishable from a estimated transformation of an independent parameter. As discussed in the supplementary Material (section: Identifiability of Identifiability), especially for two parameters, there is a strong dependence of optimal transformations on the level of noise. However, this is, to a large extent, compensated by the bootstrap approach and the rank transformation of the optimal transformations, see Equation (8).
| 5 RESULTS |
|---|
|
|
|---|
We apply our method to a non-linear dynamical model motivated by a modelling-project dealing with endocytosis. Our goal is to identify groups of functionally related parameters. Once such sets of non-identifiable, interdependent parameters are detected, either new experiments can be suggested to render the parameters identifiable, or parameters have to be fixed in order to improve identifiability. For the latter case, we suggest following guidelines for handling the two most frequent types of non-identifiable parameters: (1) structurally non-identifiable parameters may, per definition, be fixed at an arbitrary value in parameter space. The model's; dynamical properties are not changed or restricted by the fixation. (2) For practically non-identifiable parameters, the model's; dynamical properties depend slightly on the chosen parameter value. The choice which parameter to fix, depends on the experimental possibilities and available reference values from the literature. If there is no additional information, we suggest to fix the parameter within physiological constraints to that value which belongs to the best fit. (3) Iterate identifiability analysis and fixation of parameters until all free parameter are rendered identifiable. This is a necessary prerequisite for subsequent analysis techniques like sensitivity analysis. Further guidelines are provided in the Supplementary Material.
We proceed as follows : first, we fit the model to simulated data. In order to test identifiability under realistic experimental conditions, we add observational noise. Second, we apply our algorithm for identifiability analysis and interpret the result. Consider the following non-linear model derived from the biological system in Figure 4:
|
|
10-5; k3 = 0.10; k4 = 0.25; k5 = 0.15; k6 = 0.075; k7 = 0.01; Bmax = 1000; kD = 100; EPO (t = 0) = 3000; EPOR (t = 0) = 1000. All simulations were conducted with our developed Systems Biology Multi-Experiment Fitting Matlab Toolbox PottersWheel (www.potterswheel.de). The model was fitted 500 times to simulated data, each fit started at the true parameter values, disturbed like follows: |
|
|
|
|
Note that in general it is possible for three functionally related parameters to lie on a 2D manifold, i.e. a curved surface in space. This would require to fix two parameters in order to completely resolve non-identifiability. For an invertible optimal transformation
, however, it is also possible to reduce the number of non-identifiable parameters, not the non-identifiablity itself, by expressing one of the parameters through the other ones | 6 CONCLUSION |
|---|
|
|
|---|
We presented a non-parametric data-based algorithm for identifiability testing. Its major characteristics, generality and robustness, render it a valuable tool for identifiability analysis of non-linear dynamical models. The use of the bootstrap method reduces computational time tremendously.
Non-identifiabilities arise due to the structure of a model and the observation function. The proposed method is capable of identifying structural subunits causing non-identifiabilities by detecting groups of functionally related parameters. It reveals which parameters are uniquely determinable and which are not. Thus, it also provides a more realistic picture of what can be inferred from a model. In an example, the variance of estimated parameter values could be reduced from around 81 to 1%, after fixing one parameter as suggested by our approach.
The ability to apply the proposed approach to simulated data renders it directly applicable for Experimental Design, i.e. the quest for optimal experimental conditions, e.g. observables to maximize the information while minimizing experimental effort. For this purpose, we may test the identifiability of a model at all possible combinations of realistic input–output scenarios. Among the identifiable structures we then choose the input–output combination which is easiest to realize experimentally.
| ACKNOWLEDGEMENTS |
|---|
|
|
|---|
The authors thank Martin Peifer for his support in the analytical calculations, as well as Verena Becker, Marcel Schilling, and Ursula Klingmüller from the German Cancer Research Center in Heidelberg for providing the model discussed. S.H., C.K. and T.M. gratefully acknowledge financial support by the grant BMBF 0313074D & FP6 EU-grant COSBICS LSHG-CT-2004-0512060.
Conflict of Interest: none declared.
| FOOTNOTES |
|---|
Associate Editor: Limsoon Wong
Received on April 4, 2007; revised on July 10, 2007; accepted on July 10, 2007
| REFERENCES |
|---|
|
|
|---|
Audoly S, et al. Global identifiability of nonlinear models of biological systems. IEEE Trans. Biomed. Eng. (2001) 48:55–65.[CrossRef][Web of Science][Medline]
Breiman L, Friedman J. Estimating optimal transformations for multiple regression and correlation. J. Am. Stat. Assoc. (1985) 80:580–598.[CrossRef][Web of Science]
Buja A. Remarks on functional canonical variates, alternating least squares methods and ace. Ann. Stat. (1990) 18:1032–1069.
Faller D, et al. Simulation methods for optimal experimental design in systems biology. Simulation (2003) 79:717–725.[Abstract]
Godfrey KR, DiStefano JJ. Identifiabiliy of Model Parameters in Identifiability of Parametric Models, Chapter 1. (1987) Oxford: Pergamon Press. 1–20.
Jacquez J, Greif P. Numerical parameter identifiability and estimability: Integrating identifiability, estimability, and optimal sampling design. Math. Biosci (1985) 77:201–227.[CrossRef][Web of Science]
Lecourtier Y, et al. Volterra and Generating Power Series Approach to Identifiability Testing in Identifiability of Parametric Models, Chapter 5. (1987) Oxford: Pergamon Press. 50–66.
Ljung L, Glad T. On global identifiability for arbitrary model parametrizations. Automatica (1994) 30:265–276.[CrossRef][Web of Science]
Pohjanpalo H. System identifiability bases on the power series expansion of the solution. Math. Biosci (1978) 41:21–33.[CrossRef][Web of Science]
Quinn G, Keough M. Experimental Design and Data Analysis for Biologists (2002) Cambridge, UK: Cambridge University Press.
Saccomani M, et al. Parameter identifiability of nonlinear systems: the role of initial conditions. Automatica (2003) 39:619–632.[CrossRef][Web of Science]
Seber G, Wild C. Nonlinear Regression (1988) New York: Wiley.
Timmer J, et al. Parametric, nonparametric and parametric modelling of a chaotic time series. Phy. Lett. A (2000) 274:123–134.[CrossRef]
Vajda S, et al. Similarity transformation approach to identifiability analysis of nonlinear compartmental models. Math. Biosci. (1989) 93:217–248.[CrossRef][Web of Science][Medline]
Voss H, Kurths J. Reconstruction of non-linear time delay models from data by the use of optimal transformation. Phy. Lett. A (1997) 234:166–344.
Wang D, Murphy M. Identifying nonlinear relationships in regression using the ACE algorithm. J. Appl. Stat. (2005) 32:243–258.[CrossRef][Web of Science]
This article has been cited by other articles:
![]() |
A. Raue, C. Kreutz, T. Maiwald, J. Bachmann, M. Schilling, U. Klingmuller, and J. Timmer Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood Bioinformatics, August 1, 2009; 25(15): 1923 - 1929. [Abstract] [Full Text] [PDF] |
||||
![]() |
O. Kotte and M. Heinemann A divide-and-conquer approach to analyze underdetermined biochemical models Bioinformatics, February 15, 2009; 25(4): 519 - 525. [Abstract] [Full Text] [PDF] |
||||
![]() |
B. Finkenstadt, E. A. Heron, M. Komorowski, K. Edwards, S. Tang, C. V. Harper, J. R. E. Davis, M. R. H. White, A. J. Millar, and D. A. Rand Reconstruction of transcriptional dynamics from gene reporter data using differential equations Bioinformatics, December 15, 2008; 24(24): 2901 - 2907. [Abstract] [Full Text] [PDF] |
||||
![]() |
T. Maiwald and J. Timmer Dynamical modeling and multi-experiment fitting with PottersWheel Bioinformatics, September 15, 2008; 24(18): 2037 - 2043. [Abstract] [Full Text] [PDF] |
||||
![]() |
H. Schmidt, M. F. Madsen, S. Dano, and G. Cedersund Complexity reduction of biochemical rate expressions Bioinformatics, March 15, 2008; 24(6): 848 - 854. [Abstract] [Full Text] [PDF] |
||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||





) and specificity (^) is calculated. Same is done for T2 (solid lines) with T1 = 0.01. Lines at 100% are separated for clarity. We see that analytically determined threshold values, T1 = 0.01 and T2 = 0.07, have optimal specificity and sensitivity. (Inset) power of proposed algorithm. The percentage of recovered functional relations is plotted versus the mean contribution. To yield comparable results, we tested a set of standard functions fi with
f(p3) of the second right-hand side term to the left-hand side is defined like follows: M



2 values. In order to render the constrained structure identifiable, one of the three parameters has to be fixed.
