DMPy: a Python package for automated mathematical model construction of large-scale metabolic systems

Constraint based metabolic models

which is a set of linear equations that can be solved for an optimal vector v when the system is optimized to, for example, maximize biomass production (although there are a variety of other functions, as shown in [7]). This method is commonly referred to as Flux Balance Analysis (FBA) [4, 6, 12]. Notable examples highlighting the utility of FBA approaches can be found in [3, 13–15] whilst [16] have recently used FBA to obtain a genome-scale model of the human metabolic network, which has seen use in applications such as the discovery of anticancer drugs [17, 18]. Over recent years, a number of extensions to the FBA method have been developed. However, systematic analysis of these methods has shown that there is not one single best-performing method that provides a reasonable match between simulated and measured reaction fluxes [2, 7, 8].

One notable development of FBA is dynamic FBA (dFBA) that aims to match time-dependent changes in system outputs with internal flux dynamics [19–21]. In principle, dFBA is a set of FBA computations conducted independently at several time-points. This results in a set of dynamically changing optimal flux vectors that cover the analysed time-period. However, such an approach can be problematic as the solution of FBA is often non-unique, i.e. there are multiple optimal flux vectors that can solve Eq. 2 for a given set of objective functions [22, 23]. Consequently, discontinuities in the flux vector can be observed when v is plotted against time, which is suggestive of a jump from one optimal flux vector to another [20, 21]. One method of solving this issue is to constrain the flux optimisation such that the time-dependent change from one flux vector to another is not allowed to be large, thus limiting the search space for the optimisation routine [19, 20, 24]. The addition of constraints such as these enforces the reaction fluxes of a network to change continuously through time, in a similar manner to trajectories obtained from Eq. 1, but it is not clear how the trajectories depend on reaction rates or concentrations.

Dynamic metabolic models

with u _l representing the enzyme concentration for reaction R _l and \(\mathbf {k} =\allowbreak \{k^{cat+},\allowbreak k^{cat-},\allowbreak k^{M}_{a},\allowbreak k^{M}_{b},\allowbreak \boldsymbol {\alpha },\allowbreak \boldsymbol {\beta }\}\), where k ^cat are the catalytic rates, k ^M are the Michaelis constants and α _i (or β _j ) the reaction stoichiometry associated with concentration a _i of reactant A _i (or b _j of product B _j ) [9]. What one should notice immediately is that, for even a simple reaction, a large number of kinetic rates need to be known or accurately estimated to ensure that Eq. 3 matches observed data (at least 4 rates for a reaction with a single reactant and product). Thus, whereas the pitfall of dFBA is that it requires a large number of system constraints to provide an accurate solution, the downside to kinetic approaches is that a vast amount of data is required to accurately parameterise models of large metabolic networks.

Parameter balancing

Fortunately, an increasing number of experimentally-measured parameters required for the use of Michaelis-Menten approximations are becoming available in online databases (for example BRENDA [29], SABIO-RK [30] and eQuilibrator [31]). Parenthetically, even if only in vitro estimates of kinetic rates are available, in certain cases the relationship to their in vivo value has been shown [32]. Furthermore, the Parameter Balancing (referred to as PB from hereon) method has been developed to utilise Bayesian inference techniques and include constraints on thermodynamic relationships between different parameters [33, 34]. Thus, one could obtain either measured or realistic estimates for a number of the parameters required to construct a kinetic model and simulate changes in metabolite concentration. Examples of such steps can be found in [1, 11], whilst [28] have shown that relatively small models constructed with fluxes given by Eq. 3 together with measured or estimated parameter values provides better matches to data than other tested functions.

Model reduction

The examples of [1, 11] show two different desired cases of modelling metabolism. In the first instance, [1] construct a detailed model describing the central metabolic pathways of L. lactis whilst, in the second, [11] produce a genome-scale model of yeast metabolism. In principle, the conclusions drawn from larger models should be consistent with those of smaller, more detailed systems and vice versa when created in the same species. Thus one important consideration is that of model reduction. Based on current methods, model reduction could occur at one of three stages: pruning unimportant reactions directly from the genome-scale network [35], grouping subsets of reactions into a single effective reaction on the subsystem of interest [36], or by assessing the kinetic rates of the system using time-scale separation arguments fixing a subset of system components to constant values [37]. Notably, the methods of [35, 36] require only the reaction network and no information about the dynamics of the pathway in order to reduce the system size.

DMPy

Here, inspired by the workflows of [5, 11], we present an automated pipeline that can translate a (genome-scale) reaction network into a dynamic reaction equation model. Given the network as an input, our pipeline automatically integrates measured kinetic rates from different sources (both online databases and measured or estimated values) with the PB technique, optionally reduces model size and, finally, translates the network into ordinary differential equations using Eq. 3. Importantly, this method differs from other (semi-automatic) methods of constructing ODE models (such as CellNetAnalyzer and COPASI [38, 39]) as parameters are obtained from readily-available experimental measurements rather than including an optimisation step to time-series data that can prove very difficult for large-scale systems. By generating simulated data of both the L. lactis central metabolic pathway and randomly generated reaction networks, we go on to analyse the accuracy of our pipeline and determine how many kinetic rates must be measured experimentally to obtain an accurate model. We also show the effect of flux-based model reduction techniques on the resulting system dynamics, extensions to larger networks, and the utility of our pipeline by including compartmentalisation and regulation within metabolic pathways. The presented framework is intended to provide a first approximation of large-scale metabolic dynamics within which further details and computational methods can be added as more data and information come to light. In the Discussion we will highlight how our pipeline can be extended to improve the resulting models (by incorporating different model reduction and parameter optimisation strategies) if the appropriate datasets are available.

Pipeline

In this section we shall provide an overview of the computational workflow. For further implementation details and required inputs please refer to Additional file 1. All simulations and testing were performed using Python version 2.7 (Python Software Foundation, www.python.org) and MATLAB R2012b (MathWorks, Massachusetts, USA).

Figure 1 provides a pictorial overview of the pipeline and consists of three main parts. First, a genome-scale reaction network in SBML (Systems Biology Markup Language, www.sbml.org [40]) format is used as an input into DMPy. The model is parsed and relevant kinetic rates are found from a range of online databases (Steps 1-4, Fig. 1). Then, the obtained parameter estimates (both from online databases and experimental measurements) are used as an input into the PB software that generates a set of thermodynamically consistent kinetic rates (Step 5, Fig. 1). Finally, the kinetic rates are input into Michaelis-Menten functions with similar form to Eq. 3 such that the SBML reaction network can be translated into parameterised ordinary differential equations (Step 6, Fig. 1). The model can now be simulated and edited as required. To look at the effect of model reduction upon system dynamics, we have used flux-based model reduction techniques prior to inputting the genome-scale model into the pipeline [35]. In the Discussion we will explain the implication of doing this and consider alternative methods.

In the following we shall discuss each stage of the pipeline and its purpose.

Parameter gathering (steps 1-4, Fig. 1)

The creation of large-scale metabolic models can be a time-consuming process and this is, in part, due to the manual collection of kinetic rates from online databases that use different naming conventions for reactions and models. To overcome this burden we have created a subroutine that is able to convert between naming structures of models and databases to exhaustively search for all possible measurements of kinetic rates from online databases (Fig. 2). In Additional file 1 we provide the pseudocode for this routine and discuss possible extensions to further constrain the database search for particular cases. The output of the script is a table that is required for input into the PB stage of the computational workflow [34, 41].

In Step 1 of Fig. 2, the databases a user wishes to search are registered. It is noted which kinetic rates they contain and how they can convert reaction or component IDs (identifiers, see below) to these kinetic rates or other IDs. The default databases used in this study are listed in Table 1. In Step 2 (Fig. 2), the SBML model is parsed such that names and identifiers of reactions or metabolites are obtained. This list of model components (both metabolites and reactions) plus the parameters we wish to find represents a set of task objects that are inputs into Algorithm 1 (Additional file 1). Based on these set of tasks, the online databases are searched exhaustively to find all experimentally-measured rates related to the metabolite reactions of interest.

Database	Purpose	Reference
BRENDA	Michaelis constants, catalytic rates	[29]
SABIO-RK	Kinetic rates	[30]
eQuilibrator	Equilibrium constants	[31]
Rhea	Identifier translation	[42]
MetaCyc	Identifier translation	[43]
MetaNetX	Identifier translation	[44]

One issue with using online databases is that different sources use different naming conventions, e.g. BRENDA uses EC numbers, whilst eQuilibriator uses KEGG identifiers [25-27] to classify reactions. Consequently, methods of integrating these identifiers and translating from one database to another are required [29–31]. Thus in Step 3 (Fig. 2) a set of possible paths is constructed that links online data values to the rates needed to complete the model. The Rhea, MetaCyc and MetaNetX databases are used here to translate identifiers [42–44]. For example, assume we have a reaction name and we wish to obtain an equilibrium constant, k ^eq , then one possible path would be:

$$\begin{array}{*{20}l} &\mathrm{Name\ (Glucose-6-phosphate\ isomerase)}\\ &\qquad\qquad\qquad\qquad\rightarrow \mathrm{MetaCyc\ ID\ (ENZRXN-2863)} \\ &\qquad\qquad\qquad\qquad\rightarrow \mathrm{Rhea\ ID\ (11816)} \\ &\qquad\qquad\qquad\qquad\rightarrow \mathrm{KEGG\ reaction\ ID\ (R00771)} \\ &\qquad\qquad\qquad\qquad\rightarrow \text{eQuilibrator}\ k^{eq} \mathrm{ (0.361)}. \end{array} $$

Each of these pathways are exhaustively examined (Step 4, Fig. 2) until all the useful experimental measurements have been obtained. Full provenance and, if available in the original database, literature references are saved and can be manually examined if needed, for example, in the case of conflicts between identifiers.

Parameter balancing (step 5, Fig. 1)

Upon obtaining measured kinetic rates from online sources (and combining these with experimental measurements), a table of found values and their reference is automatically constructed, which is used as an input into the PB software. In order to improve the scalability of the parameter balancing algorithm and to support the parameter balancing of large genome scale models, the parameter balancing algorithm was implemented directly into the pipeline based on the original implementation by Lubitz et al. [34]. Here we shall briefly review the PB method and further details can be found in [33, 34, 41].

The principle of PB is to relate measured constants with unknown values via the Haldane relationships and Wegscheider conditions that ensure kinetic rates are thermodynamically feasible whilst reactions take place in an ideal solution [33]. These relationships are able to relate equilibrium constants, k ^eq , to catalytic rates, k^cat± (s^-1), maximal velocities, v ^max (mmol ·s^-1), Michaelis constants, k ^M (mM), enzyme concentrations u (mM), and standard chemical potentials, μ° (kJ ·mol^-1) via

$$\begin{array}{*{20}l} \ln{k_{j}^{eq}} &= -\frac{1}{RT}\sum_{i}n_{ij}\mu_{i}^{\circ}, \\ h_{j}\ln{k_{j}^{eq}} &= \ln{k_{j}^{cat+}} - \ln{k_{j}^{cat-}} + \sum_{i}h_{j}n_{ij}\ln{k_{ji}^{M}}, \\ \ln{k_{j}^{cat\pm}} &= \ln{k_{j}^{V}} \mp \frac{h_{j}}{2}\sum_{i}n_{ij}\left(\mu_{i}^{\circ}/RT + \ln{k_{ji}^{M}}\right), \\ \ln{v_{j}^{max\pm}} &= \ln{u_{j}} \,+\, \ln{k_{j}^{V}} \mp \frac{h_{j}}{2}\sum_{i}n_{ij}\left(\mu_{i}^{\circ}/RT \,+\, \ln{k_{ji}^{M}}\right), \end{array} $$

where k^cat+ is the forward catalytic rate, k^cat− is the backwards catalytic rate, \(k_{j}^{V} = \sqrt {k_{j}^{cat+}k_{j}^{cat-}}\) (s^-1) is the geometric mean of catalytic rates, h _j is the cooperativity factor for sigmoidal kinetics, n _ij is the stoichiometric coefficient of metabolite i in reaction j, R is the gas constant (J ·mol^-1 ·K^-1) and T is the temperature in Kelvin [34].

Furthermore, one can relate chemical potentials, μ (kJ ·mol^-1), and reaction affinities, R ^A , to metabolite concentrations, c (mM), using

$$\begin{array}{*{20}l} \mu_{j} &= \mu_{j}^{\circ} + RT\ln{c_{j}}, \\ R^{A}_{j} &= -\Delta_{j}G = -\sum_{i}n_{ij}\mu_{i}, \end{array} $$

(4)

where ΔG is the chemical potential through a reaction [34]. The concentrations, c, obtained from the balancing routine are used as initial conditions when simulating the resulting ordinary differential equations (see below).

Based on these relationships, any unmeasured values can be calculated directly given that a subset of the system parameters are known or can be estimated. This is the principle behind PB. By bringing together as much information as possible into prior distributions, Bayesian approaches are used to obtain the maximum likelihood estimates for kinetic rates that follow the thermodynamic relationships above. This results in posterior distributions for each kinetic rate in the metabolic system (e.g. Fig. 3, green distribution), providing a feasible range for these parameters [34]. However, in some cases where a system has not been well characterised, prior information for certain values may not be found within the databases or measured experimentally. In this case a pseudo-prior distribution is used to approximate the range of values one may expect to see experimentally (compare blue and red distributions in Fig. 3a-c as an example). Further details about the construction of our pseudo-prior distributions can be found in Additional file 1.

Test models

To highlight the utility of our pipeline, we will illustrate our analysis for multiple systems. In order to fully analyse the accuracy of the pipeline, we have used the central metabolic pathway of L. lactis [1]. The initial mathematical model of this system is already in the correct form of Eq. 7 for most of the reactions and provides ranges for kinetic rates within the system that we can include in our pipeline. To show that the pipeline can also be used for larger scale systems, we use the E. coli core and the iJO1366 genome scale models [51, 52] in combination with the optional reduction steps.

In order to show how our automated search for kinetic rates improves the number of parameters found from online databases, we also included the S. cerevisiae model iTO977 and the iJO1366 E. coli model that is well annotated relative to other species [52, 53].

Finally we have included randomly generated reaction networks in order to facilitate testing effects of including regulation and compartmentalization. See Additional file 1 for details on the generation of these models.

Pipeline analysis tests

In order to test the accuracy and reproducibility of our framework we generated a simulated dataset from a parameter balanced model of our tested reaction networks. To obtain the model we input the reaction network into our pipeline such that we have an idealized ‘gold standard’ in silico system with fixed values (using the medians of the balanced distributions) for every parameter in the network (Fig. 4). Using randomized subsets of the fixed parameter values as prior inputs into a second round of PB we simulate the availability of a subset of the parameter data. The subset of data is used to create a new model through parameter balancing and is subsequently simulated. Finally, the mean square error is calculated compared to the simulation of the ‘gold standard’ reference model. We then performed the following tests:

• how many times the posterior distributions need to be sampled before there is convergence in the mean square error between the samples and the simulated ‘gold standard’ data;

• the amount of prior distributions that need to be known to obtain minimal differences between the model obtained from the pipeline and simulated data, and;

• the effect of altering the width of the posterior distribution sampled from.

Convergence was tested by taking the standard error of deviation relative to the mean score. After an initial period of 25 samples, the simulation was marked as converged when \((\sigma _{s} / \sqrt {n}) / \mu _{s} < 0.05\), where σ _s is the standard deviation of the score, μ _s is the mean score and n is the number of simulations. If convergence was not obtained within 10000 samples, the simulation run was halted and marked as unconverged.