The systems biology simulation core algorithm
 Roland Keller†^{1},
 Alexander Dörr†^{1},
 Akito Tabira^{2},
 Akira Funahashi^{2},
 Michael J Ziller^{3},
 Richard Adams^{4},
 Nicolas Rodriguez^{5},
 Nicolas Le Novère^{6},
 Noriko Hiroi^{2},
 Hannes Planatscher^{1, 7},
 Andreas Zell^{1} and
 Andreas Dräger^{1, 8}Email author
DOI: 10.1186/17520509755
© Keller et al.; licensee BioMed Central Ltd. 2013
Received: 14 December 2012
Accepted: 18 June 2013
Published: 5 July 2013
Abstract
Background
With the increasing availability of high dimensional time course data for metabolites, genes, and fluxes, the mathematical description of dynamical systems has become an essential aspect of research in systems biology. Models are often encoded in formats such as SBML, whose structure is very complex and difficult to evaluate due to many special cases.
Results
This article describes an efficient algorithm to solve SBML models that are interpreted in terms of ordinary differential equations. We begin our consideration with a formal representation of the mathematical form of the models and explain all parts of the algorithm in detail, including several preprocessing steps. We provide a flexible reference implementation as part of the Systems Biology Simulation Core Library, a communitydriven project providing a large collection of numerical solvers and a sophisticated interface hierarchy for the definition of custom differential equation systems. To demonstrate the capabilities of the new algorithm, it has been tested with the entire SBML Test Suite and all models of BioModels Database.
Conclusions
The formal description of the mathematics behind the SBML format facilitates the implementation of the algorithm within specifically tailored programs. The reference implementation can be used as a simulation backend for Java™based programs. Source code, binaries, and documentation can be freely obtained under the terms of the LGPL version 3 from http://simulationcore.sourceforge.net. Feature requests, bug reports, contributions, or any further discussion can be directed to the mailing list simulationcoredevelopment@lists.sourceforge.net.
Keywords
Systems biology Biological networks Mathematical modeling Simulation Algorithms Ordinary differential equation systems Numerical integration Software engineeringBackground
As part of the movement towards quantitative biology, the modeling, simulation, and computer analysis of biological networks have become integral parts of modern biological research[1]. Ambitious national and international research projects such as the Virtual Liver Network[2] strive to derive even organwide models of biological systems that include all kinds of processes taking place at several levels of detail. Largescale efforts like this require intensive collaboration between various research groups, including experimenters, modelers, and bioinformaticians. The exchange, storage, interoperability, and the possibility to combine models have been recognized as key aspects of this endeavor[3–6].
XMLbased standard description formats such as the Systems Biology Markup Language (SBML)[7, 8] and CellML[9, 10] enable encoding of quantitative biological network models. To facilitate sharing and reuse of the models, online databases such as BioModels Database[11] and the CellML model repository[12] provide large collections of published models. Software libraries for reading and manipulating the content of these formats are also available[13–15] as well as enduser programs supporting these model description languages.
The models encoded in these formats can be interpreted in terms of several modeling frameworks, including, but not limited to, differential equation systems, with additional structures such as discrete events and algebraic equations. The diversity of modeling approaches and experimental data often requires customized software solutions for very specific tasks. For efficient analysis, simulation, and calibration (e.g., the estimation of parameter values) of biological network models a multiplepurpose and efficient numerical solver library is prerequisite. Although the language specifications of SBML[16–22] and CellML[23] describe the semantics of models in these formats and their interpretation, the algorithmic implementation is still not straightforward.
The SBML community offers standardized and manually derived benchmark tests[24] in order to evaluate the quality of simulation results, because it has been recognized that in many cases different solver implementations lead to divergent results[25]. The availability of this test suite and the currently much larger variety of supporting software for SBML^{a} in comparison to CellML are the reasons that in this work we focus on the simulation of models encoded in the SBML format.
We address the question of how to precisely interpret these models in terms of ordinary differential equation systems. Furthermore, we show how to adapt existing numerical integration routines in order to simulate these models. To this end, we derive a new algorithm for the accurate interpretation and simulation of all currently existing levels and versions of SBML. To demonstrate the usefulness of the algorithm, we introduce an exhaustive reference implementation in Java™. The algorithm described in this paper is, however, not limited to any particular programming language.
It is also important to note that the interpretation of these models must be strictly separated from the numerical method that solves the implied differential equation system. In this way, a similar approach would also be possible for other systems biology community formats. In particular, the architecture of the reference implementation described herein has been ab ovo designed with the aim to be complemented by a CellML module.
As the result, we present the Systems Biology Simulation Core Library, a platformindependent, welltested generic opensource library. The library is completely decoupled from any graphical user interface and can therefore easily be integrated into thirdparty programs. It comprises several ordinary differential equation (ODE) solvers and an interpreter for SBML models. It is the first simulation library based on JSBML[15].
Furthermore, the Systems Biology Simulation Core Library contains classes to both export simulation configurations to the Simulation Experiment Description Markup Language (SEDML)[26], and facilitate the reuse and reproduction of these experiments by executing SEDML files.
Results and discussion
In order to derive an algorithm for the interpretation of SBML models in a differential equation framework, it is first necessary to take a closer look at the mathematical equations implied by this data format. Based on this general description, we will then discuss all necessary steps to deduce an algorithm that takes all special cases for the various levels and versions of SBML into account.
A formal representation of models in systems biology
with t_{0}≡0 and${t}_{T}\in {\mathbb{R}}_{+}$. The vector$\overrightarrow{Q}$ of quantities contains the sizes of the compartments$\overrightarrow{C}$, amounts (or concentrations) of reacting species$\overrightarrow{S}$, and the values of all global model parameters$\overrightarrow{P}$. It should be noted that these models may contain local parameters$\overrightarrow{P}$ that influence the reactions’ velocities, but which are not part of the global parameter vector$\overrightarrow{P}$, and hence also not part of$\overrightarrow{Q}$.
All vector function terms may involve a delay function, i.e., an expression of the form delay(x,τ) with τ>0. It is therefore possible to address values of x computed in the earlier integration step at time t−τ, turning equation (3) into a delay differential equation (DDE). Note that x can be an arbitrarily complex expression.
In the general case of equation (3), not all species’ amounts can be computed by integrating the transformation$\mathbf{N}\overrightarrow{\nu}$: the change of some model quantities may be given in the form of rate rules by function$\overrightarrow{g}(\overrightarrow{Q},t)$. Species whose amounts are determined by rate rules must not participate in any reaction and hence only have zerovalued corresponding entries in the stoichiometric matrix N. Thereby, the rate rule function$\overrightarrow{g}(\overrightarrow{Q},t)$ directly gives the rate of change of these quantities, and returns 0 for all others.
In addition, SBML introduces the concept of events${\overrightarrow{f}}_{E}(\overrightarrow{Q},t)$ and assignment rules$\overrightarrow{r}(\overrightarrow{Q},t)$. An event can directly manipulate the value of several quantities, for instance, reduce the size of a compartment to a certain portion of its current size, as soon as a trigger condition becomes satisfied. An assignment rule also influences the absolute value of a subject quantity.
A further concept in SBML is that of algebraic rules, which are equations that must evaluate to zero at all times during the simulation of the model. These rules can be solved to determine the values of quantities whose values are not determined by any other construct. In this way, conservation relations or other complex interrelations can be expressed in a very convenient way. With the help of bipartite matching[29] and a subsequent conversion it is possible to turn algebraic rules into assignment rules and hence include these into the term$\overrightarrow{r}(\overrightarrow{Q},t)$. Such a transformation, however, requires symbolic computation and is thus a complicated endeavor.
When the system under study operates at multiple time scales, i.e., it contains a fast and a slow subsystem, a separation of the system is necessary, leading to differential algebraic equations (DAEs). Some species can be declared to operate at the system’s boundaries, assuming a constant pool of their amounts or concentrations. Care must also be taken with respect to the units of the species, because under certain conditions division or multiplication with the sizes of their surrounding compartments becomes necessary in order to ensure the consistent interpretation of the models. For all these reasons, solving equation (3) is much more complicated than computing the solution of the simple equation (2) alone.
From the perspective of software engineering, a strict separation of the interpretation of the model and the numerical treatment of the differential equation system is necessary to ensure that regular numerical methods can be used to solve equation (3). In order to efficiently compute this solution, multiple preprocessing steps are required, such as the conversion of algebraic rules into assignment rules, or avoiding repeated recomputation of intermediate results. The next sections will give a detailed explanation of the necessary steps to solve these systems and how to efficiently perform their numerical integration with standard numerical solvers.
Initialization
After the creation of this graph, the initial assignments and the assignment rules (including transformed algebraic rules) are processed and initial values defined by these constructs are computed.
Solving algebraic rules
Event handling
An event in SBML is a list of assignments that is executed depending on whether a trigger condition switches from false to true. In addition, SBML enables modellers to define a delay which may postpone the actual execution of the event’s assignments to a later point in time. With the release of SBML Level 3 Version 1, the processing of events has been raised to an even higher level of complexity: in earlier versions it was sufficient to determine, when an event triggers and when its assignments are to be executed. In Level 3 Version 1 only a few new language elements have been added, but these have a significant impact on how to handle events: for example, the order, in which events have been processed, used to be at programmer’s discretion in SBML Level 2, but in Level 3 Version 1 it is given by the event’s priority element. Coordinating the sequence, in which events are to be executed, has now become the crucial part of event handling. Furthermore, there exists the option to cancel an event during the time since its trigger has been activated and the actual time when the scheduler picks the event for execution. Events that can be cancelled after the activation of their triggers are called nonpersistent.
The interpretation of events is the most time consuming step of the integration procedure. This is why efficient and clearly organized data structures are required to ensure high performance of the algorithm.
Time step adaptation considering events and the calculation of derivatives
After that, the events and the assignment rules are processed at the new point in time t_{τ−1}+h. If the previous step causes a change in$\overrightarrow{Q}$, the adaptive step size is decreased by setting h to h/10 and the calculation is repeated until either the minimum step size is reached or the processing of events and assignment rules does not change$\overrightarrow{Q}$ anymore. Hence, the time at which an event takes place is precisely determined.
A reference implementation of the algorithm
Due to the strict separation between numerical differential equation solvers, and the definition of the actual differential equation system, it is possible to implement support for other community standards, such as CellML[9].
In order to support the standard Minimum Information About a Simulation Experiment (MIASE)[33], the library also provides an interpreter of Simulation Experiment Description Markup Language (SEDML) files[26]. These files allow users to store the details of a simulation, including the selection and all settings of the numerical method, hence facilitating the creation of reproducible results. A simulation experiment can also be directly started by passing a SEDML file to the interpreter in this library. Each solver has a method to directly access its corresponding Kinetic Simulation Algorithm Ontology (KiSAO) term[34] to facilitate the execution of SEDML files.
Many interfaces, abstract classes, and an exhaustive source code documentation in the form of JavaDoc facilitate the customization of the library. For testing purposes, the library contains a sample program that benchmarks its SBML interpreter against the entire SBML Test Suite version 2.3.2[24].
Benchmark and application to published models
Simulation of the models from the SBML Test Suite using the Rosenbrock solver
Level  Version  Models  Correct simulations  Total running time (in s) 

1  2  252  252  2.9 
2  1  885  885  6.8 
2  2  1,041  1,041  6.8 
2  3  1,041  1,041  6.8 
2  4  1,043  1,043  6.8 
3  1  1,077  1,077  38.5 
The total simulation time for all models in SBML Level 3 Version 1 is significantly higher than for the models in other SBML levels and versions. This can be explained by the fact that the test suite contains some models of this version whose evaluation requires a timeconsuming processing of a large number of events. In particular, the simulation of model No. 966 of the SBML Test Suite, which is only provided in SBML Level 3 Version 1, takes 20s because it contains 23 events to be processed. Two events fire every 10^{−2} time units within the simulation time period of 1,000 time units. These events must therefore be evaluated thousandfold within the specified time interval. The evaluation of this model accounts for over 50% of the total simulation time for the models in SBML Level 3 Version 1.
An implementation of an SBML solver that passes the test suite should in principle also be capable of computing the solution of all models from BioModels Database, a resource that contains a collection of published and curated models. This online database currently provides neither reference data for the models, nor any settings for the numerical computation (such as step size, end time etc.). However, it offers precomputed plots of the time courses for the vast majority of models. Therefore, while it cannot be directly used as a benchmark test, it can help checking that a solver implementation supports all features of many published models and that the algorithm always successfully terminates. The Systems Biology Simulation Core Library solves all curated models from BioModels Database (release 23, October 2012) without raising any errors, see Methods for details. These results suggest the reliability of the simulation algorithm described in this work.
In the following, we select two models that exhibit diverse features from this repository to illustrate the capabilities of this library: BioModels Database model No. 206 by Wolf et al.[35] and BioModels Database model No. 390 by Arnold and Nikoloski[36].
Comparison to existing solver implementations for SBML
In order to benchmark our software, we chose similar tools exhibiting the following features from the SBML software matrix[39]:

The last updated version was released after the final release of the specification for SBML Level 3 Version 1 Core, i.e., October 6^{th} 2010.

Support for SBML Level 3.

Opensource software

No dependency on commercial products that are not freely available (e.g., MATLAB™ or Mathematica™)
Comparison of SBMLcapable simulators
Program  Version  Difficult SBML  Fully SBML test  SEDML  Programming  GUI  API  Platform  Comments  

elements  Suite compliant  language  access  
Fast  Algebraic  Events  
reactions  rules  
BioUML  0.9.4  ✓  ✓  ✓  ✓  In α  Java  ✓  JavaScript  Independent  
version  
COPASI  4.9.45  –  –  (✓)  –  –  C++ (with  ✓  ✓  Windows, Mac OS X,  
multiple bindings)  Linux, Solaris  
iBioSim  2.4.5  ✓  ✓  ✓  ✓  In α  Java, C  ✓  (✓)  Windows, Mac OS X,  
version  Linux (Fedora 17)  
JSim  2.10  –  ✓  –  –  –  Java  ✓  ✓  Windows, Mac OS X,  
Linux  
LibSBMLSim  1.1.0  ✓  ✓  ✓  (✓)  –  C (with multiple)  –  ✓  Windows, Mac OS X,  
bindings)  Linux, Free BSD  
Simulation  1.3  ✓  ✓  ✓  ✓  ✓  Java  –  ✓  Independent  
core library  
VCell  5.0  ✓  –  ✓  –  –  Java frontend,  ✓  –  Independent  Internet 
C/C++ server  connection  
backend  required 
Limitations and perspective
The modifications done to the Rosenbrock solver enable a precise timing of events during simulation. However, this precise timing can lead to a noticeable increase in runtime when events are triggered in very small intervals, e.g., every 10^{−3} time units. This behavior can, for example, be observed in BioModels Database model No. 408[48] (a model with three events). When the precise timing of events is not of utmost importance, a solver other than Rosenbrock can be chosen. Furthermore, there are plans to improve the runtime behavior of the Rosenbrock solver for the simulation of models containing events.
When dealing with stiff problems, Rosenbrock’s method is a good choice, because it is has been designed for stiff pODE. However, our experiments show, that the Rosenbrock solver can be inefficient for nonstiff problems in comparison to other solvers. This issue can lead to an increased runtime regarding large models such as model No. 235 of the BioModels Database, which contains 622 species that participate in 778 reactions, distributed accross three compartments[49]. In some cases, tuning the relative and absolute tolerance can help, but depending on the system’s structure, Rosenbrock’s method is sometimes stretched to its limits. The RungeKuttaFehlberg method[50] (KiSAO term 86), which is included in iBioSim, shows also an increase in runtime concerning this model.
The performance of the RungeKuttaFehlberg and Rosenbrock methods show, however, that simpler ODE solvers can have more difficulties with some biological models than more advanced solvers, such as CVODE from SUNDIALS[51] that can adapt to both nonstiff and stiff problems. The SUNDIALS library, which is incorporated into BioUML, can handle complicated pODE significantly better, but since it is not available under the LGPL and no opensource Java version of these solvers can currently be obtained, we disregarded its use.
Algebraic rules constitute an important problem for any implementation of the SBML standard. The unbound variable of each such equation can be efficiently identified[29], whereas the transformation of an algebraic rule into an assignment rule includes symbolic computation and is very difficult to implement. In some cases, such a transformation is not even possible. Alternatively, the current value of the free variable in an algebraic equation could, for instance, be identified using nested intervals. However, this approach consumes a significantly higher runtime, because the nested intervals would have to be recomputed at every time step, whereas the transformation approach considers every algebraic rule only once (during the initialization).
Since Level 3, SBML entails one further aspect: it is now possible to add additional features to the model by declaring specialized extension packages. The algorithm discussed in this paper describes the core functionality of SBML. The extension packages are very diverse, reaching the graphical representation[53], the description of qualitative networks, such as Petri nets[54], and many more. It is therefore necessary to separately derive and implement algorithms for the interpretation of individual SBML packages.
The agenda for the further development of the opensource project, the Systems Biology Simulation Core Library, includes the implementation of SBML extension packages, support for CellML, and the incorporation of additional numerical solvers. Contributions from the community are welcome.
Conclusions
The aim of this work is to derive a formal description of the mathematics behind SBML together with an algorithm that efficiently solves it in terms of an ordinary differential equation framework. As an important design feature, the algorithm can be combined with existing numerical solvers in a plugin fashion. The Rosenbrock solver embodies a universal approach for simulation that can deal with stiff problems and precisely solve models containing arbitrary SBML elements. The description in this paper is intended to facilitate the implementation of the algorithm within specifically tailored programs.
Our tests indicate that at the moment only two other programs pass the entire test suite for all SBML levels and versions: BioUML, which is a workbench for modelling, simulation, and parameter fitting, and iBioSim. The reference implementation of the algorithm introduced in this work, the Systems Biology Simulation Core Library, is therefore the only API simulation library exhibiting this capability.
The Systems Biology Simulation Core Library is an efficient Java tool for the simulation of differential equation systems used in systems biology. It can be easily integrated into larger customized applications. For instance, CellDesigner[55] has already been using it since version 4.2 as one of its integral simulation libraries. The standalone application SBMLsimulator[56] provides a convenient graphical user interface for the simulation of SBML models and uses it as a computational backend. The abstract class structure of the library supports the integration of further model formats, such as CellML, in addition to its SBML implementation. To this end, it is only necessary to implement a suitable interpreter class.
By including support for the emerging standard SEDML, we hope to facilitate the exchange, archival and reproduction of simulation experiments performed using the Systems Biology Simulation Core Library.
Methods
Implementation
All the solver classes are derived from the abstract class AbstractDESSolver (Figure 6). Several solvers of the Apache Commons Math library (version 3.0) are integrated with the help of wrapper classes[32]. Numerical methods and the actual differential equation systems are strictly separated. The class MultiTable stores the results of a simulation within its Block data structures.
The abstract description of differential equation systems, with the help of several distinct interfaces, makes it possible to decouple them from a particular type of biological network. It is therefore possible to pass an instance of an interpreter for a respective model description format to any available solver. The interpretation of SBML models is split between evaluation of events and rules, computation of stoichiometric information, and computation of the current values for all model components (such as species and compartments).
For a given state of the ODE system, the class SBMLinterpreter, responsible for the evaluation of models encoded in SBML, returns the current set of timederivatives of the variables. It is connected to an efficient MathML interpreter of the expressions contained in kinetic laws, rules and events (ASTNodeInterpreter). The nodes of the syntax graph for those expressions depend on the current state of the ODE system. If the state has changed, the values of the nodes have to be recalculated (see Results).
An important aspect in the interpretation of SBML models is the determination of the exact time at which an event occurs because this influences the precision of the system’s variables. To this end, we adjusted an implementation of the Rosenbrock solver[57], an integrator with an adaptive step size, to a very precise timing of the events. In addition to events, rules are also treated during integration. Basically, rules are treated like events that occur at every given point in time and are therefore processed in the same manner. For every object of the type AlgebraicRule, a new AssignmentRule object is generated by means of the preceding bipartite matching. They represent only temporary rules, that are incorporated in the simulation process but do not influence the model in the SBML file.
In the SBMLinterpreter, events are represented via an array containing one instance of EventInProgress for every event in the model. Thereby, the distinction between events with and without delays is made. Both types of events can be triggered multiple times before being executed. If no delay is defined, the assignments of the event are usually executed at the same point in time when the event has been triggered. However, when such an event is cancelled by other events, all of its assignments are also cancelled before execution. An event with delay can produce multiple further assignments within the time frame between the trigger time and the actual execution time. In order to deal with delayed events, the class SBMLEventInProgressWithDelay keeps track of this via a list containing the points in time, at which the respective event has to be executed. When events are triggered more than once before execution, they have to be sorted in ascending order by their delay. This is neccesary, because in this case the delay of the very same event may vary.
When the SBMLinterpreter is processing events with priority, the events with the highest priority are stored in a list until one of them is selected for execution. Technically, the method of choice for the organization of such priority queues would apply a binary max heap data structure instead. The root of the heap represents the largest value in the heap. After its extraction, the heap property is restored so that the next largest value is moved to the root. However, as stated in Results, the execution of one event can influence the priority of the remaining events. It can possibly happen that many priorities simultaneously change, whereby the standard method to restore the max heap characteristic after extraction is not sufficient anymore. For this reason, we disregarded the use of more complex data structures for the current implementation.
Since SBML Level 2 Version 1, it has also become possible to create userdefined functions. These function definition objects contain lambda calculus including an optional list of arguments together with the actual mathematical expression of the function. During the initialization phase, function definitions are also incorporated into the abstract syntax graph (Figure 1). For each function definition, its arguments defined in its lambda expression are mapped to their corresponding nodes in the abstract syntax graph. The evaluation of a syntax graph node with a userdefined function consists of several steps. The arguments are evaluated and then passed to their corresponding node in the graph via the mapping established before. After this step the nodes representing arguments have a specific value attached to them. Finally, the complete abstract syntax graph can be evaluated. Care must be taken, because several function definitions may have arguments with identical identifiers. All possible naming conflicts must be preempted.
As part of the calculation of reaction velocities, the StoichiometryMath construct allows a dynamic change of a reaction’s stoichiometry over the course of the simulation. Since SBML Level 3 Version 1, the stoichiometry of a reaction can be directly altered, because it is now possible to address the identifier of a SpeciesReference as the target variable within rules or events. The SBMLinterpreter class flags reactions with changing stoichiometry during initialization and evaluates the corresponding abstract syntax graph anew if the stoichiometry is needed for calculation.
The constraints introduce assumptions about a model’s behavior. Similar to the trigger of events, the abstract syntax graph of each constraint is evaluated at every time step. In case of a violation the SBMLinterpreter generates an instance of ConstraintEvent that is then processed by the corresponding ConstraintListener class. The user is informed about the constraint upon its violation via the standard Java Logger. The output message includes the point in simulation time and the message of the constraint. In addition, more advanced userdefined implementations of ConstraintListener can be added to the SBMLinterpreter, for instance, to notify a GUI about violations or display the associate messages in a more userfriendly way.
SEDML support is enabled by inclusion of the jlibsedml library[58] in the binary download. Clients of the Systems Biology Simulation Core Library can choose to use the jlibsedml API directly, or access SEDML support via facade classes in the org.simulator.sedml package that do not require direct dependencies on jlibsedml in their code.
Default settings and configuration
The standard preferences for simulating an SBML model consist of the Rosenbrock solver with an absolute tolerance of 10^{−12} and a relative tolerance of 10^{−6}. On the basis of our experiments, this setup can handle most of the problems without further tuning. The Rosenbrock solver with its adaptive step size is the most effective solver in this library for stiff pODE. Nevertheless, the user has the possibility the choose another solver for integration. According to the SBML specifications, a model has to be simulated starting at time point 0.0. Since this library is not limited to SBML, the solvers also accept arbitrary start times. The user has also the possibility to specify the end of the simulation. Modifying the relative and absolute tolerance can increase the accuracy of the results or decrease computation time.
Simulation of models from BioModels Database
All 424 curated models from BioModels Database[11] (release 23, October 2012) have been simulated with identical settings, as suggested by Bergmann et al.[25]: time interval [0,10], the Rosenbrock solver, 10^{−6} as relative and 10^{−12} as absolute tolerance, and a step size of 0.01 time units. For the models No. 234[59] and No. 339[60] from BioModels Database the absolute tolerance had to be set to 10^{−10} in order to achieve the necessary accuracy and to avoid that the algorithm surpasses its minimal step size. On a sample basis, individual models have been selected and manually compared to the precomputed plots provided by BioModels Database in order to check the correctness of the simulation results.
Simulation of the SBML test suite
The models from SBML Test Suite version 2.3.2[24] were first simulated with the Rosenbrock solver, 10^{−6} as relative and 10^{−12} as absolute tolerance. For six models (No. 863, 882, 893, 994, 1109, and 1121) we had to set the relative tolerance to 10^{−8} in order to simulate as accurately as desired. For three other models (No. 872, 987, 1052) the relative tolerance even had to be set to 10^{−12} and the absolute tolerance to 10^{−14}.
Hard and software configuration
For all runtime tests, an Intel ^{ Ⓡ } Core™ i5 CPU with 3.33GHz and 4GB RAM was used with Microsoft ^{ Ⓡ } Windows ^{ Ⓡ } 7 (Version 6.1.7600) as operating system and Java Virtual Machine version 1.6.0_25.
The Systems Biology Simulation Core Library was also successfully tested under Linux (Ubuntu version 10.4) and Mac OS X (versions 10.6.8 and 10.8.2).
Availability and requirements
The current version of Systems Biology Simulation Core Library is available at the project’s homepage. The entire project, including source code and documentation, several versions of jar files containing only binaries, binaries together with source code, can be downloaded, optionally also as a version including all required thirdparty libraries.
Project name: Systems Biology Simulation Core Library
Project homepage: http://simulationcore.sourceforge.net
Operating systems: Platform independent, i.e., for all systems for which a JVM is available.
Programming language: Java™
Other requirements: Java Runtime Environment (JRE) 1.6 or above
License: GNU Lesser General Public License (LGPL) version 3
Endnote
^{a} More than 250 available programs now support the SBML data format (April 19^{th} 2013).
Notes
Abbreviations
 ADP:

Adenosine diphosphate
 API:

Application programing interface
 ATP:

Adenosine5’triphosphate
 DHAP:

Dihydroxyacetone phosphate
 DAE:

Differential algebraic equation
 DDE:

Delay differential equation
 F1:

6BP: Fructose 1,6bisphosphate
 GA3P:

Glyceraldehyde 3phosphate
 GUI:

Graphical user interface
 JAR:

Java archive file
 JDK:

Java development kit
 JRE:

Java runtime environment
 JVM:

Java virtual machine
 KiSAO:

Kinetic simulation algorithm Ontology
 MIASE:

Minimum information about a simulation experiment
 LGPL:

GNU lesser general public license
 ODE:

Ordinary differential equation
 RuBisCO:

Ribulose1,5bisphosphate carboxylase oxygenase
 NAD+:

Nicotinamide adenine dinucleotide
 SBML:

Systems biology markup language.
Declarations
Acknowledgements
The authors are grateful to B. Kotcon, S. Mesuro, D. Rozenfeld, A. Yodpinyanee, A. Perez, E. Doi, R. Mehlinger, S. Ehrlich, M. Hunt, G. Tucker, P. Scherpelz, A. Becker, E. Harley, and C. Moore, Harvey Mudd College, USA, for providing a Java implementation of Rosenbrock’s method, and to M. T. Cooling, University of Auckland, New Zealand, for fruitful discussion. The authors thank D. M. Wouamba, P. Stevens, M. Zwießele, M. Kronfeld, and A. Schröder for source code contribution and fruitful discussion.
This work was funded by the Federal Ministry of Education and Research (BMBF, Germany) as part of the Virtual Liver Network (grant number 0315756). The Japan Society for the Promotion of Science and the Ministry of Education, Culture, Sports, Science and Technology of Japan supported this work by GrantsinAid for Scientific Research on Innovative Areas (KAKENHI), grant number 23136513. In addition, this work was funded by the UK Biotechnology and Biological Sciences Research Council, grant number BB/D019621/1. We acknowledge support by the German Research Foundation (DFG) and the Open Access Publishing Fund of the University of Tuebingen.
Authors’ Affiliations
References
 Macilwain C: Systems biology: evolving into the mainstream. Cell. 2011, 144 (6): 839841. 10.1016/j.cell.2011.02.044. [http://dx.doi.org/10.1016/j.cell.2011.02.044] 10.1016/j.cell.2011.02.044PubMedView Article
 Holzhutter HG, Drasdo D, Preusser T, Lippert J, Henney AM: The virtual liver: a multidisciplinary, multilevel challenge for systems biology. Wiley Interdiscip Rev Syst Biol Med. 2012, 4 (3): 221235. 10.1002/wsbm.1158. [http://dx.doi.org/10.1002/wsbm.1158] 10.1002/wsbm.1158PubMedView Article
 Schulz M, Uhlendorf J, Klipp E, Liebermeister W: SBMLmerge, a system for combining biochemical network models. Genome Inform Ser. 2006, 17: 6271.
 Klipp E, Liebermeister W, Helbig A, Kowald A, Schaber J: Systems biology standards—the community speaks. Nat Biotechnol. 2007, 25 (4): 390391. 10.1038/nbt0407390. [http://dx.doi.org/10.1038/nbt0407390] 10.1038/nbt0407390PubMedView Article
 Chelliah V, Endler L, Juty N, Laibe C, Li C, Rodriguez N, Le Novere N: Data integration and semantic enrichment of systems biology models and simulations. Data Integration in the Life Sciences Volume 5647 of Lecture Notes in Computer Science. Edited by: Paton NW, Missier P, Hedeler C. 2009, Berlin, Heidelberg: Springer, 515. [http://dx.doi.org/10.1007/9783642028793_2]
 Liebermeister W, Krause F, Uhlendorf J, Lubitz T, Klipp E: SemanticSBML: a tool for annotating, checking, and merging of biochemical models in SBML format. 3rd International Biocuration Conference. 2009, Nature Publishing Group, [http://dx.doi.org/10.1038/npre.2009.3093.1]
 Hucka M, Finney A, Sauro HM, Bolouri H, Doyle JC, Kitano H, Arkin AP, Bornstein BJ, Bray D, CornishBowden A, Cuellar AA, Dronov S, Gilles ED, Ginkel M, Gor V, Goryanin II, Hedley WJ, Hodgman TC, Hofmeyr JHS, Hunter PJ, Juty NS, Kasberger JL, Kremling A, Kummer U, Le NovereN, Loew LM, Lucio D, Mendes P, Minch E, Mjolsness ED, Nakayama Y, Nelson MR, Nielsen PF, Wagner JM, Forum S, Sakurada T: The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics. 2003, 19 (4): 524531. 10.1093/bioinformatics/btg015. [http://bioinformatics.oxfordjournals.org/cgi/content/abstract/19/4/524] 10.1093/bioinformatics/btg015PubMedView Article
 The Systems Biology Markup Language. [http://sbml.org]
 Lloyd CM, Halstead MDB, Nielsen PF: CellML: its future, present and past. Prog Biophys Mol Bio. 2004, 85 (2–3): 433450. [http://dx.doi.org/10.1016/j.pbiomolbio.2004.01.004]View Article
 The CellML project. [http://cellml.org]
 Li C, Donizelli M, Rodriguez N, Dharuri H, Endler L, Chelliah V, Li L, He E, Henry A, Stefan MI, Snoep JL, Hucka M, Le Novere N, Laibe C: BioModels Database: An enhanced, curated and annotated resource for published quantitative kinetic models. BMC Syst Biol. 2010, 4: 9210.1186/17520509492. [http://dx.doi.org/10.1186/17520509492] 10.1186/17520509492PubMedPubMed CentralView Article
 Lloyd CM, Lawson JR, Hunter PJ, Nielsen PF: The CellML Model Repository. Bioinformatics. 2008, 24 (18): 21222123. 10.1093/bioinformatics/btn390. [http://dx.doi.org/10.1093/bioinformatics/btn390] 10.1093/bioinformatics/btn390PubMedView Article
 Bornstein BJ, Keating SM, Jouraku A, Hucka M: LibSBML: an API library for SBML. Bioinformatics. 2008, 24 (6): 880881. 10.1093/bioinformatics/btn051. [http://dx.doi.org/10.1093/bioinformatics/btn051] 10.1093/bioinformatics/btn051PubMedPubMed CentralView Article
 Miller AK, Marsh J, Reeve A, Garny A, Britten R, Halstead M, Cooper J, Nickerson DP, Nielsen PF: An overview of the CellML API and its implementation. BMC Bioinformatics. 2010, 11: 17810.1186/1471210511178. [http://dx.doi.org/10.1186/1471210511178] 10.1186/1471210511178PubMedPubMed CentralView Article
 Drager A, Rodriguez N, Dumousseau M, Dorr A, Wrzodek C, Le Novere N, Zell A, Hucka M: JSBML: a flexible Java library for working with SBML. Bioinformatics. 2011, 27 (15): 21672168. 10.1093/bioinformatics/btr361. [http://bioinformatics.oxfordjournals.org/content/27/15/2167] 10.1093/bioinformatics/btr361PubMedPubMed CentralView Article
 Hucka M, Finney A, Sauro H, Bolouri H: Systems Biology Markup Language (SBML) level 1: structures and facilities for basic model definitions. Tech. rep., Systems Biology Workbench Development Group JST ERATO Kitano Symbiotic Systems Project Control and Dynamical Systems. 2001, CA USA: MC 10781, California Institute of Technology, Pasadena
 Hucka M, Finney A, Sauro H, Bolouri H: Systems Biology Markup Language (SBML) level 1: structures and facilities for basic model definitions. Tech. Rep. 2, Systems Biology Workbench Development Group JST ERATO Kitano Symbiotic Systems Project Control and Dynamical Systems, MC 10781. 2003, CA USA: California Institute of Technology, Pasadena
 Finney A, Hucka M: Systems Biology Markup Language (SBML) level 2: structures and facilities for model definitions. Tech. rep., Systems Biology Workbench Development Group JST ERATO Kitano Symbiotic Systems Project Control and Dynamical Systems, MC 10781. 2003, California Institute of Technology
 Finney A, Hucka M, Le Novere N: Systems Biology Markup Language (SBML) level 2: structures and facilities for model definitions. Tech. rep. 2006
 Hucka M, Finney AM, Hoops S, Keating SM, Le Novere N: Systems Biology Markup Language (SBML) level 2: structures and facilities for model definitions. Tech. rep. 2007
 Hucka M, Finney A, Hoops S, Keating SM, Le Novere N: Systems biology markup language (SBML) Level 2: structures and facilities for model definitions. Tech. rep. Nature Precedings. 2008, [http://dx.doi.org/10.1038/npre.2008.2715.1]
 Hucka M, Bergmann FT, Hoops S, Keating SM, Sahle S, Schaff JC, Smith L, Wilkinson DJ: The Systems Biology Markup Language (SBML): language specification for level 3 version 1 core. Tech. rep. Nature Precedings. 2010, [http://precedings.nature.com/documents/4959/version/1]
 Cuellar A, Nielsen P, Halstead M, Bullivant D, Nickerson D, Hedley W, Nelson M, Lloyd C: CellML 1.1 Specification. Tech. rep., Bioengineering Institute. 2006, University of Auckland, [http://www.cellml.org/specifications/cellml_1.1/]
 SBML Test Suite. [http://sbml.org/Software/SBML_Test_Suite]
 Bergmann FT, Sauro HM: Comparing simulation results of SBML capable simulators. Bioinformatics. 2008, 24 (17): 19631965. 10.1093/bioinformatics/btn319. [http://bioinformatics.oxfordjournals.org/cgi/content/abstract/24/17/1963] 10.1093/bioinformatics/btn319PubMedPubMed CentralView Article
 Waltemath D, Adams R, Bergmann FT, Hucka M, Kolpakov F, Miller AK, Moraru II, Nickerson D, Sahle S, Snoep JL, Le Novere N: Reproducible computational biology experiments with SEDML–the Simulation Experiment Description Markup Language. BMC Syst Biol. 2011, 5: 19810.1186/175205095198. [http://dx.doi.org/10.1186/175205095198] 10.1186/175205095198PubMedPubMed CentralView Article
 Liebermeister W, Klipp E: Bringing metabolic networks to life: convenience rate law and thermodynamic constraints. Theor Biol Med Model. 2006, 3 (42): 41[http://dx.doi.org/10.1186/17424682341]PubMedPubMed CentralView Article
 Liebermeister W, Uhlendorf J, Klipp E: Modular rate laws for enzymatic reactions: thermodynamics, elasticities and implementation. Bioinformatics. 2010, 26 (12): 15281534. 10.1093/bioinformatics/btq141. [http://dx.doi.org/10.1093/bioinformatics/btq141] 10.1093/bioinformatics/btq141PubMedView Article
 Hopcroft JE, Karp RM: An n5/2 algorithm for maximum matchings in bipartite graphs. SIAM J Comput. 1973, 2 (4): 225231. 10.1137/0202019. [http://dx.doi.org/10.1137/0202019] 10.1137/0202019View Article
 Press WH, Teukolsky SA, Vetterling WT, Flannery BP: Numerical Recipes in FORTRAN; The Art of Scientific Computing. New York: Cambridge University Press. 1993
 Drager A, Hassis N, Supper J, Schroder A, Zell A: SBMLsqueezer: a CellDesigner plugin to generate kinetic rate equations for biochemical networks. BMC Syst Biol. 2008, 2: 3910.1186/17520509239. [http://www.biomedcentral.com/17520509/2/39] 10.1186/17520509239PubMedPubMed CentralView Article
 Commons Math: The Apache Commons Mathematics Library. [http://commons.apache.org/proper/commonsmath/]
 Waltemath D, Adams R, Beard DA, Bergmann FT, Bhalla US, Britten R, Chelliah V, Cooling MT, Cooper J, Crampin EJ, Garny A, Hoops S, Hucka M, Hunter P, Klipp E, Laibe C, Miller AK, Moraru I, Nickerson D, Nielsen P, Nikolski M, Sahle S, Sauro HM, Schmidt H, Snoep JL, Tolle D, Wolkenhauer O, Le Novere N: Minimum Information About a Simulation Experiment (MIASE). PLoS Comput Biol. 2011, 7 (4): e100112210.1371/journal.pcbi.1001122. [http://dx.doi.org/10.1371/journal.pcbi.1001122] 10.1371/journal.pcbi.1001122PubMedPubMed CentralView Article
 Courtot M, Juty N, Knupfer C, Waltemath D, Zhukova A, Drager A, Dumontier M, Finney A, Golebiewski M, Hastings J, Hoops S, Keating S, Kell DB, Kerrien S, Lawson J, Lister A, Lu J, Machne R, Mendes P, Pocock M, Rodriguez N, Villeger A, Wilkinson DJ, Wimalaratne S, Laibe C, Hucka M, Le Novere N: Controlled vocabularies and semantics in systems biology. Mol Syst Biol. 2011, 7: 543[http://dx.doi.org/10.1038/msb.2011.77]PubMedPubMed CentralView Article
 Wolf J, Passarge J, Somsen OJG, Snoep JL, Heinrich R, Westerhoff HV: Transduction of intracellular and intercellular dynamics in yeast glycolytic oscillations. Biophys J. 2000, 78 (3): 11451153. 10.1016/S00063495(00)766720. [http://dx.doi.org/10.1016/S00063495(00)766720] 10.1016/S00063495(00)766720PubMedPubMed CentralView Article
 Arnold A, Nikoloski Z: A quantitative comparison of CalvinBenson cycle models. Trends Plant Sci. 2011, 16 (12): 676683. 10.1016/j.tplants.2011.09.004. [http://dx.doi.org/10.1016/j.tplants.2011.09.004] 10.1016/j.tplants.2011.09.004PubMedView Article
 Le Novere N, Bornsteinm BJ, Broicher A, Courtot M, Donizelli M, Dharuri H, Li L, Sauro H, Schilstra M, Shapiro B, Snoep JL, Hucka M: BioModels Database: a free, centralized database of curated, published, quantitative kinetic models of biochemical and cellular systems. Nucleic Acids Res. 2006, 34: D689D691. 10.1093/nar/gkj092. [http://nar.oxfordjournals.org/cgi/content/full/34/suppl_1/D689] 10.1093/nar/gkj092PubMedPubMed CentralView Article
 Hairer E, Norsett SP, Wanner G: Solving Ordinary Differential Equations. 1 Nonstiff Problems. 2000, Berlin: Springer
 SBML Software Matrix (October 8th 2012). [http://sbml.org/SBML_Software_Guide/SBML_Software_Matrix]
 Kolpakov FA, Tolstykh NI, Valeev TF, Kiselev IN, Kutumova EO, Ryabova A, Yevshin IS, Kel AE: BioUML–open source plugin based platform for bioinformatics: invitation to collaboration. Moscow Conference on Computational Molecular Biology. 2011, Department of Bioengineering and Bioinformatics of MV Lomonosov Moscow State University, 172173. [http://mccmb.genebee.msu.ru/2011/mccmb11.pdf]
 Hoops S, Sahle S, Gauges R, Lee C, Pahle J, Simus N, Singhal M, Xu L, Mendes P, Kummer U: COPASI–a COmplex PAthway SImulator. Bioinformatics. 2006, 22 (24): 30673074. 10.1093/bioinformatics/btl485. [http://dx.doi.org/10.1093/bioinformatics/btl485] 10.1093/bioinformatics/btl485PubMedView Article
 Myers CJ, Barker N, Jones K, Kuwahara H, Madsen C, Nguyen NPD: iBioSim: a tool for the analysis and design of genetic circuits. Bioinformatics. 2009, 25 (21): 28482849. 10.1093/bioinformatics/btp457. [http://dx.doi.org/10.1093/bioinformatics/btp457] 10.1093/bioinformatics/btp457PubMedView Article
 Raymond GM, Butterworth E, Bassingthwaighte JB: JSIM: Free software package for teaching physiological modeling and research. Exp Biol. 2003, 280: 102107.
 Takizawa H, Nakamura K, Tabira A, Chikahara Y, Matsui T, Hiroi N, Funahashi A: LibSBMLSim: A reference implementation of fully functional SBML simulator. Bioinformatics. 2013, [http://dx.doi.org/10.1093/bioinformatics/btt157]
 Moraru II, Schaff JC, Slepchenko BM, Blinov ML, Morgan F, Lakshminarayana A, Gao F, Li Y, Loew LM: Virtual Cell modelling and simulation software environment. IET Syst Biol. 2008, 2 (5): 352362. 10.1049/ietsyb:20080102. [http://dx.doi.org/10.1049/ietsyb:20080102] 10.1049/ietsyb:20080102PubMedPubMed CentralView Article
 Resasco DC, Gao F, Morgan F, Novak IL, Schaff JC, Slepchenko BM: Virtual Cell: computational tools for modeling in cell biology. Wiley Interdiscip Rev: Syst Biol Med. 2012, 4 (2): 129140. 10.1002/wsbm.165. [http://dx.doi.org/10.1002/wsbm.165] 10.1002/wsbm.165
 SBML Test Suite Database—Test results for SBMLcompatible software systems. [http://sbml.org/Facilities/Database/Simulator]
 Hettling H, van Beek: JHGM: Analyzing the functional properties of the creatine kinase system with multiscale ‘sloppy’ modeling. PLoS Comput Biol. 2011, 7 (8): e100213010.1371/journal.pcbi.1002130. [http://dx.doi.org/10.1371/journal.pcbi.1002130] 10.1371/journal.pcbi.1002130PubMedPubMed CentralView Article
 Kuhn C, Wierling C, Kühn A, Klipp E, Panopoulou G, Lehrach H, Poustka AJ: Monte Carlo analysis of an ODE model of the Sea Urchin Endomesoderm network. BMC Syst Biol. 2009, 3 (3): 83[http://dx.doi.org/10.1186/17520509383]PubMedPubMed CentralView Article
 Fehlberg E: Klassische RungeKuttaFormeln vierter und niedrigerer Ordnung mit SchrittweitenKontrolle und ihre Anwendung auf Wärmeleitungsprobleme. Computing. 1970, 6 (1–2): 6171. [http://dx.doi.org/10.1007/BF02241732]View Article
 Hindmarsh AC, Brown PN, Grant KE, Lee SL, Serban R, Shumaker DE, Woodward CS: SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. ACM T Math Software. 2005, 31 (3): 363396. 10.1145/1089014.1089020. [https://computation.llnl.gov/casc/sundials/documentation/documentation.html] 10.1145/1089014.1089020View Article
 Madsen C, Myers CJ, Patterson T, Roehner N, Stevens JT, Winstead C: Design and test of genetic circuits using iBioSim. Design Test Comput, IEEE. 2012, 29 (3): 3239.View Article
 Gauges R, Rost U, Sahle S, Wegner K: A model diagram layout extension for SBML. Bioinformatics. 2006, 22 (15): 18791885. 10.1093/bioinformatics/btl195. [http://bioinformatics.oxfordjournals.org/cgi/content/abstract/22/15/1879] 10.1093/bioinformatics/btl195PubMedView Article
 Chaouiya C, Keating SM, Berenguier D, Naldi A, Thieffry D, van Iersel MP, Helicar T: Qualitative models. Tech. rep. 2013, [http://sbml.org/images/4/40/SBMLL3qualspecification20130415.pdf]
 Funahashi A, Tanimura N, Morohashi M, Kitano H: CellDesigner: a process diagram editor for generegulatory and biochemical networks. BioSilico. 2003, 1 (5): 159162. 10.1016/S14785382(03)023709. [http://www.sciencedirect.com/science/article/B75GS4BS08JD5/2/5531c80ca62a425f55d224b8a0d3f702] 10.1016/S14785382(03)023709View Article
 SBMLsimulator—An efficient Java solver implementation for SBML. [http://www.cogsys.cs.unituebingen.de/software/SBMLsimulator]
 Kotcon B, Mesuro S, Rozenfeld D, Yodpinyanee A: Final Report for Community of Ordinary Differential Equations Educators. Harvey Mudd College Joint Computer Science and Mathematics Clinic. 2011, Claremont CA: 301 Platt Boulevard, 9171191711. [http://www.math.hmc.edu/clinic/projects/2010/]
 Java Resources for SEDML. [http://jlibsedml.sourceforge.net]
 Tham LS, Wang L, Soo RA, Lee SC, Lee HS, Yong WP, Goh BC, Holford NHG: A pharmacodynamic model for the time course of tumor shrinkage by gemcitabine + carboplatin in nonsmall cell lung cancer patients. Clin Cancer Res. 2008, 14 (13): 42134218. 10.1158/10780432.CCR074754. [http://dx.doi.org/10.1158/10780432.CCR074754] 10.1158/10780432.CCR074754PubMedView Article
 Wajima T, Isbister GK, Duffull SB: A comprehensive model for the humoral coagulation network in humans. Clin Pharmacol Ther. 2009, 86 (3): 290298. 10.1038/clpt.2009.87. [http://dx.doi.org/10.1038/clpt.2009.87] 10.1038/clpt.2009.87PubMedView Article
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.