Performance of objective functions and optimisation procedures for parameter estimation in system biology models

Abstract

Mathematical modelling of signalling pathways aids experimental investigation in system and synthetic biology. Ever increasing data availability prompts the development of large dynamic models with numerous parameters. In this paper, we investigate how the number of unknown parameters affects the convergence of three frequently used optimisation algorithms and four objective functions. We compare objective functions that use data-driven normalisation of the simulations with those that use scaling factors. The data-driven normalisation of the simulation approach implies that simulations are normalised in the same way as the data, making both directly comparable. The scaling factor approach, which is commonly used for parameter estimation in dynamic systems, introduces scaling factors that multiply the simulations to convert them to the scale of the data. Here we show that the scaling factor approach increases, compared to data-driven normalisation of the simulations, the degree of practical non-identifiability, defined as the number of directions in the parameter space, along which parameters are not identifiable. Further, the results indicate that data-driven normalisation of the simulations greatly improve the speed of convergence of all tested algorithms when the overall number of unknown parameters is relatively large (74 parameters in our test problems). Data-driven normalisation of the simulations also markedly improve the performance of the non-gradient-based algorithm tested even when the number of unknown parameters is relatively small (10 parameters in our test problems). As the models and the unknown parameters increase in size, the data-driven normalisation of the simulation approach can be the preferred option, because it does not aggravate non-identifiability and allows for obtaining parameter estimates in a reasonable amount of time.

Fig. 1

Relations between parameter estimates and non-identifiability. a Clustergram visualising the relations between the parameter estimates. b Scatter plot illustrating that the space occupied by the estimates is a low-dimensional manifold: here a 1D curve in 10D space; shown is a projection in 3D (blue dots), and 2D (grey dots). c Principal component analysis of the parameter estimates. Bars illustrate the number of principal components required to explain the variability of the parameter estimates. The colours indicate how much variance is explained by each principal component (PC). The data were brought onto the same scale by normalising each parameter with respect to the best estimate from GLSDC-DNS-LS (see Methods). The number of PCs required indicates the dimensionality of the estimated parameter space (linear approximation), thus indicating the degree of non-identifiability

Fig. 2

Minima reached by LevMar SE and LevMar FD after the optimisation terminated. Boxplots show the median, 25th and 75th percentile, and extreme points (dots) outside 1.5 times the interquartile range (whiskers) of 96 independent runs. a STYX-1-10 optimisation problem. b EGF/HRG-8-10 optimisation problem. DNS normalisation of the simulations, LL log-likelihood, LS least squares, SF scaling factors. To facilitate the comparison, we always report the log-likelihood values on the y axis, even when LL was optimised

Fig. 3

Convergence speed of LevMar SE and LevMar FD. Boxplots (n = 96) as in Fig. 1. a Required number of function evaluations to terminate and b computation time to terminate for the STYX-1-10 problem. c Required number of function evaluations to terminate and d computation time to terminate for the EGF/HRG-8-10 problem. DNS normalisation of the simulations, LL log-likelihood, LS least squares, SF scaling factors

Fig. 4

Convergence of the three optimisation algorithms for a, b the STYX-1-10 problem, c, d the EGF/HRG-8-10 problem and e, f the EGF/HRG-8-74 problem. The plots show the median (thick line), and 25th and 75th percentiles (thin lines) of the objective function minima (least squares values) over the computation time from independent 96 runs. a, c, e Normalisation of the simulations (DNSs). b, d, f Scaling factors (SFs)

Fig. 5

Minima reached by different algorithm-objective function combinations after a set time for a, b the STYX-1-10 problem, c, d the EGF/HRG-8-10 problem and e the EGF/HRG-8-74 problem. Times are a after 3 min, b after 8 min, c after 10 min, d after 2 h and e after 24 h of optimising. Boxplots as in Fig 1. DNS data-driven normalisation of the simulations, LL log-likelihood, LS least squares, SF scaling factors. To facilitate the comparison, we always report the LS values on the y axis, even when LL was optimised