ERNEST: a toolbox for chemical reaction network theory

Abstract

Summary: ERNEST Reaction Network Equilibria Study Toolbox is a MATLAB package which, by checking various different criteria on the structure of a chemical reaction network, can exclude the multistationarity of the corresponding reaction system. The results obtained are independent of the rate constants of the reactions, and can be used for model discrimination.

Availability and Implementation: The software, implemented in MATLAB, is available under the GNU GPL free software license from http://people.sissa.it/∼altafini/papers/SoAl09/. It requires the MATLAB Optimization Toolbox.

Contact:  [email protected]

1 INTRODUCTION

Chemical reaction network (CRN) theory has been developed since the early 70s to study the dynamical evolution of the concentrations of the chemical species involved in a known set of reactions (Feinberg, 1987). Under the assumptions of well mixedness and constant temperature, a system of ordinary differential equations for the species concentrations can be derived from the chemical reactions.

Important aspects in the study of this type of dynamical systems are the existence of equilibria, their number and stability. In particular, the capacity of a reaction system to exhibit more than one equilibrium is not only of interest for chemistry or chemical engineering, but has become a major topic also for biologists, in particular in the study of the switch-like behaviour observed during intracellular signaling (Kaufman and Thomas, 1987; Markevich et al., 2004; Pomerening et al., 2003) or cell differentiation (Yan et al., 2009) processes.

The main difficulty which often arises in the verification of the multistability property is the poor knowledge of the rate functions of the reactions, from the type of kinetics (like mass action, Michaelis–Menten or Hill) to the value of the parameters (called rate constants) involved in each reaction kinetics. The power of CRN theory lies in the fact that its results are based on the network structure alone, and so are independent of the values of the constants and in some cases also of the type of kinetics.

In fact, Feinberg and colleagues found a number of conditions (Ellison, 1998; Feinberg, 1977, 1987, 1988), centred on the concept of network deficiency, under which the dynamical system for a CRN with mass action kinetics does not admit multiple equilibria, regardless of the rate constants. Other different conditions for monostationarity, verifiable on a graph representing the CRN, were then proved (Craciun and Feinberg, 2005, 2006; Schlosser and Feinberg, 1994). These latest conditions, based on system injectivity, were later extended [(Banaji et al., 2007); M.Banaji and g.Craciun, Graph-theoretic criteria for injectivity and unique equilibria in general chemical reaction systems, submitted for publication] also to other kinetics, provided that the system is non-autocatalytic.

An interesting application of these theoretical results to biology is model discrimination (Conradi et al., 2005, 2007; Ellison and Feinberg, 2000): if a biological process is known to be multistable, and for it there are multiple candidate reaction models, it is possible to eliminate some of them by proving that they are always monostationary.

For small CRNs, the various conditions can be verified manually, while for larger networks the only so far available software tool was the Chemical Reaction Network Toolbox, version 1.1a (Ellison and Feinberg, 2000), a closed source DOS program which implements only the criteria based on deficiency theory. We decided to implement a new, open and complete toolbox for all the previously mentioned approaches, in a modern environment with integrated support for the system of differential equations and linear programming like MATLAB.

2 FEATURES AND IMPLEMENTATION

The ERNEST toolbox is structured as a set of MATLAB functions and classes. The analysis is performed by the main function model_analysis, which needs as input a structure specifying the species and reactions of a CRN. The format of this structure is simply the one defined by the TranslateSBML function from libSBML (Bornstein et al., 2008), which imports SBML files in MATLAB. So a SBML model, after the standard import, can be directly analysed by our toolbox, but all the extra information potentially contained in the file, like compartments, constraints, reaction modifiers and kinetic laws, will be ignored.

All the criteria implemented by ERNEST aim to verify conditions on the CRN structure which are sufficient for monostationarity of the relative dynamical system, i.e. to rule out the possibility of multiple equilibria regardless of the rate constants and the initial concentrations.

The model_analysis function operates in the following way:

• (1) calculates complexes, stoichiometric matrix and rank, linkage classes, strong-linkage classes, network reversibility, weak reversibility and deficiency;

• (2a) if the Deficiency Zero theorem (Feinberg, 1977) is applicable, prints out the relative response;

• (2b.1) otherwise calculates terminal strong-linkage classes and deficiencies of the linkage classes;

• (2b.2a) if the Deficiency One theorem (Feinberg, 1987) is applicable, prints out the relative response for mass action kinetics;

• (2b.2b) otherwise verifies that the network is regular, and in case apply the Deficiency One algorithm (Feinberg, 1988) for mass action kinetics; if this algorithm verifies that the network admits multiple positive equilibria within a stoichiometric compatibility class, it also prints out an example of reaction rate constants and two equilibria for the corresponding mass action system;

• (3a) if the CRN is autocatalytic, then calculates the Species–Reaction graph, finds all its cycles and tries to verify the two conditions for monostationarity with mass action kinetics of Craciun and Feinberg (2006);

• (3b) otherwise verifies if the stoichiometric matrix is strongly sign determined (SSD) and in case exclude multiple non-degenerate equilibria Banaji et al. (2007) for general kinetics.

Obviously the runtime of this function increases with the number of species and reactions of the network. The most complex part is the Deficiency One algorithm (point 2b.2b), which involves the solution of linear programming problems with additional sign constraints, but for medium/big reaction networks this code is usually not executed since their deficiency is typically greater than one.

3 CONCLUSIONS

We verified the correctness of our toolbox by successfully reproducing the results for all the examples of Craciun and Feinberg (2006), Craciun et al. (2006) and Conradi et al. (2005), plus others selected from Feinberg's papers.

ERNEST can be quite useful if applied for model discrimination, as in the examples cited above. It has several advantages over the Chemical Reaction Network Toolbox since it verifies more criteria, it is applicable also to kinetics not of mass action type, it can be applied to SBML models, it is multiplatform and open source.

Two possible extensions of the toolbox features are the implementation of the advanced deficiency theory (Ellison, 1998) (which generalizes the Deficiency One Algorithm to CRNs of deficiency greater than one), and the verification of some sufficient conditions for multistationarity, like those in Section 4 of Craciun and Feinberg (2005) or maybe others inspired by Soulé (2003).

Funding: illycaffè, Trieste.

Conflict of Interest: none declared.

REFERENCES

Author notes