BACKGROUND
a systems level understanding of highly complex biological systems requires an integration of experimental techniques and computational research  <cit> .
current molecular biology techniques can generate high-throughput quantitative data that support in silico research using mathematical models  <cit> .
these models can be used to simulate and study the dynamic interactions among the components of cellular systems as well as the systems' responses to external perturbations and signals.
such tools offer enormous potential for understanding cellular functions at the organism level  <cit> .
mathematical models also serve as test beds for generating hypotheses and designing experiments to test them  <cit> .
furthermore, they provide bases for model-based product and process design applications.
useful insights and predictions have been obtained for several biological systems from computational modeling and analysis.
a few examples include the metabolic network analysis of escherichia coli growth on glucose and acetate  <cit> , the map kinase signaling pathways  <cit>  and caspase function in apoptosis  <cit> , and bifurcation analysis of cell cycle in saccharomyces cerevesiae  <cit> .
an iterative process for model development and the testing of hypotheses has been proposed by many researchers in the field and was recently highlighted by kitano  <cit> .
a qualitative approach of this process is described in  <cit> .
in addition, rabitz and co-workers  <cit>  have recently developed an iterative method for closed loop parameter identification in biochemical reaction networks.
a global inversion algorithm was used to extract the parameter estimates that minimize the differences between model prediction and experimental data.
unfortunately, global search methods typically have high computational requirements, and thus, do not scale very well with the system size.
in this work, a quantitative model identification is developed that effciently obtains parameter estimates and facilitates scalability to very large network sizes.
a proposed strategy is described in section  <dig>  with an emphasis on the modeling element.
the modeling strategy is decomposed into three main steps:  determining the connectivity of the biological network and the interactions of the sub-components,  formulating the kinetics of inter-conversion among the subcomponents, and  estimating the parameters in the rate equations.
to the authors' knowledge, this work represents the first documented example of multiple iterations for model refinement using such a framework in systems biology.
the parameter estimation from experimental data remains the bottleneck in the model development  <cit> .
banga and coworkers  <cit>  have compared several advanced deterministic and stochastic global optimization methods for parameter identification from available experimental data.
it was observed that the traditional gradient-based optimization methods often failed to arrive at the global optimal solutions.
deterministic methods  <cit>  can achieve global optimality for certain classes of problems, but there is no guarantee of convergence in finite time  <cit> .
stochastic strategies  <cit>  can locate the parameter region containing the global solution with relatively better effciency, but global optimality is not guaranteed.
furthermore, both methods suffer from the large computational burden required, even for moderately sized problems.
moreover, the validity of model with the estimated parameters over the entire operating space remains to be determined.
parameter identifiability tests should be performed prior to the estimation process to ensure that the parameter estimation problem is well-posed.
further, the identifiability tests assist in selection of optimal measurements.
several researchers  <cit>  have developed methods to determine whether a parameter is "identifiable a priori", i.e., identifiable from a given experiment design using the available measurements.
a similar concept known as "practical identifiability" is concerned with the achievable accuracy of the parameter estimates.
the confidence interval for the model parameters are determined using the fisher information matrix   <cit> .
doyle and coworkers  <cit>  have performed model identifiability studies for a gene regulatory network using gene expression data, in which the identifiability of the parameters was found to be strongly dependent on the driving function.
the final step in the iterative model development process is the design of "optimal experiments" that would provide rich experimental data for improving the parameter estimates.
experiments can also be designed for discrimination among competing model structures that translates to selection between multiple proposed mechanisms of cellular function.
asprey and macchietto  <cit>  have developed a strategy of optimal experiment design for model structural identifiability.
the strategy was used to identify the kinetics of the reactions in the fermentation of saccharomyces cerevesiae.
kremling and co-workers  <cit>  propose several strategies for model discrimination to identify the correct reaction mechanism of a test metabolic network.
banga and coworkers  <cit>  have formulated the optimal design problem, using a scalar function of the fisher information matrix  as the performance index, for parameter estimation of nonlinear dynamic systems.
in this work, an iterative procedure for model identification is proposed and applied to the caspase-dependent apoptosis system.
an optimal measurement set is determined using the fisher information matrix .
the parameter estimation from partial measurements is decoupled into two parts.
first, the available measurements are used to estimate the profiles of all unmeasured concentrations and reaction rates using a state regulator problem  formulation.
in the second part the concentration and rate estimates are used to determine the model parameter values.
the srp formulation in this work is an extension of the dynamic flux balance analysis  approach developed by doyle and coworkers  <cit> .
finally, the model-based experiment design uses parameter identifiability and d-optimality criteria to obtain the optimal experimental procedure that would generate the most informative data for model refinement in the next iteration.
RESULTS
the iterative scheme for model identification is shown in figure  <dig>  the optimal set of measurements is determined a priori.
for an effcient model identification, a significant fraction of the unknown model parameters should be identifiable.
in the case of poor identifiability, a higher number of measurements would be motivated.
also, the model complexity could be reduced to decrease the number of parameters; but this does not guarantee identifiability of the reduced number of parameters.
alternatively, a richer protocol  might yield improved identifiability.
in this work, selection of the a priori optimal measurement set is restricted to the "preliminary" experiment design that may be suboptimal.
the model connectivity and reaction mechanisms are developed from existing biological knowledge and are assumed to be known.
the network connectivity, along with the partial measurements , is used in the state regulator problem  to obtain estimates of all system unknowns .
here it is important to note that the kinetics of the reaction rates are not used in the srp algorithm.
next, the full estimates of the concentrations and the reaction rates are used for estimating the parameters in the kinetic model.
this decouples the model identification into two parts such that the parameters involved in the kinetic equation of each reaction are independently determined as opposed to simultaneously estimating full model parameters from limited measurements.
next, the model invalidation test, which is a critical step in model development and the last "quality control" step before the desired application  <cit> , is performed.
invalidity of the model could be determined by comparing model predictions with experimental data that is not used in the srp algorithm.
further, model invalidity can occur if the model predictions conflict with documented biological knowledge of the system.
in case of an invalid model, the model parameters are refined in subsequent iterations using the information obtained from the optimal experiment or by expanding the measurement set.
the process of model identification is repeated in an iterative manner until an "acceptable" model is obtained.
system
the mathematical model considered in this paper has the following structure:
= ax + br + c
where
r = f,
x represents the protein/metabolite/gene concentrations, r the reaction rates, and p the model parameters.
this is a very general nonlinear state space model in the variable x. the matrices a and c describe degradation and auto-generation respectively, whereas the matrix b represents the stoichiometry of the biological network.
the kinetics among the proteins/metabolites/genes interactions show up in the reaction rates r. the aforementioned model is a continuous time invariant affine system from which a discrete version can be derived by standard techniques using a zero-order hold  <cit> .
the resulting discrete model equation is represented as:
where
the discrete version of the model is used in the srp estimation algorithm.
theory
step 1: determination of measurement set
the optimal measurement set consists of species whose concentration measurements would have maximum benefit for model identification, e.g, parameter identifiability and accuracy.
in this work, the measurement set is determined such that the model parameters can be estimated accurately.
the assessment of parameter identifiability in a model is crucial prior to parameter estimation from experimental data  <cit> .
identifiability is closely linked with parametric sensitivity analysis through the fisher information matrix   <cit> .
the unidentifiable parameters are determined using the orthogonal procedure proposed by macauley and coworkers  <cit> .
here, a scaled sensitivity coeffcient matrix () shown below is computed:
where {p <dig>  …, pk } are the model parameters, {η <dig>  …, ηm } are the response variables which include all possible measurable quantities, {t <dig>  t <dig>  …, tn} are the sampling times for the measurements, and  is the "initial" parameter value that is either the guess values of the parameter or the value obtained from literature.
the orthogonal method is a geometric based approach where the number of identifiable parameters correlates with the rank of the orthogonalization of the scaled sensitivity matrix.
the parameters corresponding to the columns of orthogonalized sensitivity matrix are deemed unindentifiable if the norms are smaller than a given tolerance.
the details of the orthogonal method are not included here for the sake of brevity.
the next step is to obtain a measurement set that maximizes the expected accuracy in the identifiable parameters .
the fisher information matrix  along with the cramer-rao theorem are used to determine a measurement set such that the estimated parameters have minimum variance.
a detailed description of the procedure and its theoretical foundation can be found in  <cit> .
assuming that the measurement errors are additive and gaussian, the fim is given by  <cit> :
fim = jtwj
where w is the inverse of the measurement error covariance matrix and j denotes the sensitivity coefficient matrix for the measured response variables:
the quantity {p <dig>  …, pr} denotes the identifiable parameter vector and {, …, } denotes the measured response variable vector.
the cramer-rao inequality establishes a lower bound on the variance of the identifiable parameters given by:
σ <dig> ≥ ii
the 95% confidence interval  for a parameter is given by:
ci =  ±  <dig> σ
in equation  the lower bound of the variance is used.
symmetry of the confidence region about the nominal value is assumed.
this results in the following definition of the percentage deviation from the nominal value:
the optimal measurement set is chosen such that the sum of the percentage error  for all the identifiable parameters is minimized.
in this work, the optimal set is determined by a brute-force search over all combinations of measurement sets subject to restrictions that may be imposed by the system.
doyle and co-workers have developed effcient rational algorithms to determine the optimal measurement set with minimum computational burden  <cit> .
the confidence intervals for non-identifiable parameters are infinitely large and hence are eliminated from the analysis.
identifiability for these parameters can be obtained only by a change in the experimental design or by the selection of an alternative model structure.
step 2: state estimation algorithm
generally, it is not possible to measure all time-varying components in a metabolic or signaling network.
however, there are several techniques from systems engineering to estimate the behavior of unmeasured components given partial measurements of other system constituents.
bastin and dochain have used an adaptive nonlinear observer for estimation of specific growth rate and biomass concentration  <cit> .
given accurate models, extended kalman filters  have had success in several biological applications  <cit> .
artificial neural networks  have also found applications where dynamic models are not available  <cit> .
in this work, an extension of dynamic flux balance analysis   <cit>  is developed to estimate unmeasured concentration and reaction rate trajectories given partial measurement sets.
the premise of this approach is straightforward: cellular processes have evolved regulatory structures that optimally use cellular resources.
this premise translates into two postulates;  network flows are managed to minimize internal accumulation and  networks are managed to minimize the number of edges carrying flux at any given time.
these two requirements are analogous to a classic problem in automatic control, namely, the state regulator problem .
the srp based estimator uses the measurement set selected from step  <dig> to estimate unknown concentration and reaction rate trajectories via a constrained convex programming problem.
the srp estimator constrained by the key measurements captures the optimal cellular behavior of the system.
estimates of the reaction rates at time step k and protein/gene/metabolite concentrations at time step k +  <dig> are determined by the srp and the discrete mass balance equations.
the srp must be solved at each sampling interval to obtain estimates of the unknown rates and concentrations.
consider the model ; the discrete mass balance equations over a p step horizon are given by:
where the matrices k+ <dig>  k, x and c are defined as
and the matrix r is given by
the srp estimator is a quadratic program  with two cost terms, the cost of intermediate accumulation and the cost of operating a network reaction.
the srp penalizes intracellular metabolite or protein accumulation, but does not explicitly forbid it.
moreover, because reactions introduce an additional cost, the srp only utilizes those reactions required to satisfy the mass balances and thermodynamic constraints.
formally, the estimation problem is given by:
subject to:
k +  <dig> ≥  <dig>
αr  ≤ k ≤ βr
the srp problem is subject to non-negativity constraints , flux-directionality constraints  and constraints imposed by the measurement set.
specifically, constraint of equation  <dig> forces state estimates belonging to the measurement set to equal the corresponding measured value.
the quantities  denote measurements that may have been corrupted with noise.
specifies the tolerance around the measurement within which the estimate is constrained to lie .
the term Ξx defines the measurement set.
the matrix wx defines the cost of intermediate accumulation whereas the matrix wr represents the reaction cost:
for the caspase system considered in this work, the wx and wr are taken to be the order of magnitude of the inverse of the maximum value of the corresponding state or rate.
this requires only approximate information regarding ranges of the protein concentrations and identification of the slow versus the fast reactions:
step 3: parameter estimation
the estimates of the concentration profiles and the reaction rates allow efficient determination of the parameter values by decoupling the full parameter estimation into multiple sets.
each set consists of parameters associated with one reaction rate.
the parameters are obtained by minimizing the difference between the estimates of each reaction rate and that predicted by the kinetics r , which is a function of the concentrations.
in case of the first iteration the minimization follows:
where ri are the individual reaction rates and nr is the total number of reactions.
the kinetic parameters associated with a reaction rate equation are determined independently from those with other reactions, i.e., the parameter estimation is decoupled with respect to each reaction.
for subsequent iterations, the bayesian estimation formulation in  <cit>  is used.
in this formulation, in addition to the difference between the estimates of the reaction rates and the model predictions, the deviations of parameter values from those obtained after the previous iteration are minimized.
the formulation can be represented as:
in the above equations,  is the estimate of the ith reaction rate and  is the estimates of the concentrations obtained from the srp algorithm, ri is the predicted rate of the ith  reaction from the kinetics in equation  <dig>  p are the parameters associated with the ith  reaction rate, p <dig> are the parameter values obtained from the previous iteration, and vε and vp are the variances of the estimates of the reaction rates and the prior parameter estimates.
the parameter variances are determined using the fisher information matrix .
the variances of the non-identifiable parameters are infinite and penalty for deviations for these parameters are not considered in equation  <dig>  the variance for the estimates of the reaction rates can be determined from the expected noise in the measurements from which the estimates are obtained.
step 4: model invalidation tests
given the iterative nature of this framework, a termination criterion must be established.
poolla et al.
<cit>  have shown that for certain experimental data, it is not possible to confirm whether the model is really valid; however, one can conclude whether the model is not contradicted by the given data.
model validation tests are usually based on the difference between the simulated and measured output and some statistics about these differences.
typical statistics for the model errors include maximum absolute value, mean value and variance  <cit> .
in this work, model invalidity is tested by determining the model prediction errors using the estimated parameters.
this error is calculated as:
to implement this test, experimental data that was not used in the srp algorithm is required.
the statistic used is the maximum and mean value of the errors for the measured states.
when the prediction errors are below a certain desired value, the iterations are terminated.
step 5: model-based optimal experiment design
the optimal experiment design determines the optimal experiment to be performed for the next iteration such that there is maximum information content in the measurements.
this would maximize the accuracy of the estimated parameters.
the model-based optimal experiment design uses the fisher information matrix as a measure of the amount of information contained in a given set of measurements about the model parameters  <cit> .
the optimization searches through the space of experimental conditions or some parameterizations of the experimental protocol.
for example, an optimal ligand input can be parameterized into a time series profile such that the optimization variables are the levels of ligand at different times .
naturally, the optimization will be restricted by the limitations in the experimental conditions and apparatus.
there exist several fim-based optimality measures that quantify the overall informativeness of the measurements  <cit> .
among these, parameter identifiability and the d-optimality are the most widely used measures.
for accurate model identification it is critical that maximum number of the parameters be estimated accurately.
thus, maximizing the number of identifiable parameters is the primary criterion proposed for determining the next experimental design.
the orthogonal procedure proposed by mcauley and co-workers  <cit>  is used to determine the number of identifiable parameters.
there can be multiple experimental designs with the same maximum number of identifiable parameters.
the selection among these is done so as to maximize the informativeness of measurement data.
for this purpose, the d-optimality criteria is proposed.
the use of d-optimality translates to minimizing the confidence interval of all the identifiable parameter estimates.
the optimal experiment design criterion is shown as follows:
where e denote the feasible experimental conditions , fim is given in equation  <dig> and r denotes the number of identifiable parameters.
thus, the identifiability is maximized in the sense that the hyper-dimensional confidence interval is minimized.
for experimental designs with the maximum number of identifiable parameters, the one with the highest determinant of the fisher information matrix is selected as the optimal design for the next experiment.
a case study
the proposed iterative model identification is applied to the function of caspase- <dig> and caspase- <dig> in apoptosis.
caspase enzymes are at the core of the cell's suicide machinery.
these enzymes are activated either by an external signal or by stress, and activated enzymes will then dismantle the cells.
varner and co-workers have developed a model for the caspase function in apoptosis  <cit> .
the model describes the key elements of receptor-mediated and stress-induced caspase activations.
the model consists of  <dig> states  and  <dig> reactions with  <dig> parameters .
the in silico experiments in this study use the varner model as the "actual" system.
measurements are assumed to be obtained from this "actual" system corrupted with up to 10% noise.
the iteration starts using a "initial" parameter set, generated by perturbing the parameter values of the varner model  by 70–100%.
the external and stress signals that activate the caspase system are considered as the manipulated variables in step  <dig>  model refinement is performed either by determination of an optimal experiment or by optimal refinement of the measurement set.
the performance of each of these criteria for model refinement requires to be tested.
therefore, in the first iteration, a "preliminary" suboptimal experiment with a suboptimal measurement set is considered.
moreover, in most cases of model identification, it is expected that preliminary experimental data is available.
this data is usually obtained from a suboptimal experiment design and does not include the optimal set of measurements.
the second iteration is performed in two ways; one by improving the measurement set and second by improving the experiment design.
this tests the performance of both refinement criteria.
the sequence of events is shown in figure  <dig>  it should be noted that it would be best to use optimal experiment design along with the optimal measurement set.
however, it may not always be possible to do so due to feasibility issues specific to the particular system.
hence, this approach is not considered in this work.
model identification under less constrained conditions using a similar framework is included in  <cit> .
after the first iteration, it is observed that there is a significant improvement in the predictions with the estimated parameters for both the protein concentrations and the reaction rates, as shown in figures  <dig> and  <dig>  the high errors with the "initial" parameters demonstrate that there is no bias in the results based on the starting guess values for the parameters and that there is indeed an improvement in prediction of both the protein concentrations and the reaction rates.
however, the improvement is not suffcient as observed from the invalidation test .
this warrants a second iteration.
in general, the model predictions improve with the second iteration, as shown again in figures  <dig> and  <dig>  however, it is observed that the predictions are better for the case with the optimal experiment design, in spite of a suboptimal measurement set.
this is due to the fact that model performance depends strongly on the accuracy of the estimated parameters.
using the suboptimal experiment, only  <dig> of the  <dig> parameters were identifiable.
the optimal measurement set simply improved the confidence in these  <dig> parameters.
on the other hand, the optimal experiment increased parameter identifiability to  <dig> parameters even with the suboptimal measurement set.
these result indicates that performance of model identification is strongly linked with parameter identifiability.
discussion
in the proposed algorithm, the network topology and the mechanism of interactions in the pathway are assumed to be known.
several approaches have been proposed in literature to determine the network connectivity from experimental data  <cit> .
in the case of unknown connectivity, these approaches should be used prior to the proposed model identification.
the mismatch between the model and the actual network can appear in two different aspects of the algorithm.
first, the srp step can fail because there exist no feasible state and flux estimates that satisfy the measurement constraints.
such scenario arises mainly due to the network topology mismatch, and has low probability to occur due to the large degree of freedoms in a typical biological system.
alternatively, a well-designed model validation step catches the mismatch between the model and the real system, e.g., an independent measurement set contradicts the model prediction.
also, if multiple models are proposed for a particular biological process, model discrimination methods  <cit>  can be used to identify the correct model structure.
the effect of incorrect connectivity and/or mechanism would depend on the degree of mismatch and is case dependent.
the fact that cellular processes are carried out in an optimal manner lends tremendous promise to the success of this approach.
in case the assumed optimality does not represent the in vivo behavior, the estimates from srp may be inaccurate.
however, the measurements  can attenuate this problem by restricting part of the estimates to match the observations.
the parameter estimates may also be inaccurate if the real system deviates considerably from the assumed optimal behavior and this deviation is not captured by the measurements.
the model identification framework has applicability to all systems that could be represented in the form shown in equations  <dig> and  <dig>  all biological processes are complex, interconnected networks.
a feature common to these processes is that they have a fixed connectivity.
the proposed algorithm for model development could be applied to metabolic networks, signaling processes and gene networks.
the computation burden for solving the srp is minimal as it involves only a quadratic programming problem.
the parameter estimation is also not computationally intensive due to the decoupling facilitated by srp.
a global optimization algorithm can also be used in the parameter estimation instead of the gradient search method to avoid convergence to local minima.
however, the computation burden of parameter estimation will increase.
to avoid local minima and high computation cost, the first few iterates  can utilize a global optimization method, while the remaining iterations can implement a gradient search algorithm.
due to the iterative nature of the approach, errors in parameter estimates can be tolerated as the corrections will be made in the next iteration with the optimal experiment.
the optimal measurement selection is performed by a brute force search in this work.
for very large systems the computation burden for this process grows exponentially.
computationally efficient algorithms for optimal measurement selection  <cit>  can be used.
efficient measurement selection algorithms and the decoupling of parameter estimation for individual reaction rates into separate optimization problems result in good scalability properties of the proposed algorithm for large scale systems.
one limitation of the approach is that it is dependent on the weights in equation  <dig> for the minimization of the cellular resources.
the choice of weights used in this work has provided accurate results but it requires information of the order of magnitude of the concentrations and rates of the system under study  <cit> .
efficient schemes for determining the weights for metabolic networks has been developed by varner and co-workers  <cit> .
CONCLUSIONS
an iterative methodology for model identification from experimental data is developed in this paper.
identifiability tests are performed for an optimal measurement set selection for a given experimental design.
the optimal measurements represent maximum information such that the model identification process is maximally benefitted.
the model identification process is decoupled into two parts.
in the first, the measurements are used to estimate all the unmeasured quantities of the system.
this is achieved using the state regulator problem  formulation which is based on the assumption that the cell is an optimal strategist and uses its resources in an optimal manner.
the srp algorithm developed in this work has shown promising results.
the average errors in the estimates for a significant fraction of the unmeasured responses is less than 10%.
the accuracy of the estimates obtained by the srp decreases with decrease in the information content due to suboptimal measurement set.
in the second part, the full state and rate estimates are used to determine the model parameters.
the decoupling relaxes considerable computation burden compared to estimating all model parameters simultaneously from the limited measurements.
in the final step, a model-based experiment design determines the optimal experimental procedure that generates the most informative measurements for the next iteration.
a strong dependence is observed between parameter identifiability and model performance.
thus, it is critical that the experiment design and measurement set be chosen such that maximum number of parameters are identifiable.
tools developed for quantitative analysis of the dynamics of cellular pathways have tremendous potential in improving the predictive capabilities of biological systems especially in cases where experimental data is available but the kinetic parameters of the pathway reactions are unknown.
the model developing tools are used for a host of applications and systems analysis.
the measurement selection algorithm presented in this work is freely available as part of a model analysis and development toolkit, biosens  <cit> .

 a systems level understanding of highly complex biological systems requires an integration of experimental techniques and computational research. current molecular biology techniques can generate high-throughput quantitative data that support in silico research using mathematical models. these models can be used to simulate and study the dynamic interactions among the components of cellular systems as well as the systems' responses to external perturbations and signals. such tools offer enormous potential for understanding cellular functions at the organism level. mathematical models also serve as test beds for generating hypotheses and designing experiments to test them. furthermore, they provide bases for model-based product and process design applications. useful insights and predictions have been obtained for several biological systems from computational modeling and analysis. a few examples include the metabolic network analysis of escherichia coli growth on glucose and acetate, the map kinase signaling pathways and caspase function in apoptosis, and bifurcation analysis of cell cycle in saccharomyces cerevesiae. an iterative process for model development and the testing of hypotheses has been proposed by many researchers in the field and was recently highlighted by kitano. a qualitative approach of this process is described in. in addition, rabitz and co-workers have recently developed an iterative method for closed loop parameter identification in biochemical reaction networks. a global inversion algorithm was used to extract the parameter estimates that minimize the differences between model prediction and experimental data. unfortunately, global search methods typically have high computational requirements, and thus, do not scale very well with the system size. in this work, a quantitative model identification is developed that effciently obtains parameter estimates and facilitates scalability to very large network sizes. a proposed strategy is described in section with an emphasis on the modeling element. the modeling strategy is decomposed into three main steps: determining the connectivity of the biological network and the interactions of the sub-components, formulating the kinetics of inter-conversion among the subcomponents, and estimating the parameters in the rate equations. to the authors' knowledge, this work represents the first documented example of multiple iterations for model refinement using such a framework in systems biology. the parameter estimation from experimental data remains the bottleneck in the model development. banga and coworkers have compared several advanced deterministic and stochastic global optimization methods for parameter identification from available experimental data. it was observed that the traditional gradient-based optimization methods often failed to arrive at the global optimal solutions. deterministic methods can achieve global optimality for certain classes of problems, but there is no guarantee of convergence in finite time. stochastic strategies can locate the parameter region containing the global solution with relatively better effciency, but global optimality is not guaranteed. furthermore, both methods suffer from the large computational burden required, even for moderately sized problems. moreover, the validity of model with the estimated parameters over the entire operating space remains to be determined. parameter identifiability tests should be performed prior to the estimation process to ensure that the parameter estimation problem is well-posed. further, the identifiability tests assist in selection of optimal measurements. several researchers have developed methods to determine whether a parameter is "identifiable a priori", i.e., identifiable from a given experiment design using the available measurements. a similar concept known as "practical identifiability" is concerned with the achievable accuracy of the parameter estimates. the confidence interval for the model parameters are determined using the fisher information matrix. doyle and coworkers have performed model identifiability studies for a gene regulatory network using gene expression data, in which the identifiability of the parameters was found to be strongly dependent on the driving function. the final step in the iterative model development process is the design of "optimal experiments" that would provide rich experimental data for improving the parameter estimates. experiments can also be designed for discrimination among competing model structures that translates to selection between multiple proposed mechanisms of cellular function. asprey and macchietto have developed a strategy of optimal experiment design for model structural identifiability. the strategy was used to identify the kinetics of the reactions in the fermentation of saccharomyces cerevesiae. kremling and co-workers propose several strategies for model discrimination to identify the correct reaction mechanism of a test metabolic network. banga and coworkers have formulated the optimal design problem, using a scalar function of the fisher information matrix as the performance index, for parameter estimation of nonlinear dynamic systems. in this work, an iterative procedure for model identification is proposed and applied to the caspase-dependent apoptosis system. an optimal measurement set is determined using the fisher information matrix. the parameter estimation from partial measurements is decoupled into two parts. first, the available measurements are used to estimate the profiles of all unmeasured concentrations and reaction rates using a state regulator problem formulation. in the second part the concentration and rate estimates are used to determine the model parameter values. the srp formulation in this work is an extension of the dynamic flux balance analysis approach developed by doyle and coworkers. finally, the model-based experiment design uses parameter identifiability and d-optimality criteria to obtain the optimal experimental procedure that would generate the most informative data for model refinement in the next iteration.  the iterative scheme for model identification is shown in figure the optimal set of measurements is determined a priori. for an effcient model identification, a significant fraction of the unknown model parameters should be identifiable. in the case of poor identifiability, a higher number of measurements would be motivated. also, the model complexity could be reduced to decrease the number of parameters; but this does not guarantee identifiability of the reduced number of parameters. alternatively, a richer protocol might yield improved identifiability. in this work, selection of the a priori optimal measurement set is restricted to the "preliminary" experiment design that may be suboptimal. the model connectivity and reaction mechanisms are developed from existing biological knowledge and are assumed to be known. the network connectivity, along with the partial measurements, is used in the state regulator problem to obtain estimates of all system unknowns. here it is important to note that the kinetics of the reaction rates are not used in the srp algorithm. next, the full estimates of the concentrations and the reaction rates are used for estimating the parameters in the kinetic model. this decouples the model identification into two parts such that the parameters involved in the kinetic equation of each reaction are independently determined as opposed to simultaneously estimating full model parameters from limited measurements. next, the model invalidation test, which is a critical step in model development and the last "quality control" step before the desired application, is performed. invalidity of the model could be determined by comparing model predictions with experimental data that is not used in the srp algorithm. further, model invalidity can occur if the model predictions conflict with documented biological knowledge of the system. in case of an invalid model, the model parameters are refined in subsequent iterations using the information obtained from the optimal experiment or by expanding the measurement set. the process of model identification is repeated in an iterative manner until an "acceptable" model is obtained. system the mathematical model considered in this paper has the following structure: = ax + br + c     where r = f,     x represents the protein/metabolite/gene concentrations, r the reaction rates, and p the model parameters. this is a very general nonlinear state space model in the variable x. the matrices a and c describe degradation and auto-generation respectively, whereas the matrix b represents the stoichiometry of the biological network. the kinetics among the proteins/metabolites/genes interactions show up in the reaction rates r. the aforementioned model is a continuous time invariant affine system from which a discrete version can be derived by standard techniques using a zero-order hold. the resulting discrete model equation is represented as: where the discrete version of the model is used in the srp estimation algorithm. theory step 1: determination of measurement set the optimal measurement set consists of species whose concentration measurements would have maximum benefit for model identification, e.g, parameter identifiability and accuracy. in this work, the measurement set is determined such that the model parameters can be estimated accurately. the assessment of parameter identifiability in a model is crucial prior to parameter estimation from experimental data. identifiability is closely linked with parametric sensitivity analysis through the fisher information matrix. the unidentifiable parameters are determined using the orthogonal procedure proposed by macauley and coworkers. here, a scaled sensitivity coeffcient matrix () shown below is computed: where {p …, pk } are the model parameters, {η …, ηm } are the response variables which include all possible measurable quantities, {t t …, tn} are the sampling times for the measurements, and is the "initial" parameter value that is either the guess values of the parameter or the value obtained from literature. the orthogonal method is a geometric based approach where the number of identifiable parameters correlates with the rank of the orthogonalization of the scaled sensitivity matrix. the parameters corresponding to the columns of orthogonalized sensitivity matrix are deemed unindentifiable if the norms are smaller than a given tolerance. the details of the orthogonal method are not included here for the sake of brevity. the next step is to obtain a measurement set that maximizes the expected accuracy in the identifiable parameters. the fisher information matrix along with the cramer-rao theorem are used to determine a measurement set such that the estimated parameters have minimum variance. a detailed description of the procedure and its theoretical foundation can be found in. assuming that the measurement errors are additive and gaussian, the fim is given by: fim = jtwj     where w is the inverse of the measurement error covariance matrix and j denotes the sensitivity coefficient matrix for the measured response variables: the quantity {p …, pr} denotes the identifiable parameter vector and {, …, } denotes the measured response variable vector. the cramer-rao inequality establishes a lower bound on the variance of the identifiable parameters given by: σ ≥ ii     the 95% confidence interval for a parameter is given by: ci = ± σ     in equation the lower bound of the variance is used. symmetry of the confidence region about the nominal value is assumed. this results in the following definition of the percentage deviation from the nominal value: the optimal measurement set is chosen such that the sum of the percentage error for all the identifiable parameters is minimized. in this work, the optimal set is determined by a brute-force search over all combinations of measurement sets subject to restrictions that may be imposed by the system. doyle and co-workers have developed effcient rational algorithms to determine the optimal measurement set with minimum computational burden. the confidence intervals for non-identifiable parameters are infinitely large and hence are eliminated from the analysis. identifiability for these parameters can be obtained only by a change in the experimental design or by the selection of an alternative model structure. step 2: state estimation algorithm generally, it is not possible to measure all time-varying components in a metabolic or signaling network. however, there are several techniques from systems engineering to estimate the behavior of unmeasured components given partial measurements of other system constituents. bastin and dochain have used an adaptive nonlinear observer for estimation of specific growth rate and biomass concentration. given accurate models, extended kalman filters have had success in several biological applications. artificial neural networks have also found applications where dynamic models are not available. in this work, an extension of dynamic flux balance analysis is developed to estimate unmeasured concentration and reaction rate trajectories given partial measurement sets. the premise of this approach is straightforward: cellular processes have evolved regulatory structures that optimally use cellular resources. this premise translates into two postulates; network flows are managed to minimize internal accumulation and networks are managed to minimize the number of edges carrying flux at any given time. these two requirements are analogous to a classic problem in automatic control, namely, the state regulator problem. the srp based estimator uses the measurement set selected from step to estimate unknown concentration and reaction rate trajectories via a constrained convex programming problem. the srp estimator constrained by the key measurements captures the optimal cellular behavior of the system. estimates of the reaction rates at time step k and protein/gene/metabolite concentrations at time step k + are determined by the srp and the discrete mass balance equations. the srp must be solved at each sampling interval to obtain estimates of the unknown rates and concentrations. consider the model; the discrete mass balance equations over a p step horizon are given by: where the matrices k+ k, x and c are defined as and the matrix r is given by the srp estimator is a quadratic program with two cost terms, the cost of intermediate accumulation and the cost of operating a network reaction. the srp penalizes intracellular metabolite or protein accumulation, but does not explicitly forbid it. moreover, because reactions introduce an additional cost, the srp only utilizes those reactions required to satisfy the mass balances and thermodynamic constraints. formally, the estimation problem is given by: subject to: k + ≥     αr ≤ k ≤ βr     the srp problem is subject to non-negativity constraints, flux-directionality constraints and constraints imposed by the measurement set. specifically, constraint of equation forces state estimates belonging to the measurement set to equal the corresponding measured value. the quantities denote measurements that may have been corrupted with noise. specifies the tolerance around the measurement within which the estimate is constrained to lie. the term Ξx defines the measurement set. the matrix wx defines the cost of intermediate accumulation whereas the matrix wr represents the reaction cost: for the caspase system considered in this work, the wx and wr are taken to be the order of magnitude of the inverse of the maximum value of the corresponding state or rate. this requires only approximate information regarding ranges of the protein concentrations and identification of the slow versus the fast reactions: step 3: parameter estimation the estimates of the concentration profiles and the reaction rates allow efficient determination of the parameter values by decoupling the full parameter estimation into multiple sets. each set consists of parameters associated with one reaction rate. the parameters are obtained by minimizing the difference between the estimates of each reaction rate and that predicted by the kinetics r, which is a function of the concentrations. in case of the first iteration the minimization follows: where ri are the individual reaction rates and nr is the total number of reactions. the kinetic parameters associated with a reaction rate equation are determined independently from those with other reactions, i.e., the parameter estimation is decoupled with respect to each reaction. for subsequent iterations, the bayesian estimation formulation in is used. in this formulation, in addition to the difference between the estimates of the reaction rates and the model predictions, the deviations of parameter values from those obtained after the previous iteration are minimized. the formulation can be represented as: in the above equations, is the estimate of the ith reaction rate and is the estimates of the concentrations obtained from the srp algorithm, ri is the predicted rate of the ith reaction from the kinetics in equation p are the parameters associated with the ith reaction rate, p are the parameter values obtained from the previous iteration, and vε and vp are the variances of the estimates of the reaction rates and the prior parameter estimates. the parameter variances are determined using the fisher information matrix. the variances of the non-identifiable parameters are infinite and penalty for deviations for these parameters are not considered in equation the variance for the estimates of the reaction rates can be determined from the expected noise in the measurements from which the estimates are obtained. step 4: model invalidation tests given the iterative nature of this framework, a termination criterion must be established. poolla et al. have shown that for certain experimental data, it is not possible to confirm whether the model is really valid; however, one can conclude whether the model is not contradicted by the given data. model validation tests are usually based on the difference between the simulated and measured output and some statistics about these differences. typical statistics for the model errors include maximum absolute value, mean value and variance. in this work, model invalidity is tested by determining the model prediction errors using the estimated parameters. this error is calculated as: to implement this test, experimental data that was not used in the srp algorithm is required. the statistic used is the maximum and mean value of the errors for the measured states. when the prediction errors are below a certain desired value, the iterations are terminated. step 5: model-based optimal experiment design the optimal experiment design determines the optimal experiment to be performed for the next iteration such that there is maximum information content in the measurements. this would maximize the accuracy of the estimated parameters. the model-based optimal experiment design uses the fisher information matrix as a measure of the amount of information contained in a given set of measurements about the model parameters. the optimization searches through the space of experimental conditions or some parameterizations of the experimental protocol. for example, an optimal ligand input can be parameterized into a time series profile such that the optimization variables are the levels of ligand at different times. naturally, the optimization will be restricted by the limitations in the experimental conditions and apparatus. there exist several fim-based optimality measures that quantify the overall informativeness of the measurements. among these, parameter identifiability and the d-optimality are the most widely used measures. for accurate model identification it is critical that maximum number of the parameters be estimated accurately. thus, maximizing the number of identifiable parameters is the primary criterion proposed for determining the next experimental design. the orthogonal procedure proposed by mcauley and co-workers is used to determine the number of identifiable parameters. there can be multiple experimental designs with the same maximum number of identifiable parameters. the selection among these is done so as to maximize the informativeness of measurement data. for this purpose, the d-optimality criteria is proposed. the use of d-optimality translates to minimizing the confidence interval of all the identifiable parameter estimates. the optimal experiment design criterion is shown as follows: where e denote the feasible experimental conditions, fim is given in equation and r denotes the number of identifiable parameters. thus, the identifiability is maximized in the sense that the hyper-dimensional confidence interval is minimized. for experimental designs with the maximum number of identifiable parameters, the one with the highest determinant of the fisher information matrix is selected as the optimal design for the next experiment. a case study the proposed iterative model identification is applied to the function of caspase- and caspase- in apoptosis. caspase enzymes are at the core of the cell's suicide machinery. these enzymes are activated either by an external signal or by stress, and activated enzymes will then dismantle the cells. varner and co-workers have developed a model for the caspase function in apoptosis. the model describes the key elements of receptor-mediated and stress-induced caspase activations. the model consists of states and reactions with parameters. the in silico experiments in this study use the varner model as the "actual" system. measurements are assumed to be obtained from this "actual" system corrupted with up to 10% noise. the iteration starts using a "initial" parameter set, generated by perturbing the parameter values of the varner model by 70–100%. the external and stress signals that activate the caspase system are considered as the manipulated variables in step model refinement is performed either by determination of an optimal experiment or by optimal refinement of the measurement set. the performance of each of these criteria for model refinement requires to be tested. therefore, in the first iteration, a "preliminary" suboptimal experiment with a suboptimal measurement set is considered. moreover, in most cases of model identification, it is expected that preliminary experimental data is available. this data is usually obtained from a suboptimal experiment design and does not include the optimal set of measurements. the second iteration is performed in two ways; one by improving the measurement set and second by improving the experiment design. this tests the performance of both refinement criteria. the sequence of events is shown in figure it should be noted that it would be best to use optimal experiment design along with the optimal measurement set. however, it may not always be possible to do so due to feasibility issues specific to the particular system. hence, this approach is not considered in this work. model identification under less constrained conditions using a similar framework is included in. after the first iteration, it is observed that there is a significant improvement in the predictions with the estimated parameters for both the protein concentrations and the reaction rates, as shown in figures and the high errors with the "initial" parameters demonstrate that there is no bias in the results based on the starting guess values for the parameters and that there is indeed an improvement in prediction of both the protein concentrations and the reaction rates. however, the improvement is not suffcient as observed from the invalidation test. this warrants a second iteration. in general, the model predictions improve with the second iteration, as shown again in figures and however, it is observed that the predictions are better for the case with the optimal experiment design, in spite of a suboptimal measurement set. this is due to the fact that model performance depends strongly on the accuracy of the estimated parameters. using the suboptimal experiment, only of the parameters were identifiable. the optimal measurement set simply improved the confidence in these parameters. on the other hand, the optimal experiment increased parameter identifiability to parameters even with the suboptimal measurement set. these result indicates that performance of model identification is strongly linked with parameter identifiability. discussion in the proposed algorithm, the network topology and the mechanism of interactions in the pathway are assumed to be known. several approaches have been proposed in literature to determine the network connectivity from experimental data. in the case of unknown connectivity, these approaches should be used prior to the proposed model identification. the mismatch between the model and the actual network can appear in two different aspects of the algorithm. first, the srp step can fail because there exist no feasible state and flux estimates that satisfy the measurement constraints. such scenario arises mainly due to the network topology mismatch, and has low probability to occur due to the large degree of freedoms in a typical biological system. alternatively, a well-designed model validation step catches the mismatch between the model and the real system, e.g., an independent measurement set contradicts the model prediction. also, if multiple models are proposed for a particular biological process, model discrimination methods can be used to identify the correct model structure. the effect of incorrect connectivity and/or mechanism would depend on the degree of mismatch and is case dependent. the fact that cellular processes are carried out in an optimal manner lends tremendous promise to the success of this approach. in case the assumed optimality does not represent the in vivo behavior, the estimates from srp may be inaccurate. however, the measurements can attenuate this problem by restricting part of the estimates to match the observations. the parameter estimates may also be inaccurate if the real system deviates considerably from the assumed optimal behavior and this deviation is not captured by the measurements. the model identification framework has applicability to all systems that could be represented in the form shown in equations and all biological processes are complex, interconnected networks. a feature common to these processes is that they have a fixed connectivity. the proposed algorithm for model development could be applied to metabolic networks, signaling processes and gene networks. the computation burden for solving the srp is minimal as it involves only a quadratic programming problem. the parameter estimation is also not computationally intensive due to the decoupling facilitated by srp. a global optimization algorithm can also be used in the parameter estimation instead of the gradient search method to avoid convergence to local minima. however, the computation burden of parameter estimation will increase. to avoid local minima and high computation cost, the first few iterates can utilize a global optimization method, while the remaining iterations can implement a gradient search algorithm. due to the iterative nature of the approach, errors in parameter estimates can be tolerated as the corrections will be made in the next iteration with the optimal experiment. the optimal measurement selection is performed by a brute force search in this work. for very large systems the computation burden for this process grows exponentially. computationally efficient algorithms for optimal measurement selection can be used. efficient measurement selection algorithms and the decoupling of parameter estimation for individual reaction rates into separate optimization problems result in good scalability properties of the proposed algorithm for large scale systems. one limitation of the approach is that it is dependent on the weights in equation for the minimization of the cellular resources. the choice of weights used in this work has provided accurate results but it requires information of the order of magnitude of the concentrations and rates of the system under study. efficient schemes for determining the weights for metabolic networks has been developed by varner and co-workers.  an iterative methodology for model identification from experimental data is developed in this paper. identifiability tests are performed for an optimal measurement set selection for a given experimental design. the optimal measurements represent maximum information such that the model identification process is maximally benefitted. the model identification process is decoupled into two parts. in the first, the measurements are used to estimate all the unmeasured quantities of the system. this is achieved using the state regulator problem formulation which is based on the assumption that the cell is an optimal strategist and uses its resources in an optimal manner. the srp algorithm developed in this work has shown promising results. the average errors in the estimates for a significant fraction of the unmeasured responses is less than 10%. the accuracy of the estimates obtained by the srp decreases with decrease in the information content due to suboptimal measurement set. in the second part, the full state and rate estimates are used to determine the model parameters. the decoupling relaxes considerable computation burden compared to estimating all model parameters simultaneously from the limited measurements. in the final step, a model-based experiment design determines the optimal experimental procedure that generates the most informative measurements for the next iteration. a strong dependence is observed between parameter identifiability and model performance. thus, it is critical that the experiment design and measurement set be chosen such that maximum number of parameters are identifiable. tools developed for quantitative analysis of the dynamics of cellular pathways have tremendous potential in improving the predictive capabilities of biological systems especially in cases where experimental data is available but the kinetic parameters of the pathway reactions are unknown. the model developing tools are used for a host of applications and systems analysis. the measurement selection algorithm presented in this work is freely available as part of a model analysis and development toolkit, biosens. 
an optimal experiment design using the parameter identifiability and d-optimality criteria is formulated to provide "rich" experimental data for maximizing the accuracy of the parameter estimates in subsequent iterations. optimal experiment design was determined to be critical for model identification. the proposed algorithm has general application to modeling a wide range of cellular processes, which include gene regulation networks, signal transduction and metabolic networks. an iterative scheme is introduced for model identification using available system knowledge and experimental measurements. recent advances in molecular biology techniques provide an opportunity for developing detailed mathematical models of biological processes. the iterative scheme is tested on a model for the caspase function in apoptosis where it is demonstrated that model accuracy improves with each iteration. the importance of model identifiability tests for optimal measurement selection is also considered. the full system information is used for estimation of the model parameters. the scheme includes a state regulator algorithm that provides estimates of all system unknowns.
 recent advances in molecular biology techniques provide an opportunity for developing detailed mathematical models of biological processes. an iterative scheme is introduced for model identification using available system knowledge and experimental measurements.  the scheme includes a state regulator algorithm that provides estimates of all system unknowns. the full system information is used for estimation of the model parameters. an optimal experiment design using the parameter identifiability and d-optimality criteria is formulated to provide "rich" experimental data for maximizing the accuracy of the parameter estimates in subsequent iterations. the importance of model identifiability tests for optimal measurement selection is also considered. the iterative scheme is tested on a model for the caspase function in apoptosis where it is demonstrated that model accuracy improves with each iteration. optimal experiment design was determined to be critical for model identification.  the proposed algorithm has general application to modeling a wide range of cellular processes, which include gene regulation networks, signal transduction and metabolic networks. 