Systems biology tools for toxicology

Abstract

An important goal of toxicology is to understand and predict the adverse effects of drugs and other xenobiotics. For pharmaceuticals, such effects often emerge unexpectedly in man even when absent from trials in vitro and in animals. Although drugs and xenobiotics act on molecules, it is their perturbation of intracellular networks that matters. The tremendous complexity of these networks makes it difficult to understand the effects of xenobiotics on their ability to function. Because systems biology integrates data concerning molecules and their interactions into an understanding of network behaviour, it should be able to assist toxicology in this respect. This review identifies how in silico systems biology tools, such as kinetic modelling, and metabolic control, robustness and flux analyse, may indeed help understanding network-mediated toxicity. It also shows how these approaches function by implementing them vis-à-vis the glutathione network, which is important for the detoxification of reactive drug metabolites. The tools enable the appreciation of the steady state concept for the detoxification network and make it possible to simulate and then understand effects of perturbations of the macromolecules in the pathway that are counterintuitive. We review how a glutathione model has been used to explain the impact of perturbation of the pathway at various molecular sites, as would be the effect of single-nucleotide polymorphisms. We focus on how the mutations impact the levels of glutathione and of two candidate biomarkers of hepatic glutathione status. We conclude this review by sketching how the various systems biology tools may help in the various phases of drug development in the pharmaceutical industry.

Introduction

In recent years, system-level understanding has been a recurrent theme in biological science and has greatly improved our understanding of the function and regulation of many different processes. In particular, the new discipline systems biology has provided a framework for investigating the interactions between the separate parts of biological systems in order to understand its functioning. One field where this approach has the potential to add significant value is toxicology. In particular, it is now recognised that many commonly prescribed drugs cause rare and unexpected adverse effects in man. These drug-induced adverse reactions are a leading cause of human ill health, restricted drug usage, failed licensing and withdrawal of licensed drugs from clinical use. Therefore, they are of major clinical and socioeconomic importance.

The mechanisms that underlie adverse drug reactions are complex and involve multiple biological processes. The latter may include cellular uptake and efflux mediated by specific membrane transport proteins, metabolic biotransformation (often to chemically reactive metabolites), off-target interactions of the drug and/or its metabolites with tissue macromolecules, downstream interactions between cells and activation of innate and/or adaptive immune responses. Our understanding of the inter-relationships between these events is incomplete. In particular, it is unclear which of the different processes can explain why relatively few patients are susceptible to adverse drug reactions whereas the vast majority of drug-treated patients are not. As an essential first step, a rigorous understanding is required of the intracellular networks affected by drugs, their metabolites and their molecular targets, and of how these can act in concert with either causing toxicity or preventing toxicity from arising.

It is a frequent mistake to assume that systems biology is focussed on system-level or holistic understanding (Alberghina and Westerhoff 2005). It is not: physiology is the science of the whole; molecular biology is the science of the parts; systems biology is the scientific discipline that encompasses and describes relationships between the two. Systems biology therefore involves assessment of biological function at the level of network interactions. The relevant networks start at the molecular level, are genome-wide and ultimately relate to the functioning of the organism as a whole. An important feature of systems biology is that it focuses upon the emergent properties of the system. Emergent properties are functional properties not present within the individual components of the system, which only arise when system components interact. An example is the interaction between hydrogen and oxygen to make water. The resulting change in properties is unpredictable if only the individual properties of hydrogen and oxygen are known (Aderem 2005). Systems biology is confined to cases where the emergent properties are important for biological function, that is, for the maintenance of the living state, for dysfunction such as in disease and for malfunctions such as in the case of toxicity. In the context of this review, systems biology addresses much of the toxicity of adverse drug reactions (ADRs).

Once the system can be understood to the point that emergent properties can be explained, then biology can start to move from being a descriptive science to being a more predictive science (Materi and Wishart 2007). Once biology is a more predictive science, it should become possible to predict off-target effects of drugs as well as network effects of affected targets. Even a minor advance in predictive capabilities may help both the research into new drugs and their toxicity, and decisions concerning whether or not to advance new drug candidates towards the market (Westerhoff et al. 2008).

Glutathione metabolism is important for toxicology as it is vital in binding and neutralising reactive species produced in, for example, the liver via the metabolism of xenobiotics by the formation of glutathione conjugates (Ketterer et al. 1983). When the dose of reactive metabolites exceeds the capacity of this defence system to replenish glutathione through new synthesis, glutathione depletion occurs. For a variety of toxic compounds, it has been found that the depletion of glutathione precedes the covalent binding of chemically reactive metabolites to cellular macromolecules, oxidative stress and, ultimately, organelle injury with the potential for liver failure and death (Park et al. 2011). A more complete understanding of the events that lead to hepatic glutathione depletion might therefore highlight new and sensitive biomarkers and provide insights into key steps in glutathione production, thereby leading to improvements in risk assessment for new, and existing, drugs and xenobiotics.

In this review, we shall provide an outline of the available in silico systems biology tools that might be of use to toxicology. We will discuss to what extent these tools have been applied already and how they can be beneficial in toxicology. In addition, these tools are exemplified by being applied to a model of glutathione metabolism, where possible, to illustrate their utility.

Creating a model

Biological systems are highly complex due to both the number of components present and the nonlinear interactions between them (Kitano 2002). In order to understand emergent properties of the system, the interactions between the components must be studied. Components such as enzymes only ‘talk’ to the rest of the system through the rates at which they consume or produce small molecules. The reaction rates are governed by kinetics such as Michaelis–Menten and Hill equations, which express the dependence of the reaction rates on the concentrations of the small molecules in the systems; this is how the enzymes ‘listen’ to each other through the altering concentrations of the small molecules. Only rarely, the integration of the concentration of the small molecules over a period of time can be understood without computer help. An important part of systems biology is therefore collecting the experimental data of components of a system into a mathematical model (Alberghina and Westerhoff 2005) and integrating that model over time. Such models of biological systems have many uses (Hornberg et al. 2007), some of which are mentioned here. Firstly, comparisons between experimental observation and the mathematical model behaviour can link knowledge of system components to explanations of system behaviour. These new explanations might have been elusive in the past due to the focus on single molecules. Secondly, in silico experiments can be carried out on the model to test the effects of perturbations on the system and to identify the processes that control the system. These experiments may either be only feasible using a computer, or are faster and cheaper than laboratory experiments (Bakker et al. 2000). Such ‘dry’ experiments may generate new hypotheses about the system, which can then be tested experimentally. Simulation can also be used as a tool for experimental design. In an iterative process, the results of further experiments can then be used to improve the model until such a time as the model is sufficiently good as to be considered equivalent to its biological test systems in terms of ability to predict outcomes of experiments. At this point, the biological experiments become redundant.

There are several examples of modelling applied in toxicology (Krishnan and Andersen 2010). One popular type is physiologically based pharmacokinetic (PBPK) models. These models investigate variables such as the rates of absorption, distribution, excretion, and biotransformation of chemicals and their metabolites. In toxicology research and chemical risk assessment, they are used to make more accurate predictions of target tissue dose for different exposure situations (Andersen 1995). Other models simulate electrophysiological activity to assess the potential cardiotoxicity of cardiac drugs (Rodriguez et al. 2010), identify the important parameters in hepatocyte regeneration following toxicity using single cell–based spatial–temporal models (Höhme et al. 2007) or show the feasibility of examining toxicology pathways in kinetic detail (Reed et al. 2008).

For toxicology, using models in terms of rate equations and balance equations, together culminating in ordinary differential equations (ODEs), would be beneficial as it is a well-understood formalism, fast and mathematically robust. It also lends itself to thorough interpretation via metabolic control analysis (MCA), elementary mode analysis (EMA) and flux balance analysis (FBA), which are tools that can help understand how a system functions (vide infra).

In the ODE methodology, the biochemistry of the reactions is essentially translated into mathematics. The biological system and the corresponding network of chemical reactions are described in terms of a set of balance equations with reaction stoichiometries indicating which metabolite is produced or consumed in which reaction, plus a set of reaction rate equations. The ODEs resulting from the combination of these two sets are then solved using numerical methods. Most processes in biological systems are catalysed by enzymes or transporters. Examples of well-known and popular equations used to describe enzyme catalysed reactions are mass action and Michaelis–Menten kinetics (Cornish-Bowden 1995). For more information on enzyme kinetics and the equations used for this kind of modelling, please refer to ‘Appendix’.

In the past, modelling has been unattractive to experimental biologists due to the necessity to acquire extensive programming and mathematic knowledge. This should no longer be the case. A number of software packages have been put together to perform several mathematical tasks and to make modelling more user-friendly. Different software packages have been designed for users of different backgrounds as reviewed in some detail by Alves et al. (2006). Modelling simulators such as COPASI provide access to powerful systems biology tools while remaining accessible to non-expert modellers (Hoops et al. 2006). Examples of kinetic models are available in model repositories such as JWS online [http://www.jjj.bio.vu.nl/; (Olivier and Snoep 2004; Materi and Wishart 2007)] that enable in silico experimentation with the existing models without exposing the user to the detailed mathematics [cf. Biomodels (http://www.ebi.ac.uk/biomodels-main/)]. Models include those of glycolysis in human erythrocytes (Mulquiney and Kuchel 1999), glycolysis in Trypanosoma brucei (Helfert et al. 2001) and the methionine/threonine metabolism in Arabidopsis (Curien et al. 2009).

Below, a model of glutathione metabolism by Geenen et al. (in press; cf Fig. 1) will be investigated as an example. The model can be found under (http://jjj.mib.ac.uk/webMathematica/Examples/run.jsp?modelName=geenen). It was based on kinetic equations taken from a previously published model (Reed et al. 2008) and expanded by adding equations for the gamma-glutamyl cycle, ophthalmic acid synthesis and detoxification of paracetamol (acetaminophen). Some of the parameter values used for the creation of the model are uncertain, and the lack of data required for validation means that we are not confident at the ability of the model to make quantitatively correct predictions. Here, we use this model to illustrate systems biology tools and to create testable hypotheses for further analyses. The more specific conclusions of the model are of limited certainty.

The network of biochemical reactions modelled is presented in Fig. 1. Starting from methionine at the top, a long pathway leads to reduced glutathione, which can then be conjugated to a xenobiotic such as paracetamol. This then leads to the removal of the xenobiotic as the glutathione conjugate. Glutathione can also be made from 5-oxoproline, which can in turn be synthesised from glutathione through a partly extracellular route. The precursor of glutathione can release a cysteine to return to 5-oxoproline. The overall scheme is complex, and the presence of at least three reaction cycles makes it difficult to predict how changes in the activity in any of the enzymes (or the gene that encode them) will affect the level of glutathione, the rate of detoxification, or the concentrations of metabolites and potential biomarkers such as 5-oxoproline.

Fig. 1

The reaction network in the mathematical model of glutathione metabolism made by Geenen et al. (in press) for liver. The network embraces methionine catabolism, glutathione metabolism, 5-oxoproline and ophthalmic acid synthesis and glutathione-mediated detoxification pathways. Metabolites are indicated by ovals; enzymes and rates are shown in italics

Systems biology computational methods

A mathematical representation of biology is not the aim of systems biology. Rather, computational models are viewed as tools to help the human mind in understanding the behaviour of the networks of biology. The insight gained by modelling is more extensive than what is possible by experimentation alone. This is because the networks of biology are too complex for their behaviour to be predicted by the unaided human brain. In biology simplicity is rare, because living systems are essentially complex (Westerhoff et al. 2009). These models therefore can also get complex, and thus we need modelling tools to be able to understand these networks and particularly the role that specific molecules such as drugs, drug targets and drug metabolisers play in them.

Because detailed experimental analysis of the networks in humans is virtually impossible, and tissue culture cells and animal models are rarely reliably representative of the in vivo situation in humans, it might hereby seem that mathematical models are the methods of choice. This is not yet so, however. The catch is that for the models to be made, much experimental information is required; none of the models of the type required here can be made ab initio. It is a property of the highly nonlinear networks of life such that their behaviour depends strongly on the values of many parameters that are products of the evolution of the macromolecules. These parameter values can almost only be determined experimentally. To date there is no single network in biology for which a completely validated dynamic model has been made, although some pathways are perhaps getting close. A slightly but not yet very successful (Westerhoff et al. 2010) alternative method is found in flux balance analysis with objective functions: If one knows which functions the network has been optimised for in evolution, then one may arrive at complete descriptions of systems fluxes.

Steady state analysis

Once the concentrations of the chemical substances in the system have integrated so as to give zero further change in concentrations over time, a steady state is reached. The concentrations and fluxes of the model at this steady state often represent the functions of the biological system that is being modelled. Steady states therefore can contain a lot of information. The steady state of a model can be readily calculated in software such as COPASI, or even more directly, that is at the click of a button in JWSonline (http://jjj.mib.ac.uk/) (Olivier and Snoep 2004).

For the diagram of Fig. 1, Geenen et al. (in press) have collected the best present estimates of the parameter values, as well as the best known rate equations, and have assembled these into a rate plus balance equations model, which can now be integrated by JWSonline. The column entitled ‘Normal’ in Table 1 gives the steady state concentrations and flux values for the reactions, as predicted by the model. It is clear that without the model, it would have been difficult to predict these flux values, which correspond to what is implied by the best of our understanding of the glutathione pathway. On the other hand, this understanding is incomplete, and we should not put too much trust in the precise values of the fluxes and metabolites at this stage. If not for the precise values, what else can one use the model for?

Table 1 Changes in the steady state of the Geenen et al. (in press) standard model under normal conditions, when V _GSf1 is reduced from 948 to 200 or when V _OP is reduced from 846,930 to 10,000

We will first use the model to verify our understanding of the fundamentals of the modelling strategies and of certain principles of networks functioning at steady state. The flux values in Table 1 enable us to appreciate that at each ‘node’ (a node is any metabolite in the network for which the concentration may vary with time until the steady state is achieved), the total carbon flux into the node equals the total carbon flux out of it. For instance, at the node cglc (γ-glutamyl-cysteine) in Fig. 2, the flux synthesising glutamyl-cysteine is 940 μM/h (v₁₀) and the flux degrading glutamyl-cysteine is 830 μM/h through reaction 24 and 110 μM/h through reaction 11. This results in a net rate of change of zero and thus no change in the concentration of glutamyl-cysteine with time.

Fig. 2

The conservation of the glutamyl-cysteine concentration at steady state. Reactions v₁₀, v₁₁ and v₂₄ were taken from Fig. 1. The steady state reaction rates in the model are given in italics (from Table 1)

Another semi-quantitative issue that the modelling can address is what the major fluxes are in the pathway. Is the flux from methionine to glutathione indeed the dominant flux? The answer is perhaps surprising: no. According to Table 1, the production of γ-glutamyl-cysteine (cglc) from cysteine and glutamate in reaction 10 at a rate of 940 μM/h, the production of glutamate from 5-oxoproline at a rate of 930 μM/h (reaction v₂₅) and the production of 5-oxoproline from γ-glutamyl-cysteine (v₂₄) at a rate of 830 μM/h are ten times higher than the flux through reaction 8. This makes for a large cyclic flux around reactions v₂₄, v₂₅ and v₁₀ of 830 μM/h to be dominant. This may be unexpected considering that this is a so-called ‘futile cycle’, which hydrolyses ATP. This suggests that either a strong ability to regulate, which could be an effect of the futile cycle, is to be discovered here, or our molecular understanding of the pathway is still incomplete. Perhaps not all enzymes of the cycle are expressed simultaneously.

Interindividual variation in glutathione metabolism

The same model can help understand the disease 5-oxoprolinuria. 5-Oxoprolinuria is an inherited disorder of glutathione metabolism, which causes a decrease in glutathione concentration and an increase in 5-oxoproline concentration (a glutathione metabolism by-product) in the urine. The literature is unclear as to the mechanism of 5-oxoprolinuria. Some researchers have suggested that this occurs via a deficiency in glutathione synthetase (GS). They presumed that this deficiency causes a low glutathione concentration, which they expected to overstimulate γ–glutamyl-cysteine synthesis, and thus 5-oxoproline formation (Shi et al. 1996). Another possible cause of 5-oxoprolinuria is a disruption in the 5-oxoprolinase enzyme (which catalyses reaction 25 in Fig. 1), which is hypothesised to directly increase the 5-oxoproline concentration in the cell (Croal et al. 1998). We would hypothesise that in practice there is a plethora of causes of 5-oxoprolinuria, that is, any single-nucleotide polymorphism affecting an enzyme activity with control over the level of 5-oxoproline without causing embryonic lethality. But that is not the issue here.

Here, we examine the realism of the former two explications. If one decreases the V _max of GS in our model by a factor of 4.7, the concentration of reduced glutathione decreases, but by a mere third, while the concentration and the export rate of 5-oxoproline (v₃₈) increase by factors of 2 and 3, respectively. Using the scan function of JWSonline, it is possible to analyse more precisely how the model predicts the export rate v₃₈ to vary with the V _max of GS (V _GSf1). Fig. 3 shows that a reduction in the glutathione synthetase (GS) activity is predicted to increase the efflux for 5-oxoproline from the cells, confirming deficiencies in GS as possible basis for the disease.

Fig. 3

The variation of 5-oxoproline export (V ₃₈) with the V _max of glutathione synthetase (V _GSf1), as predicted by the model. Calculated using the scan function of JWSonline (scanning V _GSf1 from 1 to 1,000)

The other proposed explanation for 5-oxoprolinuria may also be examined by an in silico experiment, decreasing the 5-oxoprolinase activity (V _OP; reaction 25 in Fig. 1). In this case, no significant increase in 5-oxoproline export (v₃₈) was observed (Fig. 4). From this, we hypothesise that the more effective cause of 5-oxoprolinuria is a decrease in the glutathione synthetase, not a decrease in the 5-oxoprolinase activity.

Fig. 4

5-Oxoproline export (V ₃₈) predicted for the variation of the V _max of 5-oxoprolinase (V _OP). The scan function of JWSonline (scanning V _OP from 1 to 10,000) was used

These findings are partly counterintuitive. A decrease in the activity of GS (reaction 11) might be expected to decrease the level of extracellular glutathione and thereby the synthesis of 5-oxoproline through reactions 28, 30 and 31, with a concomitant decrease in the intracellular level and efflux of 5-oxoproline, contrary to the modelling results. Again, inspecting Fig. 1, one might expect a decrease in the V_max of reaction 25 to lead to an increase in the 5-oxoproline concentration and an increase in 5-oxoproline efflux (reaction 38), again deviating from the results of the modelling. By looking at the steady state values in Table 1, one can examine why these intuitive predictions would be incorrect (illustrated in Fig. 5). When GS activity is reduced, this should first reduce the rate of reaction 11, thereby causing the concentration of γ–glutamyl-cysteine (cglc in Table 1) to increase and the concentration of glutathione (cGSH) to decrease as confirmed by the model calculations (v₁₁ changes from 110 μM/h to 79 μM/h, cglc changes from 0.19 to 0.64 mM and GSH changes from 1.5 to 1.0 mM; cf. the third column in Table 1). The rather substantial increase in γ–glutamyl-cysteine might be expected to decrease the flux through v₁₀. Paradoxically, this flux is increased, however, by some 7 % from 940 to 1,000 μM/h. Consequently, the cause of the increase in 5-oxoproline cannot involve product inhibition of reactions 10 and 25. Rather, it must be due to the increase in the concentration of γ–glutamyl pushing more flux through reaction 24 (from 830 to 930 μM/h) increasing the rate of production of 5-oxoproline. Apparently, this increase is enough to compensate for the decrease in v₃₁ due to the decrease in v₁₁ and the increase in 5-oxoproline, and consequent increase in glutamate (cf. Table 1) now leads to the paradoxical increase in reaction 10, as mentioned above. This pushes the intracellular 5-oxoproline up, with a consequent increase in 5-oxoproline export. The above finding shows that indeed the proper cause–effect relationship in a network as complex as glutathione metabolism cannot be reliably identified by mere inspection of network topology. Modelling is required for following the logics of biology.

Fig. 5

Changes in fluxes around glutathione induced by reactions v₁₀, v₁₁, v₂₄, v₃₁, v₃₈ and v₂₅ taken from Fig. 1 in order to show the hypothesised mechanism of 5-oxoprolinuria. The changes in steady state flux for the shift in condition between normal (V _GSf1 = 948) and reduced GS (V _GSf1 = 200) are shown in the figure. Flux is in μM/h

This is illustrated further when V _OP (reaction 25) is reduced by a factor of 8.5 (final column in Table 1): this brings about an increase in the 5-oxoproline concentration by only 25 % (from 2.5 to 3.1 μM) so that reaction v₂₄ changes direction and degrades 5-oxoproline back into the glutathione cycle instead of much more of the substance being exported from the cell: the export reaction increases by a mere 20 % (from 0.57 to 0.79). We emphasise that we are here discussing the performance of what may be the best possible mathematical model of the pathway and that these conclusions are subject to experimental validation. The calculations show, however, how the model can serve as a vehicle for explanation and a method for the suggestion of key experimental measurements.

Metabolic control analysis

In order to understand emergent properties of the system, the interactions between the components must be studied. But studying the emergent properties does not always require the full dynamics of the network. If one is interested exclusively in the control of fluxes and concentrations in the network, then there is a simpler way to determine and quantify the control of emergent properties and to calculate how control emerges from the interactive properties of the components of the network: metabolic control analysis (MCA). In principle, this approach could help to determine the processes in the network that determine the toxicity of a given drug.

Although it is still the prevailing view that the flux through a pathway is necessarily set by one rate-limiting enzyme; this view has been shown to be false (Groen et al. 1982a). MCA enables one to deal with the possibility that the control of flux through the pathway is shared. The flux control coefficient ($C_{{_{i} }}^{J}$) for example is defined as the percentage change in a specified steady state metabolic flux (J) that is caused by a one per cent change in activity (v) of any enzyme alone, or more precisely:

$$ C_{{v_{i} }}^{J} = \frac{{\left( {\frac{\partial \ln J}{\partial p}} \right)_{{{\text{system}}\;{\text{steady}}\;{\text{state }}}} }}{{\left( {\frac{{\partial \ln v_{i} }}{\partial p}} \right)_{\text{all\,metabolite\,concentrations\,constant}} }} $$

In words, one here adds an inhibitor p specific for an enzyme i. One measures the relative (‘fold-change’) effect on the flux of the pathway in the intact system (i.e. δlnJ). One then compares this (by division) to the effect the inhibitor would have on the flux through the enzyme if the enzyme were studied in isolation δlnv _i . After all, in the pathway, homeostatic responses will reduce the effect on the flux. If the inhibitor were to reduce the pathway flux by only 20 % at a concentration where it would inhibit the enzyme in isolation for 50 %, the control coefficient would be only 20/50 = 0.4.

Here, the activities of all other enzymes are kept constant (Kholodenko et al. 1995), and p is a parameter that only affects enzyme i. When the control coefficient is 1, then the change in flux is proportional to the change in the activity of the reaction. When its flux control coefficient is 0, the enzyme is not limiting the flux at all. These two options already existed in the classical view of the control of metabolism, that is, an enzyme was either the rate-limiting step (C = 1) or not rate-limiting at all (C = 0). The above definition enables one to deal with the intermediate case where an enzyme carries partial flux control, that is, somewhat but not completely rate-limiting. In such cases, the flux control might be 0.3 or 0.7. This has allowed the establishment of the pattern of flux control for a number of cases both directly and experimentally (Bakker et al. 1999; Groen et al. 1982a) and by an integration of experimental data into a model (Bakker et al. 1999; Groen et al. 1982a).

As an increase in all enzyme activities by 1 % leads to a 1 % increase in flux through a pathway, it can be shown that the sum of the flux control coefficients is 1 (Hornberg et al. 2007). This is called the summation law (Kacser and Burns 1973; Heinrich and Rapoport 1974):

$$ \sum\limits_{i = 1}^{n} {C_{i}^{J} = C_{1}^{J} + C_{2}^{J} \ldots + C_{n}^{J} = 1} $$

The same principle can be applied to concentration control analysis, where the effect of a change in metabolite concentration on the flux is measured. Then, the analogy of the above sum amounts to zero.

A related parameter is the elasticity coefficient. This is a local parameter of an enzyme that shows how perturbations of a reaction parameter affect the local reaction rate. ‘Local’ means that the effect is to be measured without taking into account the changes in the environment of the enzyme that may be caused by its reduced activity. The elasticity coefficient does not take network effects into account, whereas the control coefficients do. Indeed, the perspective of elasticity differs from that around the control coefficients, as the elasticity is not directly a systemic property, although its magnitude does depend on properties of the system.

The elasticity coefficients are defined as the ratio of relative change in local rate to the relative change in one parameter (normally the concentration of an effector).

$$ \mathop \varepsilon \nolimits_{p}^{i} = \frac{{\partial v_{i} }}{\partial p} \cdot \frac{p}{{v_{i} }} = \left( {\frac{{\partial \ln v_{i} }}{\partial \ln p}} \right)_{{{\text{all}}\;{\text{other}}\;{\text{metabolite}}\;{\text{concentrations}}\;{\text{constant}}}}^{{}} $$

where v _i is the rate of the enzyme and p is the parameter we are perturbing. This means one could measure the relative (‘fold-change’) effect on the rate of a reaction in response to/compared with (by division) a relative (‘fold-change’) in the parameter of interest. Each enzyme has as many elasticity coefficients as the number of parameters that affect it directly, for example, reaction substrates, products and effectors.

A particularly useful and important feature of MCA is that it can relate the kinetic properties of the individual reactions (local properties) to (global) properties of the whole intact pathway. MCA therefore shows how local properties within the system give rise to system-level properties, that is, MCA shows the emergence. This is done through the connectivity theorems that relate the control coefficients to the elasticity coefficients.

The connectivity theorem for flux control coefficients (Kacser and Burns 1973) states that, for a common metabolite S, the sum of the products of the flux control coefficient of all (i) steps affected by S and its elasticity coefficients towards S is zero:

$$ \sum\limits_{i} {C_{i}^{J} \varepsilon_{[S]}^{i} = 0} $$

For the concentration control coefficients, two equations can apply, depending on whether or not the metabolite, the concentration of which is being controlled, is the same as the metabolite for which the elasticity coefficients are considered (Westerhoff and Kell 1984):

$$ \begin{array}{*{20}c} {\sum\limits_{i} {C_{i}^{[S]} \varepsilon_{[S]}^{i} = - 1} } \hfill \\ {\sum\limits_{i} {C_{i}^{[A]} \varepsilon_{[S]}^{i} = 0\quad {\text{for}}\,S \ne A} } \hfill \\ \end{array} $$

Control coefficients can also be applied in drug efficacy measurements by measuring the effect of an inhibitor of an enzyme on the flux through the system. A corresponding response coefficient is defined as the percentage change in a specified steady state metabolic flux that is caused by a one per cent change in concentration of a drug or by an amount of drug that causes 1 % inactivation of the enzyme when in isolation (Groen et al. 1982b; Chen and Westerhoff 1986; Bakker et al. 2002).

An example of how modelling of a system can lead to the discovery of new drug targets can be found in the metabolic control analysis in the glycolytic pathway of the parasite T. brucei. This parasite is responsible for the African sleeping sickness and lives part of its population life cycle in the mammalian blood stream. The strongest flux control resided in the glucose transport step. This led to the prediction that the glucose transporter is the best predicted drug target, rather than glyceraldehyde-3-phosphate dehydrogenase, which is the drug target that is mostly worked on (Bakker et al. 1999). The next step is here to perform differential control analysis, which compares the flux through a single metabolic pathway in a pathogen and host, finding the control coefficients in both. The steps with the largest difference in control coefficient (higher control in pathogen than in host) can be identified as having the greatest potential for specifically inhibiting flux through the pathogen metabolic pathway. Bakker et al. (2002) integrated this new approach with molecule-based drug design. T. brucei is also a model system for tumour cells in the mammalian host, where the tumour cells need to be killed by drugs that should not be toxic for the host cells. Again, one should develop drugs that impinge on steps with high flux control in the tumour cell and small flux control on, hence low toxicity for, normal host cells. Differential control analysis has also been performed for oncogenic signalling to gain information on potential drug sites (Hornberg et al. 2007).

By analysing which reactions have control or high fragility in a detoxification pathway, it is possible to predict which steps in a pathway are sensitive to drug-induced changes. Additionally, by analysing how control might change in situations such as malnutrition and genetic defects, it may be possible to predict an individual’s sensitivity to toxicity from a certain drug.

MCA has been applied to the glutathione model to demonstrate the control of methionine influx on the ability of the cell to protect itself from paracetamol toxicity (Geenen et al. in press). We can also use MCA to explain the findings we have made in this paper regarding 5-oxoprolinuria. The control coefficients of the model were analysed by using the metabolic control analysis task in COPASI. This model is available in SBML format from the authors and from JWSonline, which can be loaded into COPASI for analysis. Table 2 shows the flux control coefficients of GCS, GS and OP on 5-oxoproline export (i.e. the effect that a change in the individual V _max’s of GCS, GS or OP would have on the export rate of 5-oxoproline) under three conditions. Under normal and glutathione synthetase–limiting conditions the control is strongest in the glutathione synthetase (control of −0.71 and −0.79) showing this is the more limiting step. The negative value of these control coefficients reflects that an increase in the activity of the enzyme leads to a decrease in the flux, as expected. When V _OP is decreased, the control in 5-oxoprolinase is also decreased, again supporting the expectation that decreased 5-oxoprolinase does not lead to much of an increase in 5-oxoproline.

Table 2 The control exercised by GCS, GS and OP on 5-oxoproline export flux under three conditions: normal steady state, reduced glutathione synthetase activity (V _Gsf1 = 200 rather than 948) and reduced 5-oxoprolinase activity (V _OP = 10,000 rather than 846,930). Flux control is quantified in terms of the flux control coefficient (Groen et al. 1982a)

Hierarchical regulation analysis (RA)

Control analysis examines which steps control or limit a flux or a concentration. This gives an answer to the ‘what-if question’: ‘what is the percentage effect on the flux if an enzyme is activated by 1 %?’. If a certain enzyme has a flux control coefficient of, say, 0.7, then this does not guarantee that the system itself actually activates this enzyme by 1 % when it needs to increase the flux through the system by 0.7 %. The organism that is host to the pathway may or may not activate the enzyme when it responds to an external challenge, such as a drug. It has been the development of regulation analysis that has brought this crucial distinction between control (=limitation) and actual regulation to the fore. Hierarchical regulation analysis (ter Kuile and Westerhoff 2001) is a quantitative method for investigating whether the regulation of flux through a system is metabolic or depends on alterations in gene expression.

When an organism needs to increase the steady state flux through a pathway, it needs to increase the fluxes through all the enzymes in that pathway. For each enzyme, the organism is able to increase the flux by regulating it through gene expression, through signal transduction (leading to covalent modification and activation of the enzyme), or metabolically (i.e. by increasing the concentration of its substrate or allosteric modifiers or by decreasing the concentration of its product). The former two types of regulation have been quantified together by the so-called hierarchical regulation coefficient (ρ_h), consisting of a gene expression regulation coefficient and a signal-transduction regulation coefficient. The metabolic regulation has been quantified by the metabolic regulation coefficient (ρ_m). At steady state, the sum of ρ_h and ρ_m equals one. This means that 100 % of regulation is distributed between gene expression, signal transduction and metabolic regulation (ter Kuile and Westerhoff 2001; Westerhoff et al. 2009).

Historically, the flow of information from DNA to RNA to protein to function has suggested to some that regulation is exclusively hierarchical (i.e. that ρ_h = 1) and in fact dominated by regulation of gene expression. Others thought that only regulation at the metabolic (ρ_m) level, thus through the concentrations of substrate, products and modifiers, was of importance. Experimental regulation analysis has examined this for Trypanosoma brucei (ter Kuile and Westerhoff 2001) and Saccharomyces cerevisiae (Rossell et al. 2006; Daran-Lapujade et al. 2007) and found the regulation to be distributed between gene expression and metabolic regulation and to vary over time (Eunen et al. 2009). Both earlier groups of scientists were right and wrong.

ρ_h is measured as the ratio between the percentage change in enzyme concentration and percentage change in flux accompanying (and being partly caused by) this change.

$$ \rho {}_{h} = \frac{{\log V_{\max c1} - \log V_{\max c2} }}{{\log J_{c1} - \log J_{c2} }} $$

In words, the difference of the logarithm of the V _max between condition 1 and 2 (c ₁ and c ₂ ) is divided by the difference of the logarithm of the flux flowing through the enzyme in c ₁ and c ₂ . If a precisely proportional change in enzyme level (or V _max) accompanies the change in flux through that enzyme, then ρ _h will be one and there will be no regulation of flux through metabolite level changes. If this is not the case, then some of the flux is metabolically regulated (ρ _m ). ρ _m is given by the percentage change in the metabolites’ concentration dependence of the enzyme rate divided by the change in flux. Because the two coefficients add up to 1, and ρ _h can be measured relatively readily, ρ _m is usually calculated from the former.

$$ \rho_{m} = 1 - \rho_{h} $$

In the presence of drugs, fluxes such as the synthesis rate of glutathione may change either because of the direct effect of the drug or its metabolites or because of subsequent regulation of the fluxes by the cell, in response to the addition of the drug. The regulation can consist of a combination of three components, that is, alteration in the levels of intermediary metabolites such as ATP or glutamate, alterations in the expression levels of the participating enzymes or alterations in covalent modification of the enzymes. Regulation analysis enables one to dissect the three effects and to suggest ways to optimise the response.

Metabolic network analysis: elementary mode analysis (EMA), metabolic flux analysis (MFA) and flux balance analysis (FBA)

MCA and RA flourish by looking at specific aspects of dynamic networks only, that is, control and regulation, respectively. One may also focus on flux itself, in the sense of what fluxes are admitted by a given network model (EMA), what fluxes are actually occurring (MFA) and what fluxes correspond to the network fulfilling certain aims (FBA). The map of a metabolic network (or a mathematical model thereof) is completely defined by all the reaction stoichiometries: every reaction in a biological system has a stoichiometry, which defines the number of substrate and product molecules that are consumed and produced, respectively, per turnover of the enzyme (Westerhoff and van Dam 1987). The stoichiometric matrix is one where the rows represent the compounds of the reactions and the columns correspond to the reactions. For the network given in Fig. 6, the stoichiometric matrix is (1 −1 −1), indicating that the change in time of the concentration of metabolite X is v₁–v₂–v₃.

Fig. 6

Demonstration of how a network can be decomposed into a set of elementary modes and extreme pathways

Elementary mode analysis (EMA) studies the possible routes through a biochemical network. This requires only stoichiometric data. A flux mode is a set of fluxes through reaction steps in the network for which each metabolite is balanced in the sense that its production rate equals its net consumption rate. The sum of two flux modes is also a flux mode, and thereby, some flux modes can be described as linear combination of two other flux modes. An elementary flux mode is a minimum set of reactions that can operate at steady state (Fig. 6). The mode is termed non-decomposable or minimal in the sense that removal of any of the reactions in the mode (knock out of any of the corresponding genes) destroys the flux (Schuster et al. 1999). It can also be described as a minimal set of reactions needed to convert one reactant into one metabolic product or a small set of the latter (Schilling et al. 2000).

Another term often used in this kind of analysis is ‘extreme pathways’. Extreme pathways are a subset of elementary modes and are distinguishable in that they cannot be expressed as the combination of more than one elementary mode. Elementary mode analysis enumerates all distinct metabolic routes through a network, that is, all possible flux pathways in a generalised sense. Extreme pathways analysis on the other hand focuses on enumerating the unique and minimal set of convex basis vectors needed to describe all allowable steady state flux distributions through the metabolic network as sums (not differences) of the extreme pathways (Papin et al. 2003). In large biological networks, many, more than one, elementary modes can be present. These elementary modes may overlap, using the same reactions or metabolites, depending on the connectivity of the metabolites. The number of, and reaction stoichiometries within, elementary modes can be calculated using software such as COPASI (Hoops et al. 2006). For genome-wide networks, the number can be so large that it defeats the purpose of calculating all of them. In the example of Fig. 6, the three elementary modes are (1 1 0), that is, v₁ = v₂, = v₃ = 0, and (1 0 1), that is, v₁ = v₃, and (0 1 −1), that is, v₂ = reverse v ₃. Because the first pathway can be seen as the sum of the latter two, there are only two extreme pathways, that is, the latter two elementary modes.

Elementary modes may be associated with biological functions: the catabolic network of an organism may need to produce glutamate for protein synthesis and glutathione for detoxification. Starting from glucose and ammonia, there may be only a few elementary modes (fundamental pathways) doing this. In order to understand life, we need to improve our knowledge and understanding of the role of the elementary modes and which, how and why certain ones or certain combinations are used more than others.

EMA was performed on the glutathione model using the JaPathways software package for the calculation of elementary mode flux values (Taylor and Schwartz, manuscript in preparation), which is based on the quadratic programming approach suggested by Schwartz and Kanehisa (Schwartz and Kanehisa 2005). Steady state flux values calculated for 3 paracetamol concentrations from Geenen et al. (in press) were used as input to the JaPathways software along with the stoichiometric coefficients of the model. This resulted in a complete set of 86 elementary modes and the corresponding elementary utilisation (flux going through each elementary mode). The elementary mode with the highest utilisation was EM1 (shown in Fig. 7), consistent with our identification of a dominant flux mode shown above.

Fig. 7

The schematic of EM1 and its usage (flux going through this elementary mode at steady states of three paracetamol concentrations)

Elementary modes are related to function. By comparing the flux flowing through the elementary mode associated with methionine being used for glutathione recycling (EM29) to the flux through the elementary mode associated with methionine being used for detoxification (EM28), we could analyse and quantify how the cell deals with different concentrations of paracetamol (Fig. 8). As paracetamol is increased, the flux going through glutathione recycling (EM29) reduces as the cell is using more of its glutathione to detoxify paracetamol (EM28). By analysing the flux through an elementary mode rather than the flux through a single reaction, we are able to compare ‘functions’ of the model rather than merely reactions.

Fig. 8

Scheme of EM29 (elementary mode associated with methionine being used for glutathione recycling) and EM28 (elementary mode associated with methionine being used for detoxification) and the flux going through each elementary mode at three paracetamol concentrations

The ‘what is possible question’ is addressed by EMA. Although software is in place to calculate the flux through elementary modes, this software is at early stages of development and is considered rather specialised (Poolman et al. 2004; Schwartz and Kanehisa 2005; Schwartz and Kanehisa 2006). The ‘how’ question is related to kinetic properties of and expression levels of enzymes catalysing the steps in the elementary modes. This question is not addressed by EMA, but by kinetics, MCA and HRA. The ‘why’ question may have to do with the performance of the elementary modes in terms of efficiency. This question is addressed by FBA.

The ‘which’ questions are addressed by metabolic flux analysis (MFA). The aim of metabolic flux analysis (MFA) is to quantify the amount of flux going through the metabolism of an organism. The end result would be a flux map that shows the distribution of anabolic and catabolic fluxes over the metabolic network (Wiechert 2001). This is the combination of flux data gained from, for example, ¹³C-labelled substrate with the stoichiometric matrix of the network. This can use either of the two approaches. Method 1 integrates the ¹³C data with computer models of the network being measured. Calculated fluxes through the model are iteratively fitted to measured data, thereby minimising the difference between the observed and simulated isotope spectra, similar to a parameter-fitting procedure. Method 2 studies flux ratios and is derived by probabilistic equations that quantify the relative contribution of converging pathways to the formation of a metabolite from the NMR or mass spectra. However, this method is only appropriate for smaller data sets (10–15 fluxes) (Wiechert 2001).

An example of MFA is a study on the carbon cycle of hepatoma cells to investigate potential drug targets. Using a combination of metabolite concentration information and transient ¹³C-labelling experimental data obtained in hepatoma cell lines, a network model of glycolysis, the pentose-phosphate pathway (PPP) and the citric acid cycle (TCA) was set up. The fluxes from the model were then estimated and found to be in accordance with in vivo ¹³C-labelling data (Hofmann et al. 2008; Maier et al. 2008). The comprehensive network model of Maier et al. (2008) was extended and applied to modelling cholesterol and central carbon metabolism. This found 3-hydroxy-3-methylglutaryl-coenzyme A reductase (HMG-CoA reductase) to be a potent target for lowering cholesterol synthesis with hypolipidemic drugs. By applying MCA to this model, they found a high (0.5) control coefficient of HMG–CoA reductase over cholesterol efflux (Maier et al. 2009).

Another example of successful flux analysis is the investigation of hepatic metabolic pathways in order to provide nutritional support and non-surgical medical therapies for fulminant hepatic failure (FHF) (Arai et al. 2001). Here, a mass balance model was created to characterise changes in liver metabolism during the hypermetabolic response to injury. By combining experimental data about amino acid perturbation during FHF with the mass balance model, it was possible to estimate intrahepatic metabolic fluxes under treatment conditions. D-galactosamine treatment, for example, significantly decreased hepatic gluconeogenesis, which correlated with a reduction in amino acid entry into various points of the tricarboxylic acid (TCA) cycle.

By contrast, flux balance analysis (FBA) does not rely on the input of experimental data or any expression of kinetic information. Instead, it simulates the metabolic network by using the stoichiometric matrix to find the pattern of fluxes while optimising for certain criteria such as the maximal production of certain products or the minimisation of ATP hydrolysis (Stelling et al. 2002; Kauffman et al. 2003). This linear optimisation will ideally result in a single steady state reaction flux distribution in a metabolic network (Raman and Chandra 2009), but in practice it does not do so, as there are many parallel pathways in most networks with, for example, equal ATP stoichiometries.

The investigation of the fluxes is very important for a network understanding of a system. The ability to manipulate the fluxes through the system would give scientists more control over toxicology pathways. The prediction of the effects of complete inhibition of certain reactions can be carried out by FBA, but only when the optimisation criteria ‘used’ in the evolutionary selection for fitness is known, which is rarely the case (Simeonidis et al. 2010). The manipulation of fluxes by partial impediment or by stimulation of steps in the network will require MCA or full-out kinetic modelling.

Robustness

Robustness is a system’s ability to respond to changes in the external conditions or internal organisation while maintaining a constant behaviour and a similar steady state (Barabasi and Oltvai 2004). Metabolites in biological systems are clustered into nodes. Some of these nodes are hubs that are connected with many other nodes. The topology of a system has a high importance in robustness. While deletions of individual nodes of a system affect the system to a small degree, elimination of hubs causes major disruption as this leaves small isolated node clusters uninteracting (Albert et al. 2000). Not only topology but also gene duplication can play an important part in robustness (Daniels et al. 2008). To understand network robustness requires investigation into the functional and dynamic changes that a perturbation causes (Barabasi and Oltvai 2004), thereby viewing robustness as the systems property that it is seen to be (Barkai and Leibler 1997).

Many different mathematical approaches have been used to compare robustness of systems (Daniels et al. 2008; Grimbs et al. 2007). Westerhoff showed that the higher the robustness of a variable towards changes in enzyme activities, the lower the fragility and the lower the corresponding concentration control coefficient (Westerhoff 2007). He also showed (in preparation) that total fragility of a concentration is always conserved and equal to 0, but that total robustness is not conserved and is high when some fragilities are close to zero and high when all fragilities differ strongly from zero. Robustness can be calculated as being the inverse of the control coefficients (Koefoed et al. 2002; Swat et al. 2011).

In Geenen et al. (in press), it was shown that the concentrations of proposed extracellular biomarkers of intracellular glutathione levels were sensitive to methionine concentrations in the cell. To investigate this, the robustness coefficients were calculated for the biomarkers’ secretion flux (v₃₂ and v₃₈) when perturbing methionine import (v₃₉) and methionine entry into the methionine cycle (v₁ and v₂). The fragility coefficient was evaluated as the inverse of the flux control coefficient. The results in Table 3 show that at zero paracetamol, the robustness is fairly low, with a slightly higher robustness in 5-oxoproline secretion. This means that the 5-oxoproline secretion is less likely to change than the ophthalmic acid efflux when the flux of v₃₉, v₁ and v₂ changes, making 5-oxoproline the more robust biomarker of glutathione levels. The negative robustness merely means that an increase in the flux through the perturbed process results in a decrease in the 5-oxoproline or ophthalmic acid flux. Interestingly, at low paracetamol (20 μM), where Geenen et al. (in press) showed 5-oxoproline could behave as a biomarker for glutathione depletion, oxoproline secretion flux is more robust to changes in v₃₉, v₁ and v₂ than in the absence of paracetamol. In contrast, at high methionine concentrations, 5-oxoproline secretion flux again has a lower robustness again, and ophthalmic acid secretion flux has a high robustness with regard to v₃₉, v₁ and v₂. This suggests that 5-oxoproline might be a better biomarker at low paracetamol exposure and ophthalmic acid a better biomarker at high paracetamol concentrations.

Table 3 The robustness of ophthalmic acid and 5-oxoproline secretion vis-à-vis perturbation of methionine import (v ₃₉) and methionine entry into the methionine cycle (v₁ and v₂)

Another way to make quantitative predictions on the relative importance of various reactions in a network uses flux balance analysis. Several publications have shown that the flux distribution in organisms like Escherichia coli is not homogeneous (Almaas et al. 2004). Most reactions have quite small fluxes while a few have high flux. It is likely that perturbations affecting these high flux reactions will alter the performance of the system most.

Conclusion: how could systems biology help systems toxicology?

As mentioned above, we have applied a range of systems biology tools to the glutathione detoxification pathway to illustrate how these tools can be used to understand how toxicologically relevant pathways function. We have shown how steady state analysis of a system can allow us to investigate system behaviour. In silico experiments in which the V _max of the enzyme glutathione synthetase (GS) in the glutathione pathway was perturbed have allowed us to simulate a 5-oxoprolinuria patient and to predict the subsequent increase in 5-oxoproline production. Metabolic control analysis (MCA) has enabled us to explore the impact of individual reactions on control of flux through the glutathione pathways. In the case of 5-oxoprolinuria, this analysis showed that a decrease in the V _max of GS resulted in an increased impact of this enzyme on glutathione synthesis flux, thereby making GS a more rate-determining enzyme: this mechanism of 5-oxoprolinuria should become progressively more severe as the molecular defect is stronger. By contrast, the mechanism of a defect in 5-oxoprolinase should have a limited oxoprolinuria effect, which should not increase further as the molecular defect increases in severity. Tools like metabolic control analysis may well be of more general utility in toxicology, for example, for determining which enzymes in multistep detoxification pathways are rate determining and can be expected to exert the greatest effect when perturbed by dietary effects or single-nucleotide polymorphisms. This has been exemplified in the demonstration that the supply reactions of methionine have high control on paracetamol detoxification, thereby predicting that reduced amounts of this amino acid (e.g. following malnutrition) could lead to problems in detoxification (Geenen et al. in press). We have also applied tools such as metabolic flux analysis (MFA), flux balance analysis (FBA) and elementary mode analysis (EMA) to the glutathione pathway, in order to explore the mechanisms by which changes in discrete analytes arising from the pathway correlate with glutathione levels and therefore provide useful biomarkers. It was shown that EMA can give us information about the flux through a function, such as detoxification, rather than just the rate of the reaction. Robustness is an important systems property to analyse for toxicology, as it can predict the expected level of change as the system responds to external conditions such as a toxic effect. By analysing robustness of the glutathione depletion biomarkers, we could investigate under what condition these biomarkers were robust and therefore more likely to exhibit a reliable correlation with glutathione status. We would recommend that when biomarkers are proposed, not only their sensitivity towards changes in the phenomenon they are marking should be discussed, but also their robustness with respect to various SNPs and nutritional variations. This comes on top of our earlier recommendation that the variation of the biomarker concentration with the property that the biomarker is supposed to monitor should be monitonic (Geenen et al. in press).

When applying systems biology tools to toxicology problems, it pays to define the specific issue that needs to be addressed. To do this, it is important to take full account of the context. Within the pharmaceutical industry, toxicity caused by candidate drugs remains a leading cause of failure of progressing candidate drugs to the market. Pharmaceutical R&D comprises of distinct Discovery and Development activities (Fig. 9). Drug Discovery aims to predict the efficacy and toxicity in man during compound design and selection, and then to assess potential toxicity in man by undertaking preclinical safety studies in experimental animals. The approaches used at this early stage are focussed on reproducible and dose-dependent adverse effects, not on infrequent events that may occur only in relatively few highly susceptible individuals. During drug Development, studies are undertaken that provide accurate assessment of efficacy and toxicity exhibited by patients in clinical trials, which as molecules progress through the development pipeline typically involve progressively larger numbers of patients and longer durations of compound administration. This results in a shift in emphasis from ‘average efficacy and toxicity across the human population’ in Discovery and in early clinical trials to ‘individual response and toxicity in each patient’ in Development, that is, in later trials and after licensing. Where marked variability in efficacy or in toxicity is observed in man, the opportunity arises for personalised health care in which the relevant drug may be given only to certain patients who are considered likely to benefit (e.g. the anticancer drug Iressa is specific to patients with EGF-stimulated tumour cell growth (Wakeling et al. 2002; Ciardiello et al. 2000)), or may not be given to patients who are at high risk of exhibiting a serious adverse event (e.g. the association in Han Chinese between a genetic marker, HLA–B*1502, and carbamazepine-induced Stevens–Johnson syndrome (Chung et al. 2004)).

Fig. 9

The drug development process

In principle, systems modelling has the potential to aid this process significantly. Modelling of data from non-clinical sources may improve the accuracy of predictions of effects that occur in the clinic before the first time in man (FTIM) clinical trial. Such effects may arise from processes with significant control that exhibit inter-species variation between animals and man or inter-individual variation between humans. Figure 10 provides examples of how such considerations can impact upon the risk of drug-induced liver injury (DILI) caused by paracetamol, where the toxicity is caused by a reactive intermediate, which at normal therapeutic doses is detoxified by conjugation to glutathione.

Fig. 10

How can systems biology tools be applied to the risk of drug-induced liver injury (DILI) from paracetamol?

We can investigate the toxic effect of paracetamol in the liver in several different ways. Using pharmacokinetic modelling, it is possible to predict the concentration of paracetamol in the liver following administration of different doses and in different individuals (who vary in age, size, etc.). By modelling the metabolism of the drug, it is possible to predict how the steady state flux and concentration changes with different doses. This can increase the understanding of paracetamol toxicity and also identify doses at which its administration is safe and does not cause a depletion in glutathione and consequently protein binding. Tools like metabolic control analysis allow us to predict which steps have the highest control on important concentrations and fluxes in the intracellular network. They can therefore be applied to toxicology to find the steps in pathways that are more sensitive to toxic effects.

By using systems biology tools such as control and robustness analysis vis-à-vis toxicology pathways, one should be able to provide better understanding and prediction of toxicological risk. The available tools can assess which parts of a toxicology pathway are most fragile and to which variations the fragilities are high. This includes the variability caused by single-nucleotide polymorphisms in individuals putting them at risk of a, then predictable, toxic effect. Once such an analysis has been undertaken, part of the subsequent biomonitoring activities could focus on the genetic basis of the most controlling steps rather than requiring the analysis of all possible components of the pathway. This could help assess the potential range of toxicity for the human population, but also open up avenues for personalised medicine. Here, individual genome sequences and the frequencies of the polymorphisms in the human population might have impact, at last. A different set of biomonitoring activities would read the patients at more physiological levels, ranging from the expressed proteome, to metabolomics and personal drug tolerance history. Here, the results would be compared with results of systems biology models such as the one discussed here for glutathione.

In the case of paracetamol, it has already been shown experimentally that the methionine supply to the cell is a controlling factor for glutathione detoxification ability. This has been put into mechanistic perspective by modelling. Reduced methionine levels are linked to lower glutathione detoxification ability. Since hepatic glutathione protects against liver toxicity caused by the reactive metabolites of paracetamol, an understanding of the range and variability of the hepatic methionine status in the human population would be expected to improve our understanding of the range of doses at which paracetamol is safe in man. For a candidate drug that causes DILI via mechanisms similar to those described for paracetamol, this approach could greatly assist with point 2 in the drug Development phase goals (Fig. 9) and allow us to more accurately translate from clinical trial data to human population safety predictions. If it is known how an individual may be different from the average of the population in a controlling step, we can more accurately calculate individual risk and potentially avoid adverse toxic effects. This would then be addressing point 3 in the drug Development phase goals and act as a step towards personalised health care. In principle, a similar systems modelling approach could be used to explore and understand risk of toxicity arising via other mechanisms, provided that the key pathways responsible for toxicity and for detoxification can be defined and flux controlling steps can be identified.

One of the main impacts of the glutathione model was that the increase in our understanding of the robustness of the biomarkers 5-oxoproline and ophthalmic acid. Biomarkers like these could play an important role in tracking toxicity both during clinical trials and in patients. Modelling approaches would allow us to predict and understand the kinetics, magnitude and dynamic range of these biomarkers. Therefore, these techniques could allow us to optimise biomarker usage and help make decisions with regards to prioritisation between candidates. This would greatly assist with point 1 in the drug Development phase goals of Fig. 9 and possibly enable improved monitoring of drug safety in man in clinical trials.

We have shown that systems biology tools give more rigorous descriptions of complex biological processes and effects of perturbations. This has led and will lead to a better understanding of controlling events and regulatory processes in the system. Toxicology is a consequence of the loss of homeostatic regulatory processes, that is, of a loss of robustness (for example, in the case of toxicity due to paracetamol, cell protection is overwhelmed by high doses of reactive metabolites). Therefore, as illustrated here, the application of systems biology tools to toxicology has the unique opportunity to provide network insights into underlying mechanisms and basis of susceptibility to drug compounds. We also infer that by integrating data for individuals, such as enzyme activities or diet into a model, it may in the future be possible to help understand why certain individuals are more susceptible to adverse effect of drugs than others. It may become possible to adjust drug dosing accordingly, making drugs more effective for individuals not subject to their toxicity and safer for those who are.

In the above sections, we asserted that in principle, systems biology could be far superior over current toxicology methods: clinical trials in human, in vivo work in animal models and in vitro drug testing in cell lines. We here wish to repeat that this superiority is far from being in place. For models to attain the required level of quality and detail, much more experimental work is needed. At present, the model quality is insufficient to claim that the precise calculations in this review lead to reliable results. In systems biology, unlike in the genomics revolution 10 years ago, the necessary work flow towards understanding drug targeting and toxicity is clear. It is not known with certainty that this work flow will become successful across the lines, but it is likely to become so in some cases, such as glutathione-mediated detoxification of paracetamol. Perhaps it is time to try.

List of abbreviations

The complete names of the enzymes and metabolites indicated by acronyms in Fig. 1 are as follows.

Enzyme names and acronyms and EC numbers

 

Reaction	Abbreviations	Enzyme name	E.C. number
v[1]	mati	Methionine adenosyl transferase I	2.5.1.6
v[2]	matiii	Methionine adenosyl transferase III	2.5.1.6
v[3]	meth	Glycine N-methyltransferase	2.1.1.20
v[4]	gnmt	DNA methyltransferase	2.1.1.72
v[5]	ah	S-adenosylhomocysteine hydrolase	3.3.1.1
v[6]	bhmt	Betaine-homocysteine methyltransferase	2.1.1.5
v[7]	ms	Methionine synthase	2.1.1.13
v[8]	cbs	Cystathionine gamma-synthase	4.2.1.22
v[9]	ctgl	Cystathionase	4.4.1.1
v[10]	gcl	Glutamylcysteine synthetase	6.3.2.2
v[11]	gs	Glutathione synthetase	6.3.2.3
v[12]	gpx	Glutathione peroxidase	1.11.1.9
v[13]	gr	Glutathione reductase	1.8.1.7
v[24]	ggct	Gamma-glutamylcyclotransferase	2.3.2.4
v[25]	oxoase	5-oxoprolinase	3.5.2.9
v[26]	gcs	Glutamylcysteine synthetase	6.3.2.2
v[27]	gs	Glutathione synthetase	6.3.2.3
v[28]	ggtp	Gamma-glutamyltranspeptidase	2.3.2.2
v[29]	ap	Aminopeptidase	3.4.11.2
v[31]	ggct	Gamma-glutamylcyclotransferase	2.3.2.4
v[33]	gpx	Glutathione S-transferase	2.5.1.18
v[34]	gpx	Glutathione S-transferase	2.5.1.18
v[35]	ggtp	Gamma-glutamyltranspeptidase	2.3.2.2
v[36]	ccat	Cysteine-S-conjugate N-acetyltransferase	2.3.1.80

Names of metabolites of which the concentrations were variable (μM)

• THF—Tetrahydrofolate

• 5,10-MTHF—5-10-Methenyltetrahydrofolate

• 5-MTHF—5-Methyltetrahydrofolate

• met—Methionine

• SAM—S-adenosylmethionine

• SAH—S-adenosylhomocysteine

• hcy—Homocysteine

• cyt—Cystathionine

• ccys—Cytosolic cysteine

• bcys—Blood cysteine

• glc—γ-Glutamyl-cysteine

• cGSH—Cytosolic glutathione

• bGSH—Blood glutathione

• cGSSG—Cytosolic glutathione disulphide

• bGSSG—Blood glutathione disulphide

• cgly—Cytosolic glycine

• cglut—Cytosolic glutamate

Names of metabolites of which the concentrations were held constant

• AB—2-Aminobutyrate (Soga et al. 2006)

• bgly—Blood glycine

• bglut—Blood glutamate

• bmet—Blood methionine (varies in some experiments)

• cser—Cytosolic serine

• H₂O₂—Cellular hydrogen peroxide

• HCHO—Formaldehyde

• OPA—Ophthalmic acid, N-[N-(γ-glutamyl)-α-aminobutyryl]glycine

• Oxo—5-Oxoproline, pyroglutamic acid

• para—N-acetyl-p-aminophenol/paracetamol/acetaminophen

Appendix: Further details on kinetic modelling

Kinetic equations

This appendix contains additional information on the equations that are typically used in kinetic modelling.

For a reaction where nA + mB gives pC + qD, the mass action kinetics rate equation is:

$$ v = v_{ + } - v_{ - } = k_{ + } \cdot [A]^{n} \cdot [B]^{m} - k_{ - } \cdot [C]^{p} \cdot [D]^{q} $$

(1)

The properties between the square brackets refer to the chemical activity or concentration of the compounds in question. Mass action kinetics describes the time dependence of chemical reaction systems where the molecule numbers are much higher than 100 (Onsager 1931; Westerhoff and van Dam 1987) and all participating molecules are explicit in the reaction equation.

Whether a reaction in a network needs to be modelled as reversible or not depends on the thermodynamics, that is, on the free energies. Eq. 1 can be rewritten as (Westerhoff and van Dam 1987):

$$ v = v_{ + } - v_{ - } = k_{ + } \cdot [A]^{n} \cdot [B]^{m} \cdot (1 - \frac{{K_{\text{eq}} }}{\Upgamma }) $$

(2)

where the mass action ratio is defined by:

$$ \Upgamma = \frac{{[A]^{n} \cdot [B]^{m} }}{{[C]^{p} \cdot [D]^{q} }} $$

(3)

and the equilibrium constant equals the mass action ratio at equilibrium. The thermodynamics is given by:

$$ - \Updelta G = R \cdot T \cdot \ln \left( {\frac{\Upgamma }{{K_{\text{eq}} }}} \right) = - \Updelta G_{0}^{\prime } + {\text{RT}}\ln \Upgamma $$

(4)

ΔG ₀′ is the standard Gibbs free energy difference of the reaction. This represents the influence of the chemistry on the direction into which the reactions wants to run. Whenever this Gibbs free energy difference of reaction is very negative, the reaction is effectively irreversible, and Eq. 2 reduces to its irreversible version:

$$ v = v_{ + } = k_{ + } \cdot [A]^{n} \cdot [B]^{m} $$

(5)

In the mass action kinetics of the simplest possible reaction, that is, the ‘irreversible’ unimolecular degradation of S, the rate (v) of a reaction is therefore dependent on the reaction rate constant (k) and the substrate concentration [S] (Hofmeyr 1995).

$$ v = k\left[ S \right] $$

(6)

Several additional assumptions have to hold for Eq. 6 to represent the rate of the corresponding enzyme catalysed irreversible reaction. These assumptions include either that the reaction is diffusion limited or that the substrate concentration is low enough for the enzyme not to be saturated by its substrate (Westerhoff and van Dam 1987; Cornish-Bowden 1995). In cases where the reaction, which is being modelled, is reversible, the same assumption continues to apply, but the backward reaction has to be included in the rate equation, e.g.:

$$ v = k_{s} \left[ S \right] - k_{p} \left[ P \right] $$

(7)

When the reaction is not limited by diffusion, that is, by the collision of the substrate with the enzyme molecule, a more representative form of the rate equation for the enzyme catalysed breakdown of S in the absence of any product takes into account the maximum rate of the enzyme (V _max). This rate of reaction tails off as the substrate concentration exceeds the K _M of the enzyme.

$$ v = \frac{{V_{\max } \cdot [S]}}{{K_{M} + [S]}} $$

(8)

This rate equation was derived by Michaelis and Menten (1913) assuming the opposite of diffusion limitation, that is, that the binding of substrate S to the enzymes is fast enough almost to reach equilibrium. Haldane (1930) relaxed this to the assumption that the enzyme-substrate concentration reaches steady state. Whenever [S] refers to the free concentration of the substrate, this equation is valid at steady state, even if the total enzyme concentration is high compared with the total substrate concentration. When the total enzyme concentration is small as compared with the total substrate concentrations, [S] in the above can be replaced by the total concentration of the substrate in the system, which much simplifies the further integration of the system. In cases where this does not apply, a more accurate zero derivative potential method may be applied (Härdin et al. 2009).

The reversible form of the Michaelis–Menten equation is (Westerhoff and van Dam 1987; Cornish-Bowden 1995; Haldane 1931):

$$ v = \frac{{V_{\max + } \cdot \frac{[S]}{{K_{{M_{S} }} }} - V_{\max + } \cdot \frac{[P]}{{K_{{M_{P} }} }}}}{{1 + \frac{[S]}{{K_{{M_{S} }} }} + \frac{[P]}{{K_{{M_{P} }} }}}} = \frac{{V_{\max + } \cdot \frac{[S]}{{K_{{M_{S} }} }} \cdot \left( {1 - \frac{{K_{\text{eq}} }}{\Upgamma }} \right)}}{{1 + \frac{[S]}{{K_{{M_{S} }} }} + \frac{[P]}{{K_{{M_{P} }} }}}} $$

(9)

When reactions are far enough from equilibrium, that is, ‘irreversible’, this reduces to:

$$ v = \frac{{V_{\max + } \cdot \frac{[S]}{{K_{{M_{S} }} }}}}{{1 + \frac{[S]}{{K_{{M_{S} }} }} + \frac{[P]}{{K_{{M_{P} }} }}}} $$

(10)

This equation is almost equal to the commonly used equation Eq. 8, but carries the product inhibition effect with it. This shows that independent of whether a reaction is irreversible or not, a modeller should always reckon with the likelihood that the product inhibits the forward reaction rate, simply because it necessarily binds to the same site of the enzyme as the substrate does and competes with the latter. It is good practice to have product-inhibited reactions in a model by default. Absence of the product inhibition terms imposes strong limitations to what controls the system. The distinction between reversibility and product inhibition is important.

There are several random order and ordered order binding equations commonly used in biology, which take the form of Michaelis–Menten equations but contain more substrates and products or can model reversible reactions (Cornish-Bowden 1995).

For many reactions, the rate of reaction is not only influenced by substrates, products, enzyme concentration and turnover, but also be allosteric effectors that can change the rate of turnover of the enzyme. It is often seen that the function of these allosteric interactions is not only to couple demand and supply but also to maintain a high independence between fluxes in competing pathways (Curien et al. 2009). When the activator binds to a site that does not affect binding of substrate, the rate equation changes by adding a prefactor (Cornish-Bowden 1995; Liebermeister and Klipp 2006). For an activator, this term depends on the concentration of the activator (A) and the activation constant (K _A ) and takes the form:

$$ {\text{Activation prefactor}} = \frac{[A]}{{K_{A} + [A]}} $$

(11)

For an inhibitor, this term is dependent on the concentration of the inhibitor (I) and the inhibition constant (K _I ) and would take the form:

$$ {\text{Inhibition prefactor }} = \frac{{K_{I} }}{{K_{I} + [I]}} $$

(12)

Modelling a dynamic system requires rate equations and balance equations. There is a balance equation for each metabolite. For a metabolite S, this reads:

$$ \frac{{{\text{d}}S}}{{{\text{d}}t}} = v_{{{\text{production}}\;{\text{of}}\;S}} - v_{{{\text{consumption}}\;{\text{of}}\;S}}, $$

(13)

that is, the increase in the concentration of metabolite S is equal to the rate of its production minus the rate of its destruction. In this equation, the reaction stoichiometries play a role, that is, if a certain reaction produces two molecules of S rather than 1 per turnover, then its term in the balance equation has to be multiplied by 2.

Kinetic parameters necessary for this modelling (K _I , K _A , V _max, K _Ms, K _Mp, etc.) have been measured in quite a few enzyme studies. They can be obtained from the literature or by searching in databases such as BRENDA (http://www.brenda-enzymes.org/) and by implementing other text mining tools. Alternatively, they can be measured by well-known biochemical techniques. Rates of reactions can also be fitted directly from experimental data. However, this will not be discussed in detail here.