|Home | About | Journals | Submit | Contact Us | Français|
The organization, regulation and dynamical responses of biological systems are in many cases too complex to allow intuitive predictions and require the support of mathematical modeling for quantitative assessments and a reliable understanding of system functioning. All steps of constructing mathematical models for biological systems are challenging, but arguably the most difficult task among them is the estimation of model parameters and the identification of the structure and regulation of the underlying biological networks. Recent advancements in modern high-throughput techniques have been allowing the generation of time series data that characterize the dynamics of genomic, proteomic, metabolic, and physiological responses and enable us, at least in principle, to tackle estimation and identification tasks using “top-down” or “inverse” approaches. While the rewards of a successful inverse estimation or identification are great, the process of extracting structural and regulatory information is technically difficult. The challenges can generally be categorized into four areas, namely, issues related to the data, the model, the mathematical structure of the system, and the optimization and support algorithms.
Many recent articles have addressed inverse problems within the modeling framework of Biochemical Systems Theory (BST). BST was chosen for these tasks because of its unique structural flexibility and the fact that the structure and regulation of a biological system are mapped essentially one-to-one onto the parameters of the describing model. The proposed methods mainly focused on various optimization algorithms, but also on support techniques, including methods for circumventing the time consuming numerical integration of systems of differential equations, smoothing overly noisy data, estimating slopes of time series, reducing the complexity of the inference task, and constraining the parameter search space. Other methods targeted issues of data preprocessing, detection and amelioration of model redundancy, and model-free or model-based structure identification.
The total number of proposed methods and their applications has by now exceeded one hundred, which makes it difficult for the newcomer, as well as the expert, to gain a comprehensive overview of available algorithmic options and limitations. To facilitate the entry into the field of inverse modeling within BST and related modeling areas, the article presented here reviews the field and proposes an operational “work-flow” that guides the user through the estimation process, identifies possibly problematic steps, and suggests corresponding solutions based on the specific characteristics of the various available algorithms. The article concludes with a discussion of the present state of the art and with a description of open questions.
The task of biomathematical modeling comprises the conversion of a biological system into a simplified analogue that is easier to analyze, interrogate, predict, extrapolate, manipulate, and optimize than the biological system itself. The typical approach to mathematical model construction consists of nine phases: (1) data selection; (2) collection of information on network structure and regulation; (3) specification of assumptions and simplifications; (4) selection of a mathematical modeling framework; (5) estimation of parameter values; (6) model diagnostics; (7) model validation; (8) model refinements; and(9) model application (see Figure 1).
The first phase requires the identification and selection of data that are available and suitable to support the purposes of the modeling effort. The second phase is dedicated to collecting information regarding the structure and regulation of the system from the literature and, where feasible, from de novo experiments. This phase is confounded by the fact that the true topology of the biological network is often not fully understood and that regulatory details are in many cases incomplete, obscure, or entirely missing. Under these circumstances the task of this and the next phases includes inferences about the network structure and its regulation from biological data. After collecting all relevant information regarding the biological system, the third phase is dedicated to combining this information with additional, acceptable assumptions and simplifications that aim to fill the gaps in the available information. During this phase one also decides which components and interactions of the system should be included in the model. The results are usually visualized as diagrams with nodes denoting the components and arrows representing interactions between them. The fourth phase includes the choice of a suitable mathematical modeling framework and the formulation of symbolic equations. The process usually starts with converting the “wiring diagram” or the “network topology” obtained from the third phase into model equations. In many biological systems analyses, these form a set of ordinary differential equations that represent the dynamic changes in system variables and are governed by fluxes between them. The particular symbolic format of the fluxes depends on the chosen mathematical modeling framework and defines what questions can be asked of the model and what types of methods will be applicable. After the symbolic equations are formulated, the task of the fifth phase is to determine appropriate numerical parameter values that convert the symbolic model into a numerical model and makes the latter consistent with experimental observations. Once this parameterized, initial model is obtained, the sixth phase is dedicated to diagnostics of the model, including analyses associated with steady states, sensitivities, gains, and stability, as well as dynamic features, such as bifurcations and oscillations. In the seventh phase the validity of the model is tested, either with experimental data that had not been used for model construction in a process called cross-validation, or with information from a different biological level. For instance, a metabolic model could be tested against physiological or clinical observations. As presented so far, the modeling process appears to be quite straightforward. However, in most cases it is not linear but cyclic, requiring the return to earlier phases. Addressing the iterative nature of modeling, an eighth phase of model refinement is almost always necessary. Finally, once the model is deemed reliable and appropriate for the stated purposes, phase nine allows the modeler to reap the fruit of the laborious model design. It is now possible to make predictions, generate new, testable hypotheses, suggest the design of novel biological experiments, or manipulate and optimize the model, for instance, toward increases in yield of some organic compound in metabolic engineering, or toward the development of drug treatments in disease.
Among the nine modeling phases, the most challenging task is usually the estimation of parameter values. This estimation is not an isolated task but closely related to other phases in the modeling process. For instance, the size and complexity of the hypothesized model in the second phase have a direct impact on the difficulty of parameter estimation and also affect later analyses as well as the interpretation of results. Most importantly, the choice of the modeling framework naturally influences the degree of complexity, feasibility, and practicability of the parameter estimation task. As a simple example, an explicit linear model permits the use of linear regression methods, which are very well worked out. As soon as the model becomes nonlinear, many of these methods become inapplicable.
Because of the importance of issues related to model selection and to implications for parameter estimation, we will use Section 2 to review the rationale and special demands on mathematical models for biological pathways and to introduce some of the most prevalent and representative modeling frameworks. Generally, model selection and parameter estimation depend on knowledge about the system, the purpose of the modeling effort, and available data. If much is known about the details of the mechanisms governing the biological system, mechanistic formulations, maybe based on principles of physics, are often a good choice. By contrast, if details are lacking, it has been shown that canonical models, such as Lotka-Volterra models and models within Biochemical Systems Theory (BST) are very advantageous for the purposes of mathematical modeling in biological systems. Pertinent details of canonical models will be reviewed in Section 2.3.3.
The development of parameter estimation methods is driven by the availability of experimental data. Different types of data require distinctly different classes of methods. Conversely, the various estimation methods require different types of data. As a pertinent example, data for metabolic pathway models have traditionally characterized the kinetic properties of individual enzymes catalyzing particular steps within a metabolic pathway of interest. The generation of this information occurred hand in hand with the concepts and terminology of enzyme kinetics and the data were measured specifically to parameterize models in familiar formats such as Michaelis-Menten or Hill rate laws. The strategy of using these types of “local” descriptions of model components (one enzymatic process at a time) and merging them into a much more comprehensive, dynamical pathway model is referred to as “bottom-up” or “forward” modeling and will be discussed in Section 3.
Steady-state data may also be used in parameter estimation. These data characterize a metabolic system under conditions where all concentrations have reached constant levels and all fluxes are in balance. Specifically, this type of steady-state analysis is either based on stoichiometric analysis or on experiments that measure the responses of a biological system after a small perturbation around the steady state. Some of these methods will also be discussed in Section 3.
Recent innovations in biological technology enable us to tackle the parameter estimation task in an entirely different, more comprehensive manner, using a “top-down” or “inverse” approach. The biological tools for these purposes are geared toward generating time series data or “global” data of genes, proteins, or metabolites, sometimes under different sets of conditions, such as initial concentrations, or upon various gene knock-outs or the inhibition of specific enzymatic steps. Inverse methods are very appealing, because they provide measurements on cellular or organismal systems in a larger context. In particular, if the data are generated in vivo and with only minimal disturbance of the biological system, the insights gained are considered to be as close to reality as is presently possible and potentially much more valuable than data obtained from experiments performed in an artificial in vitro setting. Details and features of the traditional and the newly developed techniques with respect to parameter estimation will also be addressed in Section 3.
The potentially high rewards of the inverse modeling approach have motivated scientists from different backgrounds and from all over the world to dedicate considerable effort to the many challenges that must be overcome to make the approach useful. For BST models alone, about one hundred articles and numerous proceedings and book chapters describing computational methods for inverse tasks have appeared within the past decade. We address the specific challenges and requirements of inverse modeling in Section 3.4, along with different types of support algorithms. Many of the pertinent methods target the main problem of optimizing parameter values against the observed time series data. Others suggest strategies for circumventing the costly integration of differential equations, smoothing noisy data, estimating slopes, constraining the parameter search space, or reducing the complexity of the inference task. These auxiliary methods and algorithms will be reviewed in Section 4. The primary focus will be on methods applicable to models within BST, but we will also discuss related issues that are of interest to other modeling approaches.
The earlier discussion of the second modeling phase (see Figure 1) mentioned that the true topology of a biological system is sometimes not fully understood or obscure. In such a case, parameter estimation is much more complicated, because it is a priori not clear how to formulate an ill-characterized biological system mathematically. As we will discuss, the use of concept maps  and of canonical models is of great help in this situation, because it converts the task of identifying the uncertain structure and regulation of the biological system into a parameter estimation task. Generally, structure identification tasks are much more difficult than parameter estimation, which is already very challenging. Canonical models render both tasks reachable. In particular, one should note that there is no clear boundary between parameter estimation and structure identification if canonical models are used. Section 5 will introduce some of the most relevant structure identification methods, namely the determination of the Jacobian matrix, direct observations, correlation-based approaches, ‘simple-to-general’ and ‘general-to-specific’ modeling, and specially tailored time series data analysis within the framework of BST.
Among all parameter estimation or structure identification methods proposed so far, no single method has risen to the top and can be declared the clear winner in terms of efficiency, robustness and reliability for the majority of realistic case studies. However, it seems that a combination of methods may be useful in a large number of cases. To make inverse modeling more effective and translucent, we propose in Section 6 an operational “work-flow” that guides the user through the various steps of the estimation process, identifies possible problem areas, and suggests promising solutions based on the specific characteristics of the various available algorithms. An interesting consequence, and actually an advantage, of the combined approach is the general result that the optimized solution often consists of multiple, distinctly different parameter sets that are all consistent with the data and that can lead to novel hypotheses for further theoretical and experimental investigation.
Biological systems consist of many organizational layers including genetic, transcriptomic, proteomic, and metabolic levels, as well as phenomena of cell physiology, cell communication, tissue and organ function, populations, and ecology. In this review, we focus primarily on model construction at the genomic and metabolic levels, although many of the computational methods are independent of a particular application. The main reason is that the genomic and metabolic levels are currently best supported by available quantitative data. Metabolic time profiles are particularly well suited because of the stoichiometric property of metabolic pathways, which creates natural constraints on possible parameter combinations, and because its main drivers, namely metabolic concentrations and fluxes, can be measured, at least in principle. In contrast, no material flow is involved in gene-gene, gene-protein, or protein-protein interactions, and the measurable effects are seen in their consequences rather than in characteristics of the processes themselves. Therefore, gene and protein networks must often be studied with coarser methods than metabolic systems, such as graph theory and Boolean or Bayesian methods, which are applied under the simplifying assumption of binary on/off states. Nevertheless, recent methodological developments have enabled the generation of some dynamic profiles of gene networks, and these have been used for the quantitative identification of gene regulatory networks, primarily by several Japanese groups (see Section 2.3.4). Dynamic data on protein levels are still very rare, and quantitative time series responses are very difficult to obtain. Commensurate with the availability of data, we will primarily focus on the construction and analysis of metabolic pathway models but also discuss issues related to gene interaction networks.
Because the material in this review discusses numerous complementary aspects of parameter estimation and structure identification, it seems useful to summarize the structure of the review in the form of a roadmap, which is given in Table 1.
Mathematical modeling and control theory have a long history in physics and engineering. However, the demands and specific requirements in modeling biological systems are quite different and necessitate the adaptation and extension of former methods, as well as the development of novel, additional tools that are optimally suited for modeling biological phenomena. The peculiarities of biological system modeling can be generally described by five aspects (e.g., ). First, the biological processes and interactions are highly nonlinear and complex. Thus, a mathematical structure is needed that can capture nonlinearities and does not a priori exclude relevant biological phenomena, such as stable oscillations. Second, dynamic responses of biological systems are particularly interesting. Therefore, a suitable mathematical model will have to be time dependent, which almost always requires formulation as a set of differential equations. Third, real biological systems are usually composed of different levels of components and interactions with relatively large numbers, which require the ability of a mathematical framework to scale up to increasingly larger biological models. Fourth, biological systems may have stochastic features when there are only few molecules involved. Under this condition, the fundamental laws of kinetics and thermodynamics are no longer directly applicable and the biological behavior becomes difficult to predict. Thus, in addition to grasping a deterministic phenomenon, the mathematical model should also be able to capture stochastic behaviors when these dominate the process. And fifth, biological reactions rarely happen in a homogeneous environment but are often restricted to surfaces, channels, organelles or compartments. This feature is sometimes important, and therefore the ability of handling spatial processes is necessary for a comprehensive mathematical analysis. It might be added to this list that many biological phenomena evolve over distinctly different time scales and are controlled from different levels of organization. Furthermore, continuous trends are often affected by discrete events, such as the sudden activation of a gene, and by events that lie in the past and cause a delayed effect . At this point, no theoretical or computational frameworks exist to deal with all these aspects. For specific cases, ad hoc models may be developed in general purpose software like MatLab or Mathematica, or one might use hybrid methods [4,5], or agent based approaches , which however are computationally expensive.
While stochastic, spatial, and time scale effects are known to exist in biological systems, it is often possible to use approximations that greatly simplify model design and analysis. The validity of such approximations needs to be determined on a case-by-case basis. If the approximations are indeed valid, they often lead to simplified system representations based on ordinary differential equations. The generic format for such a representation may be written as
where Xi denotes the concentration or amount of a variable or variable pool and n is the total number of time dependent variables in the system. The functions and represent reactions or fluxes entering or leaving the quantity Xi. This general framework can be recast in numerous alternative ways and only becomes meaningful and specific when the functions and are mathematically defined and parameterized. In the following sections we will briefly review some particularly relevant implementations in the context of metabolic pathway analyses and dynamic models of gene regulatory systems.
Mathematical models describing metabolic pathways can be constructed with a focus either on stoichiometry or kinetics. The stoichiometric property of a pathway is typically considered time invariant, while kinetic aspects are used to capture the dynamics of a system and are driven by the state of the system and may change rather quickly. The stoichiometry of a pathway system determines the wiring diagram of the pathway and describes which fluxes enter or leave which pool and enforces that no mass is gained or lost in the process. The translation of this topological wiring diagram into a matrix equation is straightforward. The resulting stoichiometric matrix contains positive, negative or zero elements that represent which metabolites are converted into which other metabolites. The sign represents the direction of material flow and indicates whether the reaction increases or decreases the concentration of a given metabolite pool. If a metabolite and a reaction are unrelated, the corresponding element is zero. The value of each element in the matrix indicates the stoichiometric relationship and must be an integer. For instance, if one substrate molecule breaks down into two product molecules, the gain in product is coded as +2.
The stoichiometric matrix N is the core of stoichiometric models that show how metabolite concentrations change over time. A differential equation is formulated as
The main application of stoichiometric models is the determination of flux rates. Estimation methods for this purpose depend on the type of available experimental data. In most analyses, stoichiometric models are studied for metabolic systems in a steady state where, for all pools, the material flow into the pool equals the material flow out of the pool and all flux rates are constant. Under this assumption, the left-hand sides of the equations in Eq. (2) become zero and the system of differential equations turns into a set of linear algebraic equations. If the stoichiometric matrix is full rank, it is straightforward to calculate the fluxes. However, usually there are more unknown fluxes than equations, so that the system of linear equations is underdetermined. Stoichiometric analysis is most often applied to microbial systems, where it is assumed that the microbes tend to optimize their growth rate. This assumption can be formulated as an objective to maximize the availability of nutrients needed for growth. Representing this objective with a linear function transforms the underdetermined stoichiometric system into a linear programming task, which is easily solved even for large systems. Mass conservation and stoichiometry are distinctive properties of metabolic pathways and not applicable to the dynamics of gene regulatory or proteomic networks.
Flux balance analysis (FBA) inherits the properties of the stoichiometric approach but additionally imposes physical and chemical constraints to find the feasible or optimal distribution of fluxes . For instance, it is possible to account for thermodynamic limitations. The background of FBA is reviewed in Palsson  and representative developments of variations are summarized in Kauffman et al. . The modeling process in FBA consists of system identification, mass balancing, defining measurable fluxes, and optimization. System identification consists of setting up the stoichiometric equation and relevant constraints. Mass balance is the application of stoichiometry and conservation of mass. For instance, the total number of moles of carbon in the system is conserved during the time of reaction. By accounting for all material flows entering and leaving each metabolite pool in the pathway, one can determine the material distribution and also identify flows that might have been unknown or difficult to measure in experiments. As in stoichiometric analysis, the optimization step assumes an objective such as maximization of growth rate or product yield that permits solutions through linear programming. The main advantage of both the stoichiometric model and FBA is the linearity of their matrix representation at the steady state, which tremendously simplifies the analysis, since there are numerous well-established methods of linear algebra that are directly applicable. Several examples have demonstrated that FBA is capable of assessing the theoretical capabilities and operative modes of metabolic systems in the absence of kinetic information (cf. [9,13–16]).
Stoichiometric models are sometimes studied under a pseudo-steady-state (PSS) assumption in cases where the concentrations of metabolites rapidly adjust to new levels [17–19]. This PSS approximation was shown to be valid for most intracellular metabolites . Under this assumption, it is reasonable to neglect the instantaneous changes of metabolites and set the rate of change to zero.
When complete time courses of metabolites are available, the flux distribution at each time point can be determined with  or without  the PSS assumption. Distinctly different from the standard application of stoichiometric analysis, where only steady-state data are used, the dynamic changes in metabolite concentrations in the latter case are not necessary zero and can be used to gain incomparably deeper insight into the pathway at hand .
Mahadevan and coworkers  extended traditional FBA to account for dynamics and presented two different formulations: the dynamic optimization approach (DOA) and the static optimization approach (SOA). DOA involves optimization over the entire time period of interest to obtain time profiles of fluxes and metabolite levels. SOA involves dividing the batch time into several time intervals and solving the instantaneous optimization problem at the beginning of each time interval. By testing the methods in the analysis of diauxic growth in Escherichia coli, the authors concluded that SOA was computationally simpler to implement provided all of the constraints were linear, whereas DOA was more flexible and suitable for the incorporation of experimental data.
Utilization of the stoichiometric property together with dynamic changes in metabolites is a valuable option for studying flux distributions in metabolic pathways. However, the main advantage of the standard stoichiometric approach, namely linearity at the steady state, is at the same time its most severe limitation. Specifically, this approach focuses almost exclusively on the connectivity structure of the system and the flux distribution at steady state, but does not account for kinetic features, which often are necessarily nonlinear. Therefore, the predictive power of linear stoichiometric models, while successful in many cases, is also limited because regulatory signals and other nonlinear dynamic interactions cannot be included in the model without destroying its linear structure. As a compromise, Palsson and others (e.g., ) introduced binary-valued regulatory matrices that are multiplied to the stoichiometric matrix and determine whether a given flux is activated or not. However, a full account of nonlinear regulatory features requires the formulation of a pathway system as a kinetic model.
When detailed information is available about the kinetics of the metabolic reactions in the pathway, it is possible to describe its dynamics by incorporating kinetic features in the flux representations v of the general stoichiometric models in Eq. (2) . The crucial step toward combining the stoichiometric property with kinetic features is to search for appropriate functional forms to represent the flux quantities and in Eq. (1), which then translate into representations of the vector v. Approaches for this purpose can be categorized into three categories (see Figure 2): (1) mechanistically based functions (e.g., law of mass action, Michaelis-Menten rate law); (2) ad hoc approaches; and (3) different types of canonical models (e.g., BST and lin-log representations). Details of these representations are reviewed in the following sections.
Models based on the law of mass action are typically used to describe reaction networks consisting of elementary reactions. The rate of a given elementary chemical reaction is proportional to the product of concentrations of all variables reacting in the elementary process, including their moieties, and is generally formulated as the basis function
where k is the rate constant, which is always positive, and gi are kinetic orders which are nonnegative integer numbers that reflect the numbers of molecules involved in the reaction. The main advantage of models based on the law of mass action is that they can be determined directly from the elemental reactions and their stoichiometry. The drawback is that most biochemical reactions are not elemental but catalyzed by enzymes, and therefore composites of several elementary reactions . It is inconvenient and indeed infeasible to carry all these reactions along and much easier instead to develop composite rate functions. Furthermore, the underlying mechanisms of an enzyme catalyzed reaction are not always understood in sufficient detail or they are experimentally inaccessible. Therefore, mass action equations are difficult to set up and parameterize for complex pathway systems.
The Michaelis-Menten rate law (MMRL)  and its variations are among the most commonly used representations for kinetic modeling in metabolic pathway analysis. MMRL is based on the concept that a substrate and an enzyme form a transient complex which either dissolves to return the two or leads to the formation of a product and the release of the enzyme. The modeling of enzyme reactions in this approach is simplified considerably under the quasi-steady-state assumption, which states that the intermediate complex does not change appreciably over time. Even though Michaelis-Menten rate laws are often straightforward to set up, complete descriptions of more complex enzyme mechanisms may become massive if several substrates or reactions are involved . As the result, mathematical analyses become very complex and the parameter estimation requires an undue amount of experimental data [27,28]. In addition to technical issues, the model results are difficult to interpret in terms of the underlying biological system [28,29]. The estimation of parameter values in pure MMRL is easy and may even be accomplished with methods of linear regression . However, this simplicity vanishes for larger systems of MMRLs and their generalizations.
When the detailed mechanisms of a biological process are unknown or unclear, an alternative approach is to adapt a mathematical model from a different context, quasi as a black box model, to describe the biological phenomenon of interest. As a typical example, Voit and Sands  reviewed a collection of mathematical models for the uptake of nutrients by tree roots, most of which had no foundation in plant physiology but were expected to fit observed trends sufficiently well. Although such an ad hoc black-box approach might provide models that fit the observed data, it is highly arbitrary and poses several intrinsic problems. For instance, one must question the validity of such a model for other experimental data or extrapolations to un-tested experimental conditions and anticipate problems with interpretations of results, comparisons of results with those from alternative models, and extensions to more refined models. Furthermore, the estimation of model parameters may become problematic, because they do not necessarily correspond to measurable quantities. In many cases, a better option than an ad hoc model is a canonical model.
The predictive ability of stoichiometric models is limited because nonlinearities due to regulation cannot directly be accounted for. For improved models that permit nonlinearities, kinetic information of the pathway is needed. A good compromise that is capable of capturing nonlinear dynamics while keeping the mathematics relatively simple is the use of a “canonical” nonlinear model whose structure is fixed and whose individuality comes from its parameter values. In addition to their homogeneous structures, canonical models are more or less size independent and facilitate the development of customized techniques and methods for analysis, diagnostics, and parameter estimation.
Arguably the most promising canonical nonlinear models in metabolic modeling are Generalized Mass Action (GMA) and S-system structures within Biochemical Systems Theory (BST) [27,32–36]. These models are constructed by approximating fluxes with products of power-law functions, which are mathematically grounded in the well-established approximation theory of Taylor. In the S-system formalism, each equation has a particularly simple format: The change in system variables is given as one set of influxes minus one set of effluxes (cf. and in Eq. (1)), and each set is collectively written as one product of power-law functions. Thus, the generic S-system formulation reads
where Xi represents a time-dependent variable (metabolite) and n denotes the number of variables in the system. The non-negative multipliers αi and βi are rate constants which quantify the turnover rate of the production or degradation, respectively. The real numbers gij and hij are kinetic orders that reflect the strengths of the effects that the corresponding variables Xj have on a given flux term. A positive value signifies an activating or augmenting effect exerted by Xj, a negative value signifies an inhibitory effect. A kinetic order of zero implies that the corresponding variable Xj does not have any effect on a given flux. In some instances, m independent variables, which are typically constant during each mathematical experiment, may be included. They do not have their own equations but enter the power-law terms just like dependent variables, so that the products run from 1 to n+m.
In the GMA formalism, instead of aggregating all influxes and all effluxes into one term each, all influxes and effluxes are approximated individually with power-law terms such that
where the rate constants γip are non-negative and the kinetic orders fipj may have any real values as in the S-system form. Also as before, independent variables may be included. It should be noted that differences between these two formulations only exist at branch points, whereas all other steps are identical. One should also note that S-systems permit parameter estimation methods , which are not directly applicable to GMA systems, as we will discuss later.
BST models have a number of important advantages which have been discussed in detail elsewhere [27,32,33,35,36]. Four are particularly crucial here. First, these systems are rich enough in structure to capture virtually any nonlinearity including complex oscillations and chaos [38,39]. Second, symbolic BST models can be set up without mechanistic information on the underlying system, but if information is available, it can be used to simplify the symbolic representation. Third, the highly structured format facilitates mathematical and numerical analyses. These analyses include computations associated with steady states, sensitivity, stability, as well as dynamic features. Fourth, BST models are characterized by a one-to-one relationship between parameters and structural features. Thus, if structural features are known, it is explicitly clear where they will appear in the BST models. Conversely, if a parameter has been identified, its interpretation in terms of structural properties is immediate. This feature is especially crucial for structure identification and parameter estimation of metabolic models. The reasons will be discussed in detail in Section 4.4. The power-law models of BST were initially used to model metabolic pathways and gene circuits, but this formalism has also proven very beneficial for modeling other classes of biological systems, including genetic networks, immune networks, multi-level systems, and cell signaling [40–45].
An alternative canonical form is the “lin-log approximation” which was introduced by Hatzimanikatis and Bailey  and expanded by Visser and Heijnen . This form is based on taking the logarithm of each metabolite concentration and enzyme activity in relationship to a corresponding reference value. The lin-log model constitutes an extension of Metabolic Control Analysis (MCA), a theoretical framework for analyzing control and regulation in metabolic networks close to their steady state [29,48,49]. For small variations about a steady state, BST and lin-log models show similar responses. However, they differ for larger variations [50–52].
Another recently proposed canonical form is the Saturable and Cooperative Formalism (SC formalism) , which is based on Taylor approximation in a special transformation space that is defined by logarithms of power-inverses. The SC formalism exhibits improved cooperativity and saturation in comparison to other canonical formalisms. In addition, the SC formalism is expected to be accurate over a wider range around the operating point if the approximated functions are saturated. The main drawback of the SC formalism is its need for a much larger number of parameters, which on one hand bestow SC models with higher flexibility but on the other hand require increased estimation efforts.
For completeness, we should mention Lotka-Volterra (LV) models and their generalizations [54–57]. These models are canonical in the sense of the models discussed before and have found very many applications in ecology. In fact, they are very flexible in structure and, with sufficiently many variables, can model any nonlinearities that can be expressed as ordinary differential equations, just like BST models [38,58,59].
The generic LV format is
where the coefficients aij are real numbers. In other words, the ith equation contains one term that is linear in the variable itself and a product of the variable with another system variable. While this format is very useful in ecological systems analysis, it is less so in metabolic systems. For instance, it does not allow the direct formulation of an un-modulated precursor-product relationship, because the ith equation cannot contain a term of the form kXi−1.
The choice of an S-system, GMA, lin-log, or SC model depends on the available input information and the purposes of the intended model. For instance, GMA systems are natural extensions of stoichiometric models that incorporate kinetic information using a power-law approximation and are closer to biochemical intuition than S-systems. However, the GMA format, as well as the SC format, does not allow the algebraic calculation of steady states, which is important for a number of analyses. The lin-log model shares the intuitive advantage with the GMA format and, like the S-system format, allows algebraic steady-state calculations. However, it cannot represent certain nonlinear behaviors since its structure is essentially linear . Also, the lin-log approximation results in substantial errors for substrate values close to zero [50,52] and may lead to negative rates if the substrate concentration is low, whereas the BST representations become more inaccurate for very high substrate concentrations. As local approximations, all these formats perform similarly as long as the system variables do not deviate too much from some chosen state where they coincide.
The previous sections reviewed representative modeling strategies for metabolic pathway systems and suggested that power-law models are particularly useful in this context. The same power-law formalism is also beneficial for modeling gene regulatory systems, although there are other distinct alternatives (for a review, see e.g. ). Mathematically the simplest representations are linear; they typically capture the structure of the network with a connectivity graph and its associated connectivity matrix (e.g., [61–63]). The graph is modeled as static, that is, as time invariant, and allows classifications of the particular type of the network (e.g., “randomly connected,” “hub and spoke,” or “small world”), as well as analyses of the robustness to different types of perturbations (e.g., ). A different modeling strategy is based on Boolean networks, in which each gene is either turned on or off, and whose discrete dynamics is governed by rules that recursively determine the expression state of a gene during the transition from one time point to the next (e.g., ). Gene networks may also be analyzed with Bayesian methods that assess the probability of a gene being affected by other genes (e.g., [66,67]).
Of particular interest here are models that trace the full dynamics of gene networks with ODE models. While many alternatives could be used as base models, many researchers have resorted to the default of S-system models, because these do not require knowledge of the precise mechanisms that govern the regulation network. The study of gene regulatory networks with S-systems has a long history, which began with the pioneering work of Savageau . Scientists had been observing many types of gene circuitry and Savageau and his collaborators tried to answer questions such as: Is each type of circuitry organized in a manner that is optimized within its context or is the mode of regulation a coincidence? Can gene regulatory circuits be classified into types that are predictable from the organism’s surroundings (e.g. activation versus repression, positive versus negative control)? Is it possible to deduce the functionality of a given circuitry from its particular organization? Is it possible to infer general “design principles” for the modes of regulation from the analysis of alternative gene circuitries? Is it possible to implement these principles and “synthesize” gene circuits de novo in order to achieve prescribed biological responses (cf. [40,69,70])? To answer these questions, Savageau and others dedicated significant effort to studying gene regulatory circuits with techniques of canonical modeling and the method of controlled mathematical comparisons (e.g. [43–45,71–74]). A particularly intriguing result was the conceptual framework of demand theory (e.g. [75–80]), which indeed tied types of gene regulation to environmental demands and thus revealed natural design principles.
One of the models used to analyze gene circuits  has become a benchmark system for inverse methods (see Table 3). In addition, Maki et al.  proposed a network of 30 genes (see also Kimura et al.  and Kutalik et al. ), which is artificial but was inspired by realistic gene networks, for testing the performance of algorithms on larger systems. In a similar vein, Kutalik et al.  created an artificial network of 7 genes.
The collection of input information and the choice of a mathematical model framework result in a symbolic model which is typically in the form of ordinary equations, as discussed before. The next step is to assign numerical values to all parameters in the model. There is no unique recipe for this task of parameter estimation. In fact, the estimation problem is almost always complicated and continues to be the bottleneck of biomathematical modeling.
In this section, we review the classes of parameter estimation methods that are in current use, namely, forward (bottom-up) approaches, estimation from steady-state data, and inverse (top-down) modeling using time-series data. The nature of suitable data for each type of estimation is distinctly different, and so are the methods of analysis. Ideally, rich and diverse data will allow the use of several methods that complement each other or a combined strategy that has much greater potential leading to suitable models than either approach by itself [1,84–86].
Before high-throughput data were available, essentially all metabolic models were developed from “local” kinetic information of biochemical or physiological responses, which had been obtained in the traditional reductionist manner. Typically, biologists around the world worked on characterizing one particular enzyme or transport step at a time. They purified the enzyme, studied its characteristics, determined optimal temperature and pH ranges, and quantified cofactors, modulators, and secondary substrates. Isolated from these laboratory experiments, modelers converted this information into alleged mathematical rate laws. Once enough information had been collected for most steps in the pathway, the modeler attempted to merge this information into an integrative mathematical model. This type of model construction has been benefiting greatly from databases like KEGG [87,88], MetaCyc , and Brenda , which contain large collections of pathway topologies and kinetic parameters retrieved from literature sources. In many cases, the integration into a viable model does not succeed at first and requires revisiting assumptions and parameters, or even the search for additional experimental data.
If done right, this “forward” or “bottom-up” process leads to a model representation of the pathway that exhibits the same features as reality, at least qualitatively, if not quantitatively. Case studies that used this forward estimation approach in BST successfully include: the TCA cycle in Dictyostelium discoideum ; the citric acid cycle [92,93]; fermentation in Saccharomyces cerevisiae [94–96]; purine metabolism [97–99]; the Maillard-glyoxylase network with formation of advanced glycation end products ; the ferredoxin system with information from protein structure for model identification ; and sphingolipid metabolism in Saccharomyces cerevisiae [102–104]. In almost all of these cases, the strategy consisted of setting up a symbolic model, estimating local parameters, studying the integration of all individual rate laws into a comprehensive model, testing the model, and making refinements to some of the model structure and the parameter values in an iterative fashion.
While theoretically straightforward, there are several disadvantages of this approach. The main issue is that a considerable amount of local kinetic information is needed and that this information is often only available from different organisms, different species, or experimental performed under different conditions. As a consequence, more often than not the “integrated result” is not consistent with biological observations. Furthermore, this process of construction and repeated refinement is very labor intensive and requires a combination of biological and computational expertise that is still rare.
By far the most model estimations from steady-state data have been performed in the context of stoichiometric and flux balance analysis, as discussed before. Specifically, fluxes entering and leaving a metabolic system were measured under the assumption of complete metabolite and flux balancing, and internal flux rates were inferred from the assumed stoichiometric models and the maximization of growth rate per linear algebra and linear programming. Examples are plentiful (cf. [9,10]).
Using concepts of stoichiometry in a slightly different fashion, Wiechert and others developed isotopomer and cumomer methods for flux rate estimation using radioactive labeling techniques [105–108]. The technique is based on the fact that specifically labeled atoms in an input molecule, such as 1-13C-glucose, have known metabolic fates. By letting the system resume steady state after the labeled input, the distribution and positioning of label among the metabolites of the entire pathway provides clues of internal flux rates. The original steady-state estimation methods based on positional labeling were later extended to characterize the transients between input and steady state [109,110].
A distinctly different type of parameter estimation from steady-state data uses responses of a biochemical system to (infinitesimally) small perturbations around the steady state. Two types of approaches may be taken. First, parameter values can be obtained, at least in principle, by direct experimental measurements of how a variable affects the fluxes entering and leaving the metabolite pool [33,48,111]. Suppose the flux rate and metabolite concentrations in steady state of one particular biochemical process are known. One can then slightly alter the concentration of a variable systematically while keeping the other variables constant. The result of these experiments can be plotted as flux rate versus metabolite concentrations in logarithmic coordinates . In the case of power-law systems, the kinetic order of the investigated variable corresponds to the slope of the plotted line and is obtained by linear regression (e.g., [35,103]). Under ideal circumstances, sufficient experimental measurements can be collected to allow the regression analysis. However, the data usually contain noise and consist of only a few measurements, which make the regression rather susceptible to experimental uncertainties. As an alternative, parameter values can be estimated by experimental measurements of logarithmic gains (e.g., [111,112]). This approach is based on perturbing variables in the interesting portion of the pathway and recording the corresponding changes after perturbations. The information about modulation, including flux rates and concentrations, is collected to calculate the kinetic orders. Additional methods that have some similarity to these steady-state estimations are discussed in Section 5.1.1.
While steady-state measurements or simple perturbation experiments around the steady-state can effectively be used for an estimation of flux rates, incomparably more information about the system is contained in measurements of metabolite concentrations at sequential points in time that may include wider deviations from the steady state. Modern high-throughput techniques of biology are capable of producing these types of time series data, and they have begun to render distinctly different options for modelling metabolic systems possible. Because this type of estimation begins with comprehensive data at the level of intact systems and leads to inferences of parameter values at lower levels of organization, namely features of individual processes, it is called “top-down” or “inverse.” The primary experimental tools for measuring genomic time series data are microarrays and RNA-based gene silencing (e.g. [61,113]). For proteomic time series, two-dimensional gel electrophoresis and mass spectrometry could be employed, but technical challenges have kept the number of proteomic time series analyses small (e.g. ). Dynamic metabolite concentration profiles are presently obtained with methods of nuclear magnetic resonance (NMR) (e.g. [115,116]), mass spectrometry (MS) (e.g. [117,118]), and high performance liquid chromatography (HPLC) (e.g. [119,120]). In contrast to the “local” data obtained from traditional experiments, the clear advantages of using “global” data are that the information is collected within the same organism, obtained under the same experimental condition, and sometimes even in vivo. Time series data contain enormous information on the structure and regulation of the biological system they describe. However, this information is mostly implicit, and it is very challenging to extract it from these data due to the complexity and nonlinearity of biological networks. There are several distinct challenges of this approach, some of which are readily anticipated, while others are surprising and puzzling. We describe these challenges in detail in the next section.
The challenges of model identification from time series data may be the result of biological and/or technical features of the case. They generally fall into four categories, namely: data related issues, model related issues, computational issues, and mathematical issues [85,121,122]. Responding to these challenges, numerous modeling techniques for the analysis of dynamic data have been proposed. Representative examples are summarized in Table 2 and described below and in later sections.
Biological datasets usually contain noise and measurement errors, and they are seldom complete. Typical scenarios of missing data points include that the data are sparsely missing, that data collection is lacking at certain time points, or that entire time series are missing, either due to technical issues (e.g. the concentration is below the detection limit) or to the possibility that relevant metabolites may be unknown. In simple cases of relatively few missing points, standard interpolation and smoothing methods may lead to satisfactory solutions. Many methods are available to address these tasks; some will be discussed in a later section.
In other cases, data may be missing for entire time courses, yet qualitative or semi-quantitative information is available. To bridge this type of gap between specific wet experiments, biological insight and intuition, and the construction of mathematical models, one may pursue the strategy of concept map modeling, which permits the inclusion of semi-quantitative information on expected responses of the system based on biological knowledge and intuition . We will further describe the general features of concept map modeling in Section 7.
Besides missing data points or time series, other data related problems must be addressed. For instance, a pathway model may be designed in such a fashion that all mass is accounted for throughout the experiment. However, if the experimental data are noisy or incomplete, it is possible that non-negligible amounts of mass are apparently lost or even gained. This inconsistency may cause problems for the estimation of parameter values. We discuss this issue in more detail in Section 6.2.2.
Even if the time series are complete, they are almost always noisy. They are also often affected by uncertainties about the particular experimental conditions at the time of observation. For instance, temperature and pH may affect the reaction mechanism or speed, but they are not always reported. As far as possible, uncertainties should be taken into account as these can possibly influence the parameter values and thereby have an impact on the predictive accuracy of the resulting model.
Other potential problems in the dataset may include that the time series are non-informative, e.g., consist essentially of constant time profiles, or that they are collinear, as is the case when some of the variables are proportional to each other along the entire time horizon. This collinearity between time series might cause ill-conditioning in the estimation process. The problem should be diagnosed beforehand by checking the conditional number or correlation coefficient of the data, and remedied if possible by pooling variables or merging essentially constant variables with the rate constants of the term. We will return to this issue in Section 7.
At the opposite side of the spectrum of data availability, it may also happen that data have been measured but that they are difficult to include in the model. This situation occurs especially with “ubiquitous” metabolites like ATP and water, which may clearly affect the dynamics of a metabolic pathway system, but are involved in so many processes that a comprehensive model simply cannot be constructed. A possible solution is a “partial modeling” strategy that permits the mixing of well-defined components with components whose dynamics is only known in the form of time series that are observed but cannot be formulated in terms of other model components (e.g., [1,122]; see [82,123,124] for a similar strategy in a different context). As a specific example, we analyzed time series data describing glycolysis and lactate production in Lactococcus lactis . The data contained ATP and NAD+ concentrations, measured over time, but it was impossible to formulate the ATP and NAD+ dynamics as functions of the system variables, because both are involved in very many reactions, most of which were not modeled. As a solution, the better-defined components were formulated as differential equations in BST, and their dynamics included ATP and NAD+ as variables. However, ATP and NAD+ were not modeled as differential equations but as “off-line” data, which entered the system as time-dependent “forcing functions.”
The inverse problem requires a mathematical model that captures the dynamics of the data in a suitable fashion. However, there is an unlimited variety of nonlinear structures and mathematical formulations that could be potential candidates for the optimal data representation. We have introduced some of the modeling frameworks and their pros and cons in Section 2. Here we highlight the specific challenges associated with model selection in top-down modeling approaches.
It seems that there are good reasons for selecting model formalisms that are intended to represent the underlying chemical reactions mechanistically. However, this mechanistic approach is not always the best choice. The reasons are, first, that high-throughput time series data are essentially never of sufficient accuracy to discern among possible underlying reaction mechanisms. Second, the true mechanisms underlying the data are simply not known. Third, traditional kinetic rate functions, such as the Michaelis-Menten rate law, are not necessary the best choice for in vivo data [125,126]. Instead of trying out a roster of different mechanistic formulations that could potentially be appropriate, it is often more efficient to use a generic nonlinear approach in the form of a suitable approximation that is based on criteria like: the ability to capture important mathematical features of a data set, simplicity of representing the data, mathematical tractability, and interpretability of mathematical results within the biological realm. Canonical models, as they were described before, are particularly useful for this purpose.
The estimation process itself is challenging computationally, especially if the model consists of nonlinear differential equations . It is no surprise that the difficulties grow with the size and complexity of the system, which usually translate into increasing numbers of equations and variables in the model. None of the nonlinear methods are truly straightforward, and even for systems of modest size, all lead to challenging issues, such as slow algorithmic progress toward the error minimum, complicated error surfaces, lack of convergence, or convergence to local minima. Furthermore, the differential equations need to be integrated during the optimization process, unless special strategies are employed that explicitly avoid this step. The integration may be very time consuming, especially when the system is stiff, and require 95% or more of the entire estimation time . Other computational challenges include the distinction between direct and indirect effects among system variables, characterization of intermediate steps, time delays, spatial heterogeneity, and stochasticity at the level of the governing processes or the integrated system.
The actual development of algorithmic methods for extracting information from biological time series data sets is the subject of Section 4. The methods for biological inverse problems typically require a combination of techniques that include techniques for attacking the main problem of optimizing parameter values, as well as supporting algorithms, such as methods for circumventing the costly integration of differential equations, smoothing overly noisy data, constraining the parameter search space, or reducing the complexity of the inference task.
An often ignored source of problems is a (frequently unknown) mathematical redundancy in some models. Redundancy may occur in different manifestations. It is possible that different sets of parameter values, which fit the experimental data exactly equally well, are mathematically or numerically equivalent . For instance, if two parameters p and q always enter the equations in the same combination, such as (p+q), it is not possible to identify their individual values. In the context of power-law functions, collinearity between variables leads to unidentifiability of the corresponding kinetic orders [130,131]. It may also happen that solutions exhibit similar residual errors, even if they are mathematically not equivalent. One possible cause is compensation between a rate constant and the kinetic orders in the power-law terms of a particular data fit [130,132]. Mathematical redundancies and error compensation may occur within or between flux descriptions in the same or among different equations. The task of dealing with mathematical and computational redundancies has been addressed in some articles [130,132,133], but will require more work. A preliminary step in the identification of redundancies and their causes may be our recent method of Dynamic Flux Estimation , which permits unique tools of error diagnostics (see Section 7).
Despite these challenges, inverse approaches based on in vivo time series data are certainly worth pursuing, because these data are the most accurate reflections of what cells and organisms really do in environments that are as close to reality as is presently possible to measure. This level of realism is very appealing, and many researchers have worked on the development of methods that overcome some of these challenges. A representative set of methods is described in the following section.
Many of the recently developed techniques for top-down parameter estimation have been developed for BST models. Most of them are similarly applicable to other canonical models, although a few take advantage of the specific form of power laws in BST. The main algorithms and their representative references are shown in Figures 3 and and4.4. A historic listing of representative algorithms is presented in Table 3.
The central component of solving inverse problems is an efficient algorithm for determining optimal estimates. Most standard methods for this purpose naturally involve the numerical solving of the differential equations. This solution is computationally very expensive. As a specific example indicative of the problem at hand, consider a direct attempt to estimate the parameters of a five-variable S-system model from noise-free time series data with a genetic algorithm . The authors used a cluster of 1,040 CPUs, which ran for ~10 hours for each loop of the estimation program. Needing 7 loops, the entire estimation time thus was roughly 70,000 PC-hours. Analyzing this alarming situation, the distinct tasks within the optimization were clocked in detail with the result that 95% of the time spent for parameter searches involving differential equations is used for integrating the equations, while relatively little time is used to compute gradients toward the optimal estimates . If the equations are stiff, the computation time may increase to almost 100%, and even if the model is not stiff, the likelihood is high that some trial solutions during the algorithmic process could make it stiff . Therefore, even if efficient, custom-tailored integration methods are used , significant time savings can be gained by speeding up the evaluation of differential equations.
Numerical integration of the ODE system can be circumvented when the differentials are substituted with slopes that are estimated from the time series data at all measured time points. This substitution entirely eliminates the need to integrate differential equations, because the parameter estimation is subsequently executed on systems of algebraic equations. Furthermore, the equations become uncoupled so that they can be assessed in parallel or one at a time. An early implementation of this method was accomplished by manually estimating slopes from observed time series data and substituting them for the derivatives in the differential equations [35,136,137]. Voit and Almeida  implemented the slope-estimation-decoupling strategy with an artificial neural network that simultaneously permitted data smoothing and the computer-algebraic computation of slopes from the smoothed time courses.
The slope-estimation-decoupling idea has subsequently been combined with various methods such as genetic algorithms, simulated annealing, swarm methods, interval analysis, and a number of hybrid methods. Various slope estimation methods will be reviewed in detail in Section 4.2. One drawback of these methods is that it may not be easy to obtain good measurements or estimates of the slopes if the data are very noisy. However, it is still advantageous to use this approach, even if the results are coarse, since the coarse estimates may be used as initial guesses for standard nonlinear optimization methods. Other advantages of the decoupling approach are reviewed in .
In a different implementation of a similar idea, the decoupling allowed solving and fitting of one differential equation at a time instead of solving the entire system. Maki et al.  proposed this “step-by-step” strategy and Kimura et al. [82,124] introduced a similar concept called “decomposition,” which dissects a large network inference problem into many smaller sub-problems. In both methods, the variables contributing to the single differential equation being integrated are substituted with the actual observed time series data or with smoothed analogues, which are thus used as off-line inputs to the decoupled system. This approach significantly reduced the computation time. For instance, using the same artificial five-variable datasets that required 70,000 PC hours , Kimura and co-workers ran the algorithm on a single CPU, where it required only about 59 minutes to optimize each subtask.
A drawback of decoupling and decomposition approaches is that each subtask is solved independently, a procedure which does not allow the exchange of information between subtasks. For instance, the variables serving as off-line data in one equation are actually solved in another equation. Thus, if the value of one variable is updated during optimization, the information should be incorporated into optimization processes of the other subtask. This feature is especially important when there is considerable noise. Kimura et al.  proposed to solve the decomposed subtasks simultaneously using a cooperative coevolution algorithm. Since the decomposed subtasks interact with each other through their calculated time series data, the inferred model is more likely to represent the dynamics correctly.
In order to reduce the number of numerical integration steps, Matsubara et al.  proposed to use a radial basis function network (RBFN) for parameter estimation. RBFN is a type of artificial neural network that uses radial basis functions as activation functions; it has been shown to be able to approximate nonlinear time series data efficiently . In order to examine the performance of RBFN, Matsubara and co-workers proposed two schemes: one used a simple genetic algorithm (SGA) with numerical integration, and the other combined RBFN with a genetic algorithm in the input data selection phase. Both schemes were examined in metabolic pathways using Michaelis-Menten equations. While SGA improves the fitness between parameterized model and time series data and integrates every time during optimization, RBFN predicts the optimal parameter values by learning the relationship between parameters and fitness values using slopes to replace derivatives and integrates the system only once at the last step. Therefore, numerical integrations used to evaluate the fitness are reduced from many to one. The results indicated that the RBFN scheme halved the computation time and increased the success rate of the optimization task.
An alternative approach avoiding numerical integration is a modified collocation method, which converts ordinary differential equations into algebraic equations which directly adopt the measured data to approximate the dynamic profiles at sampling points. This approximation not only reduces computation time, but also decouples the equations so that parallel computation is possible for the parameter estimation. A collocation method was combined with hybrid differential evolution (HDE) to determine the global solution of an estimation task . Again, applying this type of “uncoupling” strategy in combination with other estimating methods reduced the computation time dramatically.
As a crucial part of the slope-estimation-decoupling strategy, decent estimates of the slopes are required. If the data are more or less noise-free, simple linear interpolation, splines [141–143], B-splines , the so-called three-point method , or even hand fitting  is effective. If the data are noisy, it is useful to smooth them, because the noise tends to be magnified in the slopes. Established smoothing methods again include splines, as well as different types of filters. Artificial neural networks (ANNs) have been shown to be useful in a number of applications of biochemical pathways modeling . They are so general and flexible that they are considered “universal functions” that are obtained from training the ANN [128,147–149]. The main advantage of using ANNs is that the resulting time traces can be trained to fit the data arbitrarily closely and that they have an algebraic format for which the slope can be computed straightforwardly with methods of computer algebra [146,150]. Furthermore, the universal output function provides an unlimited number of interpolated data points within the time interval of interest. Other advantages of ANN are reviewed in  and . The ANN method was shown to determine the smoothed traces very efficiently even if the data contained considerable noise, as long as the true trend was well represented. However, the interpolating function resulting from the ANN solution is a superposition of sigmoidal functions and has the tendency to lead to artifacts in the derivatives, which manifest in slight, but undesirable bias in the smoothed traces.
Alternatives to ANNs are filters, such as the popular Kalman or Savitzky-Golay filter or the Whittaker filter which was proposed over eighty years ago . Recently, Eilers presented a matrix form of this older implicit method, which was called a “perfect filter” . Vilela and co-workers configured the Whittaker-Eilers smoother  by adapting Rényi’s second-order entropy of the cross-validation error as the optimization criterion. The filter, implemented in the software AutoSmooth, can be used to extract signals and derivatives from time series with non-stationary noise structure .
To ensure that the results of a parameter estimation fall within reasonable ranges it is often useful to constrain the optimization process, including guesses for the initial values, to suitable ranges or permitted extreme values for all parameters. For instance, in BST representations, the structural features of a system are mapped onto model parameters in a unique fashion, as described in Section 2.3.3. Therefore, if the network structure and regulation are known, one may be able to decide immediately whether the kinetic order of a variable Xj is positive, negative, or zero, depending on its influence (activation, inhibition, or no effect) on variable Xi. Furthermore, rate constants are always non-negative, and kinetic orders in BST pathway models are known to be real numbers with typical values between −1 and +2.
Some other supporting techniques aiming to reduce the parameter search space include the following. Kutalik et al.  characterized a one-dimensional basin of attraction containing the true optimum with minimal error. Tucker and Moulton  proposed a method based on interval analysis which allows exhaustive searches of the entire set of parameter values with a finite number of steps. Tucker et al.  used constraint propagation to find the possible ranges of parameter values, thus significantly constraining the parameter search space.
The typical approach of modeling is to collect network information and translate the wiring diagram into a symbolic model, which only contains a limited number of parameters since the biological systems are usually sparsely connected. However, when the topology of the system is unknown or only partially known, it appears that one must initiate the search with a full symbolic model with all parameters free. When the system is relatively small, it may be feasible to exhaust all possibilities and to find the optimum. However, when the number of variables and parameters grows, all methods of parameter estimation eventually run into problems of “combinatorial explosion,” which makes the estimation process extremely difficult and the solutions problematic. This explosion can be tamed to some degree by constraining the connectivity within the system by systematically identifying the network structure or gradually “pruning” unlikely connection during optimization process. Of benefit in this context is the observation that biological systems are seldom fully connected and that indeed most nodes are only directly connected to a small number of other nodes [157,158]. These issues will be reviewed in detail in Section 5, in the context of structure identification techniques. In this section, we focus only on parameter pruning methods.
The rationale behind the pruning techniques is closely related to the characteristic of BST models. As briefly mentioned in Section 2.3.3, structure identification tasks can be translated into parameter estimation problems if the parameter values directly map to the network, as it is the case with BST representations. To recall this mapping, the kinetic orders fij, gij and hij for BST quantify the regulatory effect of variable Xj on a production or degradation term in the equation of variable Xi. If the magnitude of the corresponding kinetic orders are very small or close to zero, the connection between variable Xj and the dynamics of Xi is likely to be negligible. Therefore, these low intensity connections can be purged during optimization, which not only helps to detect a reasonable and parsimonious model of the true pathway structure, but also reduces the parameter search space for further optimization.
The simplest manner of “pruning” a possibly highly connected network is to define a threshold for the absolute value of each type of parameter, below which values are set to zero [128,159]. In addition, since the likelihood that a variable exists in both the production and degradation terms with non-zero values in the S-system model is low, the smaller of the kinetic orders is more likely to be zero and the value of the other one is adjusted accordingly .
Some authors have suggested more sophisticated methods for this pruning process. As an extension of the objective function for optimization, various articles have added to the residual error the sums of the absolute values of kinetic orders as a penalty term in the cost function. Thus, this basic pruning method for BST models penalizes all small kinetic orders, which have little effect on the system dynamics, and prevents the model from including false-positive interactions that unrealistically inflate the model [134,160]. To improve this condition further, Kimura and co-workers [82,124] introduced a penalty term that rearranged kinetic orders in ascending order based on their absolute values and eliminated those considered insignificant. Furthermore, accounting for the observation that very few factors modulate both the production and degradation of a specific variable, Noman and Iba  proposed an alternative representation of the penalty term.
No matter what kind of penalty term is chosen, pruning approaches pose an obvious challenge. Namely, the weighted coefficient in the penalty term needs to be carefully tuned since it affects the results of the structure identification task. So far there are no clear guidelines for setting suitable penalty weights. Stochastic ranking may be used to alleviate this difficulty since it aims to balance the error and penalty term in the objective function . However, this method requires an additional parameter defining the probability of the error term for comparisons in ranking. Cho et al.  proposed a distinctly different way to retain the sparseness feature in biological pathways without adding extra terms to the objective function, which they coined S-tree representation. The S-tree is a tree representation of the S-system, where the number of sub-trees corresponds to the number of ordinary differential equations in the system. Each sub-tree is divided into two parts; the left part represents the production term and right part represents the degradation term. The depth of the S-tree is always three and the root node is at depth zero. Since S-tree modeling is intrinsically suitable for representing sparse networks, an S-tree together with genetic programming has the potential to infer network topology and find parameter values in a more efficient way without any a priori knowledge or adding a penalty term. To avoid assigning a coefficient weight to the penalty term, Liu and Wang recently proposed an alternative method based on multi-objective optimization [164,165]. Instead of minimizing the residual error using a single objective function either in concentrations or slopes, they minimized the concentration error, slope error, and interaction measure simultaneously. The authors proved that the algorithm guarantees the minimum solution for the constrained problem to achieve the minimum interaction network for the inference problem. The approach avoided assigning a penalty weight for sums of magnitude of kinetic orders.
Pruning methods are also used in optimization approaches for determining parameter values, as described in the next section.
The parameter estimation task is traditionally formulated as an optimization problem that minimizes an objective function measuring a generalized distance between experimental data and model predictions. The Euclidean distance is the most commonly used and often refers to a least-squares error criterion. Other fitness evaluation methods include information based criteria [166,167]. Two objective functions are typically used for parameter estimation in metabolic models: a concentration error based objective function and a slope error based objective function (e.g. ). The concentration error based objective function is a straightforward calculation of the sum of squared distances between the metabolite measurements and the predictions. The simulation profiles are usually obtained by applying a numerical integration method to solve the differential equations like Eq. (1). The integration process can be computational costly, especially if the system is stiff (see Section 4.1). As an alternative, the slope error based objective function employs the decoupling technique as described in Section 4.1 and uses the slope information for evaluating fitness of the function. That is, it calculates the sum of squared errors between the measured slopes from the raw data (or upon smoothing) and the predicted slopes.
The most prominent methods for parameter estimation from time series data can be grouped into gradient-based methods, stochastic search algorithms, and others that do not belong to the first two groups. These optimization methods are reviewed in the following paragraphs and summarized in Figure 4.
The most natural choice for estimating parameter values is presumably a gradient based regression, and many of the commercial methods of this type have been applied to metabolic models. Among these are Gauss-Newton and Levenberg-Marquardt methods, which are included in all major software packages of the field and will not be reviewed per se [168–170]. A comparison of some of these methods in the context of pathway models can be found in . Marino and Voit  proposed a gradient based algorithm for finding the parameter values using BST models that comprises three modules in a novel fashion: model generation, parameter estimation for model fitting, and model selection. First, plausible initial models are generated in a step-by-step manner, upon decoupling and limiting connectivity (see Section 5.2 for detail). Secondly, each differential equation is fitted separately using the Levenberg-Marquardt method while replacing the other variables with raw data of smoothed traces. In the third phase, the model is compared to earlier, simpler models, and a statistical test decides whether the increased complexity of the model structure is warranted, as judged by the residual error. If the improvement is significant, a new, more complex model is generated, and this process is iterated until further advancements become insignificant.
Kutalik et al.  proposed an intriguing Newton-flow optimization method for parameter estimation in S-system models. The method starts with decoupling the differential equations and setting up an objective function for each equation. The next step is to select suitable start guesses and bounds for parameters and run a Newton method to obtain several points in the parameter space that correspond to reasonable, coarse solutions. The authors found that this space of coarse solutions contains a one-dimensional attractor. Standard regression allowed them to estimate the parameters of this attractor. Afterward, the Newton method was performed again using the initial guesses lying on the estimated attractor to find the true optimal of the parameter values. The interesting feature of this method is that most (or maybe even all) good parameter solutions seem to lie on one-dimensional manifolds within the high-dimensional parameter space. Optimization along this curve is comparatively easy. A potential problem of the method is that the original initial guesses for the parameters must lie within the basin of attraction of the one-dimensional manifold. Otherwise, each run may lead to disjoint sections of the parameter space.
Because biological systems are usually nonlinear, the problem of parameter estimation can be stated as a nonlinear programming problem (NLP) subject to nonlinear differential-algebraic constraints . Because of its nonlinear and constrained nature, this inverse problem is usually non-convex. Therefore, most of the traditional nonlinear algorithms involving gradient methods run the risk of getting trapped in local optima, depending upon the degree of system nonlinearity and the initial starting point .
Several classes of stochastic methods are available for global optimization. They include evolutionary computation, simulated annealing, adaptive stochastic methods, clustering methods, and other meta-heuristics, such as ant colony optimization and particle swarm optimization. These algorithms have been applied to parameter estimation tasks with the goal of finding global solutions, especially in the context of identifying the structures of gene regulatory networks .
Evolutionary computation (EC) techniques, also known as biologically inspired methods, include genetic algorithms, evolutionary programming, evolution strategies, genetic programming, as well as their variants. They are attractive because they have a high potential of finding (at least the approximate locations of) global optima. Genetic algorithms (GA) have been shown to be very useful and practical in parameter estimations of biological systems (e.g. [150,160,172,173]). Using a conventional simple genetic algorithm (SGA), Tominaga et al. inferred parameter values of a small network, but only with a very limited number of parameters, and the convergence rate was low . SGA typically has two problems: early convergence in the fast stage of the search and evolutionary stagnation in the last stage. Kikuchi et al.  enhanced SGA by using a more robust real coded genetic algorithm (RCGA) and improved the conventional cost function by adding a penalty term to prune unlikely connections in the investigated gene network using the S-system formalism. In addition, they employed a novel crossover method and introduced a gradual optimization strategy in the procedure. The results showed that the algorithm successfully inferred the network structure with faster convergence rate, optimization speed, and with more parameters predicted correctly, compared to the traditional GA. However, the approach turned out to be computationally very costly because of numerical integration of the entire system of differential equations (see Section 4.1).
Other modifications were made to improve the efficiency of SGA using time series data in S-system form. Examples include: a hybrid algorithm of SGA with a Modified Powell method ; a hybrid algorithm of SGA for static Boolean networks applied to an S-system with steady state and temporal data ; and a combination of RCGAs with unimodal normal distribution crossover and minimal generation gap to optimize parameters in S-systems [176–178]. Daisuke and Horton optimized an S-system model with a distributed genetic algorithm with “scale-free” properties . Ho et al.  proposed an intelligent two-stage evolutionary algorithm (iTEA), which used an intelligent to solve decomposed ODEs independently, then combined all solutions from each subtask and used an orthogonal experimental design-based simulated annealing algorithm to refine the solution.
Spieth and coworkers [181,182] proposed a memetic algorithm (MA) consisting of two parts: a local search with an evolutionary strategy (ES) for parameter estimation, and a global GA based search framework for structure identification, where the former is embedded within the latter part. They tested the algorithm on an S-system model and the results showed that MA was better suited for inferring genetic networks than a standard ES or GA. In follow-up work, they showed that feedback from the local search to the GA based search can further improve the performance of MA .
Kimura et al.  used an evolutionary algorithm called Genetic Local Search with distance independent Diversity Control (GLSDC) and combined it with a decomposition strategy to estimate S-system models of gene regulatory networks. The proposed method included an estimation technique for the initial gene expression levels and enabled the reconstruction of medium-scale genetic networks with noisy data. They also showed that the combination with a cooperative co-evolution algorithm can further improve the accuracy of prediction . Okamoto’s group also proposed evolutionary search techniques, such as the Network-Structure-Search Evolutionary Algorithm (NSS-EA) and a variant, the Grid-Oriented Genetic Algorithm Framework (GOGA Framework). They employed an S-system as the underlying mathematical model and used a GA as search engine to infer network structure [184–186].
Noman and co-workers recently incorporated their previously developed techniques into an improved memetic algorithm for inferring gene regulatory networks [161,166,187–189]. They used differential evolution (DE) along with a hill-climbing local-search method in their evolutionary algorithm. An information criterion-based fitness evaluation was introduced instead of the conventional least-squared errors approach.
Tsai and Wang  used hybrid differential evolution (HDE) for estimating a satisfactory, though not optimal solution, and then used the solution as the starting point for a gradient-based optimization method to obtain refined solutions. As described in Section 4.1, they used a modified collocation method to avoid direct numerical integration. In their recent work, they implemented HDE combined with a multiple-objective optimization approach (see Section 4.4 for review) to infer biochemical networks in S-system format .
Genetic programming (GP) has also been employed to discover the topology of metabolic pathway from time-series data (e.g. ). GP is an evolutionary algorithm that evolves mathematical expressions or computer programs. Traditionally, GP represents a mathematical expression or computer program as a tree structure, in which every tree node has an operator function and every terminal node has an operand. The general process of GP includes: initialization (randomly generate trees as individuals), evaluations (calculate fitness of each individual), selection (select individuals from the group base on probability), crossover (randomly select two individuals as parents and swap randomly chosen sub-trees of the parent trees), and mutation (such as insertion or deletion of terminal nodes) (cf.  for detail). The GP process makes mathematical expressions easy to evolve and evaluate. Therefore, in contrast to GA algorithms, which usually require defining equations before optimization, GP provides a general approach for finding arbitrary equations from time series data without specific knowledge of the equation. The ordinary GP is not always effective in finding the parameter values because the method relies mainly on the combination of randomly generated constants. Sagamoto and Iba  therefore used a least-mean-square (LMS) method along with ordinary GP to improve efficiency, using an S-system as one example. Their results showed that the fitness values improved faster in the early phase with the LMS method compared to the non-LMS method, since the former seemed to provide a better seed for the GP search.
Sugimoto and co-workers  implemented GP along with adding a penalty term to the cost function and introducing numeric mutations to the conventional procedure. They tested this method by predicting two equations of a metabolic reaction scheme regarding adenylate kinase and phosphofructokinase in Michaelis-Menten format, the equation of which is hard to derive if the underlying mechanism is not known. While their results showed that the algorithm can predict the equations with relatively simple forms, the method is already very time consuming for this relatively small system.
Kim et al.  adopted a symbolic pre-processing regression step in GP to avoid time consuming numerical integration, since the estimation of slopes for each data point in the time series can be obtained from the results of GP. Cho and co-workers  took advantage of the fact that GP has an evolving tree structure for given data and proposed S-tree based genetic programming for parameter estimation and structural identification in S-system models. As introduced in Section 4.4, this approach intrinsically accounts for the sparseness of the biological network. Therefore, even though no a priori knowledge about the network is known, the S-tree based GP can still identify the underlying structure rather efficiently without adding a penalty term in the objective function.
The discussion in the previous paragraphs indicates that a considerable number of recent proposals applied evolutionary algorithms to tackle inverse tasks with BST models. So far no comprehensive comparison among these algorithms has characterized their relative efficiency, robustness, and accuracy. However, some more limited comparisons have been presented. Moles et al.  compared some stochastic global optimization methods using the case study of a biochemical model that consisted of 36 parameters and was formulated as a set of eight ODEs. This model was formulated in Michaelis-Menten representation, which could not take advantage of the highly structured format of BST representations. Spieth et al.  compared six evolutionary algorithms in three model frameworks: linear weight matrices, S-systems, and H-systems, where one fitness function was used to evaluate the convergence of algorithms. A comprehensive comparison of evolutionary algorithms is still needed.
Simulated annealing, colony optimization, and particle swarm optimization are also stochastic optimization methods. Simulated annealing (SA), a physically inspired method, was created to simulate the heating and cooling process of metal or glass, where atoms (solutions of the parameter estimation task) are allowed to leave their current state of low energy (fitness) during heating and have a chance of finding an even lower state (better fit) during cooling . SA can behave as a global or local optimization search and automatically switches from a global to a local search when the “temperature” goes down. Gonzalez et al.  adapted SA for S-system parameter estimations from time series data. They tested the algorithm using three artificial datasets and assumed that the structure was either known or unknown and solved the entire set of ODEs per integration or upon decoupling. They also applied the algorithm to a real biological system.
Ant colony optimization (ACO) was inspired by the behavior of ants during the search for short paths between their colony and food sources, using pheromones that attract other ants, which then increase the amount of pheromone . Over time, ever shorter paths (better solutions of the estimation task) become more popular. ACO is a probabilistic technique for solving computational problems that can be reduced to finding good paths through nodes in a graph. Zuñiga et al.  adapted ACO for S-system models by treating each metabolite as a node in a graph and inferring how other nodes were connected to it. They called their algorithm for identifying network structure a “discrete ACO.” As an extension of ACO they furthermore proposed a variant of an enhanced aggregation pheromone system (eAPS) for parameter estimation tasks involving S-systems, called a “continuous ACO.” The discrete ACO starts with a fully connected graph which corresponds to a set of equations where all variables are included in every equation. Their preliminary results showed that ACO produces good results when the test systems are very small, However, although the discrete ACO was able to eliminate some nodes (metabolites) from the graph in larger systems, it had problems eliminating unlikely connections and thus still produced an unreasonably large search space for the parameter estimation task. The authors concluded that the large search space might have been the reason for the continuous ACO to get trapped in local minima during the parameter estimation phase.
Particle swarm optimization (PSO) is a stochastic, population-based evolutionary computation algorithm. The original form of the PSO algorithm, which was motivated by social-psychological principles such as bird flocking and fish schooling, was first described by Eberhart and Kennedy . In PSO, each potential solution is represented as a particle. A collection of potential solutions is called a swarm which consists of particles that fly around in a multidimensional search space. During flight, each particle adjusts its position according to its own experience and also collaborates with its neighboring particles through communication. When a particle encounters a promising solution, the area surrounding the solution is further explored by the swarm. Therefore, PSO combines local search methods with global search methods. Naval et al.  adapted and refined PSO to scan the parameter space of a BST model.
Some methods that aim to reduce the parameter search space using BST formalisms are described in Section 4.3 [83,155,156]. For linear parts of pathways, a technique of “peeling” terms  can be applied to models in BST to convert the nonlinear parameter estimation task into a series of linear regression tasks. Specifically, beginning with an equation that contains only one unknown power-law term, the differentials are substituted by slopes and the parameters of the unknown terms are estimated by linear regression. The results are fed into the equation of the subsequent metabolite, thereby making it amenable to linear regression as well, and this process is iterated to the end of the linear pathway.
As described in Section 4.5.1, the traditional gradient methods usually run the risk of getting trapped in local optima. To alleviate the problem, Polisetty et al.  proposed a branch-and-reduce algorithm to convert inverse problems involving GMA models into a convex optimization problem that is guaranteed to obtain the global solution within a predescribed parameter space. The major drawback of this method is its computational complexity.
Alternating regression (AR)  employs a decoupling technique for systems of differential equations and dissects the nonlinear parameter estimation task for S-systems into iterative steps of linear regression. The method is uniquely geared toward BST systems, because it utilizes the fact that power-law functions are linear in logarithmic space. AR is extremely fast in comparison to conventional methods and works well in many applications, if it converges. In cases where convergence is an issue, the fast speed renders it feasible to dedicate some computational effort to identifying suitable start values and search settings. AR is beneficial for the identification of system structure in S-systems as well. An extension of AR was successfully applied to S-distributions within the field of computational statistics .
As an extension of AR, eigenvector optimization (EO)  is based on a matrix formed from multiple regression equations of a decoupled S-system that is considered in logarithmic space. In contrast to AR, EO operates initially only on one term, whose parameter values are optimized completely before the complementary term is estimated. It was demonstrated that the EO algorithm converges fast and can be expected to converge in most cases, without necessarily requiring knowledge of the network structure. EO is easily extended to the optimization of network topologies with stoichiometric precursor-product constraints among equations.
Another recently proposed approach to metabolic systems estimation, called Dynamic Flux Estimation (DFE), consists of two distinct phases and applies particularly well to GMA models . It corresponds to stoichiometric analysis, but does not consider systems in steady state but rather over a time horizon. The first phase attempts to establish a linear relationship among all fluxes in the system at each time point. Under certain rank conditions, these relationships form a matrix that can be solved uniquely. This first phase is entirely model-free and essentially assumption-free and includes a diagnosis of inconsistencies within the time series, and between the assumed system topology and the given data. The result of the first phase consists of numerical representations of all fluxes as they depend on the variables affecting them. The second, model-based phase addresses the mathematical formulation of these flux representations as explicit functions of the involved variables. Different from currently available methods, this phase allows quantitative diagnostics of whether the chosen mathematical representations are suitable. The two-phased approach thus permits rigorous, quantitative diagnoses of the data, the model structure, the assumptions made in the choice of flux representations, and of the various causes of residual errors. Preliminary results suggest that the DFE approach is more effective and robust than alternatives that are presently available, if sufficient suitable data are available. Its combined model-free and model-based analyses reduce compensation of error between equations and between flux terms and promise significantly improved extrapolability toward new data or experimental conditions. Its diagnostic tools pinpoint causes of inadequate fits between model and data and suggest either changes in assumptions related to model choice or the use of data as un-modeled “off-line data.” The main drawback of DFE is the requirement of rather comprehensive metabolic time series data, which however can be obtained in cases with already existing experimental methods. Furthermore, a direct application of DFE requires that the stoichiometric flux system is of full rank, which is usually not the case and requires additional “substitute information” . Other issues needing refinement are related to missing data, missing flux information, error compensation among the parameters within a given flux, and ill-characterized systems topologies (see also Section 7).
Other methods which were developed recently for inverse problems in biological systems are discussed in [139,205–210]. However, these methods have not yet been implemented specifically for BST applications.
As mentioned in Section 1, the traditional approach of modeling begins with collecting network information that is translated into the design of a stoichiometric wiring diagram, which may then be converted into a fully kinetic metabolic pathway model, if desired. The translation and conversion more or less reflect the actual biological system as long as the input information is essentially correct and complete. In reality, information on network connectivity and regulation is often only partially known and seldom fully understood. As a consequence, the model design phase is subject to uncertainties, up to a point where it seems impossible to determine even the initial wiring diagram. Some of the identification methods discussed in this section ameliorate the situation and make it possible to deduce structure and regulation from time series data, at least in principle. Before we discuss these methods, it is beneficial to revisit aspects of the investigated system and the experimental data that have a direct effect on the complex identification task.
The need for valid system identification can be described in three aspects. First, wrong hypotheses regarding variables and interactions to be included in the model tend to lead to wrong interpretations of the results. Second, overly complex model representations may provide good fits to the observed time series data used for estimation but are unlikely to perform as well when tested against new datasets, due to over-fitting. Third, the inclusion of too many components and interactions in the model will eventually result in problems caused by combinatorial explosion, which means that any computational technique will ultimately be overwhelmed by the rapidly increasing number of equations, variables, and interactions between variables in large systems.
Fortunately, biology naturally offers a counteracting and very beneficial feature: namely the likelihood that a real biochemical, genetic, or proteomic network is fully connected is very low, because most metabolites (genes, proteins) are connected only to a limited number of other metabolites (genes, proteins); in fact, the vast majority of metabolites (genes, proteins) are involved in fewer than four or five processes each [157,211,212]. To take advantage of this fact of nature, it is therefore a desirable goal to precede any estimation attempt with a concerted effort to limit the number of candidate (structural and functional) connections within a system, as far as this is objectively possible. This a priori type of limitation can very significantly reduce the parameter space that must be searched, because structure identification and parameter estimation are closely related to each other, at least if canonical models are used.
In this section, we review some of the structure identification techniques; they can be categorized into two groups. First, model-free, coarse structure identification methods can be used to prescreen the particular situation at hand. These methods are based directly on the data and involve the determination of an estimated Jacobian matrix after small perturbations about the system’s normal operating point, deductions from direct observation of the time profiles, a correlation-based approach, and a Bayesian network technique. The second group consists of model based structure identification methods, including “simple-to-general” and “general-to-specific” modeling strategies, as well as various additional methods using time series data within the BST framework. These methods and representative references are summarized in Figure 5.
Much of the information necessary for identifying network structure depends on dynamic experiments. One type of such experiments is the measurement of transient responses of the system after small perturbations about the steady state. If the perturbations are small enough, the system can be expected to behave in a roughly linear fashion. Thus, the measurements may be used to populate the Jacobian matrix of the corresponding linearization, which then reveals the connectivity of the network. Over the past two decades, several proposals have been made to obtain the Jacobian matrix from experimental observations. Chevalier and co-workers  solved the Jacobian by applying multilinear least-square fitting to perturbed data. This approach is straightforward but very sensitive to noise and missing data points, because the crucially important differencing procedure is prone to generating large errors. To avoid instabilities due to numerical differentiation, the same group suggested using an integral representation, which expressed the solution in terms of eigenvectors and eigenvalues and solved the equation using nonlinear regression . The advantage of this approach is that no differentiation is needed and hence the slopes need not be estimated. However, the drawback of this method is that the fit to a sum of exponentials with undetermined exponents is numerically somewhat problematic, and the nonlinear regression does not necessarily provide a solution which fits the data.
To overcome this difficulty, Sorribas et al.  suggested to reformulate the integral representation of the target function by reducing it to a multilinear regression problem. As the result, the eigenvalues of the Jacobian in the previous method can be easily calculated. However, the computation of eigenvalues is again rather sensitive to noise and rounding error, rendering the method not very reliable unless the multiplicities of the eigenvalues are exactly known. In order to avoid this problem, Díaz-Sierra and co-workers  proposed a variation to the previous methods, in which they directly obtained the Jacobian by expanding it in its Taylor-series without searching for eigenvalues. This methods yielded faster convergence.
All these methods are based on linear approximation, which is valid as long as the perturbation from the steady state remains relatively small. Thus, on one hand, the range of deviations needs to be small enough to yield a sufficiently accurate representation. On the other hand, the perturbation must be large enough to generate measurable responses. To alleviate this dilemma, Veflingstad et al.  suggested using the entire time course and fit the data in a piecewise linear fashion, using as an illustration example an S-system within BST. In the proposed method, the time series is subdivided into appropriate time intervals and the linearization is computed about a chosen operating point within each subset. Therefore, instead of focusing on one operating point, most reference states are different from the steady state. The results show that the piecewise approach is more likely to capture the relationship between variables in the system and can tolerate larger perturbations. The authors also showed that the collection of estimated coefficients resulting from different variations of linearization provided very strong clues about which variables were likely to be involved in a given equation and which were not. These clues reflected likely parameter ranges or likely constraints on parameter values of the true model. A major drawback is that this method does not identify parameter values per se. For instance, as discussed in the original paper , it does not allow a distinction between various combinations of gij and hij in the S-system form because only their difference is being assessed as a single parameter in the linearization. However, this information is valuable. If additional information on the Jacobian matrix and both the concentration and fluxes at steady state are known, the difference between gij and hij can be directly calculated . Also, if the difference has a magnitude that is significantly different from 0, it is likely that one of the kinetic orders is zero, because it is rare that a variable influences both production and degradation of the same variable. Therefore, if one can detect which connection may be omitted, the kinetic order can be computed straightforwardly.
Hatzimanikatis, Floudas and Bailey [218,219] indirectly contributed to the topic of structure identification per linearization by optimizing not only the production of yield in an S-system at steady state, as it has been done many times (e.g. [34,220,221]), but by also optimizing its regulatory structure. This numerical and structure optimization task led to a mixed integer linear programming (MILP) approach, for which standard software is available.
Unlike the previous methods for determining the Jacobian matrix by examining the linear properties on small amplitude perturbation near one or more operating points, the network connectivity can be deduced to some degree from direct observations on responses to perturbations of arbitrary amplitude made at different locations in the network. Vance and coworkers  proposed a strategy based on perturbing different components in a network and showed that relationships between the perturbed component and the remaining components may be deduced by observation of features in the response profile. These features include the order and size of the extreme values of the unperturbed components in response to the perturbed component, and the initial slopes of the time series at the perturbation. The former reflects the topological distances among the perturbed components and the remaining components in the network, while the latter reveals whether the components are directly affected by the perturbed variable or not. This distinction is accomplished by checking if the initial slopes are nonzero or zero upon perturbation. Vance and collaborators showed that this approach works well in some artificial networks including branching, feedback, and regulatory interactions. This method was also applied to an in vitro experiment with a glycolysis system, where the authors measured concentration changes in the reactor following impulse changes of different reaction metabolites . From the experimental time series data the authors were able to identify some of the causal connectivities among the metabolites in the reaction pathway. Even though the method performed well in the synthetic time series and with experimental data from relatively small systems, this approach may not be applicable to more complicated networks, where the interpretation of profiles and the network reconstruction must be expected to be much harder. The emerging field of causality analysis [224,225] may be helpful for this type of analysis in the future.
Some alternative approaches have been suggested for the reconstruction of chemical reaction networks. Arkin and co-workers [226,227] showed how correlations among components measured in the system may be used to infer or reconstruct a chemical reaction pathway. The approach, termed correlation metric construction (CMC), is based on the calculation and analysis of a time-lagged multivariate correlation function of time series data that are subjected to a series of random, large amplitude changes in the input concentration. The correlation information is used to construct the distance matrix and interpreted using a two-dimensional graph obtained with a projection technique called multidimensional scaling (MDS). The graph represents the connectivity and the strength of interactions among the species in the network. For instance, the shorter distances in the graph imply stronger connections while longer distances represent weaker interactions. The approach was tested experimentally on a part of an in vitro glycolysis system containing eight enzymes and fourteen metabolites . Along the same lines, Samoilov and collaborators  proposed two methods, entropy metric construction (EMC) and entropy reduction method (ERM), for the analysis of correlations between species from time series data and the inference of their underlying network.
Another approach of network structure identification is statistical in nature and uses Bayesian ideas for assessing the probability that the dynamics of one metabolite directly depends on the dynamics of another metabolite. At the core of these methods is a Bayesian network , which is a graph model that represents a set of nodes (variables) and their conditional probabilistic dependencies upon each other. Its analysis permits the explicit detection of causal associations among nodes in the system, as long as there are no structural or regulatory cycles, for instance, in the form of material recycling or feedback signals. While the latter exclusions are clearly restrictive, Sachs and coworkers  successfully used Bayesian network methods to investigate the structure of a protein-signaling network from single cell flow cytometry data. The computational methods confirmed and elucidated most of the previously reported causalities and revealed new relationships between the involved signaling proteins. Bayesian network methods have not been applied extensively to metabolic pathway systems, but more often to genomic networks, where the task was to reconstruct networks of expression traits, or networks comprised of both expression and disease traits .
The task of structure identification using model based approaches is very difficult for non-canonical models, because there are infinitely many nonlinear models that would have to be explored, unless some additional rationale could guide the model selection. Even within the comparatively limited area of metabolic modeling, the choices of combinations of rate functions would be daunting . In any case, some types of “basis functions” are required for about any strategies of inferring network structures, as described in the following sections.
As briefly mentioned in the introduction of this section, overly complex models may fit the data very well since increasing the complexity of the model naturally allows more freedom to provide a better fit to the data, for instance, in terms of the sum of squared errors. However, an over-inflated model typically does not perform well when tested on new data. This problem is known as over-fitting. One approach for restricting model complexity and to find the optimal model size is to add a penalty term to the cost function that is minimized. The optimal model can then be determined by finding the one that minimizes the aggregate cost function . The consequent problem of using this approach is how to weigh data fit against model complexity. One approach, coined “simple to general,” calls for starting with the simplest reasonable model and adding one term at a time until a minimal cost function is found (e.g. ). In the opposite direction, the “general to specific” strategy initially includes everything possible in the model and then gradually eliminates terms until the minimum in the cost function is found . Crampin and co-workers [233,234] used these two approaches of model construction to extract kinetic information from time series data. Although their results suggested that the general-to-specific algorithm outperforms the simple-to-general approach, they indicated that when the number of chemical species included in the model is large (~10), the numbers of possible elementary reactions are massive, thus making the computation difficult if not infeasible. Therefore, it is desirable to limit the size of the basic set below a reasonable upper bound. This strategy appears to be reasonable, because genomic, metabolic, and proteomic networks are generally sparsely connected . An idea similar to the simple-to-general strategy was discussed by Marino and Voit , who addressed structure identification tasks beginning with the simplest types of S-system equations (see below).
Model-based methods of structure identification are especially powerful if they are based on time series data. As discussed earlier, parameter estimation from time series data usually requires considerable computational effort, and this effort increases dramatically when the structure of the underlying system is unknown. The challenges of identifying the correct structure and regulation of a system strongly suggest using all preprocessing tools available to limit the analysis to the most likely connections in advance, reduce the search space and identify good initial guesses for the parameters.
For the identification of structure from time series data, BST models seems particularly useful, especially if not much specific information about the mechanistic processes within the biological network is available. The advantages and features of BST representations have been reviewed in Sections 2.3.3 and 4.4 and need no further description here.
In addition to the computational pruning techniques reviewed in Section 4.4, pruning can also be achieved based on biological insight. Almeida and Voit  suggested making maximal use of other a priori biological information that might be available in addition to the time series data. As a specific example, Voit and Savageau  analyzed a yeast fermentation system in several a priori possible variations that corresponded to hypotheses regarding the existence of specific processes and regulatory signals and studied the improvement in error with statistical methods.
In a more generic fashion of “inverse pruning,” and pursuing the “specific to general” strategy, Marino et al.  proposed an algorithm based on reconstructing S-system equations in a gradual progression. Using the decoupling technique and focusing on one differential equation at a time, they began with equations with constant input and simple substrate driven degradation, obtained the best possible fit, and then gradually added other variables to the equation, always checking their statistical significance. Thus starting from the minimal (and most parsimonious) model, choosing a modest connectivity index, and increasing the number of variables step by step, until a maximally allowed level of connectivity was reached, they identified small pathway systems rather efficiently.
Daisuke and Horton  also utilized the “scale-free” property of networks [235,236] to restrict the connectivity in biological systems during optimization procedure. Their results showed that the restriction increased the conversion ratio while reducing the average number of generations and reducing both false positive and false negative estimations of links in the network. Zuñiga et al.  recently proposed applying ant colony optimization (ACO) to the network inference problem using the S-system formalism. Their preliminary results showed that, starting with a fully connected network, ACO was able to recover the connectivity of the network.
Kimura et al.  proposed a function approximation approach outside BST to infer reduced Normalized Gaussian network (NGnet) models of genetic networks from time-series data. Their results showed that the method successfully inferred the genetic network structure from artificial datasets with high specificity and sensitivity. The method was also tested on random genetic networks and actual biological data. The computational time of this method was shown to be much shorter than for other inference methods.
As described in Sections 4 and 5, many methods have been developed recently that attempt to solve parameter estimation and structure identification problems through inverse modeling using the BST formalism. Most of these methods were developed to address the main problem of optimizing parameter values against observed time series data using gradient based methods, regression algorithms, or evolutionary approaches. Other methods were proposed as support algorithms, for instance, methods for avoiding the time consuming integration of differential equations, smoothing noisy data and estimating slopes, restricting the parameter search space, excluding unlikely connections within the network, or reducing the number of parameters to be estimated.
Many of the published papers used a combination of several methods to solve the inverse problem. For instance, they may have applied decoupling techniques along with various optimization algorithms, tried to reduce the number of parameters before estimating their values, or included several objective functions to constrain the solution space.
In spite of the considerable number of methods that have been proposed for inverse modeling using BST models, each method has its pros and cons and there always seem to be conditions and situations where one method works well and the other not so. In the end, there is currently no algorithm that is perfect, or even sufficiently effective, for the majority of realistic cases. Granted, each proposed algorithm showed effective solutions for particular cases and superiority in comparison to other methods. Nevertheless, most algorithms were only tested against synthetic time series data with respect to robustness and algorithmic efficacy, and it is known well that many such results fail when the same methods are applied to real experimental data. Furthermore, different authors used different benchmarking systems, so that it is hard to tell from the published results which algorithms are superior to the others, and under what conditions
The difficulty of obtaining fair comparisons is a result of the following five issues. First, as said before, different biological systems were used to demonstrate the usefulness of the algorithms. It is clear that different systems generate synthetic time series with distinctly different properties. For instance, they strongly affect the features of the data matrices that are the basis for subsequent computation. Depending on the specific data, the matrix may be ill-conditioned or exhibit collinearities between rows or columns, and this algebraic consequence has a direct effect on the efficacy, correctness and reliability of the tested algorithms. As a result, it is difficult to compare the algorithmic methods and simultaneously to discern how they are influenced by the features of the test system. Second, the numbers of time series points included for computation are different or unstated. Thus, the effects of data point inclusion or missing data on the algorithm are unclear, but they can affect the information criteria and the fitness score of the method. Third, the objective functions selected for the optimization vary and thereby prevent direct comparisons among algorithms. Fourth, the constraints on the parameter values are often different. This seemingly minor issue makes it difficult to tell if the algorithm converges since the boundaries are relatively close to the true optimum or because of the efficiency of the algorithm. Fifth, in addition to testing the methods using noise-free data, artificial “measurement errors” are introduced to examine if the algorithms can still find the correct parameter values. However, the way and extent of adding noise and the methods used for data smoothing often differ, which renders fair comparisons difficult.
A cursory comparison of parameter estimation algorithms in biochemical pathways has been published, but only two networks were considered and neither of them was implemented for BST applications .
To address the problems indicated above, del Rosario and co-workers  recently proposed a project called MADMan (Munich-Atlanta-Diliman-Manila), which aims to compare the published parameter estimation algorithms using BST formalisms in a systematically way, including the testing of the algorithms with the same variety of networks, uniform benchmarking bases, and standardized evaluation criteria. The goal of the benchmarking framework is to develop a strategy for choosing a set of candidate algorithms given a biochemical or gene regulatory network and experimental data. MADMan is an ongoing project that constitutes a huge task, which will require substantial effort and cooperation among the participating groups.
The direct comparison of various optimization algorithms and different data and settings will ultimately be the least biased strategy to determine which algorithms are better than the others. One must be aware though that speed (or lack) of convergence and unsatisfactory performance in terms of fitness are merely some of the issues that need to be analyzed for each optimization algorithm or computational software. Other features may contribute to the problem and its solutions as well, such as data related issues, model related issues, and mathematical issues, as reviewed in Section 3.4.
While the MADMan project will attempt to clarify the applicability of methods under a wide range of conditions, we propose in this section a streamlined “work-flow” strategy for estimating parameter values in BST models with currently available methods. The work-flow diagram consists of a decision process based on potential problems that are encountered with some regularity. These include issues related to the time series data, model of choice, computational efficiency, and mathematical redundancy during the inverse modeling process. The work-flow also suggests relevant diagnostic tools or corresponding solutions. It addresses the main optimization algorithms as well as other supporting methods and diagnostic techniques, along with some assumptions and educated guesses that are required to estimate all parameter values of a system of realistic size.
The ultimate goal of inverse modeling is to find a mathematical model that describes the biological phenomenon and predicts situations that had not been used for model identification or data fitting with correctness, robustness, and efficiency. These standards may not always be fulfilled simultaneously, thereby requiring compromises. For instance, algorithms that find the optimal solution may be expensive in terms of computational time, whereas some of the fast algorithms may only be able to find coarse solutions.
The judgment on algorithms is comparatively easy when synthetic time series are used for testing, since the criterion of “correctness” is then automatically given and easy to assess by checking the fitness score, testing the validity of the inferred network structure, and comparing the estimates with the true model parameters. However, in reality, the “correct” model is not known and goodness of fit cannot always guarantee the reliability and applicability of the model. For instance, the model with the smallest residual error might a priori be deemed mathematically the “best” model. However, a small error does not necessarily imply that the model is the best choice for describing the biological system. In many actual cases, the “best model” tends to have over-fitting problems (see Section 5) and may not be able to extrapolate toward untested conditions when no extra constraints are introduced. Furthermore, a solution that fits the observed time series quite well may not necessarily be unique. Other solutions may exist, with distinctly different parameters, and with fits of a similar quality. In fact, instead of aiming to find “the one best” model, one might set the goal of every inverse modeling strategy as the task of determining all models that are consistent with the data within some acceptable error. The resulting candidate set of parameters may be clustered tightly or scattered throughout the search space. In either case, the diversity of possible solutions is helpful for exploring potential model structures and proposing possible causal relationships among the network components. Furthermore, the candidate models can be assessed with respect to stability, sensitivity, gains, or other features that might shed light on the models and the investigated biological system [1,35]. Comparative simulations with the candidate models may identify one or the other model as more likely or suggest hypotheses or critical experiments that ultimately reveal to true composition of the system at hand.
The proposed flow diagram for inverse modeling is shown in Figure 6. The global time series data are entered into a matrix, which is then screened and preprocessed with diagnostic and corrective tools (Step ➀). For instance, if the variable traces have similarly shaped dynamics, they may be (approximately) collinear with each other. The calculation of the condition number or correlation coefficient can point to possible collinearities in the data matrix (Step ➁). If the time traces are collinear, one may remove the model redundancy by pooling collinear variables or ignoring a subset of them, or merge constant variables with the rate constant (Step ➂). If there is no collinearity, a symbolic mathematical model of the system can be derived based on the model of choice, without numerical specification of parameter values (Step ➃). It has been shown that S-system and GMA representations in BST are good candidates for this propose. After setting up the full model, it is advisable to search for possible simplifications. For instance, if the network topology is known, a reduced symbolic model can be formulated with some parameters set to zero, in accordance with the network diagram (Step ➄). If the system contains ubiquitous metabolites such as ATP, which are involved in dozens of reactions, partial modeling techniques may be applied. This technique retains these variables as input to the model, rather than explicitly modeling them (see Section 3.4.1 for details). The “off-line” nature of the ubiquitous variables further reduces the complexity of the symbolic model. Since fast optimization is desirable for the initial stage, even if it is coarse, it is beneficial to employ the decoupling technique, which converts the differential into algebraic equations. The decoupling step involves the measurement of slopes, which may be assessed directly or upon smoothing (Step ➅). Once the symbolic model is decoupled, the parameters of each equation can be estimated by some fast optimization algorithm (Step ➆). Alternating regression (AR) was shown to be one of the algorithms that work quite well for many S-system models. If AR converges, the resulting model is ready for further analysis and evaluation. If the initial guesses lead to lack of convergence, the algorithm is restarted with a different set of initial guesses. Another option for this step may be a collocation method. If the system topology is not known or only partially known, algorithms or techniques for inferring the network connectivity are applied. These include prior linearization of the system dynamics or sorting of parameter combinations by their empirical likelihood of inclusion in an equation (Step ➇; see Section 5 for detail). If the network topology is not known, it is also necessary to choose an optimization algorithm that does not depend critically on topological information; a suitable algorithm for BST models is eigenvector optimization (EO) with prior decoupling (Steps ➈ and ➉). Algorithms permitting ill-characterized system topologies are usually combined with pruning methods that eliminate unlikely connections between network components and reduce the number of parameters to be estimated during the process of estimation. If the fast algorithms are not able to yield acceptable fittings, some other, computationally more expensive algorithms such as genetic algorithms or evolutionary approaches are applied (Step ). If the outcome of the initial fitting is not acceptable, the optimized parameter values may be used as start values for subsequent refining algorithms. These are typically more costly and may lead to better solutions, although they are not necessarily always effective. A significant consequence and advantage of the combination of approaches is that the result often consists of multiple parameter sets that are all consistent with the data and that can lead to new, testable hypotheses and may offer guidance for further theoretical and experimental investigation (Step ). It may be possible that the algorithms are not even able to produce acceptable fits (Step ). We will discuss this situation in detail in Section 7. Once the initial models are obtained, the next step is dedicated to model analysis, including model diagnostic and cross validation as described in Section 1. If the models are deemed reliable and appropriate for the purposes of the modeling effort, they can be used for applications and for gaining a deeper understanding of the biological phenomenon; specifically, they can be used to make predictions, generate new hypotheses, or guide the design of additional biological experiments (Step ). In contrast, a model analysis that indicates lack of robustness or discrepancies between model and observations reveals potentially fundamental problems in the model. In this not so rare situation one needs to return to earlier steps of the modeling process and refine the model in an iterative manner. For instance, the modeler might need to discuss with the expert biologists how to obtain additional information or identify possible mistakes in the assumed model structure or missing reactions or signals in the pathways. It may also be useful to resample the data with jackknife or bootstrap methods and to redo the analysis in order to explore possible alternative solutions.
Most of the steps of the inverse modeling “work-flow” can be automated. For instance, it is relatively easy to check for collinearities between time series once the data matrix is ready (Step ➁). The full symbolic model of the system can be derived directly and per computer if the number of variables in the model is known ; Step ➃). The slopes of the time series can be estimated directly or upon smoothing, using various algorithms (Steps ➅ and ➈). The actual parameter optimization is possible with many algorithms based on the time series and the structure of the model (Steps ➆, ➉, and ). The network topology can be inferred giving the data matrix and fully connected model (Step ➇). These more or less automatic steps can be worked into a data pipeline, at least in principle. Other steps are not as straightforward and thus require work manual intervention. For instance, even though methods of matrix diagnostics in Step ➁ point to collinearities in the data matrix, pooling variables or reducing the redundancy in the model in Step ➂ requires some thought regarding the location of the affected variables and their relationships in the pathway. Further model reduction includes the decision of which variables to model explicitly, in cases where the model contains highly connected metabolites. These “intelligent” steps are required to reduce the network topology and the corresponding symbolic model (Step ➄). Thus, the entire process is not yet fully automated and may always require human supervision. Nonetheless, some tools like Best-Kit , Cadlive , BSTBox , and BioinformaticStation  are beginning to provide interfaces that facilitate some of the steps involved in metabolic modeling.
Once several candidate solutions are obtained, the immediate question is whether there are reasonable guidelines for choosing between them. Several scenarios are to be anticipated. If many well-fitting solutions, either found with different optimization methods or obtained using some re-sampling scheme, are clustered closely within the parameter space, the solutions are parametrically similar and the networks they represent are essentially the same or very similar in structure. In contrast, if the optimization yields distinctly different solutions, exhibiting essentially the same residual error, it is a priori difficult to decide which model is best. Recent results in one of our estimation studies showed that a single data set allowed multiple distinctly different numerical solutions, especially if constraints on kinetic orders were set loosely. This was not entirely surprising because even one-variable S-systems are flexible enough to permit different parameter sets generating very similar graphs (cf. [204,242–247]). Without additional information, each such parameter set is an equally valid solution since it fits the data essentially equally well. However, problems may arise if a “wrong” solution is used for extrapolation to new conditions, as we will discuss in Section 7. By basing the estimation on many data sets and experimentally testing the same pathway under different conditions, the problem can often be alleviated, because the use of several data sets clearly constrains the flexibility of the underlying model considerably.
In many cases, one will obtain several alternative solutions. In other cases, the opposite may be true: in spite of the many options outlined before, it is still possible that even a combination strategy cannot find an acceptable fit. Potential reasons and suggested future work are discussed in Section 7.
As mentioned in Section 3.4 and in the previous section, the challenges of inverse modeling can be classified into issues related to data, model structure, computation, and mathematical features of the representation. Most of the recent articles have acknowledged and discussed various computational issues in great detail and some have addressed data and model related issues. However, there has been little discussion of model validity and quality beyond residual errors, the conditions under which the models can be obtained, and diagnostic tools for non-convergence or for situations where models cannot even be obtained with any degree of reliability.
These open issues fall into two categories. First, even if the algorithms are able to find a set of candidate models, it is possible that none of these models is acceptable for one reason or another. For instance, it may happen that model diagnosis and simulation studies reveal that none of the models are stable, that they are all overly sensitive, or that they do not exhibit reliable predictive ability. Other problems are lacking model fit for data not used in the estimation and model failure in extrapolations. Second, it is possible that the algorithms are not even able to produce acceptable fits. In these cases, the failure is usually imputed to the computational algorithms themselves. However, the sources of the problem may lie in a combination of the alleged model structure, the particular data sets, and the computational methods, and it is advisable to extend the diagnosis beyond the algorithmic techniques.
Following are some of the issues that should be addressed to improve the validity and reliability of the estimated model, beyond residual errors and computational efficiency.
Finally, one should emphasize the need for obtaining reliable solutions within short periods of time. In some cases, only a single estimation of the system may be needed, and it may be acceptable if this estimation takes several hours. However, once the field moves toward “estimation on the fly,” solutions must be obtained within a few minutes or, preferably, within seconds. The need for fast solutions becomes especially pertinent if biologists and modelers together engage in concept map modeling, which permits the conversion of hypothesized network diagrams into numerical mathematical models . Specifically, based on the known or hypothesized connectivity and regulatory information regarding the investigated system, the biologist designs a concept map consisting of a connectivity diagram of processes comprising the system and including known or assumed regulatory features, and provides semi-quantitative information on stimuli and measured or expected responses of the system. The modeler converts this information through combined methods of forward and inverse estimation into a mathematical construct that can subsequently be used for typical model analyses and to generate and test new hypotheses. This conversion step, which includes parameter estimation from time series, needs to proceed fast in order to permit interactive work, in which the modeler runs simulations with the model and the biologist-modeler team collaboratively interprets the results and devises improved concept maps. Because this method heavily depends on the biologist’s initial intuition and hypotheses, many iterations between hypothesis formulation and diagram-to-model conversion are needed, thus demanding fast solutions that might not be absolutely precise but allow the interactive exploration of complex biological systems.
Whether bottom-up, top-down or concept map modeling is the method of choice, we hope to have conveyed that the estimation of model parameters and the identification of structure and regulation of ill-characterized biological systems is a vibrant field that will continue to offer challenges to teams of biologists, mathematicians, computer scientists and modelers throughout the foreseeable future.
The authors are grateful to Dr. Siren Veflingstad and two anonymous reviewers for critically reading the manuscript and providing constructive suggestions. This work was supported in part by a National Heart, Lung and Blood Institute Proteomics Initiative (Contract N01-HV-28181; D. Knapp, PI), a Molecular and Cellular Biosciences Grant (MCB-0517135; E.O. Voit, PI) from the National Science Foundation, a grant from the National Institutes of Health (R01 GM063265; Y.A. Hannun, PI), and an endowment from the Georgia Research Alliance. The work was also in part funded by the BioEnergy Science Center (BESC), which is a U.S. Department of Energy Bioenergy Research Center supported by the Office of Biological and Environmental Research in the DOE Office of Science. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the sponsoring institutions.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.