DBSolve Optimum is significantly improved and extended successor of well known simulation software DBSolve5 [5]. Program has been written using C++ language and compiled with Borland Builder C++ compiler. DBSolve Optimum framework consists of two parts: DBSolve Optimum Developer Environment (DODE) и DBSolve Optimum Player (DOP). DODE is designed for creation, editing, simulation, analyzing and visualization of kinetic models. DOP is focused on animation of static visual map of the kinetic model on the basis of simulation results (time series and dependence of steady state on parameter).

Figs 1 and 2 represent architecture of DODE and DOP, correspondently. DODE workflow (Fig. 1) consists of data input (IO module), creation and editing model (Model Construction module), its verification and analysis (Solver module) and output of simulation results (Output module) and their visualization (Visualization module). Work with DODE starts with creation of new model or download existing file in "SLV" or "SBML" format. "SLV" (SoLVe) is internal format of DBSolve Optimum which includes both detail information about kinetic model and description of user settings. "SBML" (Systems Biology Markup Language) format [6] is designed for storage of kinetic models and their exchange between modelers which use various software packages for model development and analysis. DODE supports SBML version L2V4 (Level 2 Version 4) and uses libSBML library for Import/Export of SBML files.

### Model construction module of DODE

Development of kinetic model starts with construction of the stoichiometric matrix of corresponding biochemical system. DODE enables user either to enter the matrix directly or download it via "RCT" (ReaCTions) file. This is ASCII (American Standard Code for Information Interchange) text file with delimiters. Each string of this file describes stoichiometry of a reaction of the kinetic model. RCT format have been realized to simplify input of stoichiometric matrixes. The description of this format is given in additional file 1: supplementary materials. On the basis of the stoichiometric matrix DODE creates template kinetic model which reaction rates are described in accordance with mass action law. ODE (Ordinary Differential Equation) system and values of model parameters and variables are presented in the form of two plain-text files. Right hand sides of ODE system including reaction rates, conservation laws and explicit functions are defined in "RHS" (Right Hand Side) file. Initial values of all variables of the model and values of all other parameters are defined in "IV" (Initial Values) file. Each line of the "IV" and "RHS" files has the following syntax: "variable = expression;". Each line of the "RHS" and "IV" should be separated and ended with a comma point delimiter. In comparison to the previous version (DBSolve 5) DBSolve Optimum allows to use conditional operators such as: if (condition) {operators} else {operators}. This possibility allows to simulate piecewise continuous functions into the right hand sides of the differential equations. At the next stage of model construction user can change the rate equations and values of parameters manually by means of "IV editor", "RHS editor" to take into account all peculiarities of the dynamics and regulations of the biochemical system.

### Model analysis and calculation modules of DODE

Analysis of the kinetic model, its fitting to experimental data and simulation of verified model is accomplished by means of following calculation modules: "Ode solver", "Implicit solver", "Explicit solver", "Bifurcation analysis", "Fitter". Operation of each calculation module is specified by a set of corresponding control parameters (see user manual for details). Before running calculation DODE parses and validates "RHS" and "IV" texts entered by user and produces byte code. This code is further used to calculate right-hand sides of the ODE system corresponding to the model. Parser program code of DODE has been generated using Flex and Bison tools (free lexical and grammar generators available at [7] and [8]). Below we briefly describe the calculation modules of DODE.

#### ODE Solver

To describe dynamics of the kinetic models DBSolve Optimum uses a popular LSODE (Livermore Solver for Ordinary Differential Equations) algorithm [9]. These methods have special subroutines for getting output for user-defined time points, which is essential for fitting algorithms. In additional file 1: supplementary materials in section C an example of calculation of time dependence of NADH (Nicotinamide Adenine Dinucleotide (reduced form)) production by α-ketoglutarate dehydrogenase reaction is presented.

#### Explicit Solver

Users may have their own particular equation which they require solving and wish to be applied to a set of experimental data. DBSolve Optimum offers the facility to encode and solve such "explicitly stated" formulae. They should be typed at the bottom of the RHS window in the section "Explicit Function" and then solved by turning to the "Explicit Solver" tabbed page, where the appropriate variable, interval of variation of parameters and an initial step can be entered. In additional file 1: supplementary materials in section C an example of calculation of Succinate thiokinase initial rate dependence on the concentration of its substrates and products is presented.

#### Implicit Solver

This method allows the user to trace the changes in the steady state of the system as a result of variation of one or more of its parameters. This procedure is very useful for determining any functional dependencies (such as overall steady-state flux, control coefficients, product concentration, some parameters of the model) against any external (substrate concentrations) or internal (enzyme concentrations) or some model parameter. It is especially useful in the case of non-linear algebraic systems which have no explicit solution or have multiple or unstable solutions. DBSolve Optimum includes a general continuation procedure, based on a tangent predictor continuation scheme [10]. A modified Newton corrector is employed which makes adaptive step sizes on the basis of estimates from the current tangent and secant vectors. This minimizes the possibility of jumping from one branch of a curve to another, and allows the user to optimize the next step size according to computed points on the curve. In DBSolve Optimum 3D plot feature for "Implicit solver" has been realized. In additional file 1: supplementary materials in section C dependence of stationary glutamate consumption rate on concentration of external glutamate calculated from kinetic model of mitochondrial Krebs cycle is presented.

#### Bifurcation Analysis

Bifurcation theory is a more systematic and general theory of non-linear systems than the standard, steady-state analysis of metabolic networks. Computation of one or two-parameter bifurcation diagrams can quickly inform the user about what is possible for, or prohibited by, a particular type of non-linear model [11]. To calculate one and two-parameter diagrams of Equilibrium, Fold, Hopf, Flip and Focus-node bifurcations DBSolve Optimum uses numerical methods similar to LOCBIF (LOCal BIFurcation) [10]. All algorithms have been rewritten in "C ++" and modified to integrate with the DBSolve Optimum object-oriented environment. The Bifurcation Analyzer uses the same numerical continuation code as the Implicit Solver, but it is expanded with routines for the evaluation of bifurcation functions and calculating eigenvectors. Bifurcations have been found at points where black rectangles are drawn on the plot. Example of application of the solver can be found in [1].

#### Fitter

This method can be used to fit a model to experimental data (thereby discovering the values with appropriate error margins of the models parameters under the conditions of the experiment).

The fitting/optimization can exploit either a zero-order [12] or first-order [13, 14] algorithm. Fitting procedure often encounters difficulties caused by multiple minima, which may be a particular problem when many parameters are fitted. The "best" fit might not be easily found; however, to check the quality of the procedure, the standard deviation and confidence intervals for every active parameter as well as an ANOVA (ANalysis Of VAriance) table are shown in the "Message window" to help users make their assessment. When fitting to experimental data, the objective residual function between theoretical and experimental points is calculated according to a least square or absolute value (modulus) formula. These are defined by the following equations:

\text{F}0=\sum {(\text{Yti}-\text{Yei})}^{\text{2}}

\text{F}0=\sum |\text{Yti}-\text{Yei}|

\text{F}0=\sum {(\text{Yti}-\text{Yei})}^{\text{2}}/{\text{Yei}}^{\text{2}}

\text{F}0=\sum |\text{Yti}-\text{Yei}|/|\text{Yei}|

where Yti and Yei are the theoretical and experimental values, respectively and F0 is discrepancy.

In additional file 1: supplementary materials in section C examples of fitting of experimental data by means of ODE Solver, Explicit Solver and Implicit Solver are presented.

### Model visualization modules

Simulation results and schematic representation of biochemical system of the model can be visualized by DBSolve Optimum in two possible ways.

#### Conventional visualization of modeling results and its implementation in DODE

Each calculation module can either save calculation results as text file of CSV (Comma Separated Value) format or plot them as a graph. The second option has been implemented in DODE on the basis of TeeChart package [15]. In comparison to DBSolve5, 3D plot feature for "ODE solver" and "Implicit solver" has been implemented in DODE. This feature allows user to create a plain text file with three columns (see more detail description in additional file 1: supplementary materials): values of time (first column), values of the parameter (second column) and values of variable. This file can be used as input for any program (for example Excel, TechPlot) to make a 3D chart.

#### Dynamic visualization of model simulation results and its implementation in DBSolve Optimum Framework

To understand behavior of the biological system as a whole we have developed special visualization technique allowing to animate simulation results (Dynamic Visualization) and implemented it in DBsolve Optimum. The main idea of the technique consists of (1) construction of visual map of the biological system (and corresponding model) and (2) animation of the visual map, i.e. reproducing dynamic changes in concentrations and fluxes by altering shapes and volumes of geometrical objects corresponding to these model variables. Graphical objects corresponding to biological entities and reactions of the model can change their volume and shape in accordance to the calculated values of corresponding concentrations and fluxes. Dynamic Visualization (animation) can be used to represent both kinetic (time dependent) response of the model to any perturbations and dependences of steady state concentrations and fluxes to any model parameters. This Dynamic Visualization technique has been implemented in DBSolve Optimum Framework in following manner. Initial construction of the static visual map (pathway) for kinetic model and generation of data file with simulation results has been implemented in DODE as Dynamic Visualization Module (DVM) (see Fig. 3 for interface of the DVM). Animation of the visual map on the basis of the simulation results has been implemented in DOP (see Fig. 4 for DOP interface). Architecture of DVM is shown in Fig. 1 (see Visualization unit) and architecture of DOP is presented by Fig. 2. Both DVM and DOP are based on OpenGL library [16] to draw graphical objects. DVM uses stoichiometric matrix of the model as input data to construct initial visual map (see Layouter module in Fig.1). At the next step user can edit and annotate the initial visual map. For example, user can draw the arrows and nodes ("Arrow animation" mode) or bars (in "Bar animation" mode) to arrange the graphic objects in desirable manner.

Output of the DVM is visual pathway map (as XML file) and calculation results produced by DODE solvers and saved as special plain text file with "PLT" (PLain Text) extension (see visualization block on Fig.1). XML (eXtensible Markup Language) file is designed to store geometry of basic graphic objects. The XSD (XML Schema Definition Language) schema of this XML file can be found on our website [17]. The PLT file is designed to store simulation results produced by ODE and/or Implicit Solvers (see below for detailed description of the DBSolve Optimum functionalities). PLT file is an ASCII delimited text file. First row of the PLT file specifies the number of metabolites of the kinetic model - N1. Second row of the file specifies number of reactions of the kinetic model - N2. Third row of the file specifies titles of the data columns written below. For example, "X[0]" - time, "Metabolite1" - concentration of metabolite1, ...., "MetaboliteN1" - concentration of metaboliteN1, "Reaction1" - flux through reaction1, ..., "ReactionN2" - - flux through reaction1. First data column of the PLT file specifies time steps t_i (in case of usage of ODE Solver as generator of numerical data) or steps in parameter change, p_i (in case of usage of Implicit Solver as generator of numerical data). Next N1 columns of the PLT file specify changes in concentrations and next N2 columns specify changes in fluxes.

Static scheme of the biochemical system (visual map saved as XML file) and simulation results (saved as PLT file) are input data for DOP. On the basis of the input data DOP is able to animate visual map and save the animation as video file. Indeed, to reproduce the animation corresponding to the selected model, the user should download in DOP the XML file with the constructed visual map and corresponding PLT file generated in DVM (see Fig. 2). Before starting with playback of the animation user should set the values of parameters of DOP (see Params editor box on Fig. 2) responsible for the visual properties of the animation.

DBSolve Optimum allows visualization of simulation results in two different ways: "Bar animation" and "Arrow animation" (see corresponding options for construction of visual map in DVM interface Fig. 3). "Bar animation" mode (Fig. 5) allows modeler to create a set of the bars corresponding to reactions and entities of the system. The height of the rectangles corresponds to the value of either reaction rate or concentration. Animation in this mode implemented in DOP as changes in heights of the bars in accordance to data provided by simulations. "Arrow animation" mode (Fig. 6) allows us to create static map consisting of nodes (biological entities) and the directed edges (reactions). Variation of concentration of i-th metabolite, *C*_{
i
} , is represented by change in circle radius, *r*_{
i
}, around the corresponding node in such a way that the concentration is directly proportional to circle radius:

{r}_{i}=\{\begin{array}{l}0,{C}_{i}<{\overline{C}}_{\mathrm{min}}\\ \frac{{C}_{i}-{\overline{C}}_{\mathrm{min}}}{{\overline{C}}_{\mathrm{max}}-{\overline{C}}_{\mathrm{min}}}\cdot {r}_{\mathrm{max}},\\ {\overline{C}}_{\mathrm{min}}\le {C}_{i}\le {\overline{C}}_{\mathrm{max}}\\ {r}_{\mathrm{max}},{C}_{i}>{\overline{C}}_{\mathrm{max}}\end{array}

where {\overline{C}}_{\mathrm{min}}, {\overline{C}}_{\mathrm{max}} are values of minimal and maximal concentration and *r*_{max} is maximal radius of the circles specified by the user. There are two options to visualize changes in reaction rates. First one visualizes the value of the reaction rate, *V*_{
i
}, as the thickness of the arrow, *h*_{
i
}, connecting substrates and products of the reaction (Fig. 6), thus the thickness of an arrow is directly proportional to the value of the reaction rate:

{h}_{i}=\{\begin{array}{l}0,{V}_{i}<{\overline{V}}_{\mathrm{min}}\\ \frac{{V}_{i}-{\overline{V}}_{\mathrm{min}}}{{\overline{V}}_{\mathrm{max}}-{\overline{V}}_{\mathrm{min}}}\cdot {h}_{\mathrm{max}},\\ {\overline{V}}_{\mathrm{min}}\le {V}_{i}\le {\overline{V}}_{\mathrm{max}}\\ {h}_{\mathrm{max}},{V}_{i}>{\overline{V}}_{\mathrm{max}}\end{array}

Where {\overline{V}}_{\mathrm{min}}, {\overline{V}}_{\mathrm{max}} are values of minimal and maximal fluxes, *h*_{max} is maximal thickness of the arrows specified by the user. The second one represents the value of the reaction rate as a circle, thus the circle radius is directly proportional to the reaction rate.

Further, clicking button "Play" allows viewing animation and clicking "Save to avi" saves the animation in AVI format.

In additional file 1: supplementary materials we present in details all functionalities of DBsolve Optimum and exemplify them in kinetic model of mitochondrial Krebs cycle.