BACKGROUND
mathematical modelling has become an indispensable tool for modern systems biology  <cit> .
simple qualitative descriptions are proving increasingly insufficient for understanding the intricate dynamical complexity of biological phenomena.
as a result, quantitative mathematical models are now routinely used in order to describe and analyze the complex dynamics generated by protein interactions  <cit> , metabolic pathways  <cit> , regulation of gene expression  <cit> , and other biochemical processes.
a successful modelling effort is necessarily an iteration between model analysis and experiments.
testing the appropriateness of a mathematical description of any physical process should be done against experimental data, but at the same time, models should inform the design of new experiments.
traditionally, experiments have been designed using heuristic approaches: experience, intuition, or simple causal analyses.
evidently, such heuristically designed experiments are not always maximally informative, a great impediment given the cost and effort involved in the development of new measurement techniques and the implementation of standard experiments.
as a result, it is becoming increasingly necessary to systematically design more rigorous and predictive experiments, in order for the iterative process involving computational modelling to result in reliable models.
to date, the majority of studies addressing experiment design for biological networks has adopted a system identification approach.
in this context, experiments are designed such that the resulting data are most informative about model structure or parameter values - see, for example,  <cit>  and  <cit> , respectively.
several groups considered statistically-orientated frameworks for optimal structure identification  <cit>  or for parameter identification  <cit> .
these approaches aim to find the weighted least squares of differences between data and model prediction and make use of the fisher information matrix and the associated notions of a-, d-, and e-optimality.
in this framework, yue et al.
<cit>  examine optimally designed parameter estimation methods that are robust to model uncertainties .
in numerous practical situations, accumulated biological knowledge about a system of interest can constrain the set of plausible model structures.
in this case, one can enumerate a finite set of network topologies, closely corresponding to concrete biological hypotheses.
experiment design in this context would aim for the efficient discrimination between these well-defined alternative models; in more concrete terms, several mathematical models, corresponding to the different network topologies, can describe the behaviour of this system, within error bounds reflecting uncertainty in the data due to the experimental environment and inaccuracies of measurements  <cit> .
discriminatory experiment design and model invalidation can then be used to differentiate between them.
this is because mathematically, one can never validate a model  <cit> .
at best, a model will be capable of explaining all the available data and can be tested against some of its predictions.
therefore narrowing down on the correct model can only be done from the other direction through invalidation, in order to systematically 'cross out' incorrect models.
this results in an iterative cycle of system modelling, experiment design and subsequent model performance analysis that systematically proposes and then invalidates models that cannot represent the behaviour of the system.
in order to optimally discriminate between candidate models, the experiments need to be carefully designed and implemented to produce new data that can be used to invalidate a seemingly good but incorrect model.
various aspects of model discriminatory experiment design have been addressed in the literature.
bardsley et al.
<cit>  investigated the problem of how measurements should be spaced in time to perform an optimally discriminating experiment between two models, and how many of them are required.
more specifically, they compared different patterns of measurement spacings .
chen and asprey  <cit>  developed statistical approaches to parameter estimation, the assessment of model fit, and model discrimination, assuming that the response variables are uncertain.
in this framework, model discrimination is based on a bayesian approach, which assigns prior 'goodness' probabilities to each model, updates these after each experiment and chooses the model with the likelihood that has become sufficiently large compared to others.
an alternative frequentist method uses repeated hypothesis tests to reject models one by one.
donckels et al.
<cit>  separated the uncertainty of the model predictions and the uncertainty of the measurements and used these to design the next experiment such that it is most informative.
as opposed to the traditional approach, here the expected information content of the newly designed experiment is also taken into account  in order to assess the uncertainties more accurately.
kreutz and timmer  <cit>  gave a review of approaches to parameter estimation and model discrimination .
they also discussed relevant classical statistical aspects of experiment design, such as randomization, replication, and confounding.
tidor and co-workers  <cit>  developed dynamic model-based controllers that drive the output along a prescribed target trajectory .
if such a control input signal achieves the required output trajectory in an experiment, then the model is more accurate than another model which gives a different output trajectory for this particular input.
in  <cit> , kremling et al.
presented three methods for optimal pairwise discriminating experiment design, and compared them on a test example.
their first method compares combinations of certain initial input levels and subsequent changes in input in order to determine which combination will lead to the largest difference in the outputs.
their second method replaces models with their linearized counterparts in order to find a sinusoidal input with a frequency that maximizes the difference between phase shifts of the two models.
their third method follows the work in  <cit> , and aims to find an input profile that brings the output responses of the two models as far apart as possible.
the distance is measured by a weighted objective function.
the weighting is set up such that if the measurement error of an output variable is large, then the difference of these outputs contributes less to the weighted objective function.
the authors concluded that the most appropriate method strongly depends on the possible ways to stimulate the system and the quality of the measurements.
in our approach to data-supported computational modelling of biological networks, we take the view that one should follow an iterative procedure that includes model identification , model discrimination  and model invalidation .
all three tasks present serious challenges, and remain important areas of research and investigation.
in this paper we address the problem of model discrimination.
specifically, we present a framework for defining and designing optimally discriminating experiments, that is, experiments that are the best  at discriminating between rival models.
there are cases when it is difficult or even impossible to distinguish between rival models due to the incomplete observability of their internal states.
tests exist to identify such cases  <cit> .
even when model discrimination is possible, it can be expected to be difficult as the starting assumption is that the rival models both describe all available data well.
our key principle is to maximize the difference between the outputs of two different models, in particular, the l2-norm of the output difference.
although similar in principle, our investigation follows a direction distinctly different and more practical from the work in  <cit> : we use deterministic models that do not take account of measurement noise directly.
instead, we try to make the outputs of the two models as distant as possible to ensure that even a noisy measurement has a good chance of discriminating between them.
we propose three approaches to achieve this goal.
in the first approach, the initial condition design for model discrimination, we find the initial state of the system which results in the most discriminating output between the two examined models  <cit> .
the second method, input design for model discrimination, assumes the possibility for external stimulation during the experiment and searches for the best such stimulus from a set of allowable stimuli.
this approach is reminiscent of but different from the second method in  <cit>  - there, the difference between phase shifts is maximized, whereas in our method the difference between amplitudes is maximized.
the third method, design of structural changes for model discrimination, combines optimal initial condition choice with optimal systemic modifications.
the latter reflects the assumption that in the experiment it is possible, for example, to up-or down-regulate the expression of certain genes, either through genetic manipulations or other techniques such as rnai technology.
the gene product may be an enzyme whose concentration is not explicitly modelled but is reflected in a chemical rate constant, or some protein which exists in  phosphorylated and dephosphorylated forms such that the sum of their concentrations is constant.
in our mathematical model this means a free choice in some parameter values within given intervals.
in all three approaches, we cast the problem in an optimization framework and use the sum of squares  technique  <cit>  for the experiment design, allowing us to treat the nonlinear system descriptions directly.
the theoretical results are demonstrated by the application of each method to a discrimination problem for two models of signal processing for chemotaxis in dictyostelium amœbæ.
RESULTS
problem formulation
in this work, we consider different models describing the same biological system by a set of ordinary differential equations using, for example, mass action, generalized mass action, or michaelis-menten kinetics.
in general, the ith model takes the form
where u is a q-dimensional vector denoting the input, xi is an ni-dimensional vector denoting the state, yi is an ℓ-dimensional vector denoting the output  and gi is matrix-valued, with size ni × q. the structure of the functions fi, gi and hi will depend on the modelling framework in use to describe the biological system, but we assume that all of them are smooth.
here, the output function represents measurements an experimenter obtains from the system, and the input function represents the stimuli or perturbations the experimenter could introduce to the system during the experiment.
for mathematical simplicity we assume that the input does not affect the output directly.
in this paper our aim is to discriminate between two models of the form , which have n <dig> and n <dig> state variables, respectively.
as these two models represent the same underlying biological system, we require that they both generate the same steady states and fit already available experimental data: our aim is to design the next experiment that will allow their discrimination.
a natural way to formulate the discrimination problem is to concatenate the two models and generate the difference between their outputs:
we call an experiment optimal if the difference between the outputs of the two models  is maximal over a set of experimental perturbations of bounded 'size'.
in technical terms, we aim to pick the best point in a set of allowable perturbations of the initial state conditions , the set of inputs u , or the set of some admissible parameter changes and the set of common initial conditions  in order to maximize the so-called l2-distance between the outputs of the two rival models:
to facilitate interpretation, we implement a change of coordinates that places the investigated steady state at zero in both models.
we assume that the outputs are identical in this common steady state, now the origin: h <dig> = h <dig>  throughout this paper it is also assumed that the examined steady state is asymptotically stable in both models in .
since experiments must be implemented in finite time, we require that the designed input u be zero after some future time.
for convenience, we sometimes relax this requirement and only assume that u is 'very small' after a certain time.
clearly, since there is only one experimental setup in reality, the input u must be identical for the two models.
in the case of linear systems, the description of the concatenated system  becomes
with
where the above matrices are of appropriate dimensions.
we assume that all eigenvalues of both a <dig> and a <dig> have negative real parts , hence they define asymptotically stable systems.
this makes a hurwitz too.
initial condition design for model discrimination
many biological experiments drive a cellular system into an informative out-of-equilibrium state , and then glean information from the patterns of return to equilibrium in the absence of an input.
in an optimization formulation, this amounts to searching for normalized initial conditions x <dig> = x <dig> for the two models one wishes to discriminate between, that maximize the output difference ||y|| <dig> - where y is defined in  - for the unforced system .
here, we assume that the two alternative model representations of the system are written in terms of the same chemical species, thus n <dig> = n <dig> = n.
linear case
if x <dig> is not required to be equal to x <dig>  then the solution can be borrowed from standard results in linear systems theory.
in particular, the optimal direction for the initial value of  can be found by the following procedure.
<dig>  find a positive semidefinite matrix, p ≥  <dig>  that solves the so-called observability lyapunov equation
the solution p is called the observability gramian  <cit> .
<dig>  find the normalized eigenvector  corresponding to the largest eigenvalue  of p, that is, for  find  such that
indeed, the direction x =  gives the maximum output energy, since the output energy is given by
however, this computation is not satisfactory since an experimentally meaningful initial condition should satisfy x <dig> = x <dig> = .
to enforce this condition, we can partition p into blocks of size n × n,
with this decomposition, the optimal initial state is the unit norm eigenvector  corresponding to the largest eigenvalue of the matrix
to see this, substitute x =  in  to get:
hence  is maximized exactly when  is the eigenvector corresponding to the largest eigenvalue of
nonlinear case
the ideas behind model discrimination in the linear case can be generalized for application to nonlinear systems.
however, we cannot explicitly compute the exact difference in the outputs of the two rival models ||y|| <dig>  our approach avoids simulations and concentrates on finding an upper bound on ||y|| <dig> using so-called storage functions  <cit>  and sum of squares algorithmic relaxations of the resulting optimization problem.
to determine an upper bound on ||y|| <dig> for system , suppose there exists a continuously differentiable function s: ℝ2n → ℝ satisfying
where d is a neighbourhood of the steady state defined by:
here we assume that d does not include states which are not physically meaningful, and the value of α will ensure this.
this implies that the system is dissipative with supply rate -hth.
suppose that the system is released from an initial state x inside the largest level set of s that fits into d, so that ||x1|| <dig> = ||x2|| <dig> = β ≤ α. in this case, integrating condition  and using , we get
if we let t → ∞, then
by the nonnegativity of s. this implies that  ≤ s), since condition  is valid within the whole region d and level sets of s are invariant.
hence we have found a way to bound , which involves constructing the function s.
it is worth noting that the result from the linear and nonlinear cases have a similar purpose.
whereas in the linear case the result is rooted in a lyapunov equality and provides optimal solution, in the general nonlinear case one has to be content with an estimate given by inequality .
a condition missing from the above construction is that the two system models should be released from the same initial state.
hence, following our discussion from the linear case, one has to construct an appropriately modified s, .
in the linear case, the desired initial conditions correspond to those that maximize the quadratic form .
that is, the optimal direction was that of the eigenvector corresponding to the largest eigenvalue of matrix r. this is also exactly the direction corresponding to the smallest semi-axis of the ellipsoid  = r for some r >  <dig>  in other words, we were looking for the infimum of γ >  <dig> for which the set
similarly to the linear case, we will use a geometric argument to achieve initial condition design in the nonlinear case.
here, () plays the role of the quadratic form , and one can now decrease γ >  <dig> from infinity until the shrinking level set  touches d' = {x ∈ ℝ2n| ||x1|| <dig> ≤ β, ||x2|| <dig> ≤ β}.
therefore, we need to solve the following optimization problem:
the sequence of results presented so far asserts that the presence of a function s with the properties delineated above provides an upper bound on the energy of the difference of the outputs of two rival models.
this information can be exploited to generate experimental initial conditions that drive the system towards this bound.
these methods, however, do not prescribe how one would go about finding such a function.
constructing a nonnegative function is in general a difficult problem.
however, recent advances in the theory of sum of squares provide a computationally tractable way to relax this problem  <cit> .
in a nutshell, instead of searching for a general nonnegative function, we can constrain our search to functions that can be parameterized as sums of squares of polynomials.
within this class, the problem can be solved through semidefinite programming, with worst-case polynomial-time algorithms .
therefore, our strategy to find a near optimal initial state for the nonlinear model discrimination is a two-step process.
first, we construct an sos function s that satisfies .
in the second step, we search for the direction in which  is maximal, that is, we solve the optimization problem  .
input design for model discrimination
a powerful approach to discriminate between two plausible models of a biological process is to design an experimental input that maximally differentiates between the dynamical behaviours of their outputs.
if this input generates qualitatively different patterns in the model outputs, then one can subject the actual physical system to this designed input and then eliminate the model which differs from this pattern.
the most general form of this optimally discriminating input problem is the following:
here we assume that we are designing one input and measuring one output.
we also assume that the input is of unit-energy.
this can be done without a loss of generality as one can scale the equations accordingly, depending on the amount of input  available and the properties of the system under study.
our goal is to maximize the difference between the two model outputs over a transient period after application of the new input.
recall that the two models describe currently available data equally well, so that for the same  input they have the same pre-stimulus steady states and the same outputs.
solving the general optimization problem in order to generate the maximally informative input is computationally challenging.
in fact, even the first order condition of optimality is a 2-variable differential equation with boundary conditions at both ends of the time interval  <cit> .
for that reason, our strategy will be based on approximating a maximally-discriminating input using a linearization of the system in , and then assessing its suitability for the nonlinear system by comparing the value it achieves to the supremum of the output difference l2-norm over the set of possible inputs.
this supremum will again be computed using an sos decomposition approach.
the benefits of this strategy reside in the fact that we can use established, simple methods to find an input that gives the maximal l2-norm output for the linearized system.
this  input can then be applied to the nonlinear system, and an assessment  made about how the realized output l2-norm compares to the optimal, maximally discriminating l2-gain.
designing an input profile using linearization
designing an input profile for optimal discrimination using linearization is more appropriately addressed in the frequency domain.
a standard result in the theory of linear systems states that in order to find the input that maximizes the output difference ||y|| <dig> given unit input ||u|| <dig> , we need to find where the bode magnitude plot, i.e., the plot of |g| versus frequency ω, peaks.
here, g = c-1b, where the matrices a, b and c are defined in  and j denotes the imaginary unit: j <dig> = - <dig>  see  <cit>  for more details.
when the frequency at which |g| peaks is ω <dig>  the corresponding input signal takes the form
for an appropriately small ε, with a being a normalizing constant to ensure that the energy of u is unit.
a straightforward generalization of this concept to multiple-input multiple-output systems exists  <cit> , which we will also use in this work.
obtaining an upper bound on the l2-gain of the system and comparing performance
to assess the near-optimal input designed using linearization, we can compare its performance in driving the difference in the output of the two rival nonlinear models to its upper bound.
this is the l2-gain of the system, for which we can again obtain an upper bound by constructing an appropriate storage function s. to do so, given the appropriately normalized system  and ε >  <dig>  we assume that the trajectories with input u with ||u|| <dig> = ε and initial condition  <dig> remain in a region d around the steady state  for all time.
if there exists a γ >  <dig>  and a continuously differentiable function s:  satisfying
then
in other terms, γ is the desired upper bound on the maximum difference in the output of the nonlinear rival models.
to see this, integrating condition  from  <dig> to t leads to
therefore, for t → ∞, if x ∈ d for the whole time, we obtain
here again, obtaining such a function s that provides the upper bound is difficult.
the task of finding this bound can be relaxed to solving an sos programme and its subsequent solution using semidefinite programming .
design of structural changes for model discrimination
a class of experiments is based on the introduction of internal changes, such as genetic and biochemical manipulations, to the system.
to mirror such experiments, we develop a methodology to pinpoint numerical changes of parameters in a system that maximize the difference between the outputs of two rival models of its internal structure.
since two such models are different, they do not necessarily have the same number of parameters.
therefore the design concentrates on the parameters that the models have in common, which we denote by pi, i =  <dig>  ..., m. we assume that their values can be chosen within closed intervals  , that is, p ∈ Π, where .
we rewrite  to underline the dependence of the model on those parameters as:
and again assume that n <dig> = n <dig> = n, u =  <dig>  f <dig> = f <dig> =  <dig> , h <dig> = h <dig> and let
the steady state of either model may change with changing parameter p. therefore the assumption f <dig> = f <dig> =  <dig> should be interpreted as a change of coordinates that shifts the steady state of each model to the origin individually for each p. we are not interested in how far the two equilibria shift per se, which is an algebraic problem, instead we are interested in the difference in their dynamic responses.
this would reflect a situation in which a change in parameters would not be reflected in a significant change in the steady-state but which could result in a substantial difference in the dynamics of the system.
as with the previous two methods, our methodology will rely on the construction of an appropriate function s that sets an upper bound on the difference between the outputs of the two models, followed by a computationally efficient formulation for the construction of this function using sos.
for the above system, suppose that there exists a function s: ℝ2n × Π → ℝ which is sufficiently smooth and satisfies
then
if the system is released from an initial state , p) ∈ d × Π where x is in a level set of s entirely contained in d, ||x1|| = ||x2|| = β ≤ α. the last inequality in  holds since s ≥  <dig>  the computational relaxation and implementation of the search for the function s is presented in the methods section .
once this function has been constructed, one can extract the optimal point  and parameter point  that maximizes the difference between the measured outputs of the two models.
a case study: signal sensing in dictyostelium discoideum
perfect adaptation is a critical feature of many cellular signalling networks - it allows a cell to respond to a stimulus, but to re-sensitize itself so that further increases in stimulus can be detected.
adaptation is commonly used in sensory and other signalling networks to expand the input range that a circuit is able to sense, to more accurately detect changes in the input, and to maintain homeostasis in the presence of perturbations.
one of the earliest examples of cellular networks exhibiting perfect adaptation is chemotaxis, which we use as a test case to illustrate our algorithms.
specifically, we use the chemotactic response in the social amœba dictyostelium discoideum.
under starvation, dictyostelium secretes cyclic amp  thus attracting other dictyostelium amœbæ to aggregate and form a multicellular slug and then a fruiting body, which produces spores.
experiments indicate that a step input of chemoattractant triggers a transient response, after which the chemosensory mechanism returns to its pre-stimulus values , indicating perfect adaptation  <cit> .
at least two different simple models can describe the adaptation mechanism observed when an amœba encounters the chemoattractant camp  <cit>  .
in both models, a chemotaxis response regulator r becomes active  through the action of an activator enzyme a when a camp ligand s appears.
however, the deactivating mechanism determined through the interaction of an inhibiting molecule i in the two models can be different.
since the sum of the concentrations of the active and inactive response regulators in the two models is constant, we can write rt = r* + r. consequently,  can be derived assuming mass action kinetics as:
with activation and deactivation rate constants kr and k-r.
in model  <dig> both molecules, a and i, are regulated by the external signal, which is proportional to camp concentration s. with rate constants ka, k-a, k-i and , the dynamics of a and i are given by:
in model  <dig> the inhibitory molecule i is activated through the indirect action of activator a instead of direct activation by sensing of ligand binding, giving:
where  is a rate constant.
the equations for a are obviously identical in both models.
the parameter values used are given in table  <dig>  simple manipulations show that the steady-state value for r* in model  <dig> is given by
while the steady-state value for r* in model  <dig> is given by
both are independent of the stimulus, explaining perfect adaptation.
the two models share the same unique steady state if , a condition we impose.
note also that a value of α =  <dig>  needs to be used in , as the equilibrium value for r* is at  <dig> .
initial condition design for model discrimination
for the initial condition discriminating design, we set the input to a basal level of s = s <dig>  and assume that all three concentrations, a, i, and r*, can be measured.
the most discriminating initial state  can be found based on the linearization of the system around its steady state using the main linear case result.
the common unit length initial state which provides the direction of the perturbation from equilibrium to maximize ||y <dig> -y2|| <dig> is then given by x <dig> = x <dig> = , where xi  are the state vectors of model  <dig> and  <dig>  respectively.
by applying the analogous results from the nonlinear case, the unit norm direction that maximizes the above function is also found to be x <dig> = x <dig> = , which illustrates that at least in this example, linearization can be capable of providing the correct information at a lower computational cost.
figure  <dig> compares the evolution of the states in models  <dig> and  <dig> from a common arbitrary initial state  and the common initial state generated by the nonlinear method  by taking the differences between the states in the two rival models.
input design for model discrimination
to discriminate between the two models based on an optimally chosen input profile, we first obtained an upper bound on the l2-gain of the difference system from input s to output  using the algorithm in section c of methods.
this bound was about  <dig> , a value also close to that determined through the linearization of the system .
in order to evaluate the performance of our algorithm, we simulated the original, nonlinear system subjected to different inputs between  <dig> and  <dig> time units including a constant, sine, cosine, sine with an exponentially increasing then decreasing multiplier, sinc, the function given by , or a square wave function .
as described earlier, for periodic inputs, the period of the input was determined by finding the frequency corresponding to the maximum amplification in the bode magnitude plot of the linearized system .
numerical estimates of maxu ||y <dig> - y2||2/||u|| <dig> with different inputs for the dictyostelium models by simulating the nonlinear system.
here maximization is over different frequencies for inputs where frequency makes sense.
output is either all states or i or r*.
perhaps the most notable outcome of the input design is that sinusoidal input perturbations generate the best l2-gains and are therefore superior to a step function for discriminating between rival chemotaxis models.
square wave stimulation is achievable in the reality of a laboratory.
this is important since step inputs are usually used in experiments, often at the exclusion of other input signals.
our studies demonstrate how more dynamic inputs, in this case an oscillating input , might be necessary to delineate subtle features of underlying network topologies.
design of structural changes for model discrimination
for the design of most informative structural and parameter perturbations, the equilibrium is dependent on the particular choice of parameters.
therefore, we illustrate in detail how to change coordinates in order to translate the equilibrium to zero.
let the state variables for model  <dig> be x <dig> = a, x <dig> = i, and x <dig> = r* and those for model  <dig> be x <dig> = a, x <dig> = i, and x <dig> = r*.
assume that in both models the parameters that can be modified before the experiment are rt, the total chemotaxis response regulator concentration and p = kr, the response regulator activation rate constant.
parameter values and the intervals of values that can be achieved are given in table  <dig>  the dynamics of the two models are given by
the steady states for these equations are given by
let , i =  <dig>   <dig>  j =  <dig>   <dig>   <dig>  we can now perform a change of coordinates to make the origin the equilibrium point of:
in order to discriminate between models  <dig> and  <dig>  we first solve the optimization programmes given in section d of methods, as explained in the section on design of structural changes for model discrimination.
we allow parameters rt and kr to vary.
we obtain that  =  =  for the initial conditions, and rt =  <dig> and kr = p =  <dig>  for the values of the parameters that have maximal discriminating power between the two models.
.
this means that we need to over-express the total number of chemotaxis response regulators and increase their rate of activation in order to see a large difference between the two models.
the optimization problems in all three cases were solved on a desktop computer.
the most challenging was the first sos programme for the last case with eight variables .
numerical methods will need to be improved in order to deal with sos programmes resulting from the analysis of more complex systems biological models.
CONCLUSIONS
in this paper we have developed methods for designing experiments to effectively discriminate between different models of a biological system.
these methods are tailored to generate maximally informative data that can be used to invalidate models of gene regulatory pathways by ruling out certain connectivities in their underlying biochemical reaction networks  <cit> .
we approached the problem in a unified framework, developing methodologies for initial conditions design , for the design of dynamic stimulus profiles, and for parameter modifications.
these types of manipulations cover a large spectrum of what is experimentally feasible, and this has largely informed our formulation of the problem and the approach to its investigation.
if the field of systems biology is to accelerate the pace of biological discovery, rigorous mathematical methods should be developed to link computational models of biological networks to experimental data in tight rounds of analysis and synthesis.
any informative model should be analyzed in light of existing data, but it should also be able to synthesize new experiments that further delineate the features of the underlying biological system.
despite many notable examples demonstrating the success of this iterative procedure, progress has been slow due to the ad hoc nature of its implementation: the iterations between the development of models and the production of data is still mostly guided by the intuition of the modellers, and no rigorous algorithms exist to render this process more systematic and less biased.
we believe that the work presented in this paper constitutes an important step in this direction.
by design, our formulation of the problem is of sufficient generality to accommodate many experiment design procedures, and is cast in a natural optimization framework.
acknowledging that optimality of experiment designs must always be balanced with biological and other practical constraints, our formalisms allow for the incorporation of limitation and constraints as dictated by the specific biological context.
for example, if demanding a sinusoidal input may be unrealistic in a laboratory setting, and an optimal input in a smaller input function space is practically required, such constraints can be added to the nonlinear optimization criteria.
we illustrated the applicability of our algorithms using two possible and widely accepted simplified models of the adaptation mechanism in dictyostelium discoideum chemotaxis.
evidently, these models do not capture the full complexity of the biological circuit responsible for chemotactic behaviour.
the models, however, illustrate the core circuit topologies that are sufficient to implement perfect adaptation in the system.
recent work investigating perfect adaptation demonstrated that despite the diversity of biochemical enzymatic networks, only a finite set of core circuits with defined topological features can execute a desired function  <cit> .
these findings highlight the possibility of distinguishing between mechanisms that implement a given biological function using simple models, empowered by model-discrimination methods such as those presented in this work.
we also applied the optimal experiment design methods described in this paper to invalidate models of the chemotaxis pathway in rhodobacter sphaeroides  <cit> .
there, the combination of a square wave profile stimulation and protein over-expression was necessary in the most challenging model discrimination problem.
this demonstrates the practical demand for sophisticated experiment design techniques.
the recipe for model discrimination that we propose involves collecting mostly time series data.
every new time point at which measurements are made increases the cost of experiments, and thus one must carefully balance the number of time points collected against the cost, and consider where along a time series to concentrate observations.
our methods naturally present a window into this question by providing the timescales at which data collection needs to be done to be maximally informative.
furthermore, if the optimal experiment is such that a differentiating dynamical phenotype only emerges several hours after a perturbation, our methods can be easily modified to balance optimality with practically measurable dynamics.
finally, many commonly used perturbations  lead to either extreme stress responses that put a cell in a modified physiological state, cell death, or quiescent states that do not have much measurable information about the underlying regulatory network.
experiments that generate less catastrophic failures of cellular networks under study, while being maximally informative, hold great promise for the study of biological networks.
finding this region in perturbation space, however, is a nontrivial task.
model-based design of experiments will undeniably be instrumental for that, ultimately leading to many important biological discoveries.

 mathematical modelling has become an indispensable tool for modern systems biology. simple qualitative descriptions are proving increasingly insufficient for understanding the intricate dynamical complexity of biological phenomena. as a result, quantitative mathematical models are now routinely used in order to describe and analyze the complex dynamics generated by protein interactions, metabolic pathways, regulation of gene expression, and other biochemical processes. a successful modelling effort is necessarily an iteration between model analysis and experiments. testing the appropriateness of a mathematical description of any physical process should be done against experimental data, but at the same time, models should inform the design of new experiments. traditionally, experiments have been designed using heuristic approaches: experience, intuition, or simple causal analyses. evidently, such heuristically designed experiments are not always maximally informative, a great impediment given the cost and effort involved in the development of new measurement techniques and the implementation of standard experiments. as a result, it is becoming increasingly necessary to systematically design more rigorous and predictive experiments, in order for the iterative process involving computational modelling to result in reliable models. to date, the majority of studies addressing experiment design for biological networks has adopted a system identification approach. in this context, experiments are designed such that the resulting data are most informative about model structure or parameter values - see, for example, and, respectively. several groups considered statistically-orientated frameworks for optimal structure identification or for parameter identification. these approaches aim to find the weighted least squares of differences between data and model prediction and make use of the fisher information matrix and the associated notions of a-, d-, and e-optimality. in this framework, yue et al. examine optimally designed parameter estimation methods that are robust to model uncertainties. in numerous practical situations, accumulated biological knowledge about a system of interest can constrain the set of plausible model structures. in this case, one can enumerate a finite set of network topologies, closely corresponding to concrete biological hypotheses. experiment design in this context would aim for the efficient discrimination between these well-defined alternative models; in more concrete terms, several mathematical models, corresponding to the different network topologies, can describe the behaviour of this system, within error bounds reflecting uncertainty in the data due to the experimental environment and inaccuracies of measurements. discriminatory experiment design and model invalidation can then be used to differentiate between them. this is because mathematically, one can never validate a model. at best, a model will be capable of explaining all the available data and can be tested against some of its predictions. therefore narrowing down on the correct model can only be done from the other direction through invalidation, in order to systematically'cross out' incorrect models. this results in an iterative cycle of system modelling, experiment design and subsequent model performance analysis that systematically proposes and then invalidates models that cannot represent the behaviour of the system. in order to optimally discriminate between candidate models, the experiments need to be carefully designed and implemented to produce new data that can be used to invalidate a seemingly good but incorrect model. various aspects of model discriminatory experiment design have been addressed in the literature. bardsley et al. investigated the problem of how measurements should be spaced in time to perform an optimally discriminating experiment between two models, and how many of them are required. more specifically, they compared different patterns of measurement spacings. chen and asprey developed statistical approaches to parameter estimation, the assessment of model fit, and model discrimination, assuming that the response variables are uncertain. in this framework, model discrimination is based on a bayesian approach, which assigns prior'goodness' probabilities to each model, updates these after each experiment and chooses the model with the likelihood that has become sufficiently large compared to others. an alternative frequentist method uses repeated hypothesis tests to reject models one by one. donckels et al. separated the uncertainty of the model predictions and the uncertainty of the measurements and used these to design the next experiment such that it is most informative. as opposed to the traditional approach, here the expected information content of the newly designed experiment is also taken into account in order to assess the uncertainties more accurately. kreutz and timmer gave a review of approaches to parameter estimation and model discrimination. they also discussed relevant classical statistical aspects of experiment design, such as randomization, replication, and confounding. tidor and co-workers developed dynamic model-based controllers that drive the output along a prescribed target trajectory. if such a control input signal achieves the required output trajectory in an experiment, then the model is more accurate than another model which gives a different output trajectory for this particular input. in, kremling et al. presented three methods for optimal pairwise discriminating experiment design, and compared them on a test example. their first method compares combinations of certain initial input levels and subsequent changes in input in order to determine which combination will lead to the largest difference in the outputs. their second method replaces models with their linearized counterparts in order to find a sinusoidal input with a frequency that maximizes the difference between phase shifts of the two models. their third method follows the work in, and aims to find an input profile that brings the output responses of the two models as far apart as possible. the distance is measured by a weighted objective function. the weighting is set up such that if the measurement error of an output variable is large, then the difference of these outputs contributes less to the weighted objective function. the authors concluded that the most appropriate method strongly depends on the possible ways to stimulate the system and the quality of the measurements. in our approach to data-supported computational modelling of biological networks, we take the view that one should follow an iterative procedure that includes model identification, model discrimination and model invalidation. all three tasks present serious challenges, and remain important areas of research and investigation. in this paper we address the problem of model discrimination. specifically, we present a framework for defining and designing optimally discriminating experiments, that is, experiments that are the best at discriminating between rival models. there are cases when it is difficult or even impossible to distinguish between rival models due to the incomplete observability of their internal states. tests exist to identify such cases. even when model discrimination is possible, it can be expected to be difficult as the starting assumption is that the rival models both describe all available data well. our key principle is to maximize the difference between the outputs of two different models, in particular, the l2-norm of the output difference. although similar in principle, our investigation follows a direction distinctly different and more practical from the work in: we use deterministic models that do not take account of measurement noise directly. instead, we try to make the outputs of the two models as distant as possible to ensure that even a noisy measurement has a good chance of discriminating between them. we propose three approaches to achieve this goal. in the first approach, the initial condition design for model discrimination, we find the initial state of the system which results in the most discriminating output between the two examined models. the second method, input design for model discrimination, assumes the possibility for external stimulation during the experiment and searches for the best such stimulus from a set of allowable stimuli. this approach is reminiscent of but different from the second method in - there, the difference between phase shifts is maximized, whereas in our method the difference between amplitudes is maximized. the third method, design of structural changes for model discrimination, combines optimal initial condition choice with optimal systemic modifications. the latter reflects the assumption that in the experiment it is possible, for example, to up-or down-regulate the expression of certain genes, either through genetic manipulations or other techniques such as rnai technology. the gene product may be an enzyme whose concentration is not explicitly modelled but is reflected in a chemical rate constant, or some protein which exists in phosphorylated and dephosphorylated forms such that the sum of their concentrations is constant. in our mathematical model this means a free choice in some parameter values within given intervals. in all three approaches, we cast the problem in an optimization framework and use the sum of squares technique for the experiment design, allowing us to treat the nonlinear system descriptions directly. the theoretical results are demonstrated by the application of each method to a discrimination problem for two models of signal processing for chemotaxis in dictyostelium amœbæ.  problem formulation in this work, we consider different models describing the same biological system by a set of ordinary differential equations using, for example, mass action, generalized mass action, or michaelis-menten kinetics. in general, the ith model takes the form where u is a q-dimensional vector denoting the input, xi is an ni-dimensional vector denoting the state, yi is an ℓ-dimensional vector denoting the output and gi is matrix-valued, with size ni × q. the structure of the functions fi, gi and hi will depend on the modelling framework in use to describe the biological system, but we assume that all of them are smooth. here, the output function represents measurements an experimenter obtains from the system, and the input function represents the stimuli or perturbations the experimenter could introduce to the system during the experiment. for mathematical simplicity we assume that the input does not affect the output directly. in this paper our aim is to discriminate between two models of the form, which have n and n state variables, respectively. as these two models represent the same underlying biological system, we require that they both generate the same steady states and fit already available experimental data: our aim is to design the next experiment that will allow their discrimination. a natural way to formulate the discrimination problem is to concatenate the two models and generate the difference between their outputs: we call an experiment optimal if the difference between the outputs of the two models is maximal over a set of experimental perturbations of bounded'size'. in technical terms, we aim to pick the best point in a set of allowable perturbations of the initial state conditions, the set of inputs u, or the set of some admissible parameter changes and the set of common initial conditions in order to maximize the so-called l2-distance between the outputs of the two rival models: to facilitate interpretation, we implement a change of coordinates that places the investigated steady state at zero in both models. we assume that the outputs are identical in this common steady state, now the origin: h = h throughout this paper it is also assumed that the examined steady state is asymptotically stable in both models in. since experiments must be implemented in finite time, we require that the designed input u be zero after some future time. for convenience, we sometimes relax this requirement and only assume that u is'very small' after a certain time. clearly, since there is only one experimental setup in reality, the input u must be identical for the two models. in the case of linear systems, the description of the concatenated system becomes with where the above matrices are of appropriate dimensions. we assume that all eigenvalues of both a and a have negative real parts, hence they define asymptotically stable systems. this makes a hurwitz too. initial condition design for model discrimination many biological experiments drive a cellular system into an informative out-of-equilibrium state, and then glean information from the patterns of return to equilibrium in the absence of an input. in an optimization formulation, this amounts to searching for normalized initial conditions x = x for the two models one wishes to discriminate between, that maximize the output difference ||y|| - where y is defined in - for the unforced system. here, we assume that the two alternative model representations of the system are written in terms of the same chemical species, thus n = n = n. linear case if x is not required to be equal to x then the solution can be borrowed from standard results in linear systems theory. in particular, the optimal direction for the initial value of can be found by the following procedure. find a positive semidefinite matrix, p ≥ that solves the so-called observability lyapunov equation the solution p is called the observability gramian. find the normalized eigenvector corresponding to the largest eigenvalue of p, that is, for find such that indeed, the direction x = gives the maximum output energy, since the output energy is given by however, this computation is not satisfactory since an experimentally meaningful initial condition should satisfy x = x =. to enforce this condition, we can partition p into blocks of size n × n, with this decomposition, the optimal initial state is the unit norm eigenvector corresponding to the largest eigenvalue of the matrix to see this, substitute x = in to get: hence is maximized exactly when is the eigenvector corresponding to the largest eigenvalue of nonlinear case the ideas behind model discrimination in the linear case can be generalized for application to nonlinear systems. however, we cannot explicitly compute the exact difference in the outputs of the two rival models ||y|| our approach avoids simulations and concentrates on finding an upper bound on ||y|| using so-called storage functions and sum of squares algorithmic relaxations of the resulting optimization problem. to determine an upper bound on ||y|| for system, suppose there exists a continuously differentiable function s: ℝ2n → ℝ satisfying where d is a neighbourhood of the steady state defined by: here we assume that d does not include states which are not physically meaningful, and the value of α will ensure this. this implies that the system is dissipative with supply rate -hth. suppose that the system is released from an initial state x inside the largest level set of s that fits into d, so that ||x1|| = ||x2|| = β ≤ α. in this case, integrating condition and using, we get if we let t → ∞, then by the nonnegativity of s. this implies that ≤ s), since condition is valid within the whole region d and level sets of s are invariant. hence we have found a way to bound, which involves constructing the function s. it is worth noting that the result from the linear and nonlinear cases have a similar purpose. whereas in the linear case the result is rooted in a lyapunov equality and provides optimal solution, in the general nonlinear case one has to be content with an estimate given by inequality. a condition missing from the above construction is that the two system models should be released from the same initial state. hence, following our discussion from the linear case, one has to construct an appropriately modified s,. in the linear case, the desired initial conditions correspond to those that maximize the quadratic form. that is, the optimal direction was that of the eigenvector corresponding to the largest eigenvalue of matrix r. this is also exactly the direction corresponding to the smallest semi-axis of the ellipsoid = r for some r > in other words, we were looking for the infimum of γ > for which the set similarly to the linear case, we will use a geometric argument to achieve initial condition design in the nonlinear case. here, () plays the role of the quadratic form, and one can now decrease γ > from infinity until the shrinking level set touches d' = {x ∈ ℝ2n| ||x1|| ≤ β, ||x2|| ≤ β}. therefore, we need to solve the following optimization problem: the sequence of results presented so far asserts that the presence of a function s with the properties delineated above provides an upper bound on the energy of the difference of the outputs of two rival models. this information can be exploited to generate experimental initial conditions that drive the system towards this bound. these methods, however, do not prescribe how one would go about finding such a function. constructing a nonnegative function is in general a difficult problem. however, recent advances in the theory of sum of squares provide a computationally tractable way to relax this problem. in a nutshell, instead of searching for a general nonnegative function, we can constrain our search to functions that can be parameterized as sums of squares of polynomials. within this class, the problem can be solved through semidefinite programming, with worst-case polynomial-time algorithms. therefore, our strategy to find a near optimal initial state for the nonlinear model discrimination is a two-step process. first, we construct an sos function s that satisfies. in the second step, we search for the direction in which is maximal, that is, we solve the optimization problem. input design for model discrimination a powerful approach to discriminate between two plausible models of a biological process is to design an experimental input that maximally differentiates between the dynamical behaviours of their outputs. if this input generates qualitatively different patterns in the model outputs, then one can subject the actual physical system to this designed input and then eliminate the model which differs from this pattern. the most general form of this optimally discriminating input problem is the following: here we assume that we are designing one input and measuring one output. we also assume that the input is of unit-energy. this can be done without a loss of generality as one can scale the equations accordingly, depending on the amount of input available and the properties of the system under study. our goal is to maximize the difference between the two model outputs over a transient period after application of the new input. recall that the two models describe currently available data equally well, so that for the same input they have the same pre-stimulus steady states and the same outputs. solving the general optimization problem in order to generate the maximally informative input is computationally challenging. in fact, even the first order condition of optimality is a 2-variable differential equation with boundary conditions at both ends of the time interval. for that reason, our strategy will be based on approximating a maximally-discriminating input using a linearization of the system in, and then assessing its suitability for the nonlinear system by comparing the value it achieves to the supremum of the output difference l2-norm over the set of possible inputs. this supremum will again be computed using an sos decomposition approach. the benefits of this strategy reside in the fact that we can use established, simple methods to find an input that gives the maximal l2-norm output for the linearized system. this input can then be applied to the nonlinear system, and an assessment made about how the realized output l2-norm compares to the optimal, maximally discriminating l2-gain. designing an input profile using linearization designing an input profile for optimal discrimination using linearization is more appropriately addressed in the frequency domain. a standard result in the theory of linear systems states that in order to find the input that maximizes the output difference ||y|| given unit input ||u||, we need to find where the bode magnitude plot, i.e., the plot of |g| versus frequency ω, peaks. here, g = c-1b, where the matrices a, b and c are defined in and j denotes the imaginary unit: j = - see for more details. when the frequency at which |g| peaks is ω the corresponding input signal takes the form for an appropriately small ε, with a being a normalizing constant to ensure that the energy of u is unit. a straightforward generalization of this concept to multiple-input multiple-output systems exists, which we will also use in this work. obtaining an upper bound on the l2-gain of the system and comparing performance to assess the near-optimal input designed using linearization, we can compare its performance in driving the difference in the output of the two rival nonlinear models to its upper bound. this is the l2-gain of the system, for which we can again obtain an upper bound by constructing an appropriate storage function s. to do so, given the appropriately normalized system and ε > we assume that the trajectories with input u with ||u|| = ε and initial condition remain in a region d around the steady state for all time. if there exists a γ > and a continuously differentiable function s: satisfying then in other terms, γ is the desired upper bound on the maximum difference in the output of the nonlinear rival models. to see this, integrating condition from to t leads to therefore, for t → ∞, if x ∈ d for the whole time, we obtain here again, obtaining such a function s that provides the upper bound is difficult. the task of finding this bound can be relaxed to solving an sos programme and its subsequent solution using semidefinite programming. design of structural changes for model discrimination a class of experiments is based on the introduction of internal changes, such as genetic and biochemical manipulations, to the system. to mirror such experiments, we develop a methodology to pinpoint numerical changes of parameters in a system that maximize the difference between the outputs of two rival models of its internal structure. since two such models are different, they do not necessarily have the same number of parameters. therefore the design concentrates on the parameters that the models have in common, which we denote by pi, i =..., m. we assume that their values can be chosen within closed intervals, that is, p ∈ Π, where. we rewrite to underline the dependence of the model on those parameters as: and again assume that n = n = n, u = f = f =, h = h and let the steady state of either model may change with changing parameter p. therefore the assumption f = f = should be interpreted as a change of coordinates that shifts the steady state of each model to the origin individually for each p. we are not interested in how far the two equilibria shift per se, which is an algebraic problem, instead we are interested in the difference in their dynamic responses. this would reflect a situation in which a change in parameters would not be reflected in a significant change in the steady-state but which could result in a substantial difference in the dynamics of the system. as with the previous two methods, our methodology will rely on the construction of an appropriate function s that sets an upper bound on the difference between the outputs of the two models, followed by a computationally efficient formulation for the construction of this function using sos. for the above system, suppose that there exists a function s: ℝ2n × Π → ℝ which is sufficiently smooth and satisfies then if the system is released from an initial state, p) ∈ d × Π where x is in a level set of s entirely contained in d, ||x1|| = ||x2|| = β ≤ α. the last inequality in holds since s ≥ the computational relaxation and implementation of the search for the function s is presented in the methods section. once this function has been constructed, one can extract the optimal point and parameter point that maximizes the difference between the measured outputs of the two models. a case study: signal sensing in dictyostelium discoideum perfect adaptation is a critical feature of many cellular signalling networks - it allows a cell to respond to a stimulus, but to re-sensitize itself so that further increases in stimulus can be detected. adaptation is commonly used in sensory and other signalling networks to expand the input range that a circuit is able to sense, to more accurately detect changes in the input, and to maintain homeostasis in the presence of perturbations. one of the earliest examples of cellular networks exhibiting perfect adaptation is chemotaxis, which we use as a test case to illustrate our algorithms. specifically, we use the chemotactic response in the social amœba dictyostelium discoideum. under starvation, dictyostelium secretes cyclic amp thus attracting other dictyostelium amœbæ to aggregate and form a multicellular slug and then a fruiting body, which produces spores. experiments indicate that a step input of chemoattractant triggers a transient response, after which the chemosensory mechanism returns to its pre-stimulus values, indicating perfect adaptation. at least two different simple models can describe the adaptation mechanism observed when an amœba encounters the chemoattractant camp. in both models, a chemotaxis response regulator r becomes active through the action of an activator enzyme a when a camp ligand s appears. however, the deactivating mechanism determined through the interaction of an inhibiting molecule i in the two models can be different. since the sum of the concentrations of the active and inactive response regulators in the two models is constant, we can write rt = r* + r. consequently, can be derived assuming mass action kinetics as: with activation and deactivation rate constants kr and k-r. in model both molecules, a and i, are regulated by the external signal, which is proportional to camp concentration s. with rate constants ka, k-a, k-i and, the dynamics of a and i are given by: in model the inhibitory molecule i is activated through the indirect action of activator a instead of direct activation by sensing of ligand binding, giving: where is a rate constant. the equations for a are obviously identical in both models. the parameter values used are given in table simple manipulations show that the steady-state value for r* in model is given by while the steady-state value for r* in model is given by both are independent of the stimulus, explaining perfect adaptation. the two models share the same unique steady state if, a condition we impose. note also that a value of α = needs to be used in, as the equilibrium value for r* is at. initial condition design for model discrimination for the initial condition discriminating design, we set the input to a basal level of s = s and assume that all three concentrations, a, i, and r*, can be measured. the most discriminating initial state can be found based on the linearization of the system around its steady state using the main linear case result. the common unit length initial state which provides the direction of the perturbation from equilibrium to maximize ||y -y2|| is then given by x = x =, where xi are the state vectors of model and respectively. by applying the analogous results from the nonlinear case, the unit norm direction that maximizes the above function is also found to be x = x =, which illustrates that at least in this example, linearization can be capable of providing the correct information at a lower computational cost. figure compares the evolution of the states in models and from a common arbitrary initial state and the common initial state generated by the nonlinear method by taking the differences between the states in the two rival models. input design for model discrimination to discriminate between the two models based on an optimally chosen input profile, we first obtained an upper bound on the l2-gain of the difference system from input s to output using the algorithm in section c of methods. this bound was about, a value also close to that determined through the linearization of the system. in order to evaluate the performance of our algorithm, we simulated the original, nonlinear system subjected to different inputs between and time units including a constant, sine, cosine, sine with an exponentially increasing then decreasing multiplier, sinc, the function given by, or a square wave function. as described earlier, for periodic inputs, the period of the input was determined by finding the frequency corresponding to the maximum amplification in the bode magnitude plot of the linearized system. numerical estimates of maxu ||y - y2||2/||u|| with different inputs for the dictyostelium models by simulating the nonlinear system. here maximization is over different frequencies for inputs where frequency makes sense. output is either all states or i or r*. perhaps the most notable outcome of the input design is that sinusoidal input perturbations generate the best l2-gains and are therefore superior to a step function for discriminating between rival chemotaxis models. square wave stimulation is achievable in the reality of a laboratory. this is important since step inputs are usually used in experiments, often at the exclusion of other input signals. our studies demonstrate how more dynamic inputs, in this case an oscillating input, might be necessary to delineate subtle features of underlying network topologies. design of structural changes for model discrimination for the design of most informative structural and parameter perturbations, the equilibrium is dependent on the particular choice of parameters. therefore, we illustrate in detail how to change coordinates in order to translate the equilibrium to zero. let the state variables for model be x = a, x = i, and x = r* and those for model be x = a, x = i, and x = r*. assume that in both models the parameters that can be modified before the experiment are rt, the total chemotaxis response regulator concentration and p = kr, the response regulator activation rate constant. parameter values and the intervals of values that can be achieved are given in table the dynamics of the two models are given by the steady states for these equations are given by let, i = j = we can now perform a change of coordinates to make the origin the equilibrium point of: in order to discriminate between models and we first solve the optimization programmes given in section d of methods, as explained in the section on design of structural changes for model discrimination. we allow parameters rt and kr to vary. we obtain that = = for the initial conditions, and rt = and kr = p = for the values of the parameters that have maximal discriminating power between the two models.. this means that we need to over-express the total number of chemotaxis response regulators and increase their rate of activation in order to see a large difference between the two models. the optimization problems in all three cases were solved on a desktop computer. the most challenging was the first sos programme for the last case with eight variables. numerical methods will need to be improved in order to deal with sos programmes resulting from the analysis of more complex systems biological models.  in this paper we have developed methods for designing experiments to effectively discriminate between different models of a biological system. these methods are tailored to generate maximally informative data that can be used to invalidate models of gene regulatory pathways by ruling out certain connectivities in their underlying biochemical reaction networks. we approached the problem in a unified framework, developing methodologies for initial conditions design, for the design of dynamic stimulus profiles, and for parameter modifications. these types of manipulations cover a large spectrum of what is experimentally feasible, and this has largely informed our formulation of the problem and the approach to its investigation. if the field of systems biology is to accelerate the pace of biological discovery, rigorous mathematical methods should be developed to link computational models of biological networks to experimental data in tight rounds of analysis and synthesis. any informative model should be analyzed in light of existing data, but it should also be able to synthesize new experiments that further delineate the features of the underlying biological system. despite many notable examples demonstrating the success of this iterative procedure, progress has been slow due to the ad hoc nature of its implementation: the iterations between the development of models and the production of data is still mostly guided by the intuition of the modellers, and no rigorous algorithms exist to render this process more systematic and less biased. we believe that the work presented in this paper constitutes an important step in this direction. by design, our formulation of the problem is of sufficient generality to accommodate many experiment design procedures, and is cast in a natural optimization framework. acknowledging that optimality of experiment designs must always be balanced with biological and other practical constraints, our formalisms allow for the incorporation of limitation and constraints as dictated by the specific biological context. for example, if demanding a sinusoidal input may be unrealistic in a laboratory setting, and an optimal input in a smaller input function space is practically required, such constraints can be added to the nonlinear optimization criteria. we illustrated the applicability of our algorithms using two possible and widely accepted simplified models of the adaptation mechanism in dictyostelium discoideum chemotaxis. evidently, these models do not capture the full complexity of the biological circuit responsible for chemotactic behaviour. the models, however, illustrate the core circuit topologies that are sufficient to implement perfect adaptation in the system. recent work investigating perfect adaptation demonstrated that despite the diversity of biochemical enzymatic networks, only a finite set of core circuits with defined topological features can execute a desired function. these findings highlight the possibility of distinguishing between mechanisms that implement a given biological function using simple models, empowered by model-discrimination methods such as those presented in this work. we also applied the optimal experiment design methods described in this paper to invalidate models of the chemotaxis pathway in rhodobacter sphaeroides. there, the combination of a square wave profile stimulation and protein over-expression was necessary in the most challenging model discrimination problem. this demonstrates the practical demand for sophisticated experiment design techniques. the recipe for model discrimination that we propose involves collecting mostly time series data. every new time point at which measurements are made increases the cost of experiments, and thus one must carefully balance the number of time points collected against the cost, and consider where along a time series to concentrate observations. our methods naturally present a window into this question by providing the timescales at which data collection needs to be done to be maximally informative. furthermore, if the optimal experiment is such that a differentiating dynamical phenotype only emerges several hours after a perturbation, our methods can be easily modified to balance optimality with practically measurable dynamics. finally, many commonly used perturbations lead to either extreme stress responses that put a cell in a modified physiological state, cell death, or quiescent states that do not have much measurable information about the underlying regulatory network. experiments that generate less catastrophic failures of cellular networks under study, while being maximally informative, hold great promise for the study of biological networks. finding this region in perturbation space, however, is a nontrivial task. model-based design of experiments will undeniably be instrumental for that, ultimately leading to many important biological discoveries. 
such experiments should be optimally designed to produce data that are most likely to invalidate incorrect model structures. the numerical implementation of each method is considered in an example, signal processing in starving dictyostelium amœbæ. a problem commonly encountered in the computational modelling of biological networks is that alternative, structurally different models of similar complexity fit a set of experimental data equally well. at this point, new experiments maximizing the difference between the measured values of alternative models should be proposed and conducted. this opens the way to model invalidation, which can be used to perfect our understanding of biochemical networks. the first approach determines the'best' initial condition that maximizes the l distance between the outputs of the rival models. in this paper we develop methodologies for the optimal design of experiments with the aim of discriminating between different mathematical models of the same biological system. our general problem formulation together with the three proposed experiment design methods give the practitioner new tools for a systems biology approach to experiment design. our third method uses optimized structural changes to achieve the same goal. in this case, more than one molecular mechanism can explain available data.
 the success of molecular systems biology hinges on the ability to use computational models to design predictive experiments, and ultimately unravel underlying biological mechanisms. a problem commonly encountered in the computational modelling of biological networks is that alternative, structurally different models of similar complexity fit a set of experimental data equally well. in this case, more than one molecular mechanism can explain available data. in order to rule out the incorrect mechanisms, one needs to invalidate incorrect models. at this point, new experiments maximizing the difference between the measured values of alternative models should be proposed and conducted. such experiments should be optimally designed to produce data that are most likely to invalidate incorrect model structures.  in this paper we develop methodologies for the optimal design of experiments with the aim of discriminating between different mathematical models of the same biological system. the first approach determines the'best' initial condition that maximizes the l distance between the outputs of the rival models. in the second approach, we maximize the l2-distance of the outputs by designing the optimal external stimulus profile of unit l2-norm. our third method uses optimized structural changes to achieve the same goal. the numerical implementation of each method is considered in an example, signal processing in starving dictyostelium amœbæ.  model-based design of experiments improves both the reliability and the efficiency of biochemical network model discrimination. this opens the way to model invalidation, which can be used to perfect our understanding of biochemical networks. our general problem formulation together with the three proposed experiment design methods give the practitioner new tools for a systems biology approach to experiment design. 