Refined elasticity sampling for Monte Carlo-based identification of stabilizing network patterns.

MOTIVATION: Structural kinetic modelling (SKM) is a framework to analyse whether a metabolic steady state remains stable under perturbation, without requiring detailed knowledge about individual rate equations. It provides a representation of the system's Jacobian matrix that depends solely on the network structure, steady state measurements, and the elasticities at the steady state. For a measured steady state, stability criteria can be derived by generating a large number of SKMs with randomly sampled elasticities and evaluating the resulting Jacobian matrices. The elasticity space can be analysed statistically in order to detect network positions that contribute significantly to the perturbation response. Here, we extend this approach by examining the kinetic feasibility of the elasticity combinations created during Monte Carlo sampling.
RESULTS: Using a set of small example systems, we show that the majority of sampled SKMs would yield negative kinetic parameters if they were translated back into kinetic models. To overcome this problem, a simple criterion is formulated that mitigates such infeasible models. After evaluating the small example pathways, the methodology was used to study two steady states of the neuronal TCA cycle and the intrinsic mechanisms responsible for their stability or instability. The findings of the statistical elasticity analysis confirm that several elasticities are jointly coordinated to control stability and that the main source for potential instabilities are mutations in the enzyme alpha-ketoglutarate dehydrogenase.

1 Introduction

Metabolic systems tend to exhibit steady states that can be measured in terms of the concentrations and fluxes of the metabolites involved. These measurements can be regarded as a phenotypic representation of all the complex interactions and regulatory mechanisms taking place in the underlying metabolic pathway. Such interactions determine the system’s response to external perturbations and are responsible, for example, for its asymptotic stability or for oscillatory trajectories around the steady state. However, determining these perturbation responses in the absence of fully specified kinetic models remains an important challenge of computational systems biology.

Structural kinetic modelling (SKM) is a framework to analyse such responses to perturbations, without requiring detailed knowledge about individual rate equations (Reznik and Segre, 2010; Reznik ; Steuer ). It provides a parameterized representation of the system’s Jacobian matrix in which the model parameters encode information about the enzyme-metabolite interactions. Stability criteria can be derived by generating a large number of structural kinetic models (SKMs) with randomly sampled parameter sets and evaluating the resulting Jacobian matrices. The parameter space can be analysed statistically in a Monte Carlo experiment in order to detect pathway positions that contribute significantly to the perturbation response. Because the sampled parameters are equivalent to the elasticities used in metabolic control analysis (MCA) (Kacser and Porteous, 1987), the results are easy to interpret biologically.

In order to perform an SKM experiment, it is not necessary to specify the kinetic rate laws in detail. However, the sampling procedure requires the selection of predefined intervals for each elasticity. The common procedure is to choose the interval (0,1] for sampling elasticities that should represent enzymatic reactions. This interval corresponds to the range of possible elasticity values for reactions following Michaelis Menten kinetics.

In this work, we demonstrate that if the assumption of Michaelis Menten kinetics holds true, only very specific combinations of elasticities allow the steady state to be maintained with non-negative kinetic parameters. Because of this, an easily applicable criterion is introduced that enables detecting elasticity combinations associated with negative kinetic parameters of enzymatic reactions.

Furthermore, we examine how the sampled elasticities can be analysed to derive quantitative thresholds for instabilities. This approach was already used in the past to examine stability conditions of the Calvin Benson cycle (Girbig ). In this work, we illustrate how these conditions affect the underlying kinetic parameters of the model. We examine how a simple example system reacts to inducing such instabilities. Furthermore, we demonstrate how the proposed methodology is applied to study destabilizing mechanisms in the neuronal citric acid (TCA) cycle.

1.1 Introduction to SKM

Given a metabolic system with m metabolites and n reactions , which can be defined by the ordinary differential equation system with stoichiometric matrix N, a steady state is defined as a point in the state space where no net changes in the concentrations can occur and all reaction rates fulfill the mass balance equation The response of the system in steady state to small perturbations depends on the asymptotic stability of the steady state. When it is asymptotically stable, a coordinated system response enables the return of concentrations and fluxes to the same values prior to the perturbation. If the steady state is unstable, such a return is not supported.

Local dynamic properties of a steady state, like stability or oscillatory behaviour, can be derived from the Jacobian matrix evaluated in the steady state (Heinrich and Schuster, 1996). Hence, only if the largest real part of the eigenvalues is negative, changes evoked by perturbations diminish over time and the steady state is asymptotically stable. Oscillatory trajectories occur if eigenvalues form a complex conjugate pair. The oscillation is damped, stable or increasing if the largest real part of the eigenvalues is negative, zero or positive, respectively.

Computation of the Jacobian matrix typically requires knowledge of all kinetic rate laws and kinetic parameters describing the reactions in the system in order to compute the partial derivatives . However, SKM enables the computation of the Jacobian matrix without relying on such knowledge. Instead, it is derived by with and elasticities .

The Jacobian matrix can thus be computed by the simple matrix multiplication where the matrix E contains the elasticities , and is the matrix of normalized stoichiometric coefficients .

Elasticities are a central concept from MCA (Fell and Sauro, 1985; Fell, 1997; Kacser and Porteous, 1987). If associated with enzyme-catalysed reactions, they indicate the amount of saturation of the enzyme with a particular metabolite. Because of this interpretation, the sampled elasticities in SKM experiments are also called ‘saturation parameters’ (Steuer ). In particular, the larger the absolute value of an elasticity, the less the enzyme catalysing the reaction is saturated with its metabolite (Wang ).

As demonstrated by Equation (4), the Jacobian matrix for a given steady state, for which experimental measurements are available, can be derived solely from a set of model parameters (elasticities) , the stoichiometric matrix N, and the steady state measurements and . Although the stoichiometry, the steady state concentrations and the fluxes are experimentally accessible, the elasticities are often unknown in practice. However, due to the normalization step, they are restricted to pre-defined intervals, from which they can be sampled in a Monte Carlo approach. This enables the creation of a large number of models followed by the exploration of the parameter space to detect regions associated with stability or instability.

1.2 Deriving elasticities for enzymatic reactions

The interval boundaries for elasticity sampling are chosen according to the type of kinetics employed by an enzyme. For example, if a reaction v(S) follows irreversible Michaelis Menten kinetics , the corresponding normalized rate law can be derived as

The derivative with respect to the normalized substrate is then given by

Evaluation of the derivative at the steady state (indicated by x = 1) provides the substrate elasticity

Reversely, the original kinetic parameters can be computed from the elasticities for a given steady state by

Repeating the procedure described earlier for the reversible Michaelis Menten equation where S and P are the substrate and product concentrations, leads to the following equations for the elasticities and :

Elasticities of the respective reverse reactions are given by

The relationship to the original kinetic parameters derived from these equations leads the following representations of the K values:

The corresponding maximum velocities are given by

2 Results

2.1 Reducing the elasticity sampling space to ensure kinetically realistic models

When performing an SKM experiment, the sampling interval (0, 1] is often chosen for enzymatic reactions. This interval represents the range of possible elasticity values for reactions following the reversible Michaelis Menten Equation (10). However, calculation of Michaelis constants from elasticities, as demonstrated in Equations (15) and (16) shows that if the assumption of reversible Michaelis Menten kinetics holds true, only very specific combinations of elasticities allow the steady state to be maintained with non-negative kinetic parameters. In particular, a closer inspection of Equations (15) and (16) shows that and are only non-negative for Because Michaelis constants are given in units of concentrations and, as such, are required to be non-negative, elasticity combinations not fulfilling this criterion are kinetically infeasible for any given steady state. We therefore define elasticity filtering as restricting the randomly chosen values to those that fulfil this equation.

2.2 Effects of elasticity filtering on local dynamic steady state properties

In order to systematically evaluate the effects of elasticity filtering, three simple artificial example models were constructed (Fig. 1). For each of these example pathways, a kinetic model with known rate laws and kinetic parameters was created. These kinetic models were used to compute steady state concentrations and fluxes for SKM (details about the pathway structure and steady state values used for SKM can be found in the Supplementary Material to this article.).

Fig. 1.

Pathway topologies underlying the small example models used for SKM analysis

Out of 10 000 SKMs sampled for the linear pathway, only 907 (9.07%) were kinetically feasible in the sense that the randomly sampled elasticities led to non-negative kinetic parameter values according to Equation (19). For the branched and cyclic pathways, only 891 (8.91%) and 944 (9.44%) feasible models remained. These low numbers show that it is important to account for kinetical feasibility when performing SKM experiments.

In order to assess whether filtering for feasible models can impact the results of SKM experiments, probabilities of the different types of dynamic steady state properties were compared before and after filtering. As shown in Figure 2, the chance to observe oscillations strongly increased in each pathway. The elasticity combinations fulfilling Equation (19) tended to have larger values than independently sampled elasticities and hence described fast perturbation responses of their associated enzymes. We can conclude that such fast responses enhanced the chance for observing oscillations in the examined systems.

Fig. 2.

Distribution of dynamic steady state properties in each pathway. The probability of observing oscillations increases after filtering for kinetically feasible kinetic parameters

2.3 Detecting conditions for stability or oscillatory trajectories by multivariate pattern search

So far, most SKM experiments have focused on the detection of individual enzymes to identify single reactions important for maintaining the stability of a steady state (Bulik ; Grimbs ). We recently extended this approach by demonstrating how SKM enables the detection of enzyme- or metabolite-ensembles that act together in an orchestrated manner to coordinate the pathway’s response to perturbations (Girbig , b). This was achieved by replacing the previously used univariate approach by supervised machine learning in order to search for multivariate patterns of elasticities associated with stability or instability.

The decision tree algorithm lends itself to this task because it preserves the original feature space, making it possible to interpret the detected patterns in a biological context (Quinlan, 2013). In order to assess the performance of decision trees in classifying stable versus unstable as well as oscillatory versus non-oscillatory models, trees were trained on balanced datasets retrieved from the small example pathways. In order to analyse the derived decision trees, we studied the corresponding rulesets produced by the C5.0 algorithm for decision trees (Quinlan, 2013). These rulesets summarized all paths through each tree in an easily readable format. After training, the generalizability of each of the obtained decision tree rulesets was assessed by the Laplace ratio where k is the ruleset index, h is the number of hits in the test samples, and e is the number errors on the test samples (Quinlan, 2013).

The smallest ruleset describing instability conditions for the branched pathway contained two criteria. A closer look at this ruleset revealed that it imposed conditions on elasticities and (Table 1). It correctly predicted instability with a Laplace value of 0.98. Both elasticities described the acceleration of reactions v2 and v6 in Figure 1b by their substrates and were restricted to low values. Consequently, the ruleset described a scenario of high saturation of the enzymes by their substrates.

Table 1.

The smallest instability inducing ruleset found for the branched pathway

Elasticity	Threshold type	Threshold value
ϵS_1v2+	≤	0.11
ϵS_4v6+	≤	0.16

We next used the associated elasticity thresholds to compute the corresponding ranges of possible kinetic parameters. In doing so, we set both elasticities to fixed values and computed the matching Michaelis constants using Equations (15) and (16). Maximum velocities were derived in a similar manner using Equations (17) and (18). If the kinetic parameters lay within these boundaries, it was guaranteed that the system would still exhibit the same steady state while fulfilling the constraints on the elasticities shown in the ruleset. As shown in Table 2, the new parameters were distinctly reduced compared with the original values in order to account for the reduced degree of saturation described by the elasticity thresholds in the ruleset.

Table 2.

New kinetic parameters for which the associated elasticities fulfil the instability conditions in Table 1

Reaction	Model	Elasticity	KM+	KM−	Vmax+
v2	New	0.1	0.41	6.57	0.62
	Old	0.44	3.00	30.00	1.00
v6	New	0.1	0.03	2.44	0.11
	Old	0.90	3.00	30.00	1.00

In order to examine the behaviour of the unstable system around the steady state, reactions v2 and v6 were updated by the new kinetic parameters in the kinetic model for the branched pathway. The Jacobian matrix of the updated model evaluated at the steady state now contained the positive eigenvalue 0.0205, which shows that the steady state indeed became unstable. In the original model, the eigenvalues with maximum real parts had been a complex pair . Not only did this show that the new kinetic parameters caused instabilities in the steady state, but the loss of the imaginary part also showed that they prevented oscillatory trajectories around it.

To investigate the global response of the system to the instability, the trajectories of the modified model were simulated starting in the neighbourhood of the now unstable state. As shown in Figure 3, the new stable steady state obtained by the system was not located in close neighbourhood to the old one but showed a strong decrease in the concentration of S1 accompanied by smaller increases in S2, S3 and S4. Computation of the rate laws of the original (unmodified) model at these concentrations confirmed that they did not describe a steady state in the original model.

Fig. 3.

Concentrations of the branched pathway model approaching a new steady state after inducing instabilities

2.4 The neuronal TCA cycle

The TCA cycle is of fundamental importance for cellular energy metabolism because it is the major source for reduced nicotinamide adenine dinucleotide (NADH) and ubiquinol (QH2), which are required for the production of adenosine triphosphate (ATP) in the mitochondrion (Fig. 4). Its reactions are tightly controlled by allosteric feedback regulators to enable an adjustment of the steady state fluxes to varying ATP demands (Nelson and Cox, 2004). In neurons, a steady supply of ATP is crucial for restoring the cellular membrane potential after triggering an action potential. Once the system obtains a functional working state that enables it to meet the ATP demand of the cell, we can therefore expect this state to be robust against perturbations from the cytosol (Koopman ). Such perturbations can arise, for example, due to neuronal activity and the resulting fast fluctuations of ATP turnover.

Fig. 4.

Network model of the TCA cycle and connected pathways. In the TCA cycle, PYR from the cytosol is incorporated into citrate (CIT) and converted to by a series of oxidation steps. The energy released in this process is utilized to form the reduced metabolites NADH and QH2. In the respiratory chain complexes CI, CII and CIV these metabolites are oxidized, driving the import of protons into the mitochondrial matrix. As a result of this process, a proton gradient is created which then serves as a driving force for ATP synthesis. The picture has been adapted from Berndt

To elucidate the mechanisms responsible for perturbation responses in this system, we investigated two steady states representing different scenarios in terms of cytosolic ATP demand. Although state 1 (reference state) was characterized by a moderate extent of neuronal activity, state 2 described a phenomenon called ‘gamma oscillations’. Gamma oscillations are defined as rhythmic brain activity with alternating epochs of enhanced and reduced neuron firing in a frequency of 30–100 Hz (Fell and Axmacher, 2011; Singer, 2013). They have been associated with cognitive processing and memory formation (Brittain and Brown, 2014; Hanslmayr and Staudigl, 2014).

A recently published kinetic model (Berndt ) served as a reference for network stoichiometry and steady-state information. The resulting SKMs covered 24 metabolites, 20 reactions and 71 elasticities. For each steady state, 20 000 SKMs were sampled and the proportions of stable and oscillatory models were determined. For both steady states, the majority of observed models was stable (99.75% for the reference state, 99.85% under gamma oscillations) and a large proportion did not display oscillatory trajectories (oscillatory models in the reference state: 37.24%; gamma oscillations: 39.66%).

The system contained three reactions with only one substrate and product each (aconitase, fumarase, proton leak). The corresponding elasticities were filtered using Equation (19). Using this criterion, 76.08% of all models could be corrected for kinetical feasibility in these reactions within each steady state. In 63.3% of the resulting models the flux control coefficients computed from the elasticities indicated negative control of pyruvate (PYR) import on pyruvate dehydrogenase (PDH). Consequently, an additional criterion was implemented during sampling that only admitted models in which PYR import exhibited positive control on the PDH. The proportions of stable and oscillatory models remained similar after filtering. Stability occurred in 99.56% of all models in the reference state (99.84% under gamma oscillations), and oscillatory trajectories occurred in 34.96% of all models in the reference state (40.06% under gamma oscillations).

Using these models, we next analysed whether multivariate analysis of the elasticity space could lead to more refined information about the system dynamics than univariate comparisons of elasticity distributions. In doing so, we compared the most distinctive elasticities determined by the univariate Kolmogorov–Smirnoff test (KS-test) with those elasticities that occurred most frequently in each stabilizing decision tree ruleset (Table 3). This test enables the comparison of two empirical distribution functions and using the test statistic , which describes the maximum distance between the distribution functions (Darling, 1957). Both methods identified elasticities associated with citrate synthase (CITS) and α-ketoglutarate dehydrogenase (AKGDH) as having strongest influence on stability. In contrast to the KS-test, however, decision trees identified two network positions that had exceptionally strong impact on controlling stability under gamma oscillations but not in the reference state (elasticities marked in red in Table 3). Both positions were associated with reactions downstream of the cycle that were involved in the respiratory chain and ATP synthesis. Because such a shift in control can be expected under increased ATP demand, this result is not a surprise in itself. However, it shows that multivariate pattern search can provide a more refined picture of stability control under varying environmental conditions.

Table 3.

The five most influential elasticities providing stabilizing sites in the TCA cycle

KS-test		Decision trees
State 1	State 2	State 1	State 2
ϵOAACITS	ϵOAACITS	ϵAKGAKGDH	ϵAKGAKGDH
ϵAKGAKGDH	ϵAKGAKGDH	ϵNADHAKGDH	ϵH^+ATP_syn
ϵSUCCOAAKGDH	ϵCITCITS	ϵSUCCOAAKGDH	ϵH^+CIII
ϵSUCCOACITS	ϵSUCCOACITS	ϵCITCITS	ϵOAACITS
ϵCITCITS	ϵICITIDH	ϵOAACITS	ϵSUCCOACITS

In order to understand possible causes for instabilities in the network, decision tree patterns associated with the unstable class were obtained. The simplest pattern contained only two conditions which both restricted substrate activation and product inhibition of AKGDH to low values. The pattern was highly reproducible between both states. This showed that a slow perturbation response of AKGDH posed a threat to the stability of the observed steady states. In order to understand the reasons for instabilities caused by the described conditions, the eigenvectors belonging to the largest real parts of the affected models were analysed. Eigenvectors of the Jacobian matrix can give an approximation of the system’s behaviour in the neighbourhood of the steady state because the trajectory of a perturbation δS in its neighbourhood can be approximated by where and λ are the eigenvectors and eigenvalues of the Jacobian matrix in the steady state. Each eigenvector component therefore describes the time-dependent concentration changes in one metabolite. Analysis of the eigenvector with largest real part indicated that perturbations led to an accumulation of AKG, accompanied by a decrease in the subsequent TCA cycle metabolites (Fig. 5).

Fig. 5.

Eigenvector belonging to the largest positive eigenvalue in unstable models of the reference state that fulfilled the simplest instability condition derived by decision trees. Each eigenvector component describes the time-dependent concentration changes in one metabolite in neighbourhood of the steady state. Concentrations that increased due to a perturbation are shown in blue; decreasing metabolites are shown in red. Instabilities arise due to accumulation of AKG and a depletion in the subsequent metabolites

3 Discussion

In this work, we demonstrated that if the assumption of Michaelis Menten kinetics holds true, only very specific combinations of elasticities allow the steady state to be maintained with non-negative kinetic parameters. Because of this, an easily applicable filtering criterion for enzymatic reactions with a single substrate and product was introduced, that enabled the detection of elasticity combinations associated with negative kinetic parameters.

Without elasticity filtering, the large abundance of SKMs with kinetically implausible elasticity combinations can bias the output of the numerical simulation. For example, analysis of a set of simple example pathways showed that even for small models with 9–11 sampled elasticities, ∼90% of the resulting SKMs contained elasticity combinations that required at least one negative value in the kinetic parameters in order to enable emergence of the observed steady state. Focusing only on those models with non-negative kinetic parameters hardly changed the numbers of stable models, but strongly affected the frequency of oscillatory trajectories around the steady state.

The idea of filtering the sampling space can be extended by taking into account further kinetic or dynamic restrictions on enzymes. For instance, using the thermodynamic principle to limit the analysis to thermodynamically feasible fluxes improves the predictive power of metabolic models (Tepper ). Furthermore, Monte Carlo approaches for modelling metabolism could be refined by incorporating such additional constraints (Murabito ).

The filtering criterion introduced in Equation (19) offers a helpful starting point when sampling elasticities for reactions that can be assumed to follow reversible Michaelis Menten kinetics. If Michaelis Menten kinetics cannot be assumed it might still be possible to adapt the filtering criterion to reflect the underlying approximated kinetics. However, when analysing complex biological systems, one has to be aware that the analytical derivation of similar criteria for more complex rate laws is not possible in the same straight-forward manner. For example, the TCA cycle model analysed contained only four reaction with a single substrate and product each. Although 76.08% of all models could be corrected for kinetical feasibility in these reactions, the challenge remains to detect similar criteria for the remaining reactions in order to refine the results of the Monte Carlo sampling. Recent improvements along this line were presented by Ivanov , who analytically derived general stability criteria for certain structural motifs of irreversible reactions with multiple substrates or products in the absence of regulatory interactions.

In addition to introducing elasticity filtering, we illustrated in this work how decision trees can derive refined multivariate elasticity patterns associated with different dynamic steady state properties. Such patterns can also capture scenarios in which a certain dynamic property only emerges for a small parameter interval. Because the Monte Carlo sampling in a high-dimensional space only allows a sparse coverage of the parameter space, one has to keep in mind that the derived rulesets do not capture all possible scenarios that might theoretically emerge. When applying the proposed strategy to a system of interest, it is therefore important to chose a sample size that is high enough to enable the emergence of reproducible patterns, and to assess the generalization error on independently sampled test sets.

4 Methods

4.1 Model construction

In order to account for the impact of pathway topology on the dynamic properties of steady states of a metabolic network, kinetic models were constructed using three different pathway structures (Fig. 1). Kinetic parameters and rate laws were chosen as uniformly as possible in order to eliminate their influence on the outcome of the analysis. K values were set to the 3 for forward reactions and to 30 for backward reactions. Vmax values were homogeneously set to 1 (Supplementary Material for details). All reactions were modelled using reversible Michaelis Menten kinetics. A positive feedback term was included to increase the range of possible dynamic behaviours. Models were implemented in Matlab using the Systems Biology Toolbox 2 for Matlab (Schmidt and Jirstrand, 2006). Steady states were computed using the SBsteadystate function of the Toolbox. SKMs were constructed using the MATLAB Toolbox for SKM (Girbig ). Elasticity filtering was conducted during Monte Carlo sampling by testing for each set of randomly sampled elasticities whether they fulfilled the criterion in Equation (19).

For TCA cycle analysis, the kinetic model by Berndt served as a reference for network stoichiometry and steady-state information. Two steady states were computed for SKM analysis which represented different degrees of workload imposed on the cell. The first state (also called reference state) represented a situation where cytosolic O2 consumption was 50% of its maximum value. This corresponded to a moderate work amount of ATP consumption by the cell (Berndt ). In contrast, the second state (gamma oscillations) represented a scenario of strong workload and was computed by setting cytosolic consumption to 90% of the maximum value.

4.2 Model analysis

For each steady state, the numbers of stable and unstable SKMs were counted in the unbalanced dataset of 104 randomly sampled models with and without filtering for biological feasibility. Oscillatory and non-oscillatory models were analysed in a similar manner.

Decision tree training was performed using the 5.0 library in R (Kuhn ) in the RULES mode. In total, 105 SKMs were created as independent test samples. Because our aim was to derive reliable conditions for stability and instability, we selected only those rulesets with Laplace ratio > 0.95 that were derived from five trees trained on 20 000 training samples each.