Sara Berthoumieux

Matteo Brilli

Daniel Kahn

Hidde De Jong

Eugenio Cinquemani
email: eugenio.cinquemani@inria.fr

On the Identifiability of Metabolic Network Models

Keywords: Systems biology, Metabolic network modeling, Parameter estimation, Structural and practical identifiability, Principal Component Analysis, Singular Value Decomposition, Escherichia coli carbon metabolism

A major problem for the identification of metabolic network models is parameter identifiability. The identifiability of a parameter consists in the possibility to unambiguously infer the parameter from the data. Identifiability problems may be due to the structure of the model, in particular implicit dependencies between the parameters, or to limitations in the quantity and quality of the available data. We address the detection and resolution of identifiability problems for a class of pseudo-linear models of metabolism, so-called linlog models. Linlog models have the advantage that parameter estimation reduces to linear or orthogonal regression, which facilitates the analysis of identifiability. We develop precise definitions of structural and practical identifiability, and clarify the fundamental relations between these concepts. In addition, we use Singular Value Decomposition (SVD) to detect identifiability problems and reduce the model to an identifiable approximation by a Principal Component Analysis (PCA) approach. The criterion is adapted to real data, which are frequently scarce, incomplete, and noisy. The test of the criterion on a model with simulated data shows that it is capable of correctly identifying the principal components of the data vector. The application to a state-of-the-art dataset on central carbon metabolism in Escherichia coli yields the surprising result that only 4 out of 31 reactions, and 37 out of 100 parameters, are identifiable. This underlines the practical importance of identifiability analysis and model reduction in the modeling of large-scale metabolic networks. Although our approach has been developed in the context of linlog models, it carries over to other pseudo-linear models and provides useful hints for the identifiability analysis of more general classes of nonlinear models of metabolism.

Introduction

Kinetic models of biochemical reaction systems usually contain a large number of parameters, many of which are difficult to measure in a direct way. This makes parameter estimation from experimental data a crucial problem for quantitative systems biology [START_REF] Ashyraliyev | Systems biology: Parameter estimation for biochemical models[END_REF][START_REF] Crampin | System identification challenges from systems biology[END_REF][START_REF] Jaqaman | Linking data to models: data regression[END_REF][START_REF] Chou | Recent developments in parameter estimation and structure identification of biochemical and genomic systems[END_REF]. For all but the simplest systems, parameter estimation of biological systems is a difficult problem. Models contain a large number of variables, whose dynamics evolve on different time-scales and are described by complex, nonlinear rate equations. Moreover, the available experimental data usually consist of noisy and incomplete measurements, obtained under different conditions and by means of heterogeneous experimental methods.

These difficulties have stimulated the use of approximate kinetic models, with rate equations following, for example, linear, piecewise-linear, linlog, or power-law kinetics [START_REF] De | Modeling and simulation of genetic regulatory systems: A literature review[END_REF][START_REF] Heijnen | Approximative kinetic formats used in metabolic network modeling[END_REF][START_REF] Chou | Recent developments in parameter estimation and structure identification of biochemical and genomic systems[END_REF]. In general, these models have fewer parameters to estimate and, in many cases, parameter estimation can be reduced to linear or orthogonal regression [START_REF] Nikerel | Model reduction and a priori kinetic parameter identifiability analysis using metabolome time series for metabolic reaction networks with linlog kinetics[END_REF][START_REF] Sands | Flux-based estimation of parameters in S-systems[END_REF]. Therefore, approximate kinetic models have been chosen and successfully used for the quantitative modeling of both metabolic networks and gene regulatory networks. In this paper, we focus on the case where the structure of the model is fixed by a priori knowledge on the network (i.e., the chemical species considered, the reactions among them and the possible regulatory interactions).

A major and often overlooked problem in parameter estimation is the identifiability of the model, that is, the problem of unambiguously reconstructing the unknown parameter values from the observed network behavior. A distinction is usually made between structural (or a priori) and practical (or a posteriori) identifiability [START_REF] Ljung | System Identification, Theory for the User[END_REF][START_REF] Walter | Identification of Parametric Models[END_REF]. Structural identifiability is an intrinsic property of the model family, guaranteeing that unique parameter reconstruction would be possible from perfect observations of the system response to an arbitrarily rich set of inputs. Practical identifiability refers to the ability of estimating unknown parameter values from the available experimental data within a prespecified degree of accuracy. In classical control theory, this concept is essentially related to the notion of persistence of excitation [START_REF] Ljung | System Identification, Theory for the User[END_REF]; unfortunately, limitations in the variety and quality of the data make this notion inapplicable for biological applications. In recent years, the topic of identifiability has gained considerable interest in the field of systems biology [START_REF] Chen | Classic and contemporary approaches to modeling biochemical reactions[END_REF]Chis et al., 2011b;[START_REF] Nikerel | Model reduction and a priori kinetic parameter identifiability analysis using metabolome time series for metabolic reaction networks with linlog kinetics[END_REF][START_REF] Raue | Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood[END_REF][START_REF] Raue | Addressing parameter identifiability by model-based experimentation[END_REF][START_REF] Gutenkunst | Universally sloppy parameter sensitivities in systems biology models[END_REF][START_REF] Nemcova | Structural identifiability of polynomial and rational systems[END_REF][START_REF] Srinath | Parameter identifiability of power-law biochemical system models[END_REF]Voit et al., 2006a) and several specialized software packages have been developed to support the modeler [START_REF] Bellu | DAISY: a new software tool to test global identifiability of biological and physiological systems[END_REF]Chis et al., 2011a;[START_REF] Maiwald | Dynamical modeling and multi-experiment fitting with Pot-tersWheel[END_REF]. Despite these efforts, however, no common agreement on definitions and links between structural and practical identifiability exist to date.

The aim of this paper is to develop methods for the analysis of parameter identifiability of kinetic models of metabolism, and for the reduction of nonidentifiable models to identifiable approximations. These methods should have a solid mathematical foundation, but at the same time be applicable to practical problems and be scalable to currently available datasets, such as those obtained by means of recent high-throughput methods in biology [START_REF] Ishii | Multiple high-throughput analyses monitor the response of E. coli to perturbations[END_REF]. While many of our definitions are of general applicability, identifiability results will be developed primarily for approximate kinetic models known as linlog models [START_REF] Visser | Dynamic simulation and metabolic re-design of a branched pathway using linlog kinetics[END_REF], whose pseudo-linear form enables us to apply tools from linear algebra and estimation theory in a straightforward manner. Similar results can be derived for many other approximate kinetic modeling formalisms in pseudo-linear form, such as the linear, loglin and generalized mass-action kinetic formats [START_REF] Delgado | Metabolic control analysis using transient metabolite concentrations[END_REF][START_REF] Hatzimanikatis | Effects of spatiotemporal variations on metabolic control: Approximate analysis using (log)linear kinetic models[END_REF][START_REF] Savageau | Biochemical Systems Analysis: A Study of Function and Design in Molecular Biology[END_REF]. Moreover, estimated parameters of approximate kinetic models provide useful hints for the identification of more detailed nonlinear models.

The present paper builds upon and elaborates earlier work, in which we developed an identification method for approximate kinetic models in the case of incomplete data [START_REF] Berthoumieux | Identification of linlog models of metabolic networks from incomplete high-throughput datasets[END_REF]. More precisely, the main contributions of the present paper are threefold. First, we precisely define the notions of structural and practical identifiability of linearized kinetic models, drawing upon the systems identification literature. This conceptual clarification allows us to develop the relations between structural and practical identifiability in a fundamental way. Second, we show how model reduction using Singular Value Decomposition (SVD) [START_REF] Jolliffe | Principal Component Analysis[END_REF] provides a suitable theoretical framework for addressing identifiability problems. We discuss several different criteria for model reduction, based on the singular values returned by the SVD analysis, and we show to which extent these criteria are appropriate for dealing with actual biological datasets, which are typically scarce, noisy and incomplete. Third, we apply the methods for identifiability analysis and model reduction to both simulated data and a published dataset concerning central metabolism in E. coli [START_REF] Ishii | Multiple high-throughput analyses monitor the response of E. coli to perturbations[END_REF]. These examples show that the mathematical tools developed in this paper are of practical utility for the estimation of parameters in metabolic network models, and beyond, from current high-throughput data sets.

For the readers'convenience, the notation adopted in the paper is summarized in Appendix A. To simplify the reading, all mathematical proofs are deferred to Appendix B.

Parameter estimation in linearized kinetic models

The dynamics of biochemical reaction networks are described by kinetic models having the form of systems of ordinary differential equations (ODEs) [START_REF] Heinrich | The Regulation of Cellular Systems[END_REF].

In this paper we focus on kinetic models of metabolism, where the rate functions describe enzyme-catalyzed reactions. This leads to models of the general form:

ẋ = N • v(x, u, e), (1)
with x(0) = x 0 ∈ R nx >0 , where x ∈ X ⊆ R nx >0 denotes the vector of internal metabolite concentrations, u ∈ U ⊆ R nu >0 the vector of external metabolite concentrations, e ∈ E ⊆ R m >0 the vector of enzyme concentrations, and v : R nx+nu+m >0 → V , with V ⊆ R m , the vector of reaction rate functions. N ∈ Z nx×m is a stoichiometry matrix.

Kinetic modeling formalisms differ in the choice of rate functions. Examples of classical functions in enzyme kinetics are the Michaelis-Menten, reversible Michaelis-Menten, and Monod-Wyman-Changeux rate laws [START_REF] Heinrich | The Regulation of Cellular Systems[END_REF]. Approximate formalisms simplify the mathematical form of the rate laws, in particular the nonlinear dependency of the reaction rates on the metabolite concentrations. They usually assume all reactions to follow the same simplified kinetic format, thus giving a uniform structure to the models. Moreover the parameter estimation problem for linearized kinetic models can be recast as multiple linear regression, as we will now show on linear-logarithmic (linlog) models.

The linlog approximation [START_REF] Heijnen | Approximative kinetic formats used in metabolic network modeling[END_REF][START_REF] Visser | Dynamic simulation and metabolic re-design of a branched pathway using linlog kinetics[END_REF] expresses the reaction rates as proportional to the enzyme concentrations and to a linear function of the logarithms of internal and external metabolite concentrations. where diag(e) is the square diagonal matrix with the elements of e on the diagonal, and the logarithm of a vector means the vector of logarithms of its elements. For conciseness, in the sequel we shall often drop the dependence of v on (x, u, e) from the notation.

v(x, u, e) = diag(e) • a + B x • ln x + B u • ln u (2) (a) v 1 = e 1 • (a 1 + B x 11 ln x 1 + B x 12 ln x 2 ) v 2 = e 2 • (a 2 + B x 21 ln x 1 + B x 22 ln x 2 ) v 3 = e 3 • (a 3 + B x 32 ln x 2 ) (b)
Example 1. Fig. 1(a) illustrates a prototype of a metabolic reaction network with negative feedback regulation. In terms of Eq. 1, we have

x = [x 1 x 2 ] T , e = [e 1 e 2 e 3 ] T , v = [v 1 v 2 v 3 ] T , and
N = 1 -1 0 0 1 -1 .
The linlog rate equations for this system are shown in Fig. 1(b). We will refer to this network as a running example illustrating the concepts introduced below.

The identification of metabolic networks in the linlog formalism amounts to estimating the generally unknown parameters a ∈ R m , B x ∈ R m×nx and B u ∈ R m×nu from experimental data. In most experiments, concentrations of enzymes and external metabolites are under partial control of the experimentalist, and the concentrations of internal metabolites and metabolic fluxes are measured after the system has relaxed to the steady-state

N • v(x, u, e) = 0.
(3)

In accordance with this, we shall assume that, from each of q ∈ N experiments, the data are noisy measurements (ṽ k , xk , ũk , ẽk ) of (v k , x k , u k , e k ), where the latter satisfy v k = v(x k , u k , e k ) and (3), with k = 1, . . . , q. Clearly the restriction to steady-state measurements limits the informativity of the data and may affect the identifiability of the models, as will be apparent in later sections.

For the purpose of parameter estimation, it is convenient to rewrite (2) in the form of a regression model:

v e T = [1 ln x T ln u T ] •   a T (B x ) T (B u ) T   .
(4)

Note that e is a vector of strictly positive elements, which enables the formulation of Eq. ( 4). By the linearity of (4), it holds that

v e = [1 ln x T ln u T ] •   a T (B x ) T (B u ) T   .
(5)

This allows (4) to be reformulated as a mean-removed model

v e - v e T = ln x -ln x ln u -ln u T • (B x ) T (B u ) T . (6)
We can now formulate our general estimation problem.

Problem 1. Given the data matrices

     ṽ1 ẽ1 -ṽ ẽ T . . . ṽq ẽq -ṽ ẽ T      W ,    ln x1 -ln x T ln ũ1 -ln ũ T . . . . . . ln xq -ln x T ln ũq -ln ũ T    Ỹ find parameters B B x B u T minimizing ||| W -Ỹ • B|||, where ||| • ||| is a convenient matrix norm on R q×m .
To this avail we consider the probabilistic measurement error model:

W = W + ε ε = ε 1 , . . . , ε m ε i ∼ N (0, Σ ε i ) (7) Ỹ = Y + η η = η 1 , . . . , η nx+nu η j ∼ N (0, Σ η j ) (8)
with Σ ε i = σ 2 i I > 0, Σ η j = ν 2 I ≥ 0, and ε i , η j mutually independent for all i = 1, . . . , m and j = 1, . . . , n x + n u .

Notice that the parameter vector a no longer appears in the regression problem, but that an estimate of it can be recovered from estimates of B = B x B u T by way of Eq. ( 5). With the assumption above that measurement noise is independent across different reactions, it makes sense to separate the problem into the independent estimation of the parameter vector B i of each reaction i, with i = 1, . . . , m.

Example 2. Let W and Y denote the noiseless versions of W and Ỹ , respectively. Consider the case where W = W + ε = Y • B + ε, i.e. the measurement error for the metabolite concentrations is negligible. Maximum likelihood estimation of B amounts to maximize the probability (density function) of W given Y as a function of B. After simple computations and thanks to the independence assumptions on ε, one finds that the maximum likelihood estimate of B is any solution of

min B 1 2 m i=1 ( Wi -Y B i ) T Σ -1 ε i ( Wi -Y B i ),
wihch can be solved separately for every column of B by solving, for i = 1, . . . , m,

min B i ( Wi -Y B i ) T Σ -1 ε i ( Wi -Y B i ) = || Wi -Y B i || 2 Σ -1 ε i . Thus, defining || • || i = || • || Σ -1 ε i and ||| • ||| : R q×m → R ≥0 : M →    M 1 . . . M m    T    Σ -1 ε 1 . . . Σ -1 εm       M 1 . . . M m    ,
we see that Problem 1 is equivalent to the maximum likelihood estimation of B, which is in turn equivalent to separate maximum likelihood estimation of the parameters B i of each reaction i.

We can now fully detail Problem 1 and express it as a series of estimation problems on individual reactions. In doing so we note that each reaction i depends only on a known subset of metabolites C(i) ⊆ {1, . . . , n x + n u }. Therefore, the entries of B i corresponding to metabolites that do not participate in reaction i can be set to zero, and the least-squares problem can be reduced accordingly. We will address two cases, formalized by two alternative problem statements. The first we consider is a standard regression problem [START_REF] Nikerel | Model reduction and a priori kinetic parameter identifiability analysis using metabolome time series for metabolic reaction networks with linlog kinetics[END_REF][START_REF] Sands | Flux-based estimation of parameters in S-systems[END_REF]. In analogy with Example 2, it amounts to assuming negligible noise for metabolite concentrations.

Problem 2. Given Y and W as in (7), solve min

B C(i),i || Wi -Y C(i) • B C(i),i || i , i = 1, . . . , m . (9)
The second case is more challenging and less commonly addressed in the literature. It corresponds to an errors-in-variables regression model [START_REF] Van Huffel | The Total Least Squares Problems: Computational Aspects and Analysis[END_REF] and accounts explicitly for noise on the relative fluxes as well as on metabolite concentrations.

Problem 3. Given Ỹ as in (8) and W as in (7), solve min

B C(i),i || Wi -ỸC(i) • B C(i),i || i , i = 1, . . . , m . (10)
From now on we will drop subscript i from || • || i , the meaning being clear from the argument of the norm.

Remark 1. A similar parameter estimation problem can be formulated for other pseudolinear modeling formalisms. Models linear in metabolite concentrations [START_REF] Delgado | Metabolic control analysis using transient metabolite concentrations[END_REF], loglin models [START_REF] Hatzimanikatis | Effects of spatiotemporal variations on metabolic control: Approximate analysis using (log)linear kinetic models[END_REF], and generalized mass-action models [START_REF] Savageau | Biochemical Systems Analysis: A Study of Function and Design in Molecular Biology[END_REF]) can be defined analogously to Eq. ( 4). This gives rise to, respectively,

v e T = [1 x T u T ] •   a T (B x ) T (B u ) T   , (11)
(v) T -ln e T = [1 ln x T ln u T ] •   a T (B x ) T (B u ) T   , ( 12
) ln v e T = [1 ln x T ln u T ] •   a T (B x ) T (B u ) T   . ( 13
)
Notice that the modifications concern the way in which reaction rates and concentrations enter into the linear equations. The translation of these equations into variants of Problems 2-3, by removing the mean, is straightforward. In each case, we obtain a linear regression problem.

Although below we illustrate the identifiability issues and reduction methods for the case of linlog models, it should be borne in mind that analogous results also apply to the other approximate kinetic modeling formalisms defined in Eq.s (11)-( 13). However, they are not applicable to formalisms for which parameter estimation cannot be turned into linear regression, such as reversible Michaelis-Menten kinetics and convenience kinetics [START_REF] Liebermeister | Bringing metabolic networks to life: Convenience rate law and thermodynamic constraints[END_REF].

Identifiability of linlog and related models

The problem of identifiability refers to the ability to unambiguously extract parameter values of a model structure from experimental data. Here we focus on linlog models and investigate the identifiability of this model class along the lines of [START_REF] Berthoumieux | Structural and practical identifiability of approximate metabolic network models[END_REF]. We shall first discuss the problem from the perspective of structural identifiability. For practical purposes, this is equivalent to answering the question whether each parameter can be uniquely reconstructed from an arbitrarily rich and errorless dataset. Structural identifiability forms the basis for studying practical identifiability, i.e. the ability to estimate parameter values from real datasets, which will be discussed further below.

The system, described by Eq.s (1)-(3), is parametrized by the parameter vector p = [a B x B u ] T ∈ P ⊆ R (1+nx+nu)×m . Let e, u, x and v take values in the sets

E ⊆ R m >0 , U ⊆ R nu >0 , X ⊆ R nx >0
and V ⊆ R m , respectively. We assume that e and u are system inputs, i.e. independent variables whose values can be fixed at will. We make the following standing assumption.

Assumption 1. For every p ∈ P , e ∈ E and u ∈ U, the solution to the system of equations

0 = Nv (14a) v = diag(e) • (a + B x • ln x + B u • ln u) (14b) is unique in x ∈ X and v ∈ V .
This guarantees that, for every admissible parametrization and system input, a steadystate exists and is unique. In order for this steady-state to be observable experimentally, we also make the assumption that it is locally asymptotically stable. In accordance with the metabolic control theory literature [START_REF] Heinrich | The Regulation of Cellular Systems[END_REF], fluxes v at steady-state are denoted by J. We write J p (e, u) to emphasize dependence on inputs and model parameters.

For varying values of p, Assumption 1 enables us to express the linlog model with parameters p as a map

M p : E × U → V × X : (e, u) → J p (e, u), x p (e, u) . ( 15
)
Assumption 1 is met, in particular, when matrix N diag(e)B x is invertible. In this case, one may write the output (J p , x p ) = M p (e, u) as an explicit function of the input (e, u),

J p (e, u) = diag(e) • (a + B x • ln x p (e, u) + B u • ln u), (16) ln x p (e, u) = -(N diag(e)B x ) -1 • N diag(e) • (a + B u • ln u). ( 17
)
As can be easily verified, this requires that the stoichiometry matrix N is full row rank, which is the case for systems with no mass conservation constraints [START_REF] Heinrich | The Regulation of Cellular Systems[END_REF].

In agreement with Section 2, where the identification problem is split into the identification of each reaction separately, we look at the identifiability of the parameters of the generic ith reaction, and say that a model is identifiable if all its reactions are.

Identifiability from a structural perspective

We adapt the definition from [START_REF] Ljung | System Identification, Theory for the User[END_REF] to our context as follows. Recall that p i is the ith column of p, i.e. the parameter vector for reaction i. Definition 1. A reaction i of model M p is structurally identifiable at p * if there exists an input set D ⊆ E × U such that, for all p ∈ P ,

(J p ) i , x p | D = (J p * ) i , x p * | D ⇒ p i = p * i . (18)
Here (J p , x p ) is seen as a function from E × U to X × V , and "| D " is its restriction to the input set D. In words, a reaction i is considered identifiable for a particular model parametrization p * if no p ∈ P with p i = p * i exists such that the predictions of M p and M p * are identical over all possible input sets D. Note that this definition is applicable to any metabolic reaction model, provided suitable definition of the parameters of the model class. In particular, it applies to the linlog form of the reaction rates as well as to any other pseudo-linear form reviewed in Section 2.

How can we apply Def. 1 to the analysis of identifiability of linlog models? The following proposition establishes a link between this definition and the uniqueness of the solution to Problems 2-3. Given the input set D = {(e 1 , u 1 ) , • • • , (e q , u q )} and a "true" parameter vector p * , let J k * and x k * denote the outputs J p * (e k , u k ) and x p * (e k , u k ), respectively, with k = 1, . . . , q. Proposition 1. A reaction i of M p is structurally identifiable at p * if and only if there exists D = {(e 1 , u 1 ) , • • • , (e q , u q )} ⊆ E × U such that the solution of the equation

W * i = Y * B * i , with W * i = J 1 * e 1 -J * e i • • • J q * e q -J * e i T , Y * = ln x 1 * -ln x * • • • ln x q * -ln x * ln u 1 -ln u • • • ln u q -ln u T , is unique in the parameters B * i = [B x * B u * ] T i . Corollary 1. A reaction i of M p is structurally identifiable at p * if and only if there exists D = {(e 1 , u 1 ) , • • • , (e q , u q )} ⊆ E × U such that Y * C(i) is full column-rank.
To complete the link with Section 2, notice that the matrices Y * and W * computed from the inputs and outputs of the model with true parameters p * coincide with the noiseless measurement matrices Y and W , respectively.

It is clear that the rank condition of Corollary 1 can be fulfilled only for a number of experiments q = |D| greater than or equal to the number of unknown parameters n i = |C(i)| of reaction i. The possibility to find q ≥ |C(i)| experiments making Y * C(i) full column rank depends on the network model and parameters themselves. Indeed, in our framework, the experimentalist can impose different enzyme concentrations e and inputs u, but the resulting metabolite concentrations are determined by the network. In other words, there is no full control of the regression matrix Y * , which impairs the design of optimal experiments for parameter regression. We show this by a simple example. where different values of a 1 ∈ R and B 2,1 ∈ R <0 (the coefficient that determines the strength of the feedback regulation) will be considered. For all values of the enzyme concentrations e i > 0, with i = 1, 2, 3, and all a 1 , B 2,1 , the equation Nv(x, e) = N diag(e)(a + B x ln x) = 0 yields a unique solution ln x = -(N diag(e)B x ) -1 N diag(e)a. This defines the unique steadystate of the system. Provided it is asymptotically stable, this gives us a steady-state of the system that can be observed experimentally. One first consideration is that different values of a 1 and B 2,1 may lead to very different properties of the matrix Y * even when this remains full rank, i.e. the system is structurally identifiable. For a 1 = 0.0297 and values of B 2,1 equal to -0.0073 (weaker feedback action) and -7.2961 (stronger feedback action), respectively, scatter plots of the steady-state solutions for ln x from 1000 randomly generated samples of e are reported in Fig. 2. Steady-state metabolite concentrations in the case of weak feedback are spread similarly in all directions, while with stronger feedback they are essentially aligned along a one-dimensional line. Here the strong feedback exerted by metabolite X 2 on the production of X 1 induces a negative correlation between their concentrations, which may result in an ill-conditioned estimation problem. In addition this strong feedback results in a near homeostasis of X 2 that may also impede identification. These points will be further developed in the next example. Second, some pathological parameterizations may give rise to a nonidentifiable model. Indeed, if a is in the span of B x , then the unique solution of N diag(e)(a + B x ln x) = 0 always corresponds to the value of ln x satisfying B x ln x = -a, independently of the value of e. Thus, no matter the number of experiments q, the rank of Y * is at most 1 and the model is not identifiable. In our example, this is the case for a 1 = -7.5491 and B 2,1 = -7.2961.

From the example above it is clear that a reaction may be nonidentifiable for specific values of the models parameters even if it is identifiable for other parameterizations. In the light of this, a generalization of Def. 1 from reaction identifiability at a parameter p * to reaction identifiability tout court can be obtained following [START_REF] Walter | Identification of Parametric Models[END_REF]. Namely, we stipulate that a reaction (and by extension, a model) is identifiable if it is identifiable almost everywhere in P , i.e. at almost every parametrization p * ∈ P of M p . Here 'almost everywhere' and 'almost every' are interpreted in terms of a suitable (e.g. Lebesgue) measure on P . Hence, the negative feedback network structure of Example 3 is identifiable in the sense of Walter and Pronzato. The identifiability criterion of Def. 1 holds except for the "rare" parameter combinations p * such that a ∈ span(B x ).

A second observation, following from the example above, is that the mathematical conditions that the system must fulfill to be declared nonidentifiable are too strong to be useful in practice. If we look at Fig. 2, we see that strong collinearities exist between the metabolite concentrations x 1 and x 2 . As a result, an unreasonably large number of experiments would be needed to resolve the effects of the two. Moreover, the definition of identifiability assumes that the measurements are not corrupted by noise, which is even less realistic. We therefore need to weaken our definition of identifiability in order to make it more suitable for applications to actual data on metabolism. While taking into account realistic assumptions on the experimental datasets, i.e., measurements available in a limited amount and affected by experimental error, this notion of identifiability should draw upon the theoretical notion of model identifiability discussed above.

Identifiability from a practical perspective

Let D be a fixed set of q inputs (external metabolites and enzyme concentrations), and let O be the set of the corresponding system outputs (fluxes and steady-state concentrations of internal metabolites determined by M p * ). Consider the problem of estimating the parameters B * i of reaction i given observations of D and O affected by measurement error. An estimator Bi of B * i is a function of the observations of D and O, well-defined for every possible (a priori unknown) value of p * (compare (Ljung, 1999, §7.4)). Since, due to noise, the observations are stochastic variables, Bi is itself a stochastic variable. Therefore, one cannot hope to estimate B * i exactly, but only within a certain degree of approximation. In this spirit, we define identifiability in terms of the existence of an estimator satisfying prespecified statistical requirements. In doing this, we restrict attention to the nonzero entries of B i , i.e., B C(i),i . Let B i ⊂ R n i be a bounded neighbourhood of the origin, and let α ∈ (0, 1).

Definition 2. For a given D ⊆ E × U, reaction i of M p is practically identifiable at p * with uncertainty region B i and confidence level 1 -α if there exists an estimator BC(i),i such that

P p * [ BC(i),i -B * C(i),i ∈ B i ] ≥ 1 -α, (19)
where P p * is the probability measure induced by M p * .1

Note that this definition is conceptually different from the one suggested by [START_REF] Raue | Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood[END_REF], where the definition of practical identifiability requires that the uncertainty on the parameter estimates (as defined via the profile likelihood) is bounded, but contrary to our definition, can be arbitrarily large. In addition, the definition in [START_REF] Raue | Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood[END_REF] is given in terms of a specific, not necessarily optimal choice of the estimator.

The point of view expressed by Def. 2 is that the experimentalist, or the modeler, sets the requirements (estimation accuracy and confidence level) that the estimates must fulfill in order to be useful, via the a priori specification of B i and α. Then, the possibility of fulfilling ( 19), i.e. the practical identifiability of the model, depends on the system itself and on the richness of the input set D. In general, the larger the D, the tighter the requirements that one can fulfill (i.e. the smaller the values of B i and α for which practical identifiability in the sense of Def. 2 holds).

From an alternative viewpoint, one may start from a given input set D, and look for the choices of α and B i that ensure satisfaction of ( 19). Here in turn, one may fix α and look for the B i that makes (19) achievable, or fix the acceptable estimation uncertainty B i and establish at what confidence level α this performance can be attained.

In all of the above cases, the natural questions that arise are how Def. 2 can be verified in practice, how this notion of identifiability depends on the structural system identifiability discussed in the previous section, and what B i may look like. To answer these questions, the relation between observations and observed quantities must be specified. We refer to the measurement model introduced in Section 2. For the sake of simplicity, we assume for this section that ν = 0, i.e. we address Problem 2. Problem 3 can be addressed with the same tools, but at the price of technical complications.

The following proposition answers the questions above.

Proposition 2. If a reaction i of M p is structurally identifiable at p * in the sense of Def. 1 then, for every α ∈ (0, 1), it is practically identifiable in the sense of Def. 2 with confidence level at least 1 -α for any uncertainty set

B i ⊇ E b Σ (α), where E b Σ (α) denotes the (1 -α)-confidence ellipsoid of a zero-mean Gaussian distribution with variance Σ = (Y T C(i) Σ -1 ε i Y C(i) ) -1 .
The proof relies on the use of minimum variance estimators, as dictated by standard results in linear estimation theory (see (Ljung, 1999, Appendix II) and also Appendix B).

From now on, estimation will be discussed and performed based on this type of estimator.

Observe that E b Σ (α), and hence the shape and size of the uncertainty regions B i for which the model is practically identifiable, depends on the choice of inputs D. In particular, the noise on W i affects the covariance matrix Σ by its statistics Σ ε i , while the contribution of the data Y C(i) is apparent. The number of data points q enters the picture in terms of the size of the matrix Σ ε i and the number of rows of Y C(i) . Typically, the larger q, the smaller B i can be for a fixed α. We argue that similar identifiability results can be derived even in cases where the noise is not Gaussian and metabolite measurements are affected by stochastic error, at the price of a more complicated characterization of B i . Finally, one may speak about identifiability of the whole model, e.g. by requiring that each reaction i is individually identifiable with a given confidence level α and uncertainty set B i ∈ R n i . Alternatively, one may require that all reactions be simultaneously identifiable with confidence level 1 -α and a suitably defined joint uncertainty set.

A discussion of practical identifiability in terms of covariance matrix of a (linearized) parameter estimation problem also appears in [START_REF] Srinath | Parameter identifiability of power-law biochemical system models[END_REF], in the context of power-law models. However, the discussion in [START_REF] Srinath | Parameter identifiability of power-law biochemical system models[END_REF] is essentially limited to one particular choice of the admissible estimation uncertainty B i , namely the one ensuring that the sign of the parameter values is estimated correctly with probability 1 -α.

A useful tool for better understanding Proposition 2 and the links between available data and practical identifiability is the Singular Value Decomposition (SVD) of a matrix [START_REF] Jolliffe | Principal Component Analysis[END_REF]. The SVD of Y C(i) is given by

Y C(i) = USV T , S = diag(s 1 , s 2 , ..., s n i ), (20)
with U ∈ R q×n i and V ∈ R n i ×n i orthonormal matrices and

s 1 ≥ • • • ≥ s n i ≥ 0 the singular values of Y C(i) .
In the presence of dependencies between the columns, there exists an index r with 1 ≤ r < n i such that s r+1 = . . . = s n i = 0, and Y C(i) is of rank r. Based on this, the covariance matrix of Proposition 2 can be written as Σ

= (Y T C(i) Σ -1 ε i Y C(i) ) -1 = V S -1 U T Σ ε i US -1 V T . Using the assumption that Σ ε i = σ 2 i I, the previous formula simplifies to Σ = V σ 2 i S -2 V T , (21)
which is (up to resorting of the entries) the SVD of Σ, with singular values σ 2 i s -2 1 , . . . , σ 2 i s -2 n i . Multiplied by a factor λ(α) fixed by α, the square roots of these values define the length of the axes of the confidence ellipsoid of Proposition 2. Now suppose that we seek parameter estimates that, with confidence 1 -α, fall within a ball B δ = {p : |p| < δ}, for some δ > 0. That is, all the entries of the parameter vector must be estimated with accuracy at least δ. From Proposition 2 and Def. 2, reaction i is practically identifiable if it is structurally identifiable and, for the given input set D, the ellipsoid E b Σ (α) associated with (21) fits into B δ , which happens if λ(α)σ i /s ℓ < δ for ℓ = 1, . . . , n i . In turn, since s 1 ≥ . . . ≥ s n i , this holds whenever λ(α)σ i /s n i ≤ δ, i.e. the smallest singular value of Y C(i) dictates the overall estimation performance.

If Y C(i) is ill-conditioned, i.e., some data vectors are nearly collinear, large discrepancies exist between its largest and its smallest singular values. Then, the condition s 1 ≫ s n i implies that, for practical identifiability, it must hold that λ(α)σ i /s 1 ≪ δ, i.e. a high accuracy in the estimation of the components of p along direction V 1 is required. Thus, achieving an even mild accuracy λ(α)σ i /s n i ≤ δ along direction V n i would generally require an unreasonable amount of experimental effort in terms of experimental replicas and/or measurement accuracy (see also [START_REF] Gutenkunst | Universally sloppy parameter sensitivities in systems biology models[END_REF]). Note, however, that high accuracy along V 1 is solely needed to ensure that the less accurate estimates of the components of p in direction V n i be acceptable. In a sense, this hard requirement is an artifact of the problem statement: If we accept that certain components are just not relevant, the remaining part of the model is identifiable in practice with good accuracy and much less experimental effort. To quantify our discussion, let us look at a numerical example.

Example 4. To illustrate the implications for parameter identifiability of a poorly conditioned data matrix, consider the estimation results from noisy and finite datasets for the two different identifiable parameterizations of the model of Example 3. As in the latter example, datapoints were simulated from random values of enzyme concentrations. Noise was added to W by drawing values from normal distributions with standard deviations proportional to the corresponding elements of W . Two different dataset sizes (q = 20 and q = 100) and two different noise levels (20% and 50% of the standard deviation) were tested, resulting in a total of 4 experimental scenarios for each model parameterization. For each scenario, 100 datasets were simulated and the corresponding estimates for reaction 2 are reported in the scatter plots of Fig. 3(a) (a 1 = 0.0297 and B 2,1 = -0.0073) and 3(b) (a 1 = 0.0297 and B 2,1 = -7.2961).

An immediate observation is that the shape of the 95%-confidence ellipse of the parameter estimates is different for the two model parameterizations. While estimation accuracy for B 1,1 and B 2,1 is comparable in the case of weaker feedback (B 2,1 = -0.0073), the shape of the uncertainty ellipse becomes very skewed in the case of stronger feedback (B 2,1 = -7.2961).

In particular, in the latter case estimation accuracy is much higher for B 1,1 than for B 2,1 regardless of the features of the dataset because of the strong homeostasis on x 2 (note the change in scale of the vertical axes of the plots in Fig. 3(a) and Fig. 3(b)). The plots also show that larger and/or less noisy datasets improve estimation performance, as expected. However, in the case of B 2,1 in Fig. 3(b), this improvement is seen to require extremely large and high-quality datasets. In other words, accurately estimating all parameters of the model demands a significant increase in experimental effort, even when most individual parameters are easy to estimate.

Following upon Example 3, it is apparent that the skewed estimation uncertainty in part (b) is related to the poor conditioning of the data matrix Y in the case of stronger feedback (the shape of the ellipsoid is determined by the ratio of the singular values of Y ). In terms of practical identifiability, assuming a modeler has set a maximum allowable uncertainty B 1 for some confidence level α, it is clear that in this case the system will not be practically identifiable (even if the model is structurally identifiable at the given p * ), unless B 1 is large enough, i.e. rather sloppy estimates are deemed acceptable.

To summarize the main points of the section, we have discussed practical identifiability as a relative concept that depends on the parameter estimation uncertainty that is deemed acceptable. If this is compatible with the quality of the data (dataset size, amount of noise) and the dataset is sufficiently diverse (more independent components than unknown parameters), then practical identifiability follows from structural identifiability. On the other hand, In each panel, scatter plots are reported for four different experimental scenarios: 20% noise and q = 20 datapoints (blue); 20% noise and q = 100 datapoints (green); 50% noise and q = 20 datapoints (red); 50% noise level and q = 100 datapoints (magenta). 95%-confidence ellipses are drawn for each scenario if the experiments are not informative enough or the network itself implies heavy correlations among the data, we detect lack of practical identifiability from the existence of nearly zero singular values in the decomposition of the regression matrix. For metabolic systems this will occur when feedback regulation results in strong homeostasis or when metabolite concentrations are correlated, which could be due in general to steady-state constraints. A particularly relevant situation concerns reactions that operate close to equilibrium. Indeed in this situation mass-action law relates metabolite concentrations as follows:

j N j,i ln x j ≈ ln K i (22)
where K i is the equilibrium constant of reaction i. This results in a quasi-dependency between the ln x j that makes the reaction nonidentifiable in practice.

In addition to providing little information for the estimation of the model parameters, the smallest components of Y C(i) have negligible effect on the solution of the steady-state equations ( 3)-( 6), i.e. in determining the system steady-state. In the next section, we will build upon this analysis to define a criterion for eliminating the components of Y C(i) associated with the smallest singular values and reduce the network model accordingly. This will make parameter estimation (regression) a well-conditioned problem for every single reaction while minimally affecting the quality of the model.

Reduction to identifiable models

It was shown in Example 4 that attempting to identify the parameters of a reaction that is not practically identifiable leads to an ill-posed estimation problem. That is, certain parameter combinations are practically indistinguishable from the data. Eliminating redundant components of the model parameters by Principal Component Analysis (PCA, see [START_REF] Jolliffe | Principal Component Analysis[END_REF]) is a way to ensure well-posedness of parameter estimation (i.e. practical identifiability) by a "minimal" approximation of the model.

The method applies as usual reaction by reaction. For any given reaction i, with i = 1, . . . , m, we compute and manipulate the SVD of the matrix ỸC(i) , so as to get transformed data and a corresponding model with a reduced number of parameters that can be estimated reliably.

Identifiability analysis and model reduction by PCA

We start by considering the case where the regression matrix is noiseless, i.e. ỸC(i) = Y C(i) , and rank Y C(i) = r, with r < n i . Notice that the latter is always the case for structurally nonidentifiable models. The extension of the method to practical identifiability (where Y C(i) is full column rank but ill-conditioned) and to noisy and incomplete data Ỹ will be discussed in the next sections.

Consider again the SVD Y C

(i) = U • diag(s 1 , s 2 , ..., s n i ) • V T , with s 1 ≥ s 2 ≥ . . . ≥ s n i ≥ 0. Since r < n i , it holds that s 1 ≥ . . . ≥ s r > 0 and s r+1 = . . . = s n i = 0. Then, Y C(i) has an (n i -r)-dimensional kernel K Y , given by K Y = range(V r+1:n i ). For any B C(i),i and any k Y ∈ K Y , it holds that Y C(i) • B C(i),i = Y C(i) • (B C(i),i + k Y ).
For the purpose of identification, this means that B C(i),i cannot be uniquely reconstructed from the data. On the other hand, range(Y

C(i) ) = range(Y C(i) V 1:r ), where Y C(i) V 1:r is full column rank. Then, for every B C(i),i there exists a unique Bi ∈ R r×1 such that Y C(i) • B C(i),i = Y C(i) V 1:r • Bi . This suggests to modify the regression problem W i = Y C(i) • B C(i),i + ε i into W i = Yi • Bi + ε i Yi = Y C(i) V 1:r (23)
which has a unique solution in Bi , i.e. Bi is identifiable. We call (23) the reduced model and Bi the reduced parameter vector.

For a fixed outcome of the noise ε i , from the unique solution B in the reduced parameter space one can infer a whole subspace of equivalent solutions in the original parameter space as

{B C(i),i = V 1:r • Bi +k Y , k Y ∈ K Y }.
Thus, in general, a fixed solution Bi does not determine uniquely any of the parameters B j,i (j being element of C(i)). However, depending on the structure of V , we may be able to isolate some parameters B j,i that can be reconstructed without ambiguity.

Proposition 3. Let index j be an element of C(i), i.e., for some ℓ = 1, . . . , n i , j = C ℓ (i). Suppose that the entries of V ℓ,r+1:n i are all zero. If Bi is the (unique) solution to (23), then B j,i = V ℓ,1:r Bi is uniquely determined.

A similar, but less general approach to separate identifiable from nonidentifiable parameters has been considered by [START_REF] Nikerel | Model reduction and a priori kinetic parameter identifiability analysis using metabolome time series for metabolic reaction networks with linlog kinetics[END_REF].

Model reduction put in practice

In a real setting, as shown in Example 4, small nonzero values of s r+1 , . . . , s n i can also make the problem of estimating B C(i),i ill-conditioned, thus preventing practical identifiability. In addition, measurement error can make certain components of the data indistinguishable from noise. The idea here is to remove the components of the parameters that are poorly determined from the data, thus ensuring smaller estimation uncertainty and hence practical identifiability in a reduced parameter space. In order to develop and explain the rationale of our method, we will first reconsider model reduction in the setting of Problem 2 where the metabolite data are assumed noiseless, and then move on to the more realistic scenario of Problem 3 where metabolite data are noisy. In the remarks concluding the section, we will briefly discuss the application of the method to datasets with missing or corrupted entries (e.g. outliers) and its biological interpretation. We will then summarize the model reduction procedure in Section 4.3.

The scenario of Problem 2. Here ỸC(i) = Y * C(i) . One may consider the rank-r approximation of the data matrix

Y * C(i) = U • diag(s 1 , s 2 , ..., s n i ) • V T ≃ U • diag(s 1 , s 2 , ..., s r , 0, . . . , 0 n i -r ) • V T = Ŷi .
Following the previous section, for Yi = Ŷi V 1:r , consider the reduced model

W i = Yi • Bi + ε i Yi = Ŷi V 1:r . ( 24
)
with unique least-squares solution for the reduced parameter Bi . With the same arguments as in Section 3.2, for Σ ε i = σ2 i I, one observes that the confidence ellipsoid associated with the estimate of Bi is determined by the matrix σ 2 i Y T i Yi -1 . In particular, the largest axis length, corresponding to the largest parameter estimation uncertainty, is proportional to σ i /s r , i.e. it has been reduced by a factor s r /s n i with respect to the original estimation problem. This analysis suggests a criterion for the choice of r based on our definition of practical identifiability. Suppose that, with a given confidence 100 • (1 -α)%, the admissible uncertainty B i is a ball of radius δ. Since the radii of the estimation error confidence ellipsoid are given by λ(α)σ i /s r ≥ . . . , ≥ λ(α)σ i /s 1 it suffices to choose r as the minimum value for which λ(α)σ i /s r ≤ δ for the reduced model to be practically identifiable. If this holds for r = n i , the full model is practically identifiable and needs no reduction.

The scenario of Problem 3. Here the noisy versions Ỹ of Y are the available data. The idea is to remove from the problem not only the components that make estimation ill-conditioned, but also those components that are detrimental in that they are dominated by noise. To do this, let us look at the empirical covariance matrix of the data Ỹ T C(i) ỸC(i) /q. From the approximation 2

Ỹ T C(i) ỸC(i) = Y T C(i) Y C(i) + η T C(i) η C(i) + Y T C(i) η C(i) + η T C(i) Y C(i) ≃ Y T C(i) Y C(i) + qΣ η C(i) , where Σ η C(i) = ν 2 I, it follows that Ỹ T C(i) ỸC(i) /q ≃ U diag(s 1 , . . . , s n i )V T T U diag(s 1 , . . . , s n i )V T /q + ν 2 I = V diag(s 2 1 , . . . , s 2 n i )/q + ν 2 I V T .
The expression after the last equality sign is clearly the SVD of Ỹ T C(i) ỸC(i) /q, with singular values s2 ℓ = s 2 ℓ /q + ν 2 , with ℓ = 1, . . . , n i , composed of the signal contribution s 2 j /q and the noise contribution ν 2 . In the light of this, to remove the components of the data dominated by noise, we compute the (noisy) singular values s2

1 ≥ s2 2 ≥ . . . ≥ s2 n i ≥ 0 from the SVD of Ỹ T C(i) ỸC(i) /q, draw estimates ŝ2 ℓ of the true (noiseless) singular values s 2 ℓ by posing ŝ2 ℓ = max(0, s2 ℓ -ν 2 ) for ℓ = 1, . . . , n i , and define what we call the "effective rank" of the data matrix as follows.

Definition 3. The effective rank of the data matrix is

r = max{ℓ : ŝ2 ℓ ≥ ν 2 , ℓ = 1, . . . , n i } (25)
According to this definition, the effective rank indicates the number of independent components that can be safely distinguished in the data in that not blurred by noise. Notice that noise, by its very nature, tends to decorrelate all matrix entries. Following on the discussion for the scenario of Problem 2, this criterion may also be seen as implementing model reduction for practical identifiability, with a choice of the uncertainty region B i depending on ν, i.e. adapted to the presence of noise on metabolite data.

An alternative approach to determine the effective rank of a data matrix, useful when ν is assumed small but not known with certainty, is to remove the components associated with the smallest singular values by setting r so that a suitable proportion θ ∈ (0, 1) of the total variance n i ℓ=1 s2

ℓ of the data is retained [START_REF] Berthoumieux | Identification of linlog models of metabolic networks from incomplete high-throughput datasets[END_REF]. This gives rise to the following definition of effective rank.

Definition 4. The 100 • θ%-variance effective rank of the data matrix is

r = min r ′ : r ′ ℓ=1 s2 ℓ ≥ θ • n i ℓ=1 s2 ℓ , r ′ = 1, . . . , n i . ( 26
)
Different from the previous definition, effective rank is intended here simply as the number of components needed to express most of the data content. When applied to data with small noise, data components dominated by noise are also small and are hence excluded from the count. For large noise levels, this reasoning no longer applies. Note that precise knowledge of the noise variance ν 2 is not required here, at the price of a rather uninformed choice of parameter θ.

In both cases, after computing the effective rank r, the original model can be replaced by the reduced model ( 24), providing us with a well-behaved model for the subsequent identification of the system.

Example 5. We have seen in Example 4 that the first reaction of the feedback model with a 1 = 0.0297 and B 2,1 = -7.2961 is not practically identifiable and we want to see which approach to the choice of r enables PCA to detect this property on a limited noisy dataset. In order to mimic available experimental data [START_REF] Ishii | Multiple high-throughput analyses monitor the response of E. coli to perturbations[END_REF], noise was added to the data matrix Y by drawing q = 30 values from a normal distribution with standard deviation of 0.4, corresponding approximately to 40% noise for metabolite concentrations. 100 datasets were generated in this way and PCA was performed on each of them. Two different criteria for the choice of the reduced model order r were tested, based on Def. 3 and Def. 4, respectively. Fig. 4(a) shows the estimates of the squared singular values and the cutoff of ν 2 proposed by Def. 3, while Fig. 4(b) shows the normalized cumulative sums of the squared singular values and the cutoff of 0.99 proposed by Def. 4. The second squared singular value is always smaller than ν 2 , so that the model is correctly and consistently found nonidentifiable with the first definition, contrary to what is found with the second definition 57 times out of 100.

The above example thus illustrates that the criterion for model reduction taking into account the noise level, when applicable, is more relevant.

Remark 2. The data matrix ỸC(i) may suffer from the lack of certain data entries, typically due to the removal of outliers or faults of the experimental machinery. A simple but wasteful option to recover a full data matrix for later use in a well-defined model reduction/parameter estimation problem is to discard those data points Ỹk,C(i) , and the corresponding flux information Wk,i , suffering from the absence of some entry. In [START_REF] Berthoumieux | Identification of linlog models of metabolic networks from incomplete high-throughput datasets[END_REF], in the context of parameter estimation, we have proposed methods compensating for the missing entries by the use of statistical priors inferred from the available data. For the sake of model reduction, which requires in our approach the SVD of the data matrix, completion of the data matrix by a suitable imputation method was suggested. Several imputation methods can be considered, still relying on statistics from the available data, such as multiple imputation or completion by the mean of the available metabolite data (see [START_REF] Berthoumieux | Identification of linlog models of metabolic networks from incomplete high-throughput datasets[END_REF] and references therein). As an appealing alternative we cite "minimal rank" SVD, which has been developed and applied in [START_REF] Brand | Incremental singular value decomposition of uncertain data with missing values[END_REF] for reduced-order modelling in computer vision.

Remark 3. The reduced model (24) is expressed in terms of parameters that are linear combinations of the original parameters. As a consequence, the results of identification may be difficult to interpret from a biological point of view. Proposition 3 suggests a way to isolate identifiable parameters in nonidentifiable reactions, and thus extract partial, but unambiguous information from the dataset. Unfortunately, the condition of the proposition is not usually verified in practice since, due to noise, the entries of the data kernel-generating matrix V r+1:n i will never be exactly 0. In order to ease the interpretation of the results, one may relax this condition as follows. Bearing in mind that V r+1:n i is composed of unit-L 2 -norm column vectors, we consider negligible all entries of V r+1:n i whose square is below a threshold ρ 2 significantly smaller than 1. As we shall see, in several cases, this allows us to clarify the biological interpretation of the estimation results. Further study of the kernel-generating matrix V r+1:n i would yield a theoretically more sophisticated criterion, but we will not pursue this discussion here.

The overall procedure for identification analysis and model reduction

Based on the discussion of the previous sections, here we summarize the procedure for obtaining a practically identifiable approximate kinetic model from noisy and incomplete datasets. The procedure is also summarized in Fig. 5. Given noisy steady-state metabolite data ln x1 , . . . , ln xq :

• Compute the data matrix Ỹ ;

• In case of missing entries, complete the matrix by a method of choice (multiple imputation, minimum rank completion etc.);

• For every reaction i = 1, . . . , m:

1. Extract from Ỹ the data submatrix ỸC(i) ; 2. Compute the SVD of the empirical data covariance matrix,

Ỹ T C(i) ỸC(i) /q = V diag(s 2 1 , . . . , s2 n i )V T ;
3. Compute the effective data rank r = max{ℓ : s2 ℓ -ν 2 ≥ ν 2 }; 4. Compute Ŷi , the data matrix obtained by discarding the

n i -r smallest compo- nents, as Ŷi = Y C(i) • V 1:r 0 n i ×(n i -r) V T , 0 n i ×(n i -r) being the n i × (n i -r) null matrix; 5. Return the reduced model W i = Yi • Bi + ε i , with Yi = Ŷi V 1:r .

Applications of the model identifiability and reduction approach

Application to a network with simulated data

In order to evaluate performance of our identifiability and model reduction procedure, we now discuss its application to a more realistic simulated network originally presented in [START_REF] Visser | Dynamic simulation and metabolic re-design of a branched pathway using linlog kinetics[END_REF] and depicted in Fig. 6(a). The network involves n x = 8 internal and n u = 3 external metabolites, participating in a total of m = 8 reactions. We developed a linlog model of the network based on the state and input vectors x and u whose entries are listed in Fig. 6 .

The values of B u were taken without changes from [START_REF] Visser | Dynamic simulation and metabolic re-design of a branched pathway using linlog kinetics[END_REF], while the values of a and B x were adapted from the same paper. The stoichiometry matrix N, given in Eq. ( 27) below, is fixed by the ordering of the input and state vector entries and the scheme in Fig. 6(a). The row rank of this matrix is equal to 6, corresponding to 2 mass conservation constraints (see also Eq. ( 28) below). Following the analysis of [START_REF] Reder | Metabolic control theory: A structural approach[END_REF], it is possible to factor the matrix into a link matrix L expressing dependencies between concentrations and a reduced-order full-row rank matrix N corresponding to stoichiometries of independent data matrix Yi = Ŷi V 1:r

Return the reduced

Reduced model for reaction i identifiable

End

Missing data

No missing data metabolites. Factorization is non-unique. In our case, one such factorization gives

Set i = i + 1 identifiable Reaction i r = n i r < n i i < m i < m Compute Ỹ (missing imputation, minimal rank,...) Complete Ỹ Start Compute SVD of Ỹ T C(i) ỸC(i) /q Compute effective data rank r i = m End Reaction i nonidentifiable i = m Set i = 1 Extract ỸC(i) Compute Ŷi ∈ R q×ni

Internal metabolites External metabolites Variable

Metabolite Input Metabolite x 1 M1 u 1 S x 2 M2 u 2 P1 x 3 M3 u 3 P2 x 4 M4 x 5 M6 x 6 A x 7 M5 x 8 AH (b)
N = 2 6 6 6 6 6 6 6 6 6 4 1 -1 0 0 0 0 0 0 0 1 -1 0 0 -1 0 0 0 0 1 -1 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 1 -1 0 -1 0 1 0 0 1 0 0 0 0 0 0 1 -1 0 0 1 0 -1 0 0 -1 0
1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 -1 -1 0 0 -1 0 0 0 0 0 -1 3 7 7 7 7 7 7 7 7 7 5 • 2 6 6 6 6 6 4 1 -1 0 0 0 0 0 0 0 1 -1 0 0 -1 0 0 0 0 1 -1 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 1 -1 0 -1 0 1 0 0 1 0 3 7 7 7 7 7 5 (27) = L • N .
With this factorization, the entries of x 1:6 (M1, M2, M3, M4, M6 and A) are treated as independent quantities, and determine the values of x 7 (M5) and x 8 (AH) via conservation of mass. That is, for some fixed constants

T 1 , T 2 ∈ R >0 ,      ẋ1:6 = N diag(e)(a + B x • ln x + B u • ln u), T 1 = x 2 + x 3 + x 6 + x 7 , T 2 = x 6 + x 8 . (28)
Notice that this reformulation allows one to compute the steady-state values of all system variables by setting the differential part to zero. In our case, the method is used to compute steady-state data as a function of enzyme and external metabolite concentrations.

To assess identification performance, we considered the scenario of [START_REF] Visser | Dynamic simulation and metabolic re-design of a branched pathway using linlog kinetics[END_REF], where the external metabolite concentrations are fixed to u = 1 0.1 0.2 T , T 1 = 0.3 and T 2 = 0.1. Since u is fixed, the parameter matrix B u is obviously nonidentifiable, and the contributions of a and B u ln u are indistinguishable. To circumvent this issue, we define the lumped constant term a ′ = a + B u ln u so as to obtain the modified linlog reaction rate model v = diag(e) • (a ′ + B x ln x), and study the identifiability and reduction of the latter in terms of a ′ and B x . Identifiability analysis and model reduction are performed in accordance with Section 4.3 on R = 100 randomly generated datasets Ỹ and W with realistic statistical properties. Each dataset shares the same steady-state values Y and W computed for q = 30 different values of enzyme concentrations, generated once as in Example 5, and each differs in the randomly generated 40% noise corrupting the measurements. We generated the results from the application of both Def. 3 and Def. 4, with θ = 0.99. The results are depicted in Fig. 7 and are also reported in Table 1.

First, we notice that for every reaction, the effective rank computed with the criterion of Def. 4 is higher than the one computed with the criterion of Def. 3. Thus, the latter criterion gives more conservative results, in the sense that, on average, fewer reactions are deemed identifiable. Except for reaction 2, application of Def. 4 returns the full size of the reaction for at least 96 out of 100 datasets. That is, in this case, the ability of the criterion to detect dependencies among data is very limited. This can be explained by the presence of a significant amount of noise on metabolite data, which tends to decorrelate the observations. The criterion based on Def. 3 detects dependencies among the data for all reactions. Indeed, the effective rank determined by this criterion is consistently smaller and differs from the results of Def. 4 by an average of about 1 for reactions 1, 5 and 8, and of about 2 for the other reactions. This can be attributed to the compensation of noise in the computation of the singular value estimates ŝ2 ℓ in (25), which relies upon and exploits the knowledge of the noise level ν. This is apparent in Fig. 7, where the blue dots representing the singular value estimates drawn from each of the 100 datasets Ỹ are correctly concentrated around the true singular values from Y . Fig. 7 also clarifies the non-integer average rank values of Table 1 (notably for reactions 4, 5 and 6). This comes from the fact that the singular values lying close to the significance cutoff value ν 2 are estimated above or below this threshold depending on the simulation run. For instance, the estimates of the second singular value of reaction 6 were smaller than ν 2 in 70 of the 100 runs.

Finally, the results for reactions 5 and 8 reveal a fundamental difference between Def.s 3 and 4. The former method leads to the conclusion that no modeling is possible (the effective ℓ -ν 2 computed from 100 noisy datasets Ỹ of size q = 30 (blue dots) of the model of Fig. 6. The dashed bars correspond to the singular values s 2 ℓ /q from the noiseless dataset Y of size q = 30. In red is the area below ν 2 ; all singular values in this area are considered negligible.

rank is estimated to be zero in most runs) since the corresponding metabolite data are dominated by noise, whereas Def. 4 provides effective rank estimates that are by construction lower-bounded by one.

From the results of this section, it is clear that identifiability analysis and model reduction in the presence of noise should be performed on the basis of the effective rank of Def. 3, which outperforms Def. 4 and returns consistent results (Table 1 and Fig. 7). The actual application of the method and the effects of using Def. 3 in place of Def. 4 on a real dataset are discussed in the next section.

Application to central carbon metabolism of E.coli

As a second example, we illustrate the application of our method to a complex network of biochemical reactions involved in carbon assimilation in the enterobacterium Escherichia coli. The network we consider gathers enzymes, metabolites and reactions that make up the bulk of central metabolism, including glycolysis, the pentose-phosphate pathway, the tricarboxylic acid cycle and anaplerotic reactions such as glyoxylate shunt and PEP-carboxylase (Fig. 8). The network has been studied for a long time from different perspectives, which makes it an ideal model system for our purpose. The structure of the E. coli carbon metabolism network is known in a rather precise way, its dynamics have been modeled by means of a variety of formalisms ( [START_REF] Bettenbrock | A quantitative approach to catabolite repression in Escherichia coli[END_REF][START_REF] Kotte | Bacterial adaptation through distributed sensing of metabolic fluxes[END_REF] and references therein), and recently a high-throughput dataset containing the required information for addressing Problem 3 has been published [START_REF] Ishii | Multiple high-throughput analyses monitor the response of E. coli to perturbations[END_REF].

We now investigate which reactions are identifiable following the criteria of Section 4, given the available experimental data and a linlog model of the network. From a methodological point of view, we are interested in analyzing the differences between the results [START_REF] Berthoumieux | Identification of linlog models of metabolic networks from incomplete high-throughput datasets[END_REF]. The map shows metabolites (bold fonts) and genes (italic). Abbreviations of metabolites are glucose (Glc), glucose 6phosphate (G6P), fructose 6-phosphate (F6P), fructose 1-6-biphosphate (FBP), dihydroxyacetone phosphate (DHAP), glyceraldehyde 3-phosphate (G3P), 3-phosphoglycerate (3PG), phosphoenolpyruvate (PEP), pyruvate (Pyr), 6-phosphogluconate (6PG), 2-keto-3-deoxy-6-phospho-gluconate (2KDPG), ribulose 5-phosphate (Ru5P), ribose 5-phosphate (R5P), xylulose 5-phosphate (X5P), sedoheptulose 7-phosphate (S7P), erythrose 4-phosphate (E4P), oxaloacetate (OAA), citrate (Cit), isocitrate (IsoCit), 2-keto-glutarate (2KG), succinate-CoA (Suc-coA), succinate (Suc), fumarate (Fum), malate (Mal), glyoxylate (Glyox), acetyl-CoA (Ac-coA), acetylphosphate (Acp) and acetate (Ace). Cofactors impacting the reactions are not shown: adenosine triphosphate (ATP), adenosine diphosphate (ADP), nicotinamide adenine dinucleotide phosphate (NADP) and its reduced form (NADPH), nicotinamide adenine dinucleotide (NAD) and its reduced form (NADH) and flavin adenine dinucleotide (FAD). The gene names are separated by a comma in the case of isoenzymes, by a colon for enzyme complexes, and by a semicolon when the enzymes catalyze reactions that have been lumped together in the model. obtained with Def. 3 and those obtained with the method of Def. 4 used in [START_REF] Berthoumieux | Identification of linlog models of metabolic networks from incomplete high-throughput datasets[END_REF]. From a biological point of view, we wish to understand how much information is actually contained in a state-of-the-art dataset for the purpose of parameter estimation.

The dataset used for identification of the network in Fig. 8 was obtained by experiments with 24 single-gene deletions that were grown at a fixed dilution rate of 0.2h -1 in a glucose-limited chemostat, and with wild-type cells at 5 different dilution rates [START_REF] Ishii | Multiple high-throughput analyses monitor the response of E. coli to perturbations[END_REF]. The authors collected data using multiple high-throughput techniques, in particular DNA microarray analysis and two-dimensional differential gel electrophoresis (2D-DIGE) for genes and proteins, capillary electrophoresis time-of-flight mass spectrometry (CE-TOFMS) for metabolites, and metabolic flux analysis. They thus obtained a steady-state dataset consisting of metabolite concentrations, mRNA and protein concentrations for the enzymes, and metabolic fluxes under 29 different experimental conditions. Therefore this dataset contains the information for setting up a parameter estimation problem as defined in Section 2.

We carried out the identifiability analysis for the linlog model developed in [START_REF] Berthoumieux | Identification of linlog models of metabolic networks from incomplete high-throughput datasets[END_REF]. This model is a translation of the reaction scheme of Fig. 8 into linlog rate equations (2). When certain metabolites could not be measured, preventing their inclusion in the model, we lumped together the reactions in which they are involved. In addition to the above model simplification, imposed by the available data, we added a phenomenological reaction to model biomass production. The resulting model has n x = 16 internal metabolites, n u = 7 external metabolites and measured cofactors, and m = 31 reactions [START_REF] Berthoumieux | Identification of linlog models of metabolic networks from incomplete high-throughput datasets[END_REF].

A complication for determining the identifiability of reactions and finding a suitable model reduction is that the dataset contains a large amount of missing data. In particular, certain metabolites could not be measured in up to 80% of the experimental conditions (28% on average for all metabolites). Following Remark 2 of Section 4, we therefore completed the dataset by means of multiple imputation, generating 100 datasets to allow the computation of statistics and test the robustness of the results.

Table 2 summarizes the results of applying the reduction method of Section 4.3 to the model and the data. For each of the reactions in the model, we computed the average of the effective rank of the 100 completed datasets. The effective rank for the individual datasets was usually found to be the same (in at least 82 of the 100 datasets), which explains that the computed average values are close to integers. Remarkably, out of the 31 reactions in the model, only 4 were found to be fully identifiable: reactions 4, 5, 14, and 31. The first three reactions involve two metabolites: fructose 1-6-biphosphate (FBP) and dihydroxyacetone phosphate (DHAP) (reaction 4), DHAP and 3-phosphoglycerate (3PG) (reaction 5) and ribose 5-phosphate (R5P) and sedoheptulose 7-phosphate (S7P) (reaction 14). The identifiability of these reactions means that the method did not detect any dependencies between these pairs of metabolite concentrations. Reaction 31 involves a single metabolic variable acetyl-coA/coA (Ac-coA/coA). Among the remaining 27 nonidentifiable reactions, two reactions are known to operate close to equilibrium (reactions 2 and 7: Pgi and Gpm:Eno, respectively), which provides a robust rationale for their nonidentifiability whatever the physiological conditions that are explored in the dataset (cf. Eq. ( 22)). In the other nonidentifiable cases the effective rank is reduced by 1 (for 16 reactions), 2 (6 reactions), 3 (2 reactions), and 6 (1 reaction). The latter case concerns the growth-rate reaction, which has 11 variables. A striking observation on Tables 2 and3 In order to isolate identifiable parameters in nonidentifiable reactions, Proposition 3 proposes a criterion that has been relaxed in Remark 3 so as to make it applicable to noisy data. The approach is based on the choice of a threshold ρ 2 for neglecting components of the kernel-generating marix V r+1:n i extracted from the data. In what follows, to set a ground for discussion, we set ρ equal to 0.15. We verified that changes of this threshold within the range ρ ∈ (0.1, 0.2) do not significantly alter the results as reported in Table 3. In this table, parameters that were diagnosed as being identifiable in more than 50% of the completed datasets are highlighted in green. From the 27 nonidentifiable reactions, no individual parameter could be unambiguously extracted in 8 cases (reactions 2, 7, 8, 12, 15, 16, 22 and 24). Of the 72 parameters in the remaining 19 reactions, 30 are identifiable in more than half of the datasets. In particular, we observe that all parameters associated to glucose (Glc), DHAP, Ac-coA/coA, 6-phosphogluconate (6PG), R5P, FAD and acetate (Ace) are identifiable in the sense that no significant dependencies with other metabolites could be detected in the experimental conditions of [START_REF] Ishii | Multiple high-throughput analyses monitor the response of E. coli to perturbations[END_REF].

The results shown in Table 3 are different from those obtained in our earlier work [START_REF] Berthoumieux | Identification of linlog models of metabolic networks from incomplete high-throughput datasets[END_REF], where we used a method based on Def. 4 with θ = 0.99 instead of Def. 3. Indeed we previously found many more reactions to be identifiable (24 out of 31) although for most cases parameter estimates turned out to be unreliable because of large confidence intervals. Reactions 2 and 7 (Pgi and Gpm:Eno), that were classified as identifiable in [START_REF] Berthoumieux | Identification of linlog models of metabolic networks from incomplete high-throughput datasets[END_REF], are found here to be nonidentifiable in agreement with the fact that they operate close to equilibrium. Therefore, the results of the identifiability analysis in [START_REF] Berthoumieux | Identification of linlog models of metabolic networks from incomplete high-throughput datasets[END_REF] appear to be overly optimistic, i.e. overestimating the number of identifiable reactions because measurement errors on metabolite concentrations are not taken into account. Notwithstanding, the results of both analyses are consistent in the sense that all 7 reactions detected as nonidentifiable by means of Def. 4 also remain nonidentifiable according to Def. 3. ] and results of identifiability analysis and model reduction on the data of [START_REF] Ishii | Multiple high-throughput analyses monitor the response of E. coli to perturbations[END_REF] for the linlog model of E. coli central carbon metabolism. SVD has been applied to Y C(i) for each reaction and singular values were discarded based on Def. 3. Nonidentifiable parameters are shown in grey and identifiable ones, as well as identifiable reactions, in green. The percentages of cases for which the parameters were found identifiable are given. To avoid complications deriving from the presence of conserved moieties, some of the metabolites are modeled as ratios of metabolite concentrations, e.g. ATP/ADP. Reaction 27, labeled µ, is a phenomenological reaction for biomass production. The last row indicates the percentage of missing data per metabolite. Abbreviations are as in Fig. 8.

Discussion

A major, but often overlooked problem in the identification of metabolic network models is the identifiability of the parameters, and hence of the model. Informally speaking, the identifiability of a model (parameter) consists in the possibility to unambiguously reconstruct the model (parameter) from the observed behavior of the system. Identifiability problems may reside in the very structure of the model, notably the occurrence of implicit dependencies between parameters. These dependencies might be due to an inappropriate model formulation, constraints on kind of experimental perturbations that can be realized, and unobserved variables. In addition, practical identifiability problems may arise from limitations on the quality and quantity of available data, in particular the fact that experimental data in biology are frequently noisy, sparse and incomplete.

We have studied identifiability issues in the context of approximate kinetic modeling formalisms, notably linear-logarithmic (linlog) models. On the theoretical side, following the classical systems identification literature [START_REF] Ljung | System Identification, Theory for the User[END_REF][START_REF] Walter | Identification of Parametric Models[END_REF], we have first precisely defined the notions of structural (a priori) identifiability (Def. 1) and practical (a posteriori) identifiability (Def. 2). The latter notion is obviously related to the former, in the sense that structural nonidentifiability entails practical nonidentifiability. However, Proposition 2 goes beyond by saying that identifiability in the theoretical sense may also imply identifiability in the practical sense, provided that the uncertainty on the parameters, as determined by the available dataset, remains within the desired accuracy bounds. Notice that practical identifiability is thus not an absolute notion, but rather conditional on the data properties and the required model precision.

A second methodological contribution of this paper is the development of theoretically sound and practically applicable methods for the detection of identifiability problems and the transformation of a nonidentifiable model to a reduced identifiable model. In particular, we have formulated criteria based on the SVD of the matrix of log-transformed and centered measurements of metabolite concentrations. These criteria define the effective rank of the data matrix, corresponding to the number of parameters that can be safely distinguished from the output. The criterion privileged in this paper (Def. 3), contrary to a criterion that we proposed in earlier work (Def. 4), takes into account the estimated variance of the noise. The flow chart in Fig. 5 gives a step-by-step procedure for identifiability analysis and model reduction.

The identifiability of models of biological systems is a topic that has been much studied in mathematical biology, and that has received renewed attention in the context of systems biology (see [START_REF] Chappell | Global identifiability of the parameters of nonlinear systems with specified inputs: A comparison of methods[END_REF]; Chis et al. (2011b); [START_REF] Cobelli | Parameter and structural identifiability concepts and ambiguities: A critical review and analysis[END_REF]; [START_REF] Raue | Addressing parameter identifiability by model-based experimentation[END_REF] for reviews). Systems of biochemical reactions, which have the general form of Eq. ( 3) at steady-state, have a number of peculiarities for identifiability analysis. When reaction rates, enzyme concentrations and metabolite concentrations are measured, the identification problem can be decomposed into subproblems for each of the individual reactions. Determining the identifiability of a model then reduces to checking the identifiability of the reactions. If in addition the reaction rates are expressed in terms of linlog or other pseudolinear equations, identification becomes a linear or orthogonal regression problem, depending on whether noise on the metabolite concentrations is taken into account or not (Problem 2 and Problem 3, respectively). Identifiability analysis then amounts to checking for linear de-pendencies in a transformed data matrix, which can be done using standard techniques from linear algebra and statistics. Related ideas starting from this general approach can be found in other recent work [START_REF] Berthoumieux | Identification of linlog models of metabolic networks from incomplete high-throughput datasets[END_REF][START_REF] Srinath | Parameter identifiability of power-law biochemical system models[END_REF][START_REF] Nikerel | Model reduction and a priori kinetic parameter identifiability analysis using metabolome time series for metabolic reaction networks with linlog kinetics[END_REF]. However, our development in the current paper is different in fundamental ways, such as the very definition of structural and practical identifiability, and the application of SVD to detect and resolve identifiability problems.

Although linlog models have been central in this paper, the results directly carry over to the other approximate formalisms mentioned in Section 2. In addition, they also bear on more general classes of nonlinear models of metabolic networks. Indeed the parameters in the mean-removed linlog model of Eq. ( 6) are proportional to elasticity coefficients that describe the sensitivities of reaction rates to changes in metabolite concentrations. If a reaction in a linlog model is nonidentifiable, this means that elasticity coefficients are not identifiable, therefore any other class of kinetic models is liable to encounter similar identifiability issues.

The approach for determining the identifiability of linlog models proposed in this paper has been tested on a network with simulated data and applied to a high-throughput data set for central carbon metabolism in E. coli. The use of simulated data has made it possible to demonstrate that, for typical sizes of the dataset and realistic noise levels, our approach is able to correctly identify the principal components of the parameter vector (Fig. 7). Surprisingly, the determination of the effective rank of the datasets for the different reactions in the E. coli metabolic network shows that only a small fraction of the reactions (4 out of 31) is fully identifiable from the data of [START_REF] Ishii | Multiple high-throughput analyses monitor the response of E. coli to perturbations[END_REF]. In addition, only 37 out of the total of 100 model parameters are individually identifiable. We note that these results are different from those reported in [START_REF] Berthoumieux | Identification of linlog models of metabolic networks from incomplete high-throughput datasets[END_REF], due to the fact that here we take into account the metabolite concentration measurement error to decide whether a parameter associated with a principal component is negligible. The low numbers of identifiable reactions and parameters agree with those obtained with power-law models on other state-of-the-art datasets [START_REF] Srinath | Parameter identifiability of power-law biochemical system models[END_REF]. This further demonstrates the importance of a preliminary identifiability analysis when estimating parameters in metabolic network models.

The rank analysis carried out to determine the identifiability of a reaction also shows how the model can be reduced if the data does not allow the parameters of the model to be unambiguously determined. This reduction step yields minimal models in the form of Eq. ( 24), that have the advantage of being adapted to the informativeness of the dataset. A disadvantage of this approach, however, is that the parameters of the reduced model may be difficult to interpret from a biological point of view, as they do not generally determine the original parameters in a unique way, but rather define a linear subspace of the parameter space. Nevertheless, in some cases it may still be possible to identify some parameters of the original model (Proposition 3 and Remark 3). This criterion was shown to be useful in practice, as it allowed to uniquely identify 30 out of 93 parameters in the nonidentifiable reactions of the E.coli metabolic network model (Table 3).

If the parameters of the original, non-reduced model cannot be uniquely determined from the data, then additional experiments may be necessary. Generally speaking, experimental conditions that explore the range of possible behaviors of the network as much as possible improve identifiability. Given that experiments are usually carried out at steady-state, especially for metabolic flux measurements, the available datasets have a sampling bias that may complicate parameter estimation. In particular, metabolic systems almost invariably contain highly evolved regulatory loops that may homeostatically buffer the concentrations of some metabolic pools [START_REF] Bennett | Absolute metabolite concentrations and implied enzyme active site occupancy in Escherichia coli[END_REF][START_REF] Ishii | Multiple high-throughput analyses monitor the response of E. coli to perturbations[END_REF]. As a consequence, a range of different growth conditions and genetic backgrounds may lead to little variation in steadystate concentrations. Moreover, some metabolic pools may change in a correlated fashion against experimental perturbations, which also hinders identifiability. Finally reactions that normally operate close to equilibrium are intrinsically nonidentifiable, unless the experimentalist achieves to measure their rates further from equilibrium. The growing availability of time-series data (e.g., Voit et al. (2006b); [START_REF] Hardiman | Topology of the global regulatory network of carbon limitation in Escherichia coli[END_REF]), although more demanding from an experimental point of view, promises to relieve this problem. From this we see that p = a B T T , with a defined as above, is different from p * but is such that J p (e, u), x p (e, u) = J p * (e, u), x p * (e, u) for all (e, u) ∈ D, which contradicts the hypothesis.

Figure 1 :

 1 Figure 1: (a) Structure of a small metabolic network with negative feedback. (b) Equations of the linlog model of the network

Example 3 .

 3 Consider the negative feedback network structure shown in Fig. 1. Let us define the

Figure 2 :

 2 Figure 2: Left: Scatter plot of steady-state metabolite concentrations for 1000 randomly generated enzyme concentrations, for two different parametrizations of the model of Fig. 1 (see the text of Example 3 for more details). Red: Simulation for a 1 = 0.0297 and B 2,1 = -0.0073 (weak feedback); Blue: Simulation for a 1 = 0.0297 and B 2,1 = -7.2961 (strong feedback). Right: Individual zooms of the two datasets, with consistent coloring.

Figure 3 :

 3 Figure3: Estimates of parameters B 1,1 , B 2,1 (first row of B x ) from metabolite concentrations, enzyme concentrations and flux measurements at steady-state, for a linlog model of the negative feedback network in Fig.1. In each panel, scatter plots are reported for four different experimental scenarios: 20% noise and q = 20 datapoints (blue); 20% noise and q = 100 datapoints (green); 50% noise and q = 20 datapoints (red); 50% noise level and q = 100 datapoints (magenta). 95%-confidence ellipses are drawn for each scenario (dashed lines). Reference parameter values are indicated by the intersection of horizontal and vertical dotted lines. Refer to the text of Example 4 for additional details. True parameter values are given in Example 3. (a) Case a 1 = 0.0297 and B 2,1 = -0.0073. (b) Case a 1 = 0.0297 and B 2,1 = -7.2961.

  Figure3: Estimates of parameters B 1,1 , B 2,1 (first row of B x ) from metabolite concentrations, enzyme concentrations and flux measurements at steady-state, for a linlog model of the negative feedback network in Fig.1. In each panel, scatter plots are reported for four different experimental scenarios: 20% noise and q = 20 datapoints (blue); 20% noise and q = 100 datapoints (green); 50% noise and q = 20 datapoints (red); 50% noise level and q = 100 datapoints (magenta). 95%-confidence ellipses are drawn for each scenario (dashed lines). Reference parameter values are indicated by the intersection of horizontal and vertical dotted lines. Refer to the text of Example 4 for additional details. True parameter values are given in Example 3. (a) Case a 1 = 0.0297 and B 2,1 = -0.0073. (b) Case a 1 = 0.0297 and B 2,1 = -7.2961.

Figure 4 :

 4 Figure 4: Identifiability analysis for the feedback model in Fig. 1, with a 1 = 0.0297 and B 2,1 = -7.2961. (a) Squared singular values for 100 data matrices Y (see text of Example 5). The blue dots are the estimates s2 ℓ -ν 2 of the squared singular values for all datasets, and the red box covers the area below the cutoff of ν 2 (Def. 3). (b) Normalized cumulative sum of the squared singular values s2 ℓ for 100 data matrices Y . The red box displays the area above the cutoff θ = 0.99 (Def. 4).

  (b). The parameter matrices of the model include 33 nonzero entries and are given by a = ˆ-31.4 4.41 0.13 0.31 0.31 0.13 -0.42 0.97 ˜T ,

Figure 5 :

 5 Figure 5: Overall procedure for identifiability analysis and model reduction.

Figure 6 :

 6 Figure 6: (a) A branched metabolic pathway with feedback from (Visser and Heijnen, 2003). All reactions are chemically reversible, the arrows represent the positive flux directions. Dashed lines represent allosteric interactions. (b) Model variables for internal and external metabolites.

Figure 7 :

 7 Figure 7: Singular value estimates ŝ2ℓ = s2 ℓ -ν 2 computed from 100 noisy datasets Ỹ of size q = 30 (blue dots) of the model of Fig.6. The dashed bars correspond to the singular values s 2 ℓ /q from the noiseless dataset Y of size q = 30. In red is the area below ν 2 ; all singular values in this area are considered negligible.

Figure 8 :

 8 Figure 8: Scheme of Escherichia coli central carbon metabolism[START_REF] Berthoumieux | Identification of linlog models of metabolic networks from incomplete high-throughput datasets[END_REF]. The map shows metabolites (bold fonts) and genes (italic). Abbreviations of metabolites are glucose (Glc), glucose 6phosphate (G6P), fructose 6-phosphate (F6P), fructose 1-6-biphosphate (FBP), dihydroxyacetone phosphate (DHAP), glyceraldehyde 3-phosphate (G3P), 3-phosphoglycerate (3PG), phosphoenolpyruvate (PEP), pyruvate (Pyr), 6-phosphogluconate (6PG), 2-keto-3-deoxy-6-phospho-gluconate (2KDPG), ribulose 5-phosphate (Ru5P), ribose 5-phosphate (R5P), xylulose 5-phosphate (X5P), sedoheptulose 7-phosphate (S7P), erythrose 4-phosphate (E4P), oxaloacetate (OAA), citrate (Cit), isocitrate (IsoCit), 2-keto-glutarate (2KG), succinate-CoA (Suc-coA), succinate (Suc), fumarate (Fum), malate (Mal), glyoxylate (Glyox), acetyl-CoA (Ac-coA), acetylphosphate (Acp) and acetate (Ace). Cofactors impacting the reactions are not shown: adenosine triphosphate (ATP), adenosine diphosphate (ADP), nicotinamide adenine dinucleotide phosphate (NADP) and its reduced form (NADPH), nicotinamide adenine dinucleotide (NAD) and its reduced form (NADH) and flavin adenine dinucleotide (FAD). The gene names are separated by a comma in the case of isoenzymes, by a colon for enzyme complexes, and by a semicolon when the enzymes catalyze reactions that have been lumped together in the model.

(

  x * • • • ln x q * -ln x * ln u 1 -ln u • • • ln u q -ln u T , is unique in the parameters B * i = [B x * B u * ] T i . Proof. (If) Assume that, for a given D ⊆ E × U, the solution of W * i = Y * B * i is unique. We need to prove that (J p ) i , x p | D = (J p * ) i , x p * | D implies p i = p * i .For simplicity, here we drop index i from subscripts. Given any two parameters p * = a * B * T T and p = a B T T , for which M p : (e, u) → (J p , x p ) and M p * : (e, u) → (J p * , x p * ), it holds by construction that W = Y B andW * = Y B * . If (J p , x p )| D = (J p * , x p * )| D , then it also holds that Y = Y * and W = W * , therefore we can write W * = Y * B.Because the solution in B of the latter is unique and one solution is B * , it follows that B = B * . To conclude that p = p * , we are left with showing that a = a * . This follows from Only if) Here the hypothesis is that, for a givenD ⊆ E × U, (J p ) i , x p | D = (J p * ) i , x p * | D implies p i = p *i , and we need to show that the solution in B i of W * i = Y * B i is unique. For simplicity, we will again drop i from the subscripts. For the sake of contradiction, assume that W * = Y * B admits distinct solutions. Since B * is a solution, all solutions are of the form B = B * + z, with z in the nontrivial kernel of Y * . For any such z we can write Y * B * = Y * (B * + z), i.e., ∀(e, u) ∈ D, (ln x * -ln x * ) T (ln u * -ln u * ) T B * = (ln x * -ln x * ) T (ln u * -ln u * ) T (B * + z). (B.1) Let p * = a * B * T T . For any (e, u) ∈ D, J * = J p * (e, u) and x * = x p * (e, u) are given by the solution of 0 = NJ p * , J p * = diag(e) • (a * + ln x T * ln u T * B * ) (B.2) (which is unique by virtue of Assumption 1). Using (B.1), term ln x T * ln u T * B * can be rewritten as ln x T * ln u T * (B * + z)ln x T * ln u T * z. Replacing this into (B.2) yields 0 = NJ * , J * = diag(e) • (a *ln x T *

Table 1 :

 1 Average effective rank computed for each reaction and with different definitions of r over 100 datasets of the model of Fig.6. The criterion of Def. 4 was computed choosing θ = 0.99.

	Reaction Number of Average effective rank
	number parameters Def. 3	Def. 4
	R1	2	1	1.98
	R2	4	2	3.51
	R3	3	1.05	3
	R4	4	2.11	4
	R5	1	0.21	1
	R6	3	1.3	2.99
	R7	4	2.02	3.96
	R8	1	0.03	1

Table 2 :

 2 is therefore the large number of nonidentifiable reactions and parameters. Average effective rank computed for the reactions in the linlog model of E. coli central carbon metabolism, using the data of[START_REF] Ishii | Multiple high-throughput analyses monitor the response of E. coli to perturbations[END_REF]. SVD has been applied to Y C(i) for each reaction and singular values were discarded based on Def. 3. Identifiable reactions are shown in green. Reaction 27, labeled µ, is a phenomenological reaction for biomass production.

			Average	Full			Average	Full
	Reaction	Enzyme	effective rank dimension	Reaction	Enzyme	effective rank dimension
	1	PtsG	3	4	17	GltA,PrpC	2.97	4
	2	Pgi	1	2	18	AcnA,AcnB	1	2
	3	PfkA,PfkB	2.85	4	19	IcdA	1	3
	4	FbaA,FbaB	2	2	20	SucA:SucB:LpdA;SucC:SucD	1	3
	5	TpiA	2	2	21	SdhA:SdhB:SdhC:SdhD	1	3
	6	GapA;Pgk	2.99	4	22	FumA,FumB,FumC	1	2
	7	GpmA,GpmB;Eno	1	2	23	Mdh	2.97	4
	8	PykA,PykF	2	4	24	Ppc;PckA	3	5
	9	AceE:AceF:LpdA	1.99	3	25	MaeB,SfcA	2	5
	10	Zwf;Pgl	1.98	3	26	AceA;AceB	1	3
	11	Gnd	2	3	27	µ	4.94	11
	12	Rpe	1	2	28	Edd;Eda	1	2
	13	RpiA,RpiB	1.99	3	29	Pta;AckA,AckB	3	6
	14	TktA	1.82	2	30	LdhA	1	2
	15	TalA,TalB	1	2	31	AdhE	1	1
	16	TktB	1.01	2

Table 3 :

 3 Parameter matrix [B x B u

Strictly speaking, a better version of Def.

would require that condition (19) holds for all p * in P . This would automatically rule out trivial definitions of BC(i),i such as BC(i),i B * C(i),i (which makes the reaction identifiable for any α and B i but cannot be built without the knowledge of B * C(i),i itself). Unfortunately, this is not a good choice in general, in that estimation uncertainty may severely depend on p * itself, as we shall see later on in Example 4. For simplicity, here we stick to Def. 2 with the understanding that any such triviality is avoided.

This holds as an equality in the sense of expectation, and can also be motivated by asymptotic arguments as q → +∞.

Acknowledgments

This work was supported by the Agence Nationale de la Recherche under project MetaGenoReg (ANR-06-BYOS-0003).

Appendices

Appendix A. Notation and terminology R, R >0 , Z and N denote the sets of real, positive real, integer and positive natural numbers, respectively. For an index n ∈ N, R n and R n >0 denote the n dimensional versions of R and R >0 . I denotes an identity matrix of dimension fixed by the context. Let M be any matrix. For two indices i and j and a vector of indices C compatible with the dimensions of M, M i denotes the ith column of M, M C denotes the submatrix of M formed by the columns of M with indices C, M ji denotes the element of M in row j and column i, and M j,C denotes the row vector formed by the elements of M in row j and columns indexed by C. [M] C,C denotes the minor formed by the rows and columns indexed by C. When convenient, notation M i:i ′ with i ′ ≥ i is used instead of M C with C = i, i + 1, . . . , i ′ . For vectors, a subscript i refers to the ith element of the vector. For a square matrix Σ, Σ > 0 (resp. Σ ≥ 0) means that Σ is positive definite (resp. semidefinite). For a vector µ of suitable dimension, ε ∼ N (µ, Σ) means that ε is a Gaussian random vector with mean µ and covariance matrix Σ. For a scalar quantity observed in Gaussian noise with standard deviation ν, a noise level of N%, with N ∈ R ≥0 , will stand for a value of ν such that ν is equal to N% of the value of that quantity (such that ≃ 99% of the noise outcomes fall within ±3 • N% of the observed quantity). For two vectors v and e of equal size, both v/e and v e indicate the vector obtained by elementwise division. Given a vector sequence v 1 , . . . , v q , v is the mean (1/q) q k=1 v k . For vectors and sets, | • | denotes vector dimension and set cardinality, respectively. For a function

Appendix B. Mathematical proofs Proposition 1. A reaction i of M p is structurally identifiable at p * if and only if there exists D = {(e 1 , u 1 ) , • • • , (e q , u q )} ⊆ E × U such that the solution of the equation

Corollary 1. A reaction i of M p is structurally identifiable at p * if and only if there exists D = {(e 1 , u 1 ) , • • • , (e q , u q )} ⊆ E × U such that Y * C(i) is full column-rank. Proof. From Proposition 1, we know that identifiability is equivalent to the uniqueness of the solution in B i of

Proposition 2. If a reaction i of M p is structurally identifiable at p * in the sense of Def. 1 then, for every α ∈ (0, 1), it is practically identifiable in the sense of Def. 2 with confidence level at least 1 -α for any uncertainty set

Proof. The definition of M p ensures that

Using the fact that Y * = Y and W * = W , given a noisy dataset Wi = W i + ε i and the errorless dataset Y , the regression problem becomes

From Section 3.1, identifiability of M p at p * for the given input set D is equivalent to

This enables us to define the minimum variance estimator of B * C(i),i , BC(i),i = Y † W i . From the linearity of the estimator in the Gaussian noise ε i , after simple calculations of first and second-order moments, one gets BC(i),i ∼ N B * C(i),i , Σ (also compare (Ljung, 1999, Appendix II)). Thus, from the definition of B i , P p * [ BC(i),i -B * C(i),i ∈ B i ] ≥ P p * [ BC(i),i -B * C(i),i ∈ E b Σ (α)] = 1 -α.


Caroline Bernard-Michel

Laurent Gardes

Stéphane Girard
email: stephane.girard@inrialpes.fr

Gilles Molinié
email: gilles.molinie@hmg.inpg.fr

Spatial analysis of extreme rainfalls in the Cévennes-Vivarais region

published or not. The documents may come   L'archive ouverte pluridisciplinaire