article: We present a nonlinear programming (NLP) framework for the scalable solution of parameter estimation problems that arise in dynamic modeling of biological systems. Such problems are computationally challenging because they often involve highly nonlinear and stiff differential equations as well as many experimental data sets and parameters. The proposed framework uses cutting-edge modeling and solution tools which are computationally efficient, robust, and easy-to-use. Specifically, our framework uses a time discretization approach that: i) avoids repetitive simulations of the dynamic model, ii) enables fully algebraic model implementations and computation of derivatives, and iii) enables the use of computationally efficient nonlinear interior point solvers that exploit sparse and structured linear algebra techniques. We demonstrate these capabilities by solving estimation problems for synthetic human gut microbiome community models. We show that an instance with 156 parameters, 144 differential equations, and 1,704 experimental data points can be solved in less than 3 minutes using our proposed framework (while an off-the-shelf simulation-based solution framework requires over 7 hours). We also create large instances to show that the proposed framework is scalable and can solve problems with up to 2,352 parameters, 2,304 differential equations, and 20,352 data points in less than 15 minutes. The proposed framework is flexible and easy-to-use, can be broadly applied to dynamic models of biological systems, and enables the implementation of sophisticated estimation techniques to quantify parameter uncertainty, to diagnose observability/uniqueness issues, to perform model selection, and to handle outliers.
  Dynamic modeling is essential for understanding the behavior of biological systems. Systems of interest in this domain include microbial communities and microbiome, gene regulatory networks, and metabolic pathways [1–3]. An important task that arises in modeling studies is validation against experimental data by using parameter estimation techniques. This task is computationally challenging because of the need to solve optimization problems constrained by differential equations. Challenges arise from the dimensionality, nonlinearity, and stiffness of the dynamic model, from the incomplete observation of the system states, from the need to estimate many parameters, and from the need to handle a large number of experimental data sets. Extensive research on solution methods for estimation problems with differential equations has been reported in the computational biology literature (see [4,5] for comprehensive reviews). These methods target maximum likelihood estimation formulations (which are derived from Bayesian principles). In these formulations, one aims to find parameters that maximize the likelihood function. The most used strategy to handle such formulations is the so-called simulation-based approach. Here, the idea is to perform repetitive simulations of the dynamic model at different trial parameter values to identify a set of parameters that maximizes the likelihood. The trial parameter values are updated using derivative-based or derivative-free search schemes [6–8]. While the simulation-based approach is intuitive, repetitive solutions of large dynamic models become computationally expensive and differential equation solvers can fail at trial parameter values that are non-physical or that trigger unstable responses. In addition, techniques to compute first and second order derivatives for derivative-based schemes (e. g. , finite differences, forward and adjoint sensitivities) involve intrusive procedures and are often limited to first-order derivatives [7]. The need for derivatives can be bypassed by using derivative-free search schemes [9,10], which are widely popular in computational biology. Such methods include simulated annealing [11,12], genetic algorithms [13,14], particle swarms [15], approximate Bayesian computation [16,17], and various other methods [18,19]. Derivative-free schemes do not scale well in the number of parameters (a larger number of trial parameter values often need to be explored compared to derivative-based schemes). Moreover, second order derivative information is needed to determine if the parameter estimates are unique/observable given available experimental data [20,21]. The uniqueness/observability test of the parameter estimates is based on curvature information of the likelihood function at the solution. Simulation-based estimation frameworks previously reported in the computational biology literature have focused on problems that usually contain less than 100 parameters [10,16,22]. To the best of our knowledge, the largest estimation problem solved using a simulation-based framework contains 3,780 data points and 1,801 parameters [7]. Such a problem was solved (to partial optimality) using a derivative-based search scheme that uses first-order derivative information (using an adjoint method) and required over 5 hours of computing time. The scalability limitations of simulation-based approaches present an important obstacle in considering models of higher fidelity, in exploiting high-throughput experimental data, in analyzing parameter uncertainties, and in implementing sophisticated techniques such as ensemble modeling. In this work, we propose a nonlinear programming (NLP) framework for solving estimation problems with embedded dynamic models [23,24]. The framework is based on a direct transcription approach wherein the dynamic model is converted into a large set of algebraic equations by applying time-discretization techniques. The algebraic equations are then embedded directly as constraints in the optimization problem (a nonlinear program-NLP). The NLPs arising from time discretization are of high dimension (easily reaching hundreds of thousands to millions of variables and constraints) but are also sparse and structured. Moreover, by transforming the dynamic model into algebraic equations, it becomes possible to use automatic differentiation techniques available in modern algebraic modeling languages to compute first and second derivatives. Exploitation of sparsity and structure, together with the availability of derivative information, enable the solution of estimation problems with complex dynamic models and efficient handling of many parameters and experimental data sets. Discretization-based estimation approaches have been widely studied in diverse fields such as chemical engineering [25–27] and aerospace engineering [28–30] (see [31] for a comprehensive review) but less so in computational biology. A major factor that has hindered wider adoption is the lack of easy-to-use computational frameworks that facilitate access to non-expert users. In this work, we demonstrate that modern modeling and solution tools can be combined to create scalable, robust, easy-to-use, and flexible frameworks. We demonstrate the benefits by solving challenging estimation problems arising in microbial community models. The proposed framework enables the implementation of higher level tasks such as observability analysis and uncertainty quantification. Uncertainty quantification (UQ) seeks to characterize parameter posterior distribution, which is necessary to obtain confidence levels/regions and parameter correlation information. Conventionally, UQ is performed by using second order derivative (Hessian) information of the likelihood function to construct an approximate parameter posterior covariance matrix [24,32] or by using a Markov-Chain-Monte-Carlo (MCMC) techniques [33–36]. The Hessian-based approach is scalable but it requires intrusive computation (cannot be automatically computed by the solver) and does not capture well the effect of nonlinearities and physical constraints [32]. In MCMC, one samples parameters from the prior parameter density and compares the associated model outputs with experimental data to decide whether to accept that sample or not. By repeating these accept/reject steps one can construct an approximate parameter posterior. MCMC is rather easy to implement (it is not intrusive and does not require solving an optimization problem) but, being simulation-based, also suffers from potential failures of the differential equation solver at non-physical parameter samples, it does not scale well with the number of parameters, and it cannot handle constraints directly (e. g. , nonnegative concentrations) [37]. In this work, we propose to overcome some of these challenges by using a randomized maximum a posteriori (rMAP) framework [36–39]. This method computes approximate (to second order) samples from the parameter posterior distribution by performing random perturbations on the experimental data and by re-solving the estimation problem. This allows exploration of the parameter space more efficiently compared to the MCMC scheme because each sample can be computed in parallel (MCMC is sequential). Moreover, the rMAP approach is non-intrusive, can capture nonlinear and physical constraints effects, and avoids potential failures of differential equation solvers. The proposed estimation framework is flexible and can easily accommodate advanced estimation formulations. To demonstrate this, we implement formulations that use different prior regularization schemes and k-max norms (the mean of a specified fraction of largest values) to mitigate large outliers [40,41]. The main contributions of this paper are summarized as follows. First, we provide an overview of parameter estimation and uncertainty quantification that leverages state-of-the-art optimization methods. Second, we demonstrate the computational capability of the nonlinear programming framework to handle large-scale parameter estimation problems for biological models. Lastly, we propose a novel problem formulation for parameter estimation using risk measures. Fig 1 shows an outline of parameter estimation framework presented in this paper. The paper is organized as follows. The methods section provides a general form for the estimation problem under study and discusses how this can be cast as a sparse NLP by using time-discretization techniques. Furthermore, we introduce basic concepts behind NLP solvers that exploit sparsity and structures at the linear algebra level. In addition, we discuss rMAP and outlier mitigation schemes. In the results section, we demonstrate that the proposed framework can handle challenging estimation problems arising in microbial community models.
  We consider estimation problems of the following form: min θ φ (θ) + ∑ k ∈ K φ k (η ¯ k, η k) (1) s. t. x ˙ k (t) = f k (x k (t), θ), k ∈ K, t ∈ [ 0, T k ] (2) x k (0) = x k 0, k ∈ K (3) η ¯ k (t) = ϕ k (x k (t), θ), k ∈ K, t ∈ T k (4) 0 ≤ h k (x k (t), θ), k ∈ K, t ∈ T k. (5) Here, K ≔ { 0,1 … K } is the set of experiments and T k ≔ { t 1, t 2, ⋯, t n k } is the set of measurement (sampling) times in experiment k ∈ K. Time is denoted as t ∈ [0, Tk], where T k ∈ R + is the duration for experiment k ∈ K. The variable vector x k: R → R n x k are the differential state time trajectories, θ ∈ R n θ are the model parameters, η ¯ k (t) ∈ R n η k are the model predicted outputs with corresponding experimental observations η k (t) ∈ R n η k at time t ∈ T k and experiment k ∈ K, and x k 0 ∈ R n k are initial conditions for experiment k ∈ K. For convenience in the notation, we define the output vectors η ¯ k = (η ¯ k (t 1), ⋯, η ¯ k (t n k) ) and the experimental output vectors η k = (η k (t 1), ⋯, η k (t n k) ) for experiment k ∈ K as well as the total output vector η ¯ = (η ¯ 1, ⋯, η ¯ K) and the total experimental output vector η = (η1, ⋯, ηK). The vector function fk (⋅) denotes the dynamic model mapping, φ (⋅) and φk (⋅) are objective function mappings, ϕk (⋅) is the state-to-output mapping, and hk (⋅) is the constraint mapping. All the mappings are assumed to be at least twice continuously differentiable with respect to all the arguments. The estimation formulation (1) – (5) captures all the features of our proposed framework. Our framework, however, can also accommodate more general features; for instance, the initial conditions (3) can be also considered as unknown variables that need to be estimated and we can define non-additive objective functions that penalize large errors. Problem (1) – (5) can be derived from Bayesian principles. To see this and introduce some useful notation, we start from Bayes theorem: p (θ | η) = p (η | θ) p (θ) p (η). (6) Here, p (θ | η) is the parameter posterior density (i. e. , the parameter density given knowledge on the outputs), p (η | θ) is the output posterior density (i. e. , the outputs given knowledge on the parameters), p (θ) is the prior density (i. e. , parameter density before knowledge of the output), and p (η) is the output marginal density. In a maximum a posteriori (MAP) formulation, the goal is to find the parameters that maximize p (θ | η). Because p (η) does not depend on θ, this can also be achieved by maximizing p (η|θ) p (θ). Maximizing p (η|θ) p (θ) is equivalent to maximizing the log-likelihood function L (θ) = log (p (θ) ) + log p (η|θ). If the outputs from the experiments k ∈ K are independent (which is usually the case), we have that p (η | θ) = ∏ k ∈ K p (η k | θ) and thus: L (θ) = log p (θ) + ∑ k ∈ K log p (η k | θ). (7) The observed outputs are random variables that are usually considered to be Gaussian and thus η k | θ ∼ N (η ¯ k, Σ k), where Σk is the covariance matrix. The prior density p (θ) is also often assumed to be Gaussian and thus θ ∼ N (θ ¯, Σ θ), where θ ¯ is the mean of the prior distribution and Σθ is its covariance. With this, we obtain: − log p (θ) = n θ 2 log 2 π + 1 2 log det Σ θ + 1 2 (θ − θ ¯) T Σ θ − 1 (θ − θ ¯) (8) − log p (η k | θ) = n η k 2 log 2 π + 1 2 log det Σ k + 1 2 (η k − η ¯ k) T Σ k − 1 (η k − η ¯ k). (9) By comparing (7) with (1) we can see that minimizing φ (θ) + ∑ k ∈ K φ k (η k, η ¯ k) is equivalent to minimizing −L (θ). Here, the dynamic model together with the state-to-output mapping defines a parameter-to-output mapping of the form η ¯ k ≔ m k (θ). In a simulation-based estimation approach, the mapping mk (θ) is computed by simulating the dynamic model (2) – (3) at a trial value θ using a differential equation solver and by evaluating the outputs at the sampling times t ∈ T k using (4). In a discretization-based approach, the mapping mk (θ) is not computed explicitly (but we use it here as a mathematical representation that is used to explain some relevant concepts). Constraints (5) restrict the parameter space to be explored. After dropping all constant terms in the likelihood function we obtain: φ (θ) = 1 2 (θ − θ ¯) T Σ θ − 1 (θ − θ ¯) (10) φ k (η ¯ k, η k) = 1 2 (η k − η ¯ k) T Σ k − 1 (η k − η ¯ k). (11) The function φ (θ) is usually known as the prior term and provides a regularization effect that stabilizes the solution of the estimation problem when the parameters cannot be uniquely infered from the available data [42–45]. This regularization term arises from the prior density p (θ) and provides a mechanism to encode knowledge on the parameters. Assuming that the prior density is Gaussian gives rise to a prior term that is defined by a weighted L2 norm. Recently, the machine learning community has also proposed the use of regularization terms that use L1 norms (e. g. , φ (θ) = ‖ θ − θ ¯ ‖ 1). The L1 norm induces sparsity in the parameters and corresponds to assuming that the prior density is Laplacian. One can also show that an L1 norm acts as an exact penalty function and implicitly induces constraints on the parameters. Similarly, one can also use the inequality constraints hk (⋅) to directly embed physical knowledge in the MAP formulation (e. g. , concentrations can only be positive). From (10) and (11) we see that φk (⋅) are squared error terms and thus the MAP problem minimizes the sum of the squared errors across all experiments k ∈ K. This approach offers limited control on large errors that might result from data outliers. Here, we propose to use a k-max norm to mitigate these issues. Our proposal is based on the observation that a k-max norm is equivalent to a conditional-value-at-risk (CVaR) norm [40,41,46]. The CVaRβ norm of a vector e = (e1, ⋯, eK) with components e k = φ k (η ¯ k, η k) is defined as the average of the β-fraction of largest elements of the vector (where β ∈ [0,1] is a parameter that defines the size of the fraction) [46]. One can show that, when β → 1, the CVaR norm is the largest fitting error and, when β → 0, the CVaR norm is the sum of fitting errors (as in the standard MAP formulation). A key computational property of the CVaR norm is that it can be formulated as a standard optimization problem. In particular, the MAP problem with a CVaR error norm can be expressed as: min θ, γ φ (θ) + K (γ + 1 (1 − β) ∑ k = 1 K [ e k − γ] +) (12) s. t. (2) − (5) (13) where [⋅]+ = max (0, ⋅) is the max function and γ is an auxiliary variable [41]. To solve the MAP problem (1) – (5) we approximate the differential equations by using a discretization scheme. This enables the use of computationally efficient NLP solvers and facilitates high-level UQ and observability monitoring tasks. The estimation problem under the MAP framework gives the values of the parameters θ* that maximize the parameter posterior density. However, a characterization of the entire posterior is necessary to assess parameter uncertainty. The posterior covariance may be approximated from the reduced Hessian at the solution of the problem w* and the covariance matrix can be used to determine ellipsoidal level sets of the posterior (confidence regions). This approach, however, might fail to capture nonlinear and constraint effects [24]. In this work, we circumvent these issues by using a randomized maximum a posteriori (rMAP) approach. Under this method, the posterior distribution is explored by using random perturbations on the experimental data (which can be easily parallelized). The rMAP framework can also deliver approximate samples from the parameter posterior distribution and implicitly captures nonlinear and constraint effects. To show this, we use the implicit mapping representation η ¯ k = m k (θ). Under this representation, the posterior density (6) can be expressed as: p (θ ∣ η) = 1 p (η) exp ( (θ − θ ¯) T Σ θ − 1 (θ − θ ¯) + ∑ k ∈ K (m k (θ) − η k) T Σ k − 1 (m k (θ) − η k) ) (38) = 1 p (η) exp (− 1 2 (m (θ) − η ^) T Σ − 1 (m (θ) − η ^) ) (39) where m (θ) ≔ (θ, m1 (θ), ⋯, mK (θ) ), Σ ≔ diag (Σθ, Σ1, Σ2, ⋯, ΣK), and η ^ = (θ ¯, η). Here, we redefine η ← η ^ to enable compact notation. Since θ* is a solution of the MAP problem, we have that: θ * = argmin θ (m (θ) − η) T Σ − 1 (m (θ) − η). (40) If the mapping mk (⋅) is continuously differentiable, we have that: m (θ) = m (θ *) + ∇ m (θ *) (θ − θ *) + O (‖ θ − θ * ‖ 2 2). (41) To enable compact notation we define η* = m (θ*) and ∇m* = ∇m (θ*). We have that θ* satisfies the stationary condition of (40): (∇ m *) T Σ − 1 (η * − η) = 0. (42) We use (41) to obtain a second-order Taylor approximation of the posterior as: p (θ ∣ η) ≈ 1 p (η) exp (− 1 2 (η * − η + ∇ m * (θ − θ *) ) T Σ − 1 (η * − η + ∇ m * (θ − θ *) ) ) = exp ( (η * − η) T Σ − 1 (η * − η) ) p (η) exp (− 1 2 (θ − θ *) T ∇ m * T Σ − 1 ∇ m * (θ − θ *) ) ) ∝ exp (− 1 2 (θ − θ *) T ∇ m * T Σ − 1 ∇ m * (θ − θ *) ) (43) This implies that the posterior is approximately represented as: θ ∣ η ∼ N (θ *, (∇ m * T Σ − 1 ∇ m *) − 1). (44) We recall that the output error is Gaussian and we can thus write η = m (θ) + ϵ with ϵ ∼ N (0, Σ). We now consider the MAP problem with randomly perturbed data: θ ˜ = argmin θ (m (θ) − (η + ϵ) ) T Σ − 1 (m (θ) − (η + ϵ) ) (45) and note that (m (θ) − (η + ϵ) ) T Σ − 1 (m (θ) − (η + ϵ) ) = (η * − η + ∇ m * (θ − θ *) − ϵ) T Σ − 1 (η * − η + ∇ m * (θ − θ *) − ϵ) + O (‖ θ − θ * ‖ 2 2) = (θ − θ *) T ∇ m * T Σ − 1 ∇ m * (θ − θ *) + 2 (θ − θ *) T ∇ m * T Σ − 1 ϵ + O (‖ θ − θ * ‖ 2 2) + C, (46) where C is a constant. Consequently, for sufficiently small ϵ, we can linearize the mapping m (⋅) to obtain an approximate solution of (45) of the form: θ ˜ ≈ θ * + (∇ m * T Σ − 1 ∇ m *) − 1 ∇ m * T Σ − 1 ϵ. (47) Here, we observe that the right-hand side of (47) is Gaussian with mean θ* and covariance: ( (∇ m * T Σ − 1 ∇ m *) − 1 ∇ m * T Σ − 1) Σ ( (∇ m * T Σ − 1 ∇ m *) − 1 ∇ m * T Σ − 1) T = (∇ m * T Σ − 1 ∇ m *) − 1 (∇ m * T Σ − 1 ∇ m *) (∇ m * T Σ − 1 ∇ m *) − 1 = (∇ m * T Σ − 1 ∇ m *) − 1. (48) We thus have that: θ ˜ ∼ N (θ *, (∇ m * T Σ − 1 ∇ m *) − 1), (49) Consequently, solving (45) provides an approximate sample from the posterior distribution p (θ ∣ η). The sampling procedure (45) is accurate up to second order. To obtain an exact sampling from the posterior, one needs to implement a rigorous MCMC scheme. The MCMC scheme removes the bias that appears in the rMAP sample density, which results from the second order approximation [35,36,36,37]. This effect is illustrated with an examples in [37]. Several works in the literature, however, report that accurate posterior densities can be obtained using an rMAP scheme [38,57,58]. We also note that the rMAP scheme implicitly captures nonlinearities and physical constraints when computing samples from the posterior (MCMC does not capture constraints). In particular, solving the perturbed problem (45) corresponds to solving the MAP problem (1) – (5) with randomly perturbed data and the MAP problem enforces constraints and handles the full nonlinear model. The computational framework provided in this work focuses on the solution of the estimation problem and seeks to highlight that these capabilities enable the implementation of advanced UQ procedures. A detailed analysis and comparison of different UQ methods for biological dynamical models is an interesting topic of future work. Systematic comparisons in other settings are provided in [35,36]. Having an algebraic representation of the estimation problem has many practical and computational advantages. In particular, one can implement the estimation problem in easy-to-use and open-source modeling languages such as JuMP [59], Plasmo. jl [60], and Pyomo [61,62]. These modeling languages are equipped with automatic differentiation techniques that compute exact first and second derivatives. Derivative information is communicated to optimization solvers without any user intervention. Modern algebraic modeling languages such as Plasmo. jl and Pyomo also allow users to convey structural information to the solvers. This is beneficial in the case of parameter estimation, where the structure can be exploited to enable parallelism and the use of high-performance computing clusters. In our framework, we use the modeling language Plasmo. jl to express multi-experiment estimation problems as graphs. Our implementation using Plasmo. jl is illustrated in Fig 2. The full Julia script is available at https: //github. com/zavalab/JuliaBox/tree/master/MicrobialPLOS. We highlight that the same script can be used to solve the estimation problem using a general NLP solver such as Ipopt on a single-processor computer or with a structure-exploiting parallel NLP solver such as PIPS-NLP on multiple parallel computing processors (this might be a multi-core computing server or a large-scale computing cluster). This allows users with limited knowledge on scientific computing to gain access to advanced high-performance computing capabilities.
  Parameter observability was checked by solving the MAP formulation for P1 (which uses the available experimental data) and by checking the inertia of the augmented system at the solution (reported by Ipopt). Here, we omitted the prior regularization term φ (⋅). We found that parameters obtained from P1 are not unique (not observable from the available data). Moreover, we found that the estimated parameter values without regularization have unrealistic (non-physical) values, see Fig 3 (a). This observation justifies the need to use prior information. We have used L1 and L2 priors with uniform standard deviations (other forms of prior can also be used as long as they can be reformulated as smooth functions). The results are presented in Fig 3 (b) and 3 (c). Unique parameter estimates were found when L1 or L2 priors were used. We also found that the L1 prior induces sparser solutions (many parameters are zero). For the remainder of the results, we use the formulation with an L2 prior. Model validation was performed by assessing the goodness of fit to the experimental data (Fig 4). We can see that the model is capable of fitting most of the data points, but there are a number of experiments where the model prediction deviates significantly from the experimental data (such experiments are highlighted with red boxes). Furthermore, we can observe outliers at single data points (highlighted with red circles). Poor fitting can be caused by either bad local minima (the optimization solver finds a local optimal solution rather than the global optimal solution) or by a structural model error (the model structure is incapable of capturing the actual behavior of the system). To avoid bad local minima, we solved the MAP problem with multiple starting points. Such an approach increases the probability to find the global optimum, but obtaining a rigorous certificate of a global minimum is computationally challenging (rigorous global optimization techniques are currently not scalable to large problems). We found that the use of multiple starting points does not improve the model fit. Consequently, we attribute fitting errors to the model structure itself. In particular, the gLV model neglects various physical and biological phenomena such as lag phase or interaction coefficients that change as a function of time [66]. To investigate structural errors, we solved the MAP problem with a variant of the gLV model. In particular, we investigated the saturable gLV model [72,73] (Model 2): d x s d t = (μ s + ∑ s ′ ∈ S α s s ′ x s ′ K s s ′ + x s ′) x s, (Model 2) where Kss′ > 0 are additional interaction parameters. The saturable model exhibits a much higher degree of nonlinearity than the gLV model and includes S2 more parameters (the number of degrees of freedom increases from S2 + S to 2S2 + S). As a result, the saturable gLV model provides more flexibility to improve model fitting. These results thus seek to illustrate that the proposed NLP framework can be used to handle challenging dynamic models. The model fitting obtained with the saturable gLV form is illustrated in Fig 5. As can be seen, significant improvements were made; in particular, the overall fitting error ∑ k ∈ K ∑ s ∈ S k ∑ t ∈ T k 1 2 (η k, s (t) − η ¯ k, s (t) σ k, s (t) ) 2, (51) was reduced by 30%. Increasing the number of degrees of freedom can cause overfitting, however, and this can make the model less predictive. Consequently, there is a trade-off between fit and predictability. To mitigate outliers, we solved the MAP problem with a CVaR norm and β = 0. 9 (to penalize the 10% largest errors). The relevant results are summarized in Fig 6. It can be observed that the fitting errors for the outliers obtained with the standard MAP formulation are reduced. The effect of the CVaR formulation is also evident when analyzing the prediction error histogram (see Fig 7). In particular, we observe that the tail of high prediction errors becomes less pronounced under the CVaR formulation. In particular, the mean of the 10% largest errors decreases by 18% (from 167. 81 to 137. 04). On the other hand, it can also be observed that the mean error increases under CVaR and that the tail of small errors shrinks. This illustrates the fundamental trade-off that usually arises in robust statistics. The behavior induced with CVaR aids estimator performance because it prevents overfitting experimental data sets. We used rMAP to assess the uncertainty of the 156 parameters estimated from P1 using the available experimental data. To do so, we draw data perturbations as ηk, s (t) ← ηk, s (t) + ϵk, s (t) with ϵ k, s (t) ∼ N (0, σ k, s (t) 2). We solved 500 MAP problems to obtain parameter samples and use this to approximate the covariance matrix for the posterior. The standard deviations are shown in Fig 8 (a) and the marginals for the posterior are shown in Fig 9. A large standard deviation indicates that the estimated parameter value is not reliable. We note that about half of the parameters can be estimated reliably while the other half exhibit significant uncertainty. This indicates that more experimental data should be obtained. From the sample covariance, we generated 95% ellipsoidal confidence regions for each pair of parameters. The correlation plots of μs against αss′ for s, s ′ ∈ S are shown in Fig 10 and the Pearson correlation coefficients are shown in Fig 8 (b). In an ideal case, the parameters should be uncorrelated because data should be sufficient to estimate each parameter reliably. Using our data set, however, we can observe strong correlations between the parameters μs and αss in Fig 10, and strong positive and negative correlations can also be found in Fig 8 (b). Furthermore, since the whole approximate distribution is obtained in the inference analysis based on rMAP framework, we can perform more sophisticated analysis on the characteristics of the distribution. In particular, one can investigate third and fourth moments (Fig 11) to examine the skewness and the kurtosis of the distribution. Such information can be used to investigate the deviation of the posterior distribution from the normal distribution. If the posteriors are stritly normally distributed, the third and fourth moments should be zero and 3σ4, respectively. However, we can observe that many posterior distributions deviate from such expectations. Thus, we can see that some of the distributions are not close to the normal distribution. We assessed the computational scalability of the estimation framework by analyzing problems with different sizes and characteristics. Problem P1 was implemented in the algebraic modeling platform JuMP and solved with the NLP solver Ipopt configured with the sparse linear solver MA57. Problem P1 with gLV model and L2 prior was solved in 134 seconds and 78 NLP iterations on a standard computing server with an Intel (R) Xeon (R) CPU E5-2698 v3 processor running at 2. 30GHz. Problem P1 with gLV model and L1 prior was solved in 219 seconds and 68 NLP iterations with the same hardware. A comparable problem requires over 7 hours to solve using a simulation-based approach implemented in Matlab and that uses finite differences to obtain first derivatives [66]. Despite the significant gains in computational performance obtained with Ipopt, its solution time scales nearly quadratically with the number of data sets. To overcome this scalability issue, we compared the performance of the serial solver Ipopt against that of the parallel solver PIPS-NLP (which uses a Schur complement decomposition to perform linear algebra operations). To test the scalability of PIPS-NLP, we generated a larger version of the estimation problem (labeled as P2). This problem is created by adding synthetic data sets. The NLP corresponding to P2 has over one million variables and constraints (but the number of parameters is the same as that of P1). This problem was implemented in Plasmo. jl. The benefit of using a parallel approach is clearly seen in Fig 12. Here, we highlight that PIPS-NLP solved P2 in less than 10 minutes and 94 NLP iterations using 16 cores while Ipopt requires around 30 minutes and 67 iterations. Furthermore, IPOPT found a different local solution and the solution from PIPS-NLP had a better objective value. Fig 12 (b) also shows that PIPS-NLP achieves nearly perfect strong scaling (speedup increases linearly with the number of cores). In rMAP-based uncertainty quantification, the main computational challenge was the repetitive solution of the optimization problems. However, such challenge can be overcome by using the existing solution information. The required number of iterations in NLP solver can be greatly reduced when a good starting point (initial guess of the solution) is available (often referred to as warm start). Since only small modifications are made to the original problem to formulate the rMAP problem, the NLP solution of rMAP problem is very similar to that of the original problem. Thus, by warm-starting the NLP with the original NLP solution, the computational efforts to solve rMAP problem can be significantly reduced. In particular, most rMAP sampling problem was solved in less than 10 NLP iterations while the original problem required 78 NLP iterations. We also assessed computational capability in estimation problems with the larger number of species in the microbial community (which increases the number of differential equations and parameters). Here, we generated synthetic data using simulations for larger communities. The generated data are summarized in Table 2. The number of the parameters and of data points scales nearly quadratically with respect to the size of the community. The computation times are shown in Fig 13. The results indicate that, by using PIPS-NLP, one can solve estimation problems with up to 48 species in less than 15 minutes and 40 NLP iterations (using 12 parallel computing cores). We highlight that, to the best of our knowledge, problem S4 is the largest estimation problem reported in computational biology literature. This problem contains 2,304 differential equations, 2,352 parameters, and 20,352 data points. The corresponding NLP contains 1. 3 million variables and constraints.
  The high computational efficiency achieved with the proposed framework can enable kinetic modeling of complex biological systems ranging from biomolecular networks to high-dimensional microbial communities [74]. Indeed, the proposed framework can be used to construct and analyze high-fidelity models of whole-cells or microbiomes [75,76]. In particular, these methods can be applied to develop predictive dynamic models of multi-gene synthetic circuits interacting with host-cell processes for accurately predicting cell growth and synthetic circuit activity [77] or kinetic models of metabolite transformations driving community dynamics. The proposed framework can also be used to handle more sophisticated models that arise in biological systems such as delay differential equations, differential and algebraic equations, and partial differential equations. Dynamic biological models with embedded metabolic flux formulations (giving rise to non-smooth behavior) can be handled with the proposed framework by using reformulations [78]. Estimation problems for stochastic differential equations (such as stochastic chemical kinetics models) cannot be handled with the proposed framework and remain a challenging class of problems [79]. The methods provided in this work will advance our capability of integrating mechanistic modeling frameworks with large-scale experimental data. Furthermore, uncertainty quantification and observability analysis can provide valuable information to guide and accelerate experimental data collection. These capabilities are also essential in diagnosing structural model errors. The proposed framework uses state-of-the-art and easy-to-use modeling and solution tools that can be broadly applied to diverse biological systems and accessible to a wide range of users. In the future, the proposed framework can be interfaced with systems biology markup languages such as SBML [80] and CellML [81] to allow broader applicability. Together, these advances can ultimately transform biology into a predictive and model-guided discipline.
summary: Constructing and validating dynamic models of biological systems spanning biomolecular networks to ecological systems is a challenging problem. Here we present a scalable computational framework to rapidly infer parameters in complex dynamic models of biological systems from large-scale experimental data. The framework was applied to infer parameters of a synthetic microbial community model from large-scale time series data. We also demonstrate that this framework can be used to analyze parameter uncertainty, to diagnose whether the experimental data are sufficient to uniquely determine the parameters, to determine the model that best describes the data, and to infer parameters in the face of data outliers.
keywords: classical mechanics
  random variables
  covariance
  systems science
  mathematics
  algebra
  computer and information sciences
  algebraic structures
  nonlinear dynamics
  differential equations
  probability theory
  physics
  inertia
  linear algebra
  biology and life sciences
  physical sciences
  species interactions
  eigenvalues
  motion
title: Scalable nonlinear programming framework for parameter estimation in dynamic biological system models