{'fraction_non_alphanumeric': 0.030608094366483627, 'fraction_numerical': 0.013219442749644091, 'mean_word_length': 4.3886301369863014, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A robust optimization approach to experimental design for model discrimination of dynamical systems

A high-ranking goal of interdisciplinary modeling approaches in the natural sciences are quantitative prediction of system dynamics and model based optimization. For this purpose, mathematical modeling, numerical simulation and scientific computing techniques are indispensable. Quantitative modeling closely combined with experimental investigations is required if the model is supposed to be used for sound mechanistic analysis and model predictions. Typically, before an appropriate model of a experimental system is found different hypothetical models might be reasonable and consistent with previous knowledge and available data. The parameters of the model up to an estimated confidence region are generally not known a priori. Therefore one has to incorporate possible parameter configurations of different models into a model discrimination algorithm. In this article we present a numerical algorithm which calculates a design of experiments which allows an optimal discrimination of different hypothetic candidate models of a given dynamic system for the most inappropriate parameter configurations within a parameter range via a worst case estimate. The design criterion comprises optimal measurement time points. The used criterion is derived from the Kullback-Leibler divergence. The underlying optimization problem can be classified as a semi infinite optimization problem which we solve in an outer approximation approach stabilized by a homotopy strategy. We present the theoretical framework as well as the numerical algorithmic realization.


Introduction
High-ranking goals of interdisciplinary modeling approaches in the natural sciences are quantitative prediction of system dynamics and model based optimization. In particular in modern systems biology a related issue is to link molecular attributes to dynamic mechanisms and functional properties at the system level in order to mechanistically understand emerging functionality. For these purposes, mathematical modeling, numerical simulation and scientific computing techniques are indispensable. Quantitative modeling closely combined with experimental investigations is required if the model is supposed to be used for sound mechanistic analysis and model predictions.
Typically, before an appropriate model of a system is found different hypothetical models might be reasonable and consistent with previous knowledge and available data. The goal of this article is to derive, develop, implement and apply a numerical algorithm which calculates in a suitable sense an optimal design of experiments which allows the best discrimination of different hypothetic candidate models in form of ordinary differential equations (ODE). The algorithmic idea is to iteratively separate the response of different models by use of variations of experimental conditions and perturbations to the system.
To discriminate a set of candidate models against a given set of experimental data likelihood ratio tests based on bootstrap methods have been described in the literature, see e.g. [22], [46] or [48]. Ranking methods like Stewart's method ( [44]) or the well known Akaike information criterion (see e.g. [15]) are popular as well in the field of biological modeling. Applications can be found for example in [24] or [10]. In contrast to these approaches our work deals with the problem of designing experiments so that statistical methods can be exploited in an optimal sense for model discrimination. This differs from the approach to find an experimental design to best estimate the parameters of a model for a given experimental system in terms of criteria characterizing the confidence regions [25,7,6].
Different approaches to design experiments for model discrimination exist. Besides optimization methods (see e.g [29], [19], [47] or [26]) a model-based feedback controller see e.g [3] and Markov chain Monte Carlo sampling methods [34] have been used to construct an appropriate design. An overview of various experimental design techniques can be found in [27]. Here, we present a robust numerical optimization algorithm which calculates the optimal design of experiments allowing the best discrimination of different candidate ODE models. An appropriate model and its parameters up to an estimated confidence region are not known a priori. Therefore one has to incorporate possible parameter configurations of different models into a model discrimination algorithm. The aim is to calculate the most discriminable response of different models for the most inappropriate parameter configurations within a parameter range via a worst case estimate. In that context inappropriate parameter configurations refer to the case when different models have calibrated parameter values such that the models have the most similar response. This worst case estimate leads to the formulation of a maxmin optimization problem. Building on our previous work [43] we present an algorithm to compute robust optimal experimental designs. For the robustification we set up an outer approximation approach stabilized by a homotopy strategy. The article is organized as follows. In Section 2.1 we give a brief overview of so called Kullback-Leibler(KL)-optimality as discussed by López-Fidalgo et al. [32] in the context of model discrimination. In Section 2.2 we derive our optimal design criterion by use of KL-optimality. In Section 3 the theoretical framework for the calculation of a robust design via solution of a maxmin optimization problem is presented. The numerical implementation is discussed in Section 3.1. The homotopy solution strategy is presented in Section 3.2. We demonstrate applications of the algorithmic framework for two test cases from biology, an allosteric metabolic enzyme model for glycolytic oscillations and a model describing signal sensing in dictyostelium discoideum, results are presented in Section 4.

KL-optimal design
In this section a model discrimination criterion based on the Kullback-Leibler (KL) divergence called KL-optimality as discussed by López-Fidalgo et al. [32] is introduced. López-Fidalgo et al. [32] demonstrate that KL-optimality is consistent with T-optimality [5] and generalized T-optimality [49] which are well known model discrimination criteria based on statistical tests.
We introduce the concept of a probability space and formally define the KLdivergence. Definition 1. The probability space is a triple Ω, F , P consisting of • a non-empty set Ω (sample space), • a σ-algebra F ⊆ P(Ω), E ∈ F is called an event • a probability measure P : F → [0, 1]. Definition 2. Two probability spaces Ω, F , P i , i = 1, 2, are called absolutely continuous with respect to each other, in symbols The Radon-Nikodym Theorem allows a representation of a probability measure via a measurable probability density function. Theorem 1. (Radon-Nikodym) Let λ be a probability measure such that λ ≡ P 1 , λ ≡ P 2 . Then λ-measurable functions f i : Ω → R, i = 1, 2, called generalized probability densities, exist which are unique up to sets of measure zero and non-negative, such that for all E ∈ F .
A proof of this theorem can be found e.g in [12].
In the following we use X for the generic variable and x for a specific value of X. If H i , i = 1, 2 is the hypothesis that X is from the statistical population with probability measure P i , the mean information for discrimination in favor of H 1 against H 2 given x ∈ E ∈ F , for P 1 is given by the Kullback-Leibler divergence.
Definition 3. The Kullback-Leibler (KL) divergence is given by with When E is the entire sample space Ω, we shorten the notation to I(P 1 : P 2 ). For discrete sets E the integral is substituted by a sum.
For details we refer to [28].
Now assume that the sample space Ω is split into two disjoint sets E 1 and We define a statistical test procedure to choose between hypotheses H 1 and H 2 by accepting H 1 if x ∈ E 1 and accepting Assuming that one of the hypotheses has to be true we treat H 2 as the null hypothesis and call E 1 the critical region. The following wrong test decisions can occur.

Definition 4.
Incorrectly accepting H 1 although H 2 is true is called the type I error. The probability that this error occurs is given by Definition 5. Incorrectly accepting H 2 although H 1 is true is called the type II error. The probability that this error occurs is given by We assume that the test is repeated n-times and denote by O n a sample of n independent observations. O 1 represents a sample of a single observation. β n is defined as the corresponding probability of an error of type II which depends on the number of independent observations and the splitting of the probability space Ω into disjoint sets E 1 and E 2 .
The following theorem demonstrates an asymptotic relation between the KLdivergence and the minimum possible probability β * n of an error of type II with respect to all possible splittings E 1 ∪E 2 = Ω with given α = Prob(x ∈ E 1 |H 2 ) = P 2 (E 1 ) [18].
A proof of this theorem is given in [18,28].
Assuming probability models for the outcome of a data measurement experiment depending on experimental design parameters ξ ∈ Ξ ⊂ R d , this theorem justifies the KL-divergence to be an appropriate objective functional for model-based computation of an optimal experimental design for discrimination between model hypotheses. For a design with the largest possible value of I the asymptotic probability β * n of encountering an error of type II becomes minimal with respect to all possible splittings E 1 ∪ E 2 = Ω with given α 0 . We indicate the dependency of the KL-divergence on the design by I(P 2 : P 1 , O 1 ; ξ). Our aim is to derive an algorithm to calculate the optimal designξ ∈ Ξ such that An extension of the case to test a simple null hypothesis against a simple alternative hypothesis to the more general case of both hypotheses being composite is generally of interest. This includes the situation to test whether given measurement data can be explained best by one out of a finite set of probability models based on measures P r1 , r 1 ∈ {1, ..., M 1 } parametrized by parameters θ r1 ∈ Θ r1 where Θ r1 ⊂ R pr 1 is the set of all possible parameter values to parametrize P r1 and M 1 is the cardinality of the set of probability models, against the hypothesis that the measurement can best be explained by another one out of a second finite set of probability models based on measures P r2 , parametrized by parameters θ r2 ∈ Θ r2 where Θ r2 ⊂ R pr 2 and r 2 ∈ {1, ..., M 2 }.
we can get a robust worst case estimate of an optimally discriminating design for the case of composite null and alternative hypothesis. In practical applications a simple strategy to sort different probability models into null and alternative hypothesis would be to first rank all models according to the existing measurements and then set the best model as null hypothesis and the others as alternative hypothesis. The development of a suitable and efficient strategy is subject to further work. It should be noted that the presented criteria is not symmetric, i.e. the null hypothesis is favored. In the case that both hypotheses might be equally reasonable we suggest to use the symmetrized version of the KL-divergence, i.e.

Derivation of the optimal experimental design criterion
In this section we derive a numerically computable optimization objective functional based on the framework of KL divergence. The derivation is motivated by the requirements of biological in vitro time series experiments modeled by kinetic ODE systems. In most situations such experiments are time and cost consuming. Therefore a central issue is to get the most information out of a single time series data measurement experiment taking place within a given fixed time span [0, T end ]. This means that in an optimal experimental design the most informative measurement time points for one measurement run have to be calculated in such a way that only one measurement at one time point can be performed. Often, an experiment cannot produce measurements in a time continuous way. Therefore we assume that there has to be a minimal time span ∆T for the separation of subsequent measurement time points. This contrasts to the usual approach to associate weights to a discrete or continuous time design scheme, see e.g. [5]. Additionally, the initial species concentrations of the participating species should be chosen in a most discriminating way. A commonly used experimental practice is to combine kinetic time series measurements with perturbation stimuli like external adding of species quantities. From the model discrimination point of view the optimal time point of perturbations and the optimal species quantities to be added should be determined. We further assume that a measurement cannot be done at the same time as a perturbation.
In the following we translate these experimental conditions into a statistical model. Given the measurement time-vector t ∈ R n + with entries t i for the n measurement time points t i , i ∈ {1, ..., n} such that t i+1 ≥ t i , the "internal" model response vectors y i rj at measurement time t i for the parametrized probability measures P rj , r j ∈ {1, ..., M j }, j ∈ {1, 2} are given by where y rj (t i−1 , t i , y i−1 rj + c i−1 , θ rj ) ∈ R mr j are the solutions of the initial value problems with initial state y rj (t i−1 ) := y i−1 rj + c i−1 at end time t i where t 0 := 0 and c 0 := 0. The vectors c i ∈ R mmax denote species quantities the system can be perturbed with at time points t i where m max is the maximum dimension of the models, i.e.
f rhs rj (·, ·) are the right hand side functions of the ODE models. y I ∈ R mmax denotes the initial species concentration vector of the entire experiment which is for all models the same, i.e. y 0 rj := y I . By m min we denote the minimal dimension of the models, i.e We do not assume that m min = m max therefore for a model with m rj < m max the "redundant" entries in c i and y I are "ignored".
Let O ti denote an observation at measurement time point t i . By assuming that the measurements at successive time points t i are independent with normally distributed error vectors ǫ i rj ∈ R m with zero mean and variance functions v rj (O i rj , t i , θ rj ) 2 we get for the regression models the model probability densities f rj (·; ·) at measurement time point t i given by We generally allow different error models. The error models might dependent on the expected observations O i rj , the time t i and possibly on parameters θ rj .
For the sake of notational simplicity we define For a full measurement run containing n measurement time points we get the probability density models However, by assuming such a model probability distribution we still allow that two measurements are separated by a time span less than ∆T .
To overcome this problem we extend the probability spaces Ω i = R o of a measurement at one measurement time point by one-element-containing sets N i toΩ whereΩ i is the disjoint union of Ω i and N i . The element of the set N i with measure P (N i ) ∈ [0, 1] represents the event "no measurement", i.e. O ti ∈ N i ⇔ "no measurement performed at time point t i ".
In order to derive measures onΩ i , i = 1, ..., n that allow for a density function representation according to the Radon-Nikodym theorem (Theorem 1), we introduce the Heaviside-function with By use of this Heaviside-function and σ-algebras F i , where F i contains the Lebesgue measurable sets on Ω i and additionally the union of these with the set N i , we define probability spaces Three cases have to be distinguished: For case two with E i ⊂ N i we set By introducing these modifications the probability models based on measures P i,rj do not depend on measurements which are performed in less than ∆T time after the previous measurement any more.
To take into account that a species concentration perturbation can only be applied if no measurement is done at the same time, the same procedure is repeated with the additional Heaviside-function The measuresP i,rj are defined in the same way as above replacing H(t i ) by For specific r 1 and r 2 inserting the probability modelsP i,r1 and respectivelỹ P i,r2 into the KL-divergence (Definition 3) using λ i :=P i,r2 and the additivity of the KL divergence for independent events one gets the following expression where c n := 0. With log(1) = 0 this simplifies to By inserting the normal distribution (16) in (28) one gets This is equivalent to with and where ⌊O⌋ k gives the k − th entry of the observation vector O. A i k reduces using the well known moments of the normal distribution to This further simplifies to Substituting A i k back into (30) we get This reduces to This criterion has to be maximized with respect to the initial concentration vector y I , the measurement time point vector t and the system perturbation vector c, thus ξ := (y I , t, c) ∈ Ξ ⊂ R d . For our optimal experimental design we generally start with a large number of measurement time points. By use of the Heaviside functions the number of measurement time points is reduced such that for t i − t i−1 < ∆T the corresponding measurement time point is "turned off".
These Heaviside-functions H(·) and H(·) can be replaced by any appropriate continuously differentiable switching functions with range space [0, 1].
It should be noted that we assume that we have the same time discretization for measurements of different species and the addition of further species quantities. This assumption is practical especially for the application to in vitro experiments performed by biologists. For introducing arbitrary generic controls we need a more general formulation of time schemes, i.e. simultaneously time schemes which are independent of each other. One is associated with the controls, others may be associated with distinct observables which might be measured indepently. The incorporation to the presented framework is subject of further work.

Solution of the max-min optimization problem
We formally state now the experimental design optimization problem P Θ : subject to min θr 1 ∈Θr 1 θr 2 ∈Θr 2 The auxiliary variable τ ∈ R is used to transform the maxmin optimization problem (8) to a maximization problem with an infinite number of inequality constraints. The remaining constraints model the feasible range of experimental setups.
To solve optimization problem (36) numerically by applying efficient derivative based algorithms we replace the Heaviside functions H(t i ) and H(c i ) in (35) by continuously differentiable approximations, parametrized hyperbolic tangent functions of the form The parameters a 1,2 characterize the width of the transition region between 0 and 1. The parameters b 1,2 determine the center of the transition region (see Figure 1). By setting the parameters in an adequate way arbitrarily close approximations of the Heaviside functions can be generated. A different approach to handle the discontinuous Heaviside functions would be to introduce binary variables and treat the resulting problem as Mixed Integer Nonlinear Programming problem. The drawback of this approach is that its solution can become very expensive. There seems to be little theoretical work in literature on Mixed Integer maxmin problems and an efficient solution strategy is not obvious in that case. In literature, problems as (36) fall into the class of semi-infinite inequality and equality constrained optimization problems (SIECP) [38]. Several methods to solve such SIECP problems are available, an overview can be found in [21,38]. We choose the method of outer approximation [41,42,38], whose origin can be traced back to cutting plane methods for convex problems [38]. This approach is beneficial in the presence of a complex inner problem, in our case the robustification against the parameters θ r1 and θ r2 . The outer approximation algorithm iteratively solves discretized finite counterparts P Θ N of the semi-infinite problem P Θ in each step N , successively refining the discretization Θ N until a sufficient approximation of the original problem P Θ is reached. For problem (36) this means that in each iteration of the outer approximation scheme Θ r1 and Θ r2 are replaced by finite approximations Θ N r1 and respectively Θ N r2 with Θ N := {( Θ N r1 , Θ N r2 )|r j ∈ {1, ..., M j }}. This relation between the semiinfinite problem and an infinite sequence of finite problems can be formalized in the theory of consistent approximations and epi-convergence [35,36,37,38].
We use a modified version of Algorithm 3.6.4 in [38] where "Step 1.", the calculation of augmenting vectorsθ N +1 is realized by with ξ N denoting a locally optimal design of the previous problem P Θ N . The algorithmic scheme is as follows: Algorithm 1.
• Step 2. Calculate approximate solution of P Θ N such that and the equality and inequality constraints in problem P Θ N are fulfilled up to ǫ N .
are Lipschitz continuous on bounded sets with respect to ξ and θ rj and Θ rj are compact any accumulation point of Algorithm 1 fulfills the generalized KKT conditions, compare Theorem 3.6.5 in [38].
To calculateθ N +1 rj in Step 1. of Algorithm 1 we use on heuristic base a simple random search approach coupled to a local optimization method, i.e. we have randomly generated P different start values in Θ rj , from which we have started the local optimization method for parameter estimation. The best value out of the P trials is chosen to augment the set Θ N rj . Of course there are more sophisticated approaches to search for a global minimum for a review see e.g. [4], but at this point an effective calculation of Step 1. of Algorithm 1 was not our primary goal. For the local parameter optimization we use the same optimization method as for Step 2. in Algorithm 1.
In our implementation we use a fixed ǫ at the desired final accuracy to solve problem P Θ N in every loop of Step 2. of Algorithm 1, i.e. ǫ N = ǫ, N > 0. In that way Step 1. of Algorithm 1 gives a worst case estimate of the KL divergence I(P r2 (θ r2 ) : P r1 (θ r1 ), O 1 ; ξ N ) for the current design ξ N with respect to θ rj . Therefore for practical application the algorithm can be stopped if the worst case estimate of KL divergence for the current design ξ N is big enough although no local optima might be achieved during optimization. As stopping criterion of Algorithm 1 we use: where δ is a small positive constant, then consider (θ N +1 This stop criterion has also been used in [42,39]. We call the distance ∆ RG given by (41), robustification gap.

Numerical solution of the problem P Θ N
We have implemented the resulting optimization problem in a multiple shooting setup (see for example [45,14,13]). In our multiple shooting setup the whole integration interval [0, T end ] is subdivided into several subintervals by introducing auxiliary multiple shooting node variables s rj,i,l , j ∈ {1, 2}, r j ∈ {1, ..., M j }, i ∈ {1, ..., n}, l ∈ {1, ..., N }, on each of which an independent initial value problem is solved. Each end point of a subinterval corresponds to one measurement time point. Matching conditions which enter the optimization problem as additional equality constraints assure continuity of the state trajectory from one subinterval to the next.
A graphical scheme of the multiple shooting setup is shown in Figure 2. Instead of evaluating the objective functional (35) by use of the values y i rj , given by the solution of the initial value problem (11), depending on the parameters θ rj , (35) is evaluated by use of the auxiliary multiple shooting node variables s rj ,i,l replacing the values y i rj with s rj,i,l respectively. The dependency of (35) on s rj ,i,l is indicated by I(s r1,·,l , s r2,·,l ).
We have implemented this problem within the interior point optimization package IPOPT [50,51], using the linear solver MA27 [23]. Usually for the solution of the KKT system within the direct multiple shooting ansatz the linear system is solved utilizing tailored structure-exploitation, e.g., condensing or high rank updates. See for example [30]. Since speed aspects are not our primary concern we rely on the sparse solver MA27 instead of developing a tailored solver for this problem class at the current stage. All derivatives up to second order, which are used for the calculations of the Hessian needed for a robust performance of IPOPT, are calculated by automatic differentiation using CppAD, [9,8]. For the solution of the differential equations within the optimization problem, which are commonly stiff in chemical and biochemical applications, we have implemented a fully variable step, variable order (order 1 to 6), Backward differentiation formulae (BDF) method, based on Nordsiek array polynomial interpolation similar to the EPISODE BDF method by Byrne and Hindmarsh [16], but with the step size selection strategy of Calvo and Rández [17]. For the generation of sensitivities we have adopted the sophisticated principles of internal numerical differentiation developed by Albersmeyer and Bock [2,1] in forward and adjoint mode. The idea of this principle is instead of calculating the sensitivities by use of the sensitivity differential equation, to directly differentiate the BDF integration scheme by automatic differentiation, which we implemented using CppAD [9,8]. According to some notes in the PhD thesis of Albersmeyer [1] we also have implemented the possibility to control the step size scheme not only by the local truncation error of the nominal trajectory but as well by the local truncation error of the sensitivities generated by the forward mode of automatic differentiation with respect to the sensitivity differential equation, which has shown by numerical experience to improve the robustness of the optimization approach. A different approach would be to use collocation, i.e. to incorporate a full discretization of the ODEs into the optimization problem, see e.g. [11]. Since the kinetic ODE systems in the focus of our applications are usually stiff, we prefer adaptive time integration.

Stabilizing homotopy method for subsequent P Θ N+1
By solving the subsequent optimization problems P Θ N +1 with an interior point code like IPOPT [50,51] initialized with primal and dual variables of the previous problem or with primal variables only, one often observes that the new solution may differ significantly from the previous. This is due to the fact that the solution of the previous problem P Θ N is infeasible for P Θ N +1 and thus the algorithm tries to find a feasible state before it proceeds to find a new optimum. This behavior is not desired in the context of an outer approximation algorithm, because convergence of the algorithm may be slowed down significantly. This circumstance originates from a jumping between vicinities of distinct local maxima of problem (36). The discretization Θ N of the robustification space may not be equally adequate for different local maxima. To overcome this problem we have implemented a heuristic homotopy method to gradually introduce the additional constraints g r1,r2 (τ, ξ) N +1 := I(s r1,·,N +1 , s r2,·,N +1 ) − τ ≥ 0, (46) of problem P Θ N +1 . We replace g r1,r2 (τ, ξ) N +1 bỹ with homotopy parameter κ ∈ [0, 1] and ρ is a constant which has to be set such thatg r1,r2 (τ, ξ; κ) N +1 are inactive for κ = 0 at the initial design ξ N . We choose ρ to be K is a save guard factor we set empirically to K = 1.4, which worked well in practice for our examples. For κ = 0 the augmented optimization problem should be easily solvable within a few iterations by performing a warm start from the solution of the previous problem. By increasing the homotopy parameter to κ = 1 the additional constraint is gradually introduced, which leads to a sequence of easily solvable subproblems whose solutions stay in the vicinity of the solution of the previous problem P Θ N . A similar homotopy strategy can be found e.g. in [40].

Numerical results
We have applied the algorithm developed in Section 2 and 3 to two example problems for which we present results in the following section. For the purpose of illustration we restrict ourself to the case that each hypothesis comprise only one model whereby we assume that only the alternative hypothesis is composite. We also assume that each species is "directly" measurable. We treat model 1 as null hypothesis and model 2 as alternative hypothesis.

Discriminating design for two models describing glycolytic oscillations
In the first test case for model discrimination we implemented the following models for glycolytic oscillations as described in [20].
Model 1 is an allosteric enzyme model with positive feedback under cooperativity and linear product sink. The differential equations for model 1 are given by Model 2 is an allosteric model with positive feedback in the absence of cooperativity and the product sink is represented by Michaelis-Menten kinetics. The differential equations for Model 2 are given by . α 1,2 denotes the species concentration of the substrate γ 1,2 that of the product. For both models the inflow parameter ν is the same and fixed to the value ν = 0.22. It represents the inflow of substrate to the experimental system, a CSTR (continuously stirred tank reactor). The parameters σ, q 1 , k s and L 1 of model 1 are regarded as known. Their values are given in Table 1, the parameters q 2 , r s , µ and L 2 of model 2 are regarded as unknown and subject to robustification. For the permitted parameter range see Table 1. For simplicity we consider the homoscedastic case with equal variances, i.e. v 1 = v 2 = σ 2 . In this case I(P 1 : P 2 , O 1 ) reduces to, For this test case the homotopy strategy as presented in Section 3.2 is only applied if the robustification gap ∆ RG < 0.1, then the successive problem P Θ N +1 is calculated by use of the homotopy strategy with 30 homotopy steps, i.e. κ h = h/30, h ∈ {1, ..., 30}. Otherwise the problem P Θ N +1 is solved without homotopy strategy. For each subsequent problem P Θ N +1 the solution of problem P Θ N is used as initial guess.
We first present a robust design without the possibility to disturb the system by adding species at later time points. The design is calculated within a fixed time window i.e. T end = 400. 100 equally spaced possible measurement points are defined in the initial state of the optimization procedure, the distance vector ∆t between the time points is subject to design and each entry is restricted to ∆t i ∈ [10 −7 ,  Table 2.
A plot of the functions α 1 , α 2 and γ 1 , γ 2 in the initial state and for the solution of problem P Θ 1 are shown in Figure 3. A plot for the same functions with the same solution design as for problem P Θ 1 after the next robustification step is shown in Figure 4. The final design is also shown in Figure 4. A plot of the  Figure 3: The model functions α 1 , α 2 and γ 1 , γ 2 are shown before the optimization procedure (left) and after the optimization procedure of problem P Θ 1 (right) for the glycolytic design setup without the possibility to disturb the system. One square represents one measurement time point.
robustification gap ∆ RG and as well for the objective value of problem P Θ N for each iteration N of Algorithm 1 are shown in Figure 5. A selection of design variables as solutions of problem P Θ N is shown in Figure 6(left).
In a second scenario we additionally allow for species perturbations. In this A plot of the functions α 1 , α 2 and γ 1 , γ 2 in the initial state and for the solution of problem P Θ 1 are shown in Figure 7. A plot for the same functions with the same solution design as for problem P Θ 1 after the next robustification step is shown in Figure 8. The final design is also shown in Figure 8. A plot of the robustification gap ∆ RG and the objective value of problem P Θ N for each iteration N of Algorithm 1 are shown in Figure 9. A selection of design variables as solutions of problem P Θ N is shown in Figure 6(right).

Discriminating design for two models describing signal sensing in dictyostelium discoideum
The second test case is the discrimination of two models describing the chemotactic response in the amoeba Dictyostelium discoideum as presented in [33] using the framework presented in Section 2 and 3. The two models describe the adaption mechanism observed when amoebae encounter the chemoattractant cAMP [31], see figure 10. For both models, a chemotaxis response regulator R gets activated (R * ) by an activator enzyme A, when a cAMP ligand S appears. But the deactivating mechanism determined by the interaction with an inhibitor molecule I differs for both models. Both models comprise mass action kinetics in form of ODE.   Figure 6: A selection of design variables as solutions of problem P Θ N for the glycolytic design setup without the possibility to disturb the system (left) and with the possibility to disturb the system (rigth) are shown.
In model 1 the activator enzyme as well as the inhibitor enzyme are regulated by the external signal , which is proportional to the cAMP concentration S. The overall model in this case is given by, where k −a , k a , k −i , k i1 , k r and k −r are the mass action rate constants and R T := R * + R is the total amount of the response regulator.
In model 2 the inhibitory molecule I is activated through the indirect action of activator A instead of direct activation by sensing ligand binding. The overall model in this case is given by, where k −a , k a , k −i , k i2 , k r and k −r are the mass action rate constants and R T := R * + R is the total amount of the response regulator. For modeling details we refer to [33]. We have extended these systems of ordinary differential equations by an additional state corresponding to the cAMP ligand S withṠ = 0. By allowing species concentration perturbations c only to the state S we can mimic a piecewise constant control of the system by the cAMP ligand S.
The experimental design parameters are the initial species concentrations of the four states namely, A I , I I , R I , S I , the measurement time points t and the species concentration perturbation c with respect to S. We discard the condition that either a measurement or a perturbation can be performed since in that setting by use of the perturbations c we mimic a piecewise constant input control S and therefore that restriction seems unnatural. Again for simplicity we consider the homoscedastic case with equal variances i.e. v 1 = v 2 = σ 2 , where I(P 1 : P 2 , O 1 ) reduces now to The parameters k −a , k a , k −i , k i1 , k r , k −r and R T are regarded as known and fixed, their values are given in Table 3. Parameter k i2 is regarded as un-  known and subject to robustification. The range of the parameter k i2 is set to k i2 ∈ [0, 2]. The optimal design is calculated within a fixed time window with T end = 100. 100 equally spaced possible measurement points are defined in the initial state of the optimization procedure. The distance vector ∆t between time points is subject to design and each entry is restricted to ∆t i ∈ [10 −7 ,  Table 4. With these design conditions we start the optimization procedure twice. First by use of the homotopy strategy for successive problems P Θ N +1 with 10 homotopy steps. Since the "discriminating power" of the experimental setup is very low in this case, i.e. the deviation between the two models is small, we plot the distance functions (S 1 − S 2 ), (A 1 − A 2 ), (I 1 − I 2 ) and (R 1 − R 2 ) for the initial state and  for the solution of problem P Θ 1 in Figure 11. A plot for the same functions with the same solution design as for problem P Θ 1 after the next robustification step is shown in Figure 12. The final design is also shown in Figure 12.
A plot of the robustification gap ∆ RG for each iteration N of Algorithm 1 is shown in Figure 13 (left). A plot of the objective value of problem P Θ N for each iteration N of Algorithm 1 is shown in Figure 14 (left). A selection of design variables as solutions of problem P Θ N is shown in Figure 15 (left).
Second we calculate the design without the homotopy strategy. We experience huge jumps in the final objective value of problem P Θ N for subsequent iterations N of Algorithm 1. This is due to the fact that the final design of the former problem P Θ N is an infeasible starting point for the successive problem P Θ N +1 in the interior point solution strategy. First the optimizer tries to force the iterates back into the feasible region and afterwards the new central path leads to a different design. For this case a plot of the robustification gap ∆ RG for each iteration N of Algorithm 1 is shown in Figure 13    and (R 1 − R 2 ) are shown before the optimization procedure (left) and after the optimization procedure of problem P Θ 1 (right) for two models describing signal sensing in dictyostelium discoideum. One square represents one measurement time point. P Θ N for each iteration N of Algorithm 1 is shown in Figure 14 (right). A selection of design variables as solutions of problem P Θ N is shown in Figure 15 (right). As one can clearly see, the homotopy strategy helps to considerably stabilize Algorithm 1.

Conclusion
We present a framework for the robust computation of optimal experimental designs for the purpose of model discrimination. The theoretical framework as well as the numerical realization by utilization of an outer approximation algorithm are worked out. A strategy for the numerical stabilization of the algorithm by use of a homotopy approach is suggested. The optimization procedure is successfully exemplified on two biological model systems. In our examples we clearly found that the homotopy approach is significantly superior to a cold start of successive design problems P Θ N +1 . For the first test case, the discrimi- and (R 1 − R 2 ) are shown for the same solution design as for problem P Θ 1 after the next robustification step and for the final design (right) for two models describing signal sensing in dictyostelium discoideum. One square represents one measurement time point. nation of two models describing glycolytic oscillations, the outer approximation scheme completely fails to reach the desired accuracy δ without homotopy strategy. For the second test case, the discrimination of two models describing signal sensing in dictyostelium discoideum, the outer approximation scheme also fails without warmstart, however the homotopy strategy also works with only two homotopy steps (not presented in this paper). This indicates the need of a step size strategy for reasons of efficiency which will be a next step in our work.  Figure 13: The robustification gap ∆ RG is plotted versus the number of iterations N of Algorithm 1 for the setup with two models describing signal sensing in dictyostelium discoideum, with homotopy strategy (left) and without homotopy strategy (right).    Figure 15: A selection of design variables calculated as solution of problem P Θ N is plotted versus the number of iterations N of Algorithm 1 for the setup with two models describing signal sensing in dictyostelium discoideum, on the left with homotopy strategy and on the right without homotopy strategy.