JuPOETs: a constrained multiobjective optimization approach to estimate biochemical model ensembles in the Julia programming language

Background Ensemble modeling is a promising approach for obtaining robust predictions and coarse grained population behavior in deterministic mathematical models. Ensemble approaches address model uncertainty by using parameter or model families instead of single best-fit parameters or fixed model structures. Parameter ensembles can be selected based upon simulation error, along with other criteria such as diversity or steady-state performance. Simulations using parameter ensembles can estimate confidence intervals on model variables, and robustly constrain model predictions, despite having many poorly constrained parameters. Results In this software note, we present a multiobjective based technique to estimate parameter or models ensembles, the Pareto Optimal Ensemble Technique in the Julia programming language (JuPOETs). JuPOETs integrates simulated annealing with Pareto optimality to estimate ensembles on or near the optimal tradeoff surface between competing training objectives. We demonstrate JuPOETs on a suite of multiobjective problems, including test functions with parameter bounds and system constraints as well as for the identification of a proof-of-concept biochemical model with four conflicting training objectives. JuPOETs identified optimal or near optimal solutions approximately six-fold faster than a corresponding implementation in Octave for the suite of test functions. For the proof-of-concept biochemical model, JuPOETs produced an ensemble of parameters that gave both the mean of the training data for conflicting data sets, while simultaneously estimating parameter sets that performed well on each of the individual objective functions. Conclusions JuPOETs is a promising approach for the estimation of parameter and model ensembles using multiobjective optimization. JuPOETs can be adapted to solve many problem types, including mixed binary and continuous variable types, bilevel optimization problems and constrained problems without altering the base algorithm. JuPOETs is open source, available under an MIT license, and can be installed using the Julia package manager from the JuPOETs GitHub repository


Background
Ensemble modeling is a promising approach for obtaining robust predictions and coarse grained population behavior in deterministic mathematical models. It is often not possible to uniquely identify all the parameters in biochemical models, even when given extensive training data [1]. Thus, despite significant advances in standardizing biochemical model identification [2], the problem of estimating model parameters from experimental data remains challenging. Ensemble approaches address parameter uncertainty in systems biology and other fields like weather prediction [3][4][5][6] by using parameter families instead of single best-fit parameter sets. Parameter families can be selected based upon simulation error, along with other criteria such as diversity or steadystate performance. Simulations using parameter ensembles can estimate confidence intervals on model variables, and robustly constrain model predictions, despite having many poorly constrained parameters [7,8]. There are many techniques to generate parameter ensembles. Battogtokh et al., Brown et al., and later Tasseff et al. generated experimentally constrained parameter ensembles using a Metropolis-type random walk [3,5,9,10]. Liao and coworkers developed methods to generate ensembles that all approach the same steady-state, for example one determined by fluxomics measurements [11]. They have used this approach for model reduction [12], strain engineering [13,14] and to study the robustness of nonnative pathways and network failure [15]. Maranas and coworkers have also applied this method to develop a comprehensive kinetic model of bacterial central carbon metabolism, including mutant data [16]. We and others have used ensemble approaches, generated using both sampling and optimization techniques, that have robustly simulated a wide variety of signal transduction processes [9,10,[17][18][19], neutrophil trafficking in sepsis [20], patient specific coagulation behavior [21], uncertainty quantification in metabolic kinetic models [22] and to capture cell to cell variation [23]. Further, ensemble approaches have been used in synthetic biology to sample possible biocircuit configurations [24]. Thus, ensemble approaches are widely used to robustly simulate a variety of biochemical systems.
Identification of biochemical models requires significant training data perhaps taken from diverse sources. These real-world data sets often contain intrinsic conflicts resulting from, for example, the use of different cell lines, different measurement technologies, different reagent vendors or lots, uncontrollable experimental artifacts or general cross laboratory variability. Parameter ensembles that optimally balance these inherent conflicts lead to more robust model performance. Multiobjective optimization is an ensemble generation technique that naturally balances conflicts in noisy training data [25].
Multiobjective optimization has been used to identify signal transduction models [18,23], for the design of synthetic circuits [24], to design the folding behaviors of novel RNAs [26], to design bioprocesses [27], and to understand bacterial adaptation [28]. Thus, it is a widely used approach for a variety of biochemical applications. Previously, we developed the Pareto Optimal Ensemble Technique (POETs) algorithm to address the challenge of competing or conflicting training objectives. POETs, which integrates simulated annealing (SA) and multiobjective optimization through the notion of Pareto rank, estimates parameter ensembles which optimally trade-off between competing (and potentially conflicting) experimental objectives [29]. However, the previous implementation of POETs, in the Octave programming language [30], suffered from poor performance and was not configurable. For example, Octave-POETs does not accommodate user definable objective functions, bounds and problem constraints, cooling schedules, different variable types e.g., a mixture of binary and continuous design variables or custom diversity generation routines. Octave-POETs was also not well integrated into a package or source code management (SCM) system. Thus, upgrades to the approach containing new features, or bug fixes were not centrally managed.

Implementation
In this software note, we present an open-source implementation of the Pareto optimal ensemble technique in the Julia programming language (JuPOETs). JuPO-ETs takes advantage of the unique features of Julia to address many of the shortcomings of the previous implementation. Julia is a cross-platform, high-performance programming language for technical computing that has performance comparable to C but with syntax similar to MATLAB/Octave and Python [31]. Julia also offers a sophisticated compiler, distributed parallel execution, numerical accuracy, and an extensive function library. Further, the architecture of JuPOETs takes advantage of the first-class function type in Julia allowing user definable behavior for all key aspects of the algorithm, including objective functions, custom diversity generation logic, linear/non-linear parameter constraints (and parameter bounds constraints) as well as custom cooling schedules. Julia's ability to naturally call other languages such as Python or C also allows JuPO-ETs to be used with models implemented in a variety of languages across many platforms. Additionally, Julia offers a built-in package manager which is directly integrated with GitHub, a popular web-based Git repository hosting service offering distributed revision control and source code management. Thus, JuPOETs can be adapted to many problem types, including mixed binary and continuous variable types, bilevel problems and constrained problems without altering the base algorithm, as was required in the previous POETs implementation.

JuPOETs optimization problem formulation
JuPOETs solves the K−dimensional constrained multiobjective optimization problem: The quantity O j denotes the j th objective function (j = 1, 2, . . . , K), typically the sum of squared errors for the j th data set for biochemical modeling applications. The terms f (t, x(t, p),ẋ(t, p), u(t), p) denote the system of model equations (e.g., differential equations, differential algebraic equations or linear/non-linear algebraic equations) where p denotes the decision variable vector e.g., unknown model parameters (D × 1). In typical biochemical modeling applications, the model equations f(·) are a system of continuous real-valued non-linear differential equations that comprise a kinetic model, but other types of models e.g., stoichiometric models are also common. The quantity t denotes time, x (t, p) denotes the model state (with an initial state x 0 ), and u(t) denotes an input vector. The decision variables (e.g., kinetic parameters) can be subject to bounds constraints, where L and U denote the lower and upper bounds, respectively as well as C problem specific constraints g i (t, x(t, p), u(t), p) , i = 1, . . . , C. The decision variables p are typically real-valued kinetic constants, or metabolic fluxes in the case of stoichiometric models. However, other variables types e.g., binary or categorical decision variables can also be accommodated.
JuPOETs integrates simulated annealing (SA) [32] with Pareto ranking to estimate decision variables on or near the optimal tradeoff surface between competing objectives ( Fig. 1 and Algorithm 1). A tradeoff surface defines the best possible performance for every conflicting objective, such that an increase in the performance of one objective does not decrease the performance of at least one other objective. Pareto rank is a scalar measure of distance away from the optimal tradeoff surface (low rank is near the surface, while higher ranks are progressively further away). Thus, the central idea underlying POETs is a mapping between the value of the objective vector evaluated at p i+1 (decision variable guess at iteration i + 1) and the scalar Pareto rank (Fig. 1). Traditional simulated annealing uses a scalar performance value e.g., simulation error to make a probabilistic decision to keep or reject a set of decision variables; decision variables with better performance are always accepted, while those with worse performance are sometimes accepted depending upon a parameter called the temperature. On the other hand, JuPOETs makes this same decision using the Pareto rank instead of a single performance objective. The problem of estimating biochemical model parameters from experimental data is typically posed as an error minimization problem over continuous real-valued decision variables (model parameters) subject to the model equations. A parameter set p i+1 lies along the optimal tradeoff surface if no other parameter guess leads to decreased error for Fig. 1 Schematic of multiobjective parameter mapping. The performance of any given parameter set is mapped into an objective space using a ranking function which quantifies the quality of the parameters. The distance away from the optimal tradeoff surface is quantified using the Pareto ranking scheme of Fonseca and Fleming in JuPOETs Algorithm 1: Pseudo-code for the JuPOETs run-loop. The user must specify the objective function and an initial parameter guess. The user can optionally specify the neighbor, acceptance, cooling and refinement functions (or use the default implementations). The rank archive R, solution archive S and objective archive O are initialized from the initial guess. The initial guess (potentially following a single objective local refinement step) is perturbed in the neighbor function, which generates a new solution whose performance is evaluated using the user supplied objective function. The new solution and objective values are then added to the respective archives and ranked using the builtin rank function. If the new solution is accepted (based upon a probability calculated with the user supplied acceptance function) it is added to the solution and objective archive. This solution is then perturbed during the next iteration of the algorithm. However, if the solution is not accepted, it is removed from the archive and discarded. The temperature is adjusted using the user supplied cooling function after each I iterations. When JuPOETs terminates, the parameter solution archive S, objective archive O and rank archive R are retuned to the caller.
input : User specified objective function, and initial guess (D × 1). User can also specify custom neighbor, acceptance.
cooling and refinement functions or use the default functions provided. every objective. JuPOETs calculates the performance of a candidate parameter set p i+1 by calling the user defined objective function; objective takes a parameter set as an input, evaluates the model equations, and using this solution, returns the K × 1 objective vector. Candidate parameter sets are generated by the user supplied neighbor function; the default implementation of neighbor is a random perturbation, however other perturbation logic can be implemented by the user. The error vector associated with p i+1 is ranked using the builtin Pareto rank function, by comparing the error at iteration i + 1 to the error archive O i (all error vectors up to iteration i meeting a ranking criterion). Parameter sets on or near the optimal trade-off surface between the objectives have a rank equal to 0 (no other current parameter sets are better). These rank zero parameter sets define the Pareto optimal group for the ensemble, wherein Pareto optimality is defined as a parameter set not being dominated by any other sets within the ensemble. Sets with increasing non-zero rank are progressively further away from the optimal trade-off surface. Thus, a parameter set with a rank = 0 is better in a trade-off sense than rank > 0. We implemented the Fonseca and Fleming ranking scheme in the builtin rank function [33]: where rank r is the number of parameter sets that dominate (are better than) parameter set p i+1 , and O i+1 (p i+1 ) denotes the objective vector evaluated at p i+1 . We used the Pareto rank to inform the SA calculation. The parameter set p i+1 was accepted or rejected by the SA at each iteration, by calculating an acceptance probability P (p i+1 ): where T is the simulated annealing temperature; the temperature provides control over how strictly decreasing Pareto rank is enforced.
the acceptance probability moves toward one, ensuring that we explore parameter sets along the Pareto surface.
Occasionally, (depending upon T) a parameter set with a high Pareto rank is accepted by the SA allowing a more diverse search of the parameter space. However, as T is reduced as a function of iteration count (using the cooling function), the probability of accepting a highrank set decreases. Parameter sets could also be accepted by the SA but not permanently archived in S i , where S i is the solution archive. Only parameter sets with rank less than or equal to a threshold (rank ≤4 by default) are included in S i , where the archive is re-ranked and filtered after accepting every new parameter set. Parameter bounds were implemented in the neighbor function as box constraints, while problem specific constraints were implemented in objective using a penalty method: where λ denotes the penalty parameter (λ = 100 by default). However, because both the neighbor and objective functions are user defined, different constraint implementations are easily defined.
To use JuPOETs, the user specifies the neighbor, acceptance, cooling and objective functions along with an initial decision variable guess. Default implementations of the neighbor, acceptance and cooling functions can be used directly, or they can be overridden by user defined logic. However, the user must provide an implementation of the objective function and provide an initial decision variable guess. Lastly, if the user is operating JuPOETs in hybrid mode, then a refinement function pointer must also be specified. Hybrid mode temporarily switches the search from a multiobjective to a single objective problem, where the sum of the objective functions can be used to update the best (or initial) parameter guess. The specific hybrid mode search logic is up to the user; by default hybrid mode is off, and the default refinement implementation is simply a pass through function. However, we have shown previously that POETs operated in hybrid mode (where the single objective problem used a pattern search approach) had better performance that POETs alone [29]. Thus, hybrid mode is generally recommended for most applications. In addition, there are several user configurable parameters that can be adjusted to control the performance of JuPO-ETs: maximum_number_of_iterations controls the number of iterations per temperature (default 20); rank_cutoff controls the upper rank bound on the solution archive (default 5); temperature_min controls the minimum temperature after which JuPOETs returns the error and solution archives (default 0.001); show_trace controls the level of output shown to the user (default true). After the completion of the run, JuPO-ETs returns the parameter solution archive S, objective archive O and rank archive R. The parameter solution archive S contains is an D × A array, where A denotes the number of solutions in the archive when JuPOETs terminated. On the other hand, the objective archive O is an K × A array containing the performance values for each objective corresponding the columns of S. Lastly, JuPOETs returns the rank archive R which is an A × 1 array of Pareto ranks corresponding to the columns of S. One technical note, if JuPOETs is run from multiple starting locations, and the archives from each of these runs is combined into a single collective archive, the combined parameter rank archive may become invalid. In these cases, it is required to re-rank the parameter sets using the built-in rank function to produce a collective parameter ranking.

Results and discussion
JuPOETs identified optimal or nearly optimal solutions significantly faster than Octave-POETs for a suite of multiobjective algebraic test problems ( Table 1). The algebraic test problems were constrained non-linear functions with bound constraints and additional non-linear constraints on the decision variables in one case. The problems had up to three-dimensional continuous real-valued decision vectors, and each case had two objective functions. The wall-clock time for JuPOETs and Octave-POETs was measured for 10 independent trials for each of the test problems. The same cooling, neighbor, acceptance, and objective logic was employed between the implementations, and all other parameters were held constant. For each test function, the search domain was partitioned into 10 segments, where an initial parameter guess was drawn from each partition. The number of search steps for each temperate was I = 10 for all cases, and the cooling parameter was α = 0.9. On average, JuPOETs identified optimal or near optimal solutions for the suite of test problems six-fold faster (60s versus 400s) than Octave-POETs (Fig. 2). JuPOETs produced the characteristic tradeoff curves for each test problem, given both decision variable bound and problem constraints (Fig. 3). Thus, JuPOETs estimated an ensemble of solutions to constrained multiobjective algebraic test problems significantly faster than the current Octave implementation. Next, we tested JuPOETs on a proof-of-concept biochemical model identification problem.
JuPOETs estimated an ensemble of biochemical model parameters that were consistent with the mean of synthetic training data (Fig. 4). Four synthetic training data sets were generated from a prototypical biochemical network consisting of 6 metabolites and 7 reactions (Fig. 4, inset right). We considered a common case in which the same extracellular measurements of A e , B e , C e and cellmass were made on four hypothetical cell types, each having the same biological connectivity but different performance. Network dynamics were modeled using the hybrid cybernetic model with elementary modes (HCM) approach of Ramkrishna and coworkers [34]. In the HCM approach, metabolic networks are first decomposed into a set of elementary modes (EMs) (chemically balanced steady-state pathways, see [35]). Dynamic combinations of elementary modes are then used to characterize network behavior. Each elementary mode is catalyzed by a pseudo enzyme; thus, each mode has both kinetic and enzyme synthesis parameters. The proof of concept network generated 6 EMs, resulting in 13 model parameters (continuos real-valued decision variables). The synthetic training data was generated by randomly varying these parameters.
The general form of the biochemical test problem was given by: subject to model and bounds constraints. We considered four training data sets (K = 4), each of which contained time-series measurements of A e , B e , C e and cellmass. Each objective O j , j = 1, . . . , K quantified the squared difference between the simulated (x i ) and measured extracellular species abundance (y i ) in the j th data set: where, i denotes the species index and τ denotes the time index. The abundance of extracellular species i (x i ), the pseudo enzyme e l (catalyzes flux through mode l), and cellmass were governed by the model equations: Binh and Korn function 2 Fonseca and Fleming function 3 The performance of JuPOETs on the multi-objective test suite. The execution time (wall-clock) for JuPOETs and POETs implemented in Octave was measured for 10 independent trials for the suite of test problems. The number of steps per temperature I = 10, and the cooling parameter α = 0.9 for all cases. The problem domain was partitioned into 10 equal segments, an initial guess was drawn from each segment. For each of the test functions, JuPOETs estimated solutions on (rank zero solutions, black) or near (gray) the optimal tradeoff surface, subject to bounds and problem constraints The term q l (e, p, x) denotes the specific uptake/secretion rate for mode l where e denotes the pseudo enzyme vector, p denotes the unknown kinetic parameter vector (decision variables), x denotes the extracellular species vector, and c denotes the cell mass; q l (e, p, x) is the product of a kinetic term (q l ) and a control variable governing enzyme activity. Flux through each mode was catalyzed by a pseudo enzyme e l , synthesized at the regulated specific rate r E,l (p, x), and constitutively at the rate α l . The term u l denotes the cybernetic variable controlling the synthesis of enzyme l. The term β l denotes the rate constant governing non-specific enzyme degradation, and r G denotes the specific growth rate through all modes. The specific uptake/secretion rates and the specific rate of enzyme synthesis were modeled using saturation kinetics. The specific growth rate was given by: where z μl denotes the growth flux μ through mode l. The control variables u l and v l , which control the synthesis and activity of each enzyme respectively, were given by: and v l = z slql max l=1,...,L z slql (8) where z sl denotes the uptake flux of substrate s through mode l. Each unknown kinetic parameter was continuous and real-valued, and subject to bounds constraints: L ≤ p ≤ U . JuPOETs produced an ensemble of approximately dim S 13,000 parameter sets that captured the mean of the measured data sets for extracellular metabolites and cellmass ( Fig. 4a and b). JuPOETs minimized the difference between the simulated and measured values for extracellular metabolites A e , B e , C e and cellmass, where the residual for each data set was treated as a single objective (leading to four objectives). The 95% confidence estimate produced by the ensemble was consistent with the mean of the measured data, despite having significant uncertainty in the training data. JuPOETs produced a consensus estimate of the synthetic data by calculating optimal trade-offs between the training data sets (Fig. 4c).
Multiple trade-off fronts were visible in the objective plots, for example between data set 3 (O 3 ) and data set 2 (O 2 ). Thus, without a multiobjective approach, it would be challenging to capture these data sets as fitting one leads to decreased performance on the other. However, the ensemble contained parameter sets that described each data set independently (Fig. 5). Thus, JuPOETs produced an ensemble of parameters that gave the mean of the training data for conflicting data sets, while simultaneously estimating parameter sets that performed well on each individual objective function. Currently, JuPOETs does not consider parameter identifiability when constructing parameter ensembles. Although JuPOETs produces parameter estimates that give model performance similar to the training data, we do not have strict statistical confidence that the true parameter values are contained within the ensemble. However, despite this, ensembles produced by POETs can be predictive [18,23]. Thus, JuPOETs produces a collection of Cellmass measurements (points) versus time for experiment 2 and 3 were compared with ensemble simulations. The full ensemble was sorted by simultaneously selecting the top 25% of solutions for each objective with rank ≤ 1. The best fit solution for each objective (line) ± 1-standard deviation (gray region) for experiment 2 and 3 brackets the training data despite significant differences the training values between the two data sets parameters that are constrained by the performance of the model, and not by specific hypotheses regarding the individual values of the raw model parameters. Of course, knowledge of specific parameter values, or the relationship between parameter combinations, can be used to inform the search through either bounds or problem specific constraints (for example, as demonstrated in the first example problem).

Conclusions
In this software note, we presented JuPOETs, a multiobjective technique to estimate parameter ensembles in the Julia programming language. JuPOETs is open source, and available for download under an MIT license from the JuPOETs GitHub repository at https://github. com/varnerlab/POETs.jl. We demonstrated JuPOETs on a suite of algebraic test problems, and a proof-of-concept ODE based biochemical model. While JuPOETs outperformed (and was significantly more flexible) than the previous Octave implementation, there are several areas that could be explored further. First, JuPOETs should be compared with other multiobjective evolutionary algorithms (MOEAs) to determine its relative performance on test and real world problems. Many evolutionary approaches e.g., the non-dominated sorting genetic algorithm (NSGA) family of algorithms, have been adapted to solve multiobjective problems [36,37]. However, since there is a lack of open source Julia implementations of these alternative approaches, we did not benchmark the relative performance of JuPOETs in this note. One advantage that JuPOETs may have when compared to a strictly evolutionary approaches, is the inclusion of a local refinement step (hybrid mode), which temporarily reduces the problem to a single objective formulation. Previously, POETs run in hybrid mode led to better convergence on a proof-of-concept signal transduction model compared to the same approach without the hybrid refinement step [29]. Other hybrid multiobjective methods have also been shown to be more efficient than evolutionary approaches alone, for a variety of biochemical optimization problems [24,38]. Thus, there are several different algorithms that we can use to benchmark, and improve the performance of JuPOETs, after we implement them in Julia. Another strategy to improve the performance of JuPOETs is to reduce the number (or cost) of function evaluations that are required to obtain optimal or near optimal solutions. For example, in many real world parameter estimation problems, the bulk of the execution time is spent evaluating the objective functions. One strategy to improve JuPOETs performance could be to optimize surrogates [39], while another would be parallel execution of the objective functions. Currently, JuPOETs serially evaluates the objective function vector. However, parallel evaluation of the objective functions e.g., using the parallel Julia macro or other techniques, could be implemented without significantly changing the JuPOETs run loop. Taken together, JuPOETs demonstrated improved flexibility, and performance over POETs in parameter identification and ensemble generation for multiple objectives. JuPOETs has the potential for widespread use due to the flexibility of the implementation, and the high level syntax and distribution tools native to the Julia programming language.

Availability and requirements
JuPOETs is open source, available under an MIT software license. The JuPOETs source code is freely available from the JuPOETs GitHub repository at https:// github.com/varnerlab/POETs.jl. All samples used in this study are included in the sample/biochemical and sample/test_functions subdirectories of the JuPOETs GitHub repository. JuPOETs can be run on all common. Operating system environments: (Linux, Mac OS, Windows).