This model-based design of experiments (MBDOE) method determines the input magnitudes of an experimental stimuli to apply and the associated measurements that should be taken to optimally constrain the uncertain dynamics of a biological system under study. The ideal global solution for this experiment design problem is generally computationally intractable because of parametric uncertainties in the mathematical model of the biological system. Others have addressed this issue by limiting the solution to a local estimate of the model parameters. Here we present an approach that is independent of the local parameter constraint. This approach is made computationally efficient and tractable by the use of: (1) sparse grid interpolation that approximates the biological system dynamics, (2) representative parameters that uniformly represent the data-consistent dynamical space, and (3) probability weights of the represented experimentally distinguishable dynamics. Our approach identifies data-consistent representative parameters using sparse grid interpolants, constructs the optimal input sequence from a greedy search, and defines the associated optimal measurements using a scenario tree. We explore the optimality of this MBDOE algorithm using a 3-dimensional Hes1 model and a 19-dimensional T-cell receptor model. The 19-dimensional T-cell model also demonstrates the MBDOE algorithm’s scalability to higher dimensions. In both cases, the dynamical uncertainty region that bounds the trajectories of the target system states were reduced by as much as 86% and 99% respectively after completing the designed experiments in silico. Our results suggest that for resolving dynamical uncertainty, the ability to design an input sequence paired with its associated measurements is particularly important when limited by the number of measurements.
Since experiments can be expensive and time consuming, it is important that they are planned to generate useful data. Traditional design of experiments is a well established field and has led to many advances in biology and medicine. The data obtained from strategically designed experiments has facilitated the creation of mathematical models that relate experimental stimuli to measurable outcomes. These models typically describe the system’s input-output relationship but fail to capture or encode knowledge of the system’s internal mechanisms and processes. Mechanistic and semi-mechanistic mathematical models encode the current understanding of the internal processes of the biological system even though many of these internal states or species are not directly measurable. These mechanistic models can be used to support optimal experiment design that considers the current knowledge of the system interactions and practical experimental constraints. In recent literature this type of experiment design has been referred to as model-based design of experiments (MBDOE). MBDOE produces experiments meant to reduce some measure of uncertainty in the associated model while respecting cost, time and resource constraints. Most MBDOE strategies can be categorized by three types of objectives: (1) reducing model parameter uncertainty [1–7], (2) discriminating among possible models [8–13], and (3) reducing dynamical uncertainty [14–17]. This work advances current abilities to design experiments to resolve the trajectories of target states of a biological system model, thereby reducing its dynamical uncertainty. Many of the MBDOE strategies that support reduction of parameter uncertainty and model discrimination rely on linear approximations that are locally optimal to design an experiment by optimizing a criterion of the Fisher Information Matrix (FIM) [15,18–21]. Such techniques use the local sensitivities of parameters to design an optimal experiment which requires an initial estimate of the unknown parameters. Most biological system models are not well characterized, as data is limited and noisy, so initial estimates of the model parameters are inaccurate. Furthermore, biological models are typically nonlinear, so a poor initial parameter estimate is likely to result in a local minimum, which may result in sub-optimal experiment design. Alternatively, the Sigma Point method refines estimates of the characteristic values of the parameter statistics [22,23]. This method is computationally efficient and gives a better approximation of the covariance matrix when compared to the FIM-based and boot-strap methods. However, this method is not applicable when the biological system is not well characterized with existing experimental data. These local methods for MBDOE are limited to scenarios when initial estimates for the model parameters or their distributions are already fairly well known. To overcome the shortcomings of local MBDOE strategies, global approaches that consider the entire uncertain parameter range have been developed. Global optimization has been used to improve the FIM-based design [24] and local deterministic optimization of the likelihood with multi-start initial parameter estimates [25] have been previously employed. Some Bayesian MBDOE methods use Monte Carlo generated estimates of the parameter confidence intervals to design a global optimal experiment [16,26–30]. A comparison of the FIM-based and Bayesian-based MBDOE methods by Weber et al. [31] shows that Bayesian methods give more accurate and informative designs. However, all of these global methods are computationally prohibitive for large numbers of highly uncertain parameters and mathematical models that are computationally intensive. The limitation due to computation has been partially abrogated by methods that utilize surrogate models to approximate the mathematical model of the system. Herein, as in our previous work [14,15], we use a sparse grid interpolation tool to approximate the biological model. To find the optimal design a search has to be performed over all of the experimental factors that define the feasible experimental perturbations and potential measurements. A complete exploration of the full design space over the uncertain parameter space for an optimal experiment design is impractical due to the combinatorial explosion of possible experimental factors evaluated for all possible model parameterizations. To solve the problem, most MBDOE strategies constrain the design space to a subset of the possible experiments. This reduces the design space by assuming a small predefined set of possible measurement points, measurable species and/or input levels [9,31–33]. Herein, we define the input to be the exogenous experimental factors, perturbations and/or stimuli that are applied to the biological system. Other designs have considered finding an optimal input sequence while the times and measurement species are defined a priori [2,10,24,34]. To consider both optimal measurements and an optimal input sequence, some strategies [16,35] evaluate the optimality of a pre-defined set of specified inputs and measurement points. All of these approaches may result in sub-optimal experiment designs if the best experiment is not among the specified experimental options. Ideally, experiment design determines the optimal solution within a feasible space defined by all the possible experimental measurements and inputs. Since different inputs can elicit dramatically different dynamics for some nonlinear systems, evaluating experiments for optimal inputs is not necessarily a trivial extension of the global techniques for identifying optimal measurements. Herein, we describe and evaluate a computationally efficient and tractable method for performing global MBDOE to define optimal experimental inputs and measurements to reduce the uncertainty in continuous nonlinear system dynamics. The method utilizes the sparse grid surrogate model technique first employed for global experimental designs presented in [14,15], which selected only optimal measurements and not optimal inputs. We proposed an earlier version [36] of our input design MBDOE algorithm that used a combination of scenario trees and sparse grid interpolation to screen the input design space for optimal inputs that will generate diverse dynamics of the uncertain nonlinear system. This paper is an extension of that work. Herein we utilize a greedy search to determine the input magnitude to be applied using a method that is similar to that in [37] which determined near-optimal measurements in polynomial time. Our MBDOE algorithm is made even more computationally efficient and robust than our previous sparse grid based strategies due to the uniform sampling over the dynamical space via careful selection of representative parameters, the optimization criteria for optimal experiment design for a target system, and the repeated update of probability weights over the dynamical space using predicted data. Our new MBDOE method develops a strategic experiment design that yields highly informative data that can be used to constrain the system dynamics to dynamical uncertainty regions. In the methods section, we describe the experiment design problem and define the methods used to solve the problem as sequential optimizations to first select an optimal input sequence and, subsequently, specify the associated optimal measurement pairs (species and time points). In the results section, we explore the optimality of the derived experiment design solution using a small 3-dimensional Hes1 oscillatory model and demonstrate the scalability of this MBDOE strategy with a 19-dimensional T-cell receptor model. The global applicability, support of hybrid experiment designs, and limitations of this MBDOE approach are discussed to conclude the paper.
This MBDOE algorithm is designed to resolve dynamical uncertainties associated with the biological system by reducing the predicted variance in the system dynamics. The system is described by nonlinear ordinary differential equations of the general form: x ˙ = f (x, u, θ, θ ′, t), x (t 0) = x 0 (1) where f is a deterministic smooth function of x ∈ ℝ n x, the states of the system, θ ∈ ℝ n p, the known model parameters, θ ′ ∈ ℝ n p ′, the unknown model parameters, and u ∈ ℝ n u the exogenous inputs. Using control vector parameterization (CVP) [38], we represent u = [u (τ1), …, u (τN) ] ∈ ℝnu∗N, as vector-valued function of the system inputs that can change value at each discrete time point, τj, where j = 1, …, N. There is uncertainty in the dynamics of the system which is generated by the uncertainty in the values of the unknown model parameters, θ′ ∈ Ω where Ω is a compact set. It is common in biological systems that not all states of this system are experimentally measurable. All feasible measurements for this system are modeled by: y = h (x, u, θ, θ ′, t), (2) where h is a smooth vector-valued function that relates to the system’s internal states described in Eq (1) and y ∈ ℝ n m. We abbreviate Eq (2) by the notation y (u, θ′, t). The user defines the number of measurements that they would like to have designed, K, as well as the set of possible times for measurement, 𝕋, by either specifying discrete time points or defining a time resolution, δt, between TI and TF, the initial and final model simulation times, respectively. Our MBDOE algorithm designs experiments that reduce dynamical uncertainty in a target system. This target system is a subset of the system states termed, target states, xT, when stimulated by a target input, uT. This approach employs the concept of a target system also used to evaluate an alternative MBDOE strategy [16]. We consider dynamical uncertainty in the measurable states, y, and the target states, xT, generated by the uncertainty in the unknown parameters, θ′. As a result of limited and noisy data, often biological system models will fit the data with a wide range of parameters values and associated dynamics. Therefore, for an optimal design, we need to consider all the possible data-consistent dynamical scenarios which are a function of the parameter space, Ω. The optimal experiment, D* ∈ 𝔻, defines the optimal input stimulation and optimal measurement pairs that resolve the target state dynamics. This design is defined by a piecewise input sequence, u* (Eq (3) ), and associated measurements M* = { (mk, tk): k = 1,2, …K} ∈ 𝕄, which defines K measurement pairs by the index of a measurement species mk and its corresponding measurement time tk. u * = { u 1 * τ 1 ≤ t ≤ τ 2 ⋮ ⋮ u j * τ j ≤ t ≤ τ j + 1 ⋮ ⋮ u N * τ N ≤ t ≤ T F (3) We choose D* as the experimental design that maximizes a measure of information gained from an experiment D * = argmin u ∈ ℝ n u * N, M ∈ 𝕄 γ (u). (4) The target state dynamical uncertainty (TDU), γ (u), is quantified by the sum of the maximum variance in each target state that results when analyzed assuming the input sequence, u, has been applied and measurements, M, exist to constrain the uncertain parameter space: γ (u) = ∑ i max (Var (x T i (t) ) ∣ (u, M) ) (5) where the maximum variance of the ith target state, max (Var (xT i (t) ∣ (u, M) ), is chosen across the simulated time t, from the initial, TI to the final time, TF. The ideal solution for the above optimization problem simultaneously solves for the optimal input vector, u*, together with the optimal measurements, M*. This is generally a computationally intractable problem due to the combinatorial explosion of all possible selections for the experiment design complicated by the expense of evaluating the design with the model for all possible uncertain parameter values. Our approach proposes a computationally efficient method to approximate both the optimal input vector and associated optimal measurement points by breaking the problem into smaller computationally feasible optimization problems. The first optimization problem solves for the input vector using a greedy method by minimizing the value of TDU at each iteration assuming a single best measurement is taken: u * = argmin u j ∈ ℝ n u γ (u), j = 1, ⋯ N. (6) The optimal input vector is then used in the second optimization to determine the multiple optimal measurements points. This combined solution of input and measurements approximates the optimal experiment design defined in Eq (4). A flow-chart of the MBDOE algorithm is given in Fig 1 to display the sequence of events in determining the optimal experiment design: (a) identify representative parameters that maintain the diversity of simulated dynamics that fit the existing data, (b) determine the optimal input vector that minimizes TDU, and (c) specify associated multiple measurement pairs given an optimal input vector. The first step of this process screens the uncertain parameter space, Ω, to identify the space of acceptable parameters, ΩA. The acceptable parameters are those that support model simulations to fit the available data as: Ω A = { θ ′ ∈ Ω ∣ log 10 (1 + ∑ θ ′ ∈ Ω (y ˜ i (u, θ ′, t) - y ^ i (u, t) σ i (t) ) 2) ≤ T A } (7) where y ^ i (u, t) is the mean of the experimental data for the ith model output at the time point, t, collected with the applied input u, y ˜ i (u, θ ′, t) is the corresponding model dynamics simulated with a parameter set, θ′, and σi (t) is the standard deviation of the data. In this work, we approximate y ˜ i (u, θ ′, t) using a sparse grid interpolation tool over the uncertain space Ω for computational efficiency. TA is the threshold for acceptability. This weighted least squares function has been used before to define acceptable parameters in [14,15]. The parameter space, Ω, is sampled using Latin Hypercube Sampling (LHS). If the initial screen produces fewer acceptable parameters than NA, the desired number of acceptable parameter vectors that will support the experiment design, focused grids are created [15] to improve the resolution of the grid interpolant in the acceptable regions of the uncertain parameter space. LHS is also used to sample the focused grids to find more acceptable parameters until the cardinality (ΩA) ≥ NA. In this work, we define the dynamical uncertainty region to be the region spanned by the most extreme trajectories of the acceptable parameters. We propose solving the input vector using a greedy search algorithm. The input is discretized as in control vector parameterization (CVP) by the potential admissible input times 𝓣 = {τ1⋯τN} and by potential magnitude levels. The number of possible input magnitudes is determined by the inputs bound [umin, umax] and resolution, δu: N u = u m a x - u m i n δ u + 1. (12) The algorithm for the greedy search method to select the optimal input vector is detailed in the flowchart in Fig 1 (b). The input magnitude is initially set to a base input level, u1 = ub (t) ∀t ∈ [TI TF], which is determined by the user. At each iteration, k, of the greedy search method, optimized input levels, u j *, for the admissible times, τj ∈ 𝓣k−1, are selected from all Nu possible admissible input magnitudes as described in Algorithm 1. We evaluate TDUj for each of the admissible input times τj ∈ 𝓣k−1 and associated u j * and select the input conditions, (τ j *, u j *), that minimize the value of TDU according to Eq (6) as optimal. These optimal input conditions update the current input vector as follows: u k * = { u j * τ j * ≤ t ≤ τ ′ u k - 1 * Elsewhere (13) where τ′ is an input admissible time that has previously been optimized to update the vector, u*k−1. We define 𝓢k−1, a set of all previously optimized input admissible times such that τ ′ = m i n [ t ∈ 𝓢 k − 1 ∣ t > τ j * ]. When the input vector is updated with optimal input conditions, (τ j *, u j *), the input time, τ j * is moved from the set 𝓣k−1 to the set 𝓢k−1. This process continues iteratively, and the input vector is updated until 𝓣 is empty or no input magnitude change improves the value of TDUk−1 as shown in Fig 1 (b). Algorithm 1 Algorithm for Selecting u * j Input: τj, umin, umax, δ u, 𝓢k−1, 𝓣k−1, u*k−1 Output: u j * and corresponding TDUj for p = 1: Nu do (Note: Nu = number of possible input magnitudes) 1. Define u j p = u m i n + (p − 1) δ u 2. Simulate dynamics with u = up for τj ≤ t < τ′ where τ′ = min[t ∈ 𝓢k−1∣t > τj] 3. Determine an associated measurement pair: mk ∈ nm and tk ∈ 𝕋. See Eq (15) 4. Update uncertain parameter probabilities given data for (mk, tk). See Eq (16) 5. Estimate variance of target states, xT, with target input, uT, using updated uncertain parameter probabilities 6. Calculate TDUp. See Eq (5) end for p * =argmin p∈{1, ⋯, N u } TDUp. See Eq (6) u j * = u j p * TDUj = TDUp* Our design can identify multiple informative measurements given an optimal selected input vector, u*, as shown in Fig 1 (c). We use an adaptation of the measurement scenario tree method [15] which makes multiple predictions of the data value from a given measurement point to allow the next informative point to be selected. In this work, we introduce the use of prior/posterior probability updates to minimize the bias due to outliers. The use of probabilities also enables the algorithm to work well with a smaller number of sampled uncertain parameters than previously, contributing to the computational efficiency. As in [15], a measurement scenario tree is initialized by a node which defines the first optimal measurement point. At each node, we determine the optimal sampling time point according to Eq (15) to maximize the distinguishability metric for each measurable species. Our MBDOE algorithm selects the measurement pair that minimizes the TDU according to Eq (5) as the optimal measurement pair for that node. Subsequent measurements (nodes) are chosen by considering three different scenarios (branches) from the node. These scenarios arise by predicting three possible outcomes for the measurement pair: previously these were defined by the minimum, mean, and maximum predicted measurement values. Herein, the cumulative distribution function is computed from the priors to span from the minimum to maximum predicted data values for that measurement pair. The three predicted possible outcomes are the simulated values at the measurement pairs that correspond to the 10th, 50th, and 90th percentiles. Using each of these predicted possible outcomes to approximate the data, y ^ i (u, t k), the posterior for each branch are updated for the θ′ ∈ ΩRA as in Eq (16). Given the updated branch specific priors, the measurement pair specifying the next node in each branch is determined by the optimal measurement pair selected to minimize Eq (5) assuming the estimated posterior is the prior. This process is continued until the tree structure contains at least the specified number of unique measurements. The minimum number of levels in the tree is log3 (K) +1 where K is the user-specified desired number of measurements. Evaluating all possible experiments over the design space, 𝔻, using model simulations for a large uncertain parameter space, Ω, is computationally prohibitive. The sparse grid tool offers an efficient way to approximate the system’s output dynamics using interpolating polynomials by sampling an uncertain space in a systematic way to create a surrogate model for the system. The grid approximation is then used to interpolate additional points on the uncertain space without the cost of directly simulating the model. The algorithm was coded in MATLAB (2013a). The Sparse Grid toolbox (v5. 1. 1) was obtained from http: //www. ians. uni-stuttgart. de/spinterp [40] and integrated with our model-based experiment design algorithm. We minimize the error due to inaccurate interpolants by specifying the termination conditions based on their depth or by the relative or absolute errors of 10% and 5 respectively. The max depth is set to 3 and 5 for the parameter screening and focused grids, respectively, and 5 for the input grid. We set the data-consistent acceptability threshold TA = 2, which correspond to two standard deviations. The Hes1 and T-cell receptor models for the results were numerically integrated using the stiff ode solver ode15s with the default settings. To increase computational efficiency, the sparse grid code and the mathematical models were vectorized. The code was run on four 8-core Intel Xeon 3. 4 GHz CPUs each with 16 GB of memory and running the Windows 2003 server platform with the MATLAB Parallel Computing Toolbox. Our input design MBDOE algorithm code is available upon request from the corresponding author. An upper bound on the number of model simulations used by our MBDOE algorithm can be approximated by: T e v a l = S G θ + N A ︸ (a) + m i n (N u, S G u) × N R A × N (N + 1) 2 ︸ (b) + N R A ︸ (c) (21) where SGθ and SGu are the number of model simulations (nodes) used to create sparse grid interpolants over the uncertain parameter and input space, respectively. An upper bound estimate of these values can be calculated using a Sparse Grid toolbox function, spdim (maxDepth, Dim) assuming the maximum possible depth and the dimension of the uncertain parameter or input space. The number of simulations required to build an accurate interpolant may be quite a bit lower than this approximation when the model is smooth over the uncertain parameter and/or input space. The terms of this upper bound approximation are mapped to the steps of Fig 1. The (a) term of Eq (21) is the total number of model simulations used to identify the representative parameters. The second term, (b), is the number of simulations that determine the optimal input vector. For this approximation, the sparse interpolant on the input space is only created if the estimated maximum number of nodes, SGu required is less than Nu. Furthermore, we assume a worst case scenario where the input magnitude changes at all N admissible times resulting in the fraction, N (N + 1) 2. The third term, (c), is the number of model simulations used create the scenario tree to determine the measurement pairs.
To illustrate the effectiveness of our MBDOE strategy to resolve dynamical uncertainty of a target system, we use a simple Hes1 oscillator toy model [41]. This example demonstrates how an optimal perturbation enhances the ability of the MBDOE to reduce the dynamical uncertainty. The Hes1 model was previously used to demonstrate a Bayesian design of experiments strategy [16]. This model describes the changes in the level of the Hes1 transcription factor that is important in somitogenesis. It is modeled with an ODE system that describes the changes of the levels of the Hes1 mRNA, m, the cytosolic protein, P1, and the nuclear protein, P2, as shown: d m d t = - k m + 1 1 + (P 2 P 0) h, d P 1 d t = - k P 1 + ν m - k 1 P 1, d P 2 d t = - k P 2 + k 1 P 1. (22) A description of the parameters of the model and their nominal values are listed in Table 2. The degradation rate, k, and protein transport rate, k1, are assumed to be known and their values are set to the nominal value while the rest of the parameters are considered unknown. The uncertain parameter space, Ω, is initially set to have bounds of 0. 1 and 10 times previously reported nominal parameter values. The experimental design space is 𝔻, partially defined by an input α ∈ 𝕌 that modifies the transport rate, k1, as α × k1 where 𝕌 = [0. 01,2] with a δu = 0. 05. The input magnitude can be changed at admissible times, 𝓣 = {0,2, 5,8, 10,50,100,150,200} min. In this example, ub, is set to α = 1 ∀τ ∈ 𝓣. The allowable measurement state space includes the Hes1 mRNA, m, and the total protein concentration, P1+P2. The allowable sampling time space 𝕋 = [0,300] min was defined to specify measurements with a δt = 10 min. We want the algorithm to specify a minimum of 8 distinct measurements. The experiment is designed to decrease the dynamical uncertainty on three target states, m, P1, and P2 under a target input α = 1 applied for the simulated time period. To evaluate the MBDOE algorithm in silico, we simulate a plant model for which we want to reduce the dynamical uncertainty. The plant model is initially simulated with unknown model parameters set at the nominal values shown in Table 2. We generate 6 initial measurements of both mRNA, m, and total protein concentration, P1+P2, at times t = 10,20,60 min with 10% additive Gaussian noise. Initial data identify the acceptable parameter space, ΩA, that is consistent with the data as defined by Eq (7). The resulting optimal experiment design is in Table 3. The optimal input sequence is specified by α values shown in Fig 2a and associated predicted measurements are shown on Fig 2b, superimposed on the representative dynamics for the measurable species to illustrate their connection with uncertainty. The greedy algorithm for the input selection is shown in Table 4 showing each iteration of input selection with the progressive predicted reduction in TDU values. The initial uncertainty and final uncertainty in the dynamics of the target states are superimpossed in Fig 2c to demonstrate the reduction of the region of uncertainty of the dynamics. The dynamical uncertainty regions associated with the mRNA, P1, P2 target states are reduced by 84%, 81% and 86%, respectively. In all cases, this reduced uncertainty region enclose the true dynamics of the system. The measurement pairs are shown in the scenario tree in Fig 3a. the nodes of the tree specify the measurements that maximize uncertainty in the measurable species, ξ, while minimizing the TDU, γ, under the optimal input sequence, u*. In general, the values of ξ and γ decrease as you move down the tree and the predicted dynamics are constrained. To investigate how the measurement pairs contribute to the reduction of uncertainty in the dynamics of the target system, the total, γ, and individual variances are estimated assuming each predicted measurement has been taken (moving down and from left to right along the tree, neglecting repeated measurements). These results are shown in Fig 3b and 3c. These figures show that the first five unique measurements from the tree are necessary to resolve the target system dynamics. To further explore the optimality of our MBDOE algorithm, we compare our results achieved with u* and M* to those of measurement only MBDOE. For this comparison, the measurement only MBDOE determined the optimal measurements, MT, assuming the applied input was the target input, uT. These results are summarized in Table 5 by the percentage reduction in the uncertain dynamical region for each target species and TDU. Our MBDOE algorithm outperforms all other experiment design combinations based on the TDU reduction. Although our MBDOE design is not the best in reducing the uncertain dynamical region for all target states, the results are comparable. The supplemental material contains the figures showing the uncertain dynamical regions for all measurement and input combinations in the Table. The supplemental material also contains figures that indicate the ability of our MBDOE algorithm to reduce the uncertain dynamical regions for different parameterizations of the Hes1 plant model. The computational efficiency of this MBDOE method partially arises from the use of an interpolation as a surrogate of the model to reduce numerical integration of ODEs. The sparse grid terminated with a relative accuracy of 0. 7% with an absolute tolerance of 0. 9. The algorithm required only 137 model evaluations to build an interpolant which was sampled 10,000 times to identify 3852 acceptable parameters where NA = 2000. The algorithm selected 48 representative parameters, NRA = 48, to span the distinguishable dynamics generated by the acceptable parameters. For the input design, separate 1-D sparse grids are constructed for each admissible time of input change with 9 nodes. The total number model evaluation used for the complete experiment design is 18,716 where 2137 were for finding the representative parameters, 16,531 were used for the input selection, and 48 were used for the measurement selection. Clearly the majority of the model simulations were used to determine the optimal input signal. This is exceptionally large since we allowed the input to change at 9 time points and had fine input resolution with 40 different values. For comparison purposes, the Hes1 model experiment design in [16] required more than 30,000 model simulations to determine which species should be measured to provide the most information.
In this work, we have proposed an MBDOE algorithm that extends previous work to include the design of an optimal input together with optimal measurements to reduce dynamical uncertainty in biological systems. The global nature of the algorithm overcomes the challenges posed by traditional MBDOE designs that rely on the local FIM design. Furthermore, most MBDOE techniques that consider input design only compare a predefined set of input values and/or measurements. Our computationally efficient approach enables us to avoid this limitation and search over the input space for the optimal input values at pre-specified time points. Our strategy achieves a global design by implementing computationally efficient strategies using sparse grid interpolation, probability-weighted scenario trees, and dynamical representative parameters. An interpolant of the ODE system is created by strategically sampling the uncertain parameter space and is used to evaluate additional points by LHS without simulating the ODEs. Representative parameters are used to span the dynamics of the acceptable parameter space. Together, the representative parameters and the sparse grid interpolation are used to efficiently screen the dynamical variance on the input space and facilitate the probability-weighted scenario tree to identify associated optimal measurements that minimize the uncertainty in the target system dynamics. We have confirmed that our optimal input MBDOE algorithm will specify experiments that are more informative in their ability to reduce dynamical uncertainty over MBDOE techniques that only specify optimal measurements using 3D and 19D models. The optimization of the input signal is most important when a limited number of measurements are taken. The computational expense of the MBDOE method is highly dependent on the dimension of the unknown parameter space and the feasible input space. For high dimension models, a large number of model simulations may be required to explore the large dimension uncertain parameter space to identify the acceptable and representative parameters. This computational burden can be reduced by the use of the sparse grid interpolation tool if the model output is smooth over the uncertain parameter space. For large input spaces, the majority of the model simulations in the MBDOE algorithm may be used to determine the optimal input signal. The number of model simulations increases with more flexibility in the input feasible space in terms of: (1) number of inputs that could be applied, (2) number of allowable changes in the inputs magnitudes (input resolutions), and (3) the number of admissible input times. In our examples, we specified a high degree of flexibility in the feasible input space for the lower dimension Hes1 model than the high dimension TCR model. As a result, the TCR model used 12,063 model simulations for the optimal input vector design while the Hes1 system required 16,531 model simulations. Historically, MBDOE methods have designed experiments sequentially, whereby information gained from a previous experiment is used immediately to inform the next experiment design [32,43,44]. This approach is prevalent with the local MBDOE strategies since it improves the accuracy of the parameter estimates that support the design of the next experiment. The sequential design becomes problematic when the number of uncertain parameters is large because a single experiment does not estimate all the parameters accurately [45]. A parallel design methodology improves the estimates of the parameters and reduces experimental costs by specifying multiple measurements. The parallel design in [15] uses scenario trees to specify all optimal measurements necessary to resolve the experimentally distinguishable dynamics. A hybrid design combines the advantage of the sequential design with the cost efficiency of the parallel design by specifying multiple experiments to perform and iterating with the results informing the next round of MBDOE. Our MBDOE design uniquely supports this hybrid design strategy. It uses a measurement scenario tree to determine a user-specified number of measurements as opposed to the previous scenario tree approach that indicated all measurements required to fully resolve the dynamics. Herein, the number of desired measurements controls the depth of the constructed scenario tree. (If a tree is made too deep, we believe that many of the specified measurements would be sub-optimal since the actual dynamics only lie within a fraction of the predicted tree.) Thus, if the desired dynamical resolution is not achieved by one application of our MBDOE algorithm, the design can be repeated while informed with the new additional experimental data. New experimental data restricts the acceptable parameter space further to generate a new more relevant optimal experimental design. Although our MBDOE algorithm is global, we cannot claim it is globally optimal. The algorithm makes optimal decisions within constraints. To make the problem tractable, we split the optimization problem into two sequential steps while the first specifies the optimal inputs and the second finds the associated measurement pairs. It is possible, that a better solution does exist and our algorithm misses it. Another source that may contribute to a sub-optimal result would be in the accuracy of the sparse grids. If the interpolants are not sufficiently accurate surrogates for the model responses, the results may be sub-optimal. In addition, we do not presume to have a good estimate of the number of acceptable parameter values needed to support the MBDOE since it is a function of the the biological system, the experimental setup, measurement noise, the mathematical model and its uncertainty. Hence, there is no exact formula for calculating the number of acceptable parameter vectors needed to support our experiment design algorithm. If the value is too small, the algorithm output would be sub-optimal and likely not repeatable. However, if repeated runs of our MBDOE algorithm produce similar results for increasing numbers of NA than it is likely to be sufficiently large. Thus, without putting constraints on the biological system model, NA must be sufficiently large to cover all distinguishable dynamics so that repeated application of our MBDOE algorithm derives similar optimal experiment designs. Overall, the proposed MBDOE strategy successfully extends previous MBDOE capabilities to design an optimal input with associated measurements that minimizes the uncertainty in the target system dynamics. The interpolation grid provides a computational efficient way to search for both parameter fits and optimal input over large uncertain spaces. The use of representative parameters, selected to span the dynamical space of the biological model, provides a means to sample the dynamical space uniformly enhancing the computational tractability of this approach. Furthermore, the probability weighted scenario tree designs modified from [15] supports input and multiple optimal measurement pair selection with fewer sampled parameters. Although the enhancement of the discrimination ability of the experiment is somewhat model dependent, for the examples we have presented herein, we have found that the ability to specify an optimal input in addition to optimal measurements has enhanced the ability to bound the expected target system dynamics.
Many mathematical models that have been developed for biological systems are limited because the complex systems are not well understood, the parameters are not known, and available data is limited and noisy. On the other hand, experiments to support model development are limited in terms of costs and time, feasible inputs and feasible measurements. MBDOE combines the mathematical models with experiment design to strategically design optimal experiments to obtain data that will contribute to the understanding of the systems. Our approach extends current capabilities of existing MBDOE techniques to make them more useful for scientists to resolve the trajectories of the system under study. It identifies the optimal conditions for stimuli and measurements that yield the most information about the system given the practical limitations. Exploration of the input space is not a trivial extension to MBDOE methods used for determining optimal measurements due to the nonlinear nature of many biological system models. The exploration of the system dynamics elicited by different inputs requires a computationally efficient and tractable approach. Our approach plans optimal experiments to reduce dynamical uncertainty in the output of selected target states of the biological system.
Efficient Optimization of Stimuli for Model-Based Design of Experiments to Resolve Dynamical Uncertainty