Problem statement

Given a model of a nonlinear dynamic system and a set of experimental data, the parameter estimation problem consists of finding the vector of decision variables p(unknown model parameters) that minimizes a cost function that measures the goodness of the fit of the model predictions with respect to the data, subject to a number of constraints. The output state variables that are measured experimentally are called observables.

where n _ε is the number of experiments,$n_o^{\epsilon }$ is the number of observables per experiment, and$n_s^{\epsilon ,o}$ is the number of samples per observable per experiment. The measured data will be denoted as${ym}_s^{\epsilon ,o}$ and the corresponding model predictions will be denoted as$y_s^{\epsilon ,o}$ (p). Finally, W is a scaling matrix used to balance the contributions of the observables (usually by scaling each temporal series with respect to its maximum value).

where g is the observation function, x is the vector of state variables with initial conditions x₀, f is the set of differential and algebraic equality constraints describing the system dynamics (that is, the nonlinear process model), h _eq and h _in are equality and inequality constraints that express additional requirements for the system performance, and p ^L and p ^U are lower and upper bounds for the parameter vector p.

As was mentioned in the Background section, in systems biology models this problem is often multimodal (nonconvex), due to the nonlinear and constrained nature of the system dynamics. Hence, standard local methods usually fail to obtain the global optimum. As an alternative, one may choose a multistart strategy, where a local method is used repeatedly, starting from a number of different initial guesses for the parameters. However, this approach is usually not efficient for realistic applications, and global optimization (GO) techniques, such as the ones presented in the next section, need to be used instead.

Global optimization via enhanced scatter search

Scatter search can be classified as an evolutionary optimization method which, like many other metaheuristics, arose in the context of integer optimization[26]. Several adaptations to continuous problems have been developed in recent years. The application of the method to parameter estimation in systems biology problems has provided excellent results[22].

Scatter search is a population-based algorithm which, compared with other methods like genetic algorithms, makes use of a low number of population members, called “Reference Set” (RefSet) in this context. Besides, scatter search includes the so-called “improvement method” which usually consists in a local search to speed-up the convergence to optimal solutions.

To illustrate how the method works, Figures1.(a-d) represent the contour plots of an arbitrary objective function to be optimized (e.g., a least squares function in the case of parameter estimation), together with the solutions which constitute the RefSet and their evolution.

The method starts by creating an initial population of diverse solutions within the search space (Figure1.a). From them, the best solutions in terms of quality are selected to create the initial RefSet (in red in Figure1.b). The rest of solutions in the RefSet are chosen randomly (in green in Figure1.b), although in some advanced implementations they may be chosen following a criterion of maximum diversity to cover a broader area of the search space. Once the RefSet solutions have been selected, the remaining diverse solutions generated in the first step are deleted and the so-called “subset generation method” starts. It consists in combining sets of RefSet solutions to create new ones. In our scheme, we combine every solution in the RefSet with the other RefSet by generating hyper-rectangles defined by their relative position and distance, and creating new solutions inside them. This step is illustrated in Figure1.c. Each RefSet member will create one hyper-rectangle with every other RefSet member, thus generating a number of b−1 new solutions (being b the RefSet size) per each RefSet member. After this procedure, all the RefSet members have a set of offspring solutions around them covering different distances and directions (Figure1.d). If the best offspring solution outperforms the parent solution (which is the case of Figure1.d with the best offspring solution represented as a star), the replacement (or update) is carried out. Otherwise, the same RefSet member will stay in the population in the next iteration. Although the procedure is illustrated for just one RefSet member, the same scheme applies for every solution in the RefSet. After the RefSet update, the algorithm applied the “solution combination method” again, repeating the process until a stop criterion is met.

In this work, a novel extended implementation of the enhanced scatter search metaheuristic (eSS)[24, 25], has been used. Its pseudocode is shown in Algorithm 1.

Algorithm 1

Basic pseudocode of eSS

Set parameters: dim_refset, local.n 2, balance, ndiverse

Initialize n _stuck , neval

Create set of ndiverse solutions

Generate initial RefSet with dim_refset solutions with high quality and random elements

repeat

Sort RefSet by quality [x¹,x²,…,x ^dim_refset ]so that f(x ⁱ ) ≤ f(x ^j ) where i,j∈[1,2,…,dim_refset] and i < j

if max$\left(\mathit{\text{abs}}\left(\frac{x^i-x^j}{x^j}\right)\right)\le \epsilon$ with i < j then

Replace x ^j by a random solution

end if

for i = 1 to dim_refset do

Combine x ⁱ with the rest of population members to generate a set of dim_refset new solutions,

offspring ⁱ

$x_{\text{off}}^i=$ best solution in offspring ⁱ

if$x_{\text{off}}^i$ outperforms x ⁱ then

Label x ⁱ

x_temp = x ⁱ

improvement = 1

∧ = 1

while $f\left(x_{\text{off}}^i\right)<f\left(x_{\text{temp}}\right)$ do

Create a new solution, x_new, in the hyper-rectangle defined by$\left[x_{\text{off}}^i-\frac{x_{\text{temp}}-x_{\text{off}}^i}{\wedge },x_{\text{off}}^i\right]$

$x_{\text{temp}}=x_{\text{off}}^i$

$x_{\text{off}}^i=x_{\text{new}}$

improvement = improvement + 1

if improvement = 2 then

∧ = ∧/2

improvement = 0

end if

end while

end if

end for

offspring=$\left(\begin{array}{c}{\mathbf{\text{offspring}}}^1\\ {\mathbf{\text{offspring}}}^2\\ \dots \\ {\mathbf{\text{offspring}}}^n\end{array}\right)$ with n = dim_refset

if neval ≥ local.n 2 and at least one local solution has been found then

Sort offspring solutions by quality [x_1,1,x_2,1,…,x_m,1] where m = n(n−1) and f(x_i,1) ≤ f(x_j,1) where i,j∈[1,2,…,m] and i < j

Compute the minimum distance between each element i ∈ [1,2,…,m] in offspring and all the local optima found so far, d(x _i )

Sort offspring solutions by diversity [x_1,2,x_2,2,…,x_m,2] with d(x_i,2) ≥ d(x_j,2) where i,j∈[1,2,…,m] and i < j

Choose z as the offspring member balancing quality and diversity according to balance

Perform a local search from z to obtain z^∗

if z^∗ is a new local optimum then

Update list of found local optima adding z^∗

end if

if f(z^∗) < fbest then

Update xbest, fbest

end if

neval = 0

end if

neval = neval + function evaluations of current iteration

Replace labeled population members by their corresponding$x_{\text{off}}^i$ and reset their corresponding n _stuck (i)

Add one unit to the corresponding n _stuck (j) of the not labeled population members

if any of the n _stuck values > nchange then

Replace those population members by random solutions and reset their n _stuck values

end if

until stopping criterion is met

Apply final local search over xbest and update it

The innovative mechanisms implemented in different parts of the method as well as the appropriate local search selection, make the algorithm specially suitable for solving large-scale optimization problems. These mechanisms, which address the most important question in global optimization (i.e., the balance between intensificacion, or local search, and diversification, or global search) are depicted in the following points:

Cooperative optimization

A common way of reducing the computational cost of a given algorithm is to parallelize it. This entails performing several of its calculations simultaneously in different processors. Except for some special cases, which are usually called “embarrasingly parallel”, it is not trivial to separate the problem into parallel tasks. Here we are interested in the parallelization of a metaheuristic optimization algorithm, the enhanced scatter search method described in the previous section.

Before we describe the proposed methodology, we clarify the use of a few words in order to avoid confusion. We refer to a context where there are a number of processors available for computing, which can work simultaneously (in parallel). All of the processors are assumed to be physically identical, that is, they have the same hardware characteristics. A program running in a processor is called thread; both terms are interchangeable. In our methodology we use a master processor, with the task of coordinating the remaining processors (slaves). All of the slave processors implement the same program, which means they perform similar calculations, using different data sets and different settings.

There are not many examples of research in the area of parallelization of metaheuristics exploiting cooperation in the way we present here. In[30] three different parallelizations of the scatter search were proposed and compared. In the first one, each local search was parallelized. In the second one, several local searches were carried out in parallel. In these two cases it was possible to reduce the computational time. A third option was to run the algorithm for different populations in parallel, but without any communication between threads.

More recently[31] a more ambitious approach was presented, where different kinds of metaheuristics were combined in order to profit from their complementary advantages. Information was shared between algorithms and it was dynamically tuned, so the most succesful algorithms were allowed to share more information. For a detailed review of the state of the art in parallelization of metaheuristics, see[32].

Sharing information can be seen as a way of cooperating. Such cooperation produces more than just speed-up: it can change the systemic properties of an algorithm and therefore its macroscopic behaviour. This was acknowledged in[33], where four parameters that may affect this behaviour were identified: (i) information gathered and made available for sharing; (ii) identification of programs that share information; (iii) the frequency of information sharing among the sequential programs; and (iv) the number of concurrent programs.

Each of the η threads has a fixed degree of “aggressiveness”. “Conservative” threads make emphasis on diversification (global search) and are used for increasing the probabilities of finding a feasible solution, even if the parameter space is “rugged” or weakly structured. “Aggressive” threads make emphasis on intensification (local search) and speed up the calculations in “smoother” areas. Please note that both conservative and agressive threads should be able to locate the globally optimal solution of the class of problems considered here, but their relative efficiency will vary depending on the region of the parameter space where the individual threads are operating. In this way, we achieve a synergistic gain. Communication, which takes place at fixed time intervals, enables each thread to benefit from the knowledge gathered by the others. This knowledge may include not only information about the best solution found so far, but also about the sets of diverse parameter vectors that may be worth trying for improving the solution. Thus this strategy has several degrees of freedom that have to be fixed: the time between communication (τ), the number of threads (η), and the strategy adopted by each thread (Ψ). These adjustments should be chosen carefully depending on the particular problem we want to solve. In the following lines we give some guidelines for doing this.

Settings (strategy) for each thread, Ψ

For cooperation to be efficient, every thread running in parallel must perform different calculations. Due to the random nature of metaheuristic algorithms, this can be achieved simply by choosing different random initializations, or “seeds”, for each thread. This results first in different sets of initial guesses, and subsequently in the exploration of different points of the search space. However, this is not the only difference that may be made between threads: we can also select a different behaviour (more or less aggressive) for each of the threads. An “aggressive” thread performs frequent local searches, trying to refine the solution very quickly, and keeps a small reference set of solutions. This means that it will perform well in problems with parameter spaces that have a smooth shape. “Conservative” threads, on the other hand, have a large reference set and perform local searches only sporadically. This means that they spend more time combining parameter vectors and exploring the different regions of the parameter space. Thus, they are more robust, and hence appropriate, for problems with rugged parameter spaces. Since the exact nature of the problem to consider is always unknown, it is advisable to choose a range of settings that yields conservative, aggressive, and intermediate threads. For the eSS algorithm this can be obtained by tuning the following settings:

• Number of elements in the Reference Set (dim_refset).

• Minimum number of iterations of the eSS algorithm between two local searches (local.n2).

• Balance between intensification and diversification in the selection of initial points for the local searches (balance).

• Number of diverse solutions initially generated (ndiverse).

All these settings have qualitatively the same influence on the algorithm’s behaviour: if they are large, they lead to conservative threads; if they are small, to aggressive threads. The four parameters that define the degree of aggressiveness of each of the η threads can be stored in an array Ψ_η×4. By default, we recommend to use a broad spectrum of aggressiveness using the default values 0.5·n _par < dim _refset < 20·n _par ; 0 < local.n 2 < 100; 0 < balance < 0.5; and 5·n _par < ndiverse < 20·n _par .

Figure2 shows the cooperative loop, and Algorithm 2 summarizes its pseudo-code. The text marked as SLAVE, PROCESSOR j corresponds to the code implemented in parallel by each of the slave processors; the rest of the computations are carried out by the master processor.

Algorithm 2

Basic pseudocode of the cooperative strategy

MASTER, PROCESSOR 0: Initialization

Set parameters that are common for all threads: τ, η, n _iters

Set array of aggressiveness parameters: Ψ_η×4 = {dim_refset_η×1,ndiverse_η×1,local.n2_η×1,balance_η×1}

Initialize global reference set array: Globa l _ref = []

for i = 1 to n _iters do

for j = 1 to η do

SLAVE, PROCESSOR j (concurrent processing): run optimization according to Algorithm 1

(Thread j uses aggressiveness settings vector Ψ_{j,:} and initial points Globa l _ref )

end for

if CPUtime ≥ τ then

for j = 1 to η do

Read results of thread j

if ref _j ∉ Global _ref then

Global _ref = [Global _ref ,ref _j ]

end if

end for

end if

end for

Final solution = best solution in Global _ref

Model 1: E. coli CCM

where u is an enzyme level (mmol),$k_A^M$,$k_B^M$ and$k_C^M$ are reactant constants (mM), and k^± are turnover rates (s⁻¹). The CM rate law resembles the convenience kinetics, with a slight difference for molecularities m^± ≠ 1 (in convenience kinetics, the denominator term for product C would read$1+\frac{c}{k_C^M}+{\frac{c}{k_C^M}}^2$).

The resulting model has 918 parameters: internal metabolites concentrations; activation, inhibition, and Michaelis-Menten constants; and velocity and equilibrium constants. Their values are given in the supplementary material, see Additional file2. They are estimated by least squares minimization of the difference between the model predictions and artificially generated data regarding internal metabolic concentrations (49) and fluxes (66), yielding a state vector with 115 elements.

Multistart procedures of local methods are carried out in order to estimate the nonconvexity of the problem and compare the results with those obtained with global optimization. The upper and lower bounds allowed for the parameters are, respectively, 10 times larger and 10 times smaller than the nominal values, which are shown in the Additional file2. Several local methods are tested: Direct Hill Climbing (DHC), and the MathWorks functions fminsearch (for unconstrained minimization) and fmincon (for constrained minimization). A common computation time of ten days is imposed to every one of them. Histograms of the solutions are shown in Figure3; the best value of the objective function found is f _best = 2.94·10¹, which is achieved with the DHC method. It can be noticed that for the same computation time, some algorithms perform more local searches than others. For example, with fmincon the calculation of the gradient is very costly due to the large number of decision variables, while DHC avoids these calculations and is thus capable of completing more local searches in the same time.

Suspecting that a more sophisticated method may be able to find a better solution, we test our enhanced Scatter Search (eSS) algorithm. Since the multistart calculations were carried out in a single processor, with a total running time of 10 days, it is fair to compare the results with those obtained by 10 processors in one day. We thus select 10 different settings for the eSS algorithm and launch the corresponding 10 optimizations. Since eSS provides the possibility of launching local searches, we select the local method that yields best results with the multistart procedure, DHC. After one day, the best objective function value achieved is f _best = 7.86, which represents a reduction of 73% with respect to what was obtained with a multistart of local searches.

To justify our choice of the eSS algorithm among the existing metaheuristics, we compare it with another state of the art metaheuristic such as Differential Evolution (DE). To ensure a fair comparison we set the same computation time as for the eSS and local multistart procedures: we carry out 10 DE optimizations, each with a time limit of one day. The convergence curves of the DE and eSS algorithms are shown in Figure4, which plots, for each method, the best solution found by any of its 10 threads at every time instant (note that this is different from showing 10 different convergence curves for each method). Both algorithms obtain similar results in the first few hours; after that time, however, eSS achieves fast and significant improvements whereas DE progresses more slowly. This is due to the lack of local searches, which in this problem are instrumental in refining the solution. Hence, we see that for calibrating this model it is advisable to combine the intelligent exploration of the parameter space provided by metaheuristics with the refinement of the solution provided by local searches. This supports the choice of the eSS algorithm.

The next step is to improve the eSS results by means of the parallel cooperative strategy presented in the previous section, CeSS. With this aim, the number of threads is fixed as η = 10; each of them implement the eSS algorithm with different options (dim_refset, local.n2, and balance; the fourth option ndiverse is common to all threads). In this way each of the 10 threads has a different degree of aggressiveness. Figure5 shows the convergence curves, which plot the logarithm of the objective function value against the computation time. The performance of the 10 individual, non cooperative threads (black line) is compared to the performance of the 10 cooperative threads (green line), which exchange information with τ = 2.53 (approximately every 12 hours in the hardware used). The results show that cooperation greatly speeds up the convergence: while the non-cooperative threads need 11 days to reach an objective function value of 1.88, the cooperative threads obtain the same result in approximately 1.5 days. Therefore cooperation reduces computation times by more than 85%. In 11 days, the cooperative threads achieve objective function values ranging between 8.05·10⁻¹ and 9.71·10⁻¹, thus improving the non-cooperative solution by more than 50%. More figures showing the influence of the tuning parameters η, τ are included as supplementary material, see Additional file2.

To give a better idea of the goodness of the fit, we test the ability of the calibrated model to reproduce the original data. Figure6 plots in logarithmic scale the percentage error in the estimated values of the observables, which are concentrations and fluxes. Similarly, Figure7 shows the errors in the estimated values of the parameters. It can be noticed that the model predictions for the observables (Figure6.a) are better than the fit between the nominal and estimated parameters (Figure7.a). That is, even though the optimization algorithm manages to fit the model predictions to the data quite nicely, there are some groups of parameters that are not estimated accurately. This fact reveals a lack of practical identifiability for the given model and data. The present work focuses on the performance of optimization algorithms; that is, in their ability to reduce the objective function value. This is different to the identifiability problem, which is not addressed here but can be surmounted by proper design of additional experiments. However, from the point of view of model calibration, the method presented here was able to successfully solve the inverse problem with better performance than other state of the art methods.

Model 2: E. coli metabolic model including enzymatic and transcriptional regulation

The model calibrated in the previous subsection was purely metabolic. Here we consider a previously published E. coli model[18] that also takes into account the enzymatic and transcriptional regulation layer. It consists of 47 ODEs and 193 parameters (affinity constants, specific activities, Hill coefficients, growth rates, expression rates, etc), of which 178 are considered unknown and need to be estimated. We have reformulated the model to use it in an optimization context, where the objective function to be minimized consists of the difference between the simulated concentration profiles obtained with the nominal and the estimated parameters. Upper and lower bounds for the parameters are fixed to values 10 times larger and 10 times smaller than the nominal, except for the Hill coefficients, which are assumed to be between 0.5 and 4. More information is included as supplementary material, see Additional file3.

This model was created to reproduce the way in which E. coli adapts to changing carbon sources. With this aim, it is subjected to a sequence of three consecutive environments, where the carbon source is first glucose, then acetate, and finally a mixture of both. Under these conditions, the 47 concentration profiles are sampled every 1000 seconds, for a total of 162 time points (45 hours). Therefore the overall number of available samples is 47 × 162 = 7614. The equations are integrated using LSODE (Livermore Solver for Ordinary Differential Equations)[36], which is suited for stiff systems. Integrating these equations is a computationally expensive task; one evaluation of the objective function takes several seconds.

It must be noted that in[18] a divide-and-conquer approach was adopted for estimating the parameter values of the rate equations. In this approach, the large-scale optimization problem is divided into a set of small subproblems, where only a few parameters are estimated at a time. This allows to reduce the complexity of the task and consequently the computation times. However, it is an ad hoc procedure that cannot be applied to most systems biology problems. Here we adopt a general purpose, large-scale global optimization approach that does not require any special structure of the model equations.

As with Model 1, the calibration procedure begins by carrying out several multistart procedures of local methods. The selected methods are DHC, fminsearch, and N2FB (which outperforms fmincon for this problem). We limit the number of evaluations in the fminsearch method to 2000, so that each local search lasts typically around five hours on our computers; without this limit it can go on for much longer without achieving substantial improvements. The maximum time allowed for each multistart is ten days for all of the methods. Notice that, as with Model 1, each of the algorithms carries out a different number of local searches, despite the same computational time. The best solution is f _best = 7.60·10², found by the DHC method; an histogram of the results is shown in Figure8.

Next, the enhanced Scatter Search (eSS) algorithm is tested. For this particular problem it was decided to use eSS with DHC as the local solver because, once again, it offers the best balance between quality of the solution and computational time. After selecting 10 different settings and letting the corresponding 10 optimizations run for 24 hours, the best objective function value achieved is f _best = 2.37·10² and the worst is f _best = 1.48·10³. Thus the eSS algorithm clearly outperforms any local method.

Again, we compare the performance of the eSS algorithm with the differential evolution metaheuristic (DE). We carry out 10 DE optimizations, each with a time limit of 24 hours, and do the same comparison of convergence curves performed with Model 1. As can be seen in Figure9, eSS also outperforms DE for this particular problem.

Then we test the cooperative optimization method, CeSS. With this aim, we launch η=10 cooperative threads, each of them implementing the eSS algorithm with different options (dim_refset, local.n2, balance); the interval between information sharing is τ = 2.15, which in the hardware used corresponds to 24 hours. Figure10 plots the convergence curves showing the best value found at every instant by the 10 cooperative threads compared to that found by 10 non-cooperative threads. As with Model 1, it can be seen that cooperation speeds up convergence. After 17 days CeSS reaches an objective function value of 2.51, while eSS obtains 2.58·10¹. The solution retrieved by CeSS represents a very good match between data and model predictions, as shown in Figure11. Figure12 plots in logarithmic scale the percentage error in the estimated values of the observables, while Figure13 shows the errors in the estimated values of the parameters; these plots reveal practical identifiability issues similar to those detected in Model 1. In this case, identifiability issues were actually expected, given that (i) some parameter values in the originally published model[18] are very different from those in the version available in the BioModels database[37] (model ID: 244), but (ii) despite the differences, both parameter sets yield similar simulation results. This is a clear sign of lack of identifiability. However, once again, the new method is capable of successfully solving the calibration problem, which is the main objective of this work. The practical identifiability problems are a consequence of lack of information in the considered data set and can be surmounted by proper design of additional experiments[4].

Finally, we would like to comment on the default values recommended for the search parameters of the cooperative strategy. The main potential disadvantage of using such default settings for a new problem is that the method could perform worse in two extreme situations: easy (mildly non-convex) problems, where a simple multi-start local method could perform very well, and extremely hard (highly non-convex problems) where there is no structure at all and local searches do not help. Since for each new problem we do not have any a priori information about the topology of the search space, we simply recommend a broad spectrum of values for the aggressiveness of the threads, so even in the extreme situations the cooperative strategy would be able to solve the problem. Of course, any a posteriori information about the topology, after a few runs with the default settings, can be exploited by tuning ψ. This reasoning applies for the case of new arbitrary problems, but we should add that for the particular large-scale parameter estimation problems considered above, the default recommended values have shown a good performance.

Summary

Typically, systems biology models have a large number of parameters and are highly nonlinear. Furthermore, the information content of the available experimental data is frequently scarce. As a consequence, it is hard to find the set of parameter values that gives the best fit to experimental data. This task, known as model calibration, is commonly formulated as an optimization problem, which in these cases must be solved with global optimization techniques.

Here we have presented a new method for solving this problem, called Cooperative enhanced Scatter Search (CeSS), which is specifically designed to profit from a parallel environment with several processors. Each of the processors executes a program (or “thread”) which implements the enhanced Scatter Search (eSS) algorithm[24, 25]. This is a state of the art metaheuristic capable of competing succesfully with other global optimization methods. The key feature of the strategy presented here is the cooperation between threads, which means that they exchange information at certain fixed instants. This information consists of the reference set of solutions they have found so far; it depends on the best solution found, the shape of the solution space, and the settings of each algorithm.

We have tested this strategy with two large-scale E. coli models of different sizes and characteristics. The first one models the central carbon metabolism (CCM) using a recently proposed common modular rate law[11]. The second is a previously published model that combines metabolism with a transcriptional regulatory layer[18]. In both cases the presented method clearly outperformed several state of the art techniques. The performance and capabilities of the method have also been evaluated using benchmark problems of large-scale global optimization, with excellent results.

Insights

The results for the cooperative strategy presented here shows how cooperation of individual parallel search threads modifies the systemic properties of the individual algorithms, improving its performance and outperforming other competitive methods. Importantly, it is a general purpose framework that can be easily applied to any model calibration problem. Furthermore, the underlying cooperative parallel mechanism is relatively simple and can be extended to incorporate other global and local search solvers and specific structural information for particular classes of problems.