Introduction

Development of experimental technologies in molecular biology led to accumulation of huge volumes of data relating to various levels of life organization. However, the data alone cannot be used to reconstruct the full organization of biological systems. Therefore, the interests of bioinformatics are now focused on the problems of data processing, including the problems of integration and systematization of primary experimental data and the problems of knowledge production based on mathematics and modern information technologies. The challenge of systems biology is construction of mathematical models to describe dynamic behavior of biological systems based on experimental data. Such problems involve studying a large volume of data and require software for their processing and interpreting.

The standard tools for working with biological data include access to biological databases, formalized description of biological systems, as well as visualization, simulation, parameter fitting and analysis of ODE models representing these systems. The BioUML software is an integrated environment that was developed to span all of these capabilities. Here we present optimization tools of this software intended for multi-experiment training of the models created using BioUML notation or imported in the SBML format [Hucka et al, 2003]. These tools are available both in the desktop and web editions of BioUML.

Optimization problem in BioUML

The general nonlinear optimization problem [Runarsson and Yao, 2000] can be formulated as follows: find a minimum of the objective function $\varphi \left(x\right)$, where $x$ lies in the intersection of the $\mathcal{N}$-dimensional search space

and the admissible region $ℱ\subseteq {ℝ}^{\mathcal{N}}$ defined by a set of equality and/or inequality constraints on $x$. Since the equality $g_s\left(x\right)=0$ can be replaced by two inequalities $g_s\left(x\right)\le 0$ and $–g_s\left(x\right)\le 0$, the admissible region can be defined without loss of generality as

In order to get solution situated inside $ℱ$, we minimize the penalty function

$\psi \left(x\right)={\displaystyle \sum }_{i=1}^s\max {\left\{0,g_s\left(x\right)\right\}}^2$.

The problem could be solved by different optimization methods. We implemented the following of them in the BioUML software:

Table 1 shows the generic scheme of the optimization process for these methods. SRES, MOCell, PSO and glbSolve run a predefined number of iterations $N_{it}$ considering a sequence of sets (populations) $P^i$, $i=0,\dots ,N_{it}-1$, of potential solutions (guesses). In the case of the first three methods, the size $s\in {\mathbb{N}}^{+}$ of the population is fixed, whereas in glbSolve the initial population $P^0$ consists of one guess, while the size $s_{k+1}$ of the population $P^{k+1}$ is found during the iteration with the number $k=0,\dots ,N_{it}-1$. The method ASA considers sequentially generated guesses $x^k\in \text{Ω}$, $k\in {\mathbb{N}}^{+}$, and stops if distance between $x^k$ and $x^{k+1}$ defined as Euclidean norm $x^k-x^{k+1}=\sqrt{{\displaystyle \sum }_{i=1}^{\mathcal{N}}{\left(x_i^k-x_i^{k+1}\right)}^2}$ becomes less than a predefined accuracy $\epsilon$.

All methods, excepting glbSolve, are stochastic and seek global minimum of the function $\varphi$ taking into account the admissible region $ℱ$. Thus, a guess $x\in \text{Ω}$ is more preferable than a guess $y\in \text{Ω}$ at some iteration of methods, if $\psi \left(x\right)=0$ and $\psi \left(y\right)\ne 0$ or $\psi \left(x\right)<\psi \left(y\right)$. The method glbSolve is suited to solve only the problems with $\text{Ω}\subseteq ℱ$. Values of the function $\psi$ are calculated but do not affect on the generation of potential solutions.

Implementing the optimization scheme in BioUML, we designed the OptimizationProblem interface (fig. 1) comprising the following procedures:

An abstract class OptimizationMethod provides the number of subclasses representing implementation of the foregoing methods. These subclasses involve search of optimal parameters by calling a procedure getSolution depending on the settings of optimization problem.

Figure 1. The class diagram representing implementation of the optimization process in BioUML.

Application of non-linear optimization to systems biology

We assume that a mathematical model of some biological process consists of a set of chemical species $\mathcal{S}=\left\{S_1,\dots ,S_m\right\}$ associated with variables $C\left(t\right)=\left(C_1\left(t\right),\dots ,C_m\left(t\right)\right)$ representing their concentrations, and a set of biochemical reactions $ℛ=\left\{R_1,\dots ,R_n\right\}$ with rates $v\left(t\right)=\left(v_1\left(t\right),\dots ,v_n\left(t\right)\right)$ depending on a set of kinetic constants $K$. Reaction rates are modeled by standard laws of chemical kinetics. A Cauchy problem for ordinary differential equations representing a linear combination of reaction rates is used to describe the model behavior over time:

Here $N$ is a stoichiometric matrix of $n$ by $m$. We say that $C^{ss}$ is a steady state of the system (1) if

$N\cdot v\left(C^{ss},K,t\right)=0$, $\underset{t\to \infty }{\lim }{C}_i\left(t\right)=C_i^{ss}$.

Identification of parameters $K$ and initial concentrations $C^0$ is based on experimental data represented by a set of points $C_i^{exp}\left(t_{ij}\right)$ defining dynamics of variables $C_1\left(t\right),\dots ,C_l\left(t\right)$, $l\le m$, at given times $t_{ij}$, $j=1,\dots ,r_i$, where $r_i$ is the number of such points for the concentration $C_i\left(t\right)$, $i=1,\dots l$. The problem of parameter identification consists in minimization of the function of deviations defined as the normalized sum of squares [Hoops et al, 2006]:

where normalization factors ${\omega }_{\min }/{\omega }_i$ with ${\omega }_{\min }=\underset{i}{\min }{\omega }_i$ are used to make all concentration trajectories have similar importance. The weights ${\omega }_i$ are calculated by one of the formulas on experimentally measured concentrations: ${\omega }_i^{sq}=\sqrt{r_i^{-1}\cdot {\displaystyle \sum }_j{\left(C_i^{exp}\left(t_{ij}\right)\right)}^2}$ (mean square value), ${\omega }_i^{mean}=\left|r_i^{-1}\cdot {\displaystyle \sum }_j{C}_i^{exp}\left(t_{ij}\right)\right|$ (mean value) and ${\omega }_i^{st}=\sqrt{{\omega }_i^{sq}\cdot {\omega }_i^{sq}-{\omega }_i^{mean}\cdot {\omega }_i^{mean}}$ (standard deviation).

If we want to consider additional constrains

holding for concentrations $C\left(t\right)$ and parameters $K$ for some period of time $t\in \left[t_s^{start},t_s^{end}\right]$, the penalty function is defined as

This function assumes summation of values $g_s\left(C,K\right)$ in the nodes of grid defined by an ODE solver to find a numerical solution of the system (1).

In the particular case, experimental data could be represented by steady state values of species concentrations. Then functions $\varphi$ and $\psi$ have the simpler forms:

$\varphi \left(C^0,K\right)={\displaystyle \sum }_{i=1}^l\frac{{\omega }_{\min }}{{\omega }_i}\cdot {\left(C_i^{ss}-C_i^{exp\_ss}\right)}^2$, $\psi \left(C^0,K\right)={\displaystyle \sum }_{s=1}^p{\left(\max \left\{0,g_s\left(C^{ss},K\right)\right\}\right)}^2$

where $C_i^{exp\_ss}$ and $C_i^{ss}$, $i=1,\dots ,l$, denote experimental and simulated steady state values.

Typically, researchers want to perform evaluation of model parameters using experimental data obtained with different experimental conditions, i.e. different initial concentrations $C^{01},\dots ,C^{0k}$ of species. In such case, we will consider the functions

Implementation of the parameter estimation process in BioUML

Initiation of the parameter estimation process requires definition of many details including specification of the search space, the admissible region, settings of numerical methods to solve ODE system and optimization problem, links to the files with experimental data and model description, etc. In order to structure this information, we designed an appropriate hierarchy of classes (fig. 2) taking into account the following rules:

The main class Optimization in fig. 2 comprises definition of the optimization method and parameters including the model parameters to fit, parametric constraints holding for given time intervals, and experiments. The following fields keep information about each experiment:

Figure 2. The class diagram of the optimization plug-in in BioUML.

An object of the OptimizationMethod class contains a set of parameters (methodParameters) including links to the file with the processed model (diagramPath) and the directory to save results when the parameter estimation will be completed (resultPath). The estimation process begins immediately after the object control gets command doRun. Goodness of the current fit is defined through the object optimizationProblem providing correct simulation of the model. Firstly, evaluation of functions (2) and (4) is performed separately for all initial states of the model by calling a method testGoodnessOfFit for each object of the class SingleExperimentParameterEstimation associated with the certain experiment. Then the total values of these functions are calculated in the class ParameterEstimationProblem by the formulas (5).

The parameter estimation process is optimized using the following technologies:

Fig. 3 shows a graphical user interface of optimization methods in the BioUML workbench (desktop edition). The upper left panel includes a list of methods. Description of the selected method is provided below. The upper right panel defines the search space. Under it you can find the tab panel with settings of optimization problem. In fig. 3 the selected tab contains method parameters and fields displaying intermediate values of objective/penalty functions and the number of passed evaluations. The next tab includes description of experimental data and specifies settings for all experiment fields listed above.

Figure 3. The user interface of the optimization plug-in provided by the BioUML software.

Web edition optimization

BioUML web edition (http://ie.biouml.org/bioumlweb/) is a web application providing access to BioUML tools and data via the Internet. The user can manipulate the data stored on the BioUML server and run analyses through the web browser. The web interface is a set of HTML pages with interactive JavaScript content. Ajax technology is used to communicate with the server and process user activities without page reloading.

The web edition provides a set of optimization tools, like in the workbench edition. The user can set up the boundaries and initial values of optimization parameters, select and adjust an optimization method, manage experimental data, set up constraints. After all parameters are adjusted, optimization process can be launched. Optimization results can be saved and then viewed as a table of fitted parameter values. Graphical representation of the optimization process is also available.

Analysis of the methods convergence

Stochastic methods (SRES, MOCell, PSO and ASA) for global optimization rely on probabilistic approaches and have weak theoretical guarantees of convergence to the global optimum. However, they can locate its vicinity with relative efficiency [Moles et al, 2003]. In contrast, deterministic method glbSolve guarantees global optimality, if the objective function is continuous or at least continuous in the neighborhood of a global optimum [Björkman and Holmström, 1999]. However, it can not solve general global optimization problems with certainty in finite time.

To analyze convergence rate of the implemented methods, we considered a reaction chain (fig. 4, table 2) extracted from the model by L. Neumann et al. [Neumann et al, 2010] and representing activation of caspase-8 triggered by the receptor CD95 (APO-1/Fas).

Figure 4. The test model of caspase-8 activation.

We performed estimation of parameters using the search space defined as

$0\le k_1\le 1$, $0\le k_5\le 0.1$, $0\le k_i\le {10}^{-3}$, $i=2,3,4,$

where upper bounds were chosen based on the order of magnitude of parameter values proposed by authors of the original model. Initial values of variables were fixed according to [Neumann et al, 2010]:

$C_{CD95L}\left(0\right)=113.220$, $C_{CD95R:FADD}\left(0\right)=91.266$, $C_{pro8}\left(0\right)=64.477$.

Estimation was based on the experimental data obtained by Neumann et al. for procaspase-8 and its cleaved products p43/p41 and caspase-8 (table 3).

Fig. 5 shows dependence of the objective function mean values on the number of considered guesses for 100 runs of the optimization process. Statistics of the best, mean and worst values of $\varphi$, as well as the best guesses found by the methods after consideration of 10⁷ guesses are listed in the tables 4 and 5 correspondently.

As can be seen from these tables, the best result was obtained by the particle swarm optimization (PSO) and the cellular genetic algorithm (MOCell). Methods SRES, MOCell and PSO found similar solutions. Methods ASA and glbSolve found dissimilar values of parameters $k_1$ and $k_2$ resulting in lower efficiency compared to the first three methods.

For comparison, some test cases considered in the study by Moles et al. [Moles et al, 2003] resulted in superiority of SRES. However, the authors did not explore such methods as genetic algorithms and particle swarm optimization.

Figure 5. Dynamics of mean values of the objective function for 100 runs of the optimization process. The best value obtained by the particle swarm optimization is marked by the red line.

Comparison of parameter estimation features of BioUML with other software

We compared optimization tools of BioUML with the following software applications assigned for analysis of biochemical networks and supporting the procedure of model fitting:

Details of the comparison are given in table 6.

Further, we compared computation speed of the optimization methods implemented in BioUML and COPASI. For this purpose, we considered a series of test cases. A brief description of the models used in these test cases is provided below. For more details, including specification of experimental data and fitting parameters, see Additional file 1 of the supplementary materials.

In the first test case, we analyzed three models of CD95-induced caspase-8 activation constructed on the basis of the model by Neumann et al. [Neumann et al, 2010] with varying degrees of detail (fig. 6, A-C). The second test case was proposed by Mendes et al. [Mendes et al, 2009] and corresponded to the model of the MAP kinase cascade (fig. 7) developed by Kholodenko et al. [Kholodenko, 2000]. Finally, we tested the model of Bagci et al. [Bagci et al, 2006] representing the mitochondria-depended apoptosis resulting from the cooperative formation of heptameric apoptosome complex and activation of caspase-9 and caspase-3 (fig. 8).

Figure 6. The models of caspase-8 activation constructed based on the model by Neumann et al. with the varying number of species: 7 (A), 13 (B), and 18 (C).

Figure 7. The model of the MAP kinase cascade constructed by Kholodenko et al.

Figure 8. The model of the mitochondria-depended apoptosis proposed by Bagci et al.

Analyzing the number of the objective function evaluations per second for these test cases, we found that BioUML showed a better result than COPASI (fig. 9).

Figure 9. The number of the objective function evaluations per second for different test cases in COPASI and BioUML: the model by Neumann et al. including 7 (A), 13 (B), and 18 (C) species; the model by Kholodenko et al. consisting of 8 species (D); the model of Bagci et al. consisting of 32 species (E).

Discussion

In this paper, we considered parameter estimation tools of the BioUML software. These tools can be applied to biological systems characterized with a set of ODEs. The fitting process is based on experimental time course or steady state measurements, and assumes minimization of the function of error between these measurements and the corresponding model prediction. We implemented several stochastic and deterministic global optimization methods as new plug-in for BioUML. None of these methods is effective for all cases. Nevertheless, on the basis of our observations, we concluded that adaptive simulated annealing can be used when it is necessary to quickly find the vicinity of the solution. In the case when adequacy of solution is more preferable than the rate of convergence, it is better to use such methods as MOCell, PSO and SRES.

Parameter fitting is an important part of the quantitative biological modeling. However, if the model includes more elements than are necessary to approximate experimental data with the given accuracy, we face the problem of overfitting [Hawkins, 2004]. In this case, there is no overall best solution and it is expedient to find distribution of the parameter values which are compatible with observed experimental dynamics. For this purpose, we should run parameter estimation process many times (it is better with different optimization methods) and evaluate bounds for all parameters. Implementation of this technique in BioUML is a task for the future work.

We successfully applied our optimization plug-in for creation of the combined model of CD95 and NF-κB signaling pathways [Kutumova et al, 2013], where the problem of parameter overfitting was solved using the methodology of model reduction [Gorban et al, 2010].

References