Open Access

A Boolean probabilistic model of metabolic adaptation to oxygen in relation to iron homeostasis and oxidative stress

BMC Systems Biology20115:51

DOI: 10.1186/1752-0509-5-51

Received: 27 July 2010

Accepted: 13 April 2011

Published: 13 April 2011



In aerobically grown cells, iron homeostasis and oxidative stress are tightly linked processes implicated in a growing number of diseases. The deregulation of iron homeostasis due to gene defects or environmental stresses leads to a wide range of diseases with consequences for cellular metabolism that remain poorly understood. The modelling of iron homeostasis in relation to the main features of metabolism, energy production and oxidative stress may provide new clues to the ways in which changes in biological processes in a normal cell lead to disease.


Using a methodology based on probabilistic Boolean modelling, we constructed the first model of yeast iron homeostasis including oxygen-related reactions in the frame of central metabolism. The resulting model of 642 elements and 1007 reactions was validated by comparing simulations with a large body of experimental results (147 phenotypes and 11 metabolic flux experiments). We removed every gene, thus generating in silico mutants. The simulations of the different mutants gave rise to a remarkably accurate qualitative description of most of the experimental phenotype (overall consistency > 91.5%). A second validation involved analysing the anaerobiosis to aerobiosis transition. Therefore, we compared the simulations of our model with different levels of oxygen to experimental metabolic flux data. The simulations reproducted accurately ten out of the eleven metabolic fluxes. We show here that our probabilistic Boolean modelling strategy provides a useful description of the dynamics of a complex biological system. A clustering analysis of the simulations of all in silico mutations led to the identification of clear phenotypic profiles, thus providing new insights into some metabolic response to stress conditions. Finally, the model was also used to explore several new hypothesis in order to better understand some unexpected phenotypes in given mutants.


All these results show that this model, and the underlying modelling strategy, are powerful tools for improving our understanding of complex biological problems.


A large body of data suggests that mitochondrial abnormalities may link gene defects and/or environmental challenges to many pathologies including several neurodegenerative processes (for reviews, see [14]). Mitochondria are essential organelles serving as the main site of oxygen use within cells. The divalent reduction of oxygen by the respiratory chain is tightly coupled to ATP synthesis by the oxidative phosphorylation machinery. However, a small proportion of the electrons passing through the electron transport chain reacts with molecular oxygen in a monovalent reduction reaction [5]. This process yields the superoxide anion, which can be converted into other reactive oxygen species (ROS), such as hydrogen peroxide and the highly reactive hydroxyl radical, through enzymatic and non-enzymatic reactions [6]. Cells possess an impressive arsenal of weapons for counteracting excess ROS production, including superoxide dismutases, catalases, peroxidases and low molecular mass redox compounds, such as ascorbic acid and glutathione. However, overproduction of the superoxide anion due to the abnormal reduction of key components of the respiratory chain (e.g. ubiquinone and cytochrome bc1) or to the impairment of antioxidant defences adversely affects various cellular processes and constituents (for recent reviews see [79]). Disruptions of the respiratory chain or cellular defences are thus increasingly being implicated in acquired and inherited diseases and appear to play a key role in the aetiology of many neurodegenerative disorders, including Alzheimer's and Parkinson's diseases [10, 11].

Due to its unique redox properties and chemical reactivity, iron appears to be a key player in abnormal ROS generation, principally as a catalyst of the Fenton and Haber-Weiss reactions [12]. This essential micronutrient is the redox component of the haem and iron-sulphur cluster [FeS] cofactors of many important proteins or enzymes. Iron homeostasis is thus tightly regulated, at all levels. The deregulation of iron homeostasis due to gene defects or environmental stresses leads to a wide range of diseases, from anaemia (iron deficiency) to haemochromatosis (iron overload) [1315] with consequences for cellular metabolism that remain poorly understood. The modelling of iron homeostasis in relation to the main features of metabolism, energy production and oxidative stress may provide new clues to the ways in which changes in biological processes in a normal cell lead to disease.

In the growing field of systems biology, several attempts have been made to model cellular processes. For complex systems, these models can be classified into static [16, 17] and dynamic models [1820]. Dynamic models are generally analyzed and/or simulated with specific methods as a function of the availability (or lack) of numerical data. The approaches used include ordinary differential equations, stochastic simulation algorithms [2123], agent or rule-based approaches [24, 25], Boolean approaches (probabilistic or otherwise) [2629], or Petri nets (with hybrid extensions) [3032]. Most of these models are based on well known pathways (cell cycle control, glycolysis, signal transduction, etc.), or well studied processes relevant to cell physiology (e.g. action potential propagation in neurons).

Our understanding of the relations between oxidative metabolism and iron homeostasis is based on a large body of qualitative knowledge from heterogeneous sources, often lacking numerical data. It is therefore not possible to derive mathematical relationships based on biological knowledge for the entire system, and the model has to include uncertain knowledge. As a consequence, despite several attempts to construct models accurately describing certain aspects of iron homeostasis (at the level of the cell [33] or organism [30, 34, 35]), no formal model linking iron homeostasis to metabolism control and oxidative stress has yet been developed.

We describe here an approach based on probabilistic Boolean modelling that can deal effectively with vast amounts of heterogeneous knowledge not always associated with quantitative data. Using this approach, we were able to construct a realistic model of cell fate including oxygen, carbon, nitrogen, sulphur, phosphate and iron homeostasis. This methodology dealt well with the mixture of precise and uncertain knowledge. Despite its large size (642 elements and 1007 reactions), we were able to simulate the model and analyse its dynamics. We focused on a simple unicellular eukaryotic system, the yeast Saccharomyces cerevisiae, which has many features (genes, proteins, pathways, cell compartmentalisation) similar to those of mammalian cells. We first validated the model by simulating 198 in silico mutations resulting from the deletion of individual genes from the model (with a small number of elements kept constant in the model). An independent validation was provided by analysing the key transition from anaerobic to aerobic metabolism, by comparing in silico reactions frequencies with experimental fluxomic data [36]. The results of individual deletions were compared with experimental data for real mutants, for which detailed phenotypic analyses were available. We were able to classify the in silico mutants into groups of similar "phenotype" profiles, making it possible to identify original properties of the model. The model was used to explore several alternative hypothesis in order to better understand some unexpected phenotypes in mutants.

In this study, we focused on iron homeostasis and oxidative stress. However, we believe that the proposed modelling strategy could be useful for other systems. Typically, it should allow the building of large models with a high level of biological relevance that can cope with both a lack of numerical data and precise knowledge.

Results and Discussion

The model

Choice of the type of model

In eukaryotic cells grown under aerobic conditions, energy production may depend strongly on oxidative phosphorylation. The efficiency of this process is dependent on tight coupling between the production of reduced equivalents (mainly NADH and FADH2) through interconnected metabolic reactions and their reoxidation by the complexes of the respiratory chain. The reoxidation steps generate the electrochemical gradient required to drive ATP synthesis. The assembly of the respiratory complexes requires the coordinated synthesis of the protein components and their appropriate cofactors, mainly haem or iron-sulphur clusters. Any attempt to model these processes and their relationships must therefore include a description of the main metabolic pathways, ion transport and assimilation, the expression and regulation of genes, protein complex assembly and cofactor biosynthesis. These different levels of description involve several cellular compartments, such as the cytosol, mitochondria and nucleus. Yeast is a particularly suitable cell model for studies of this type, as most of its biochemical pathways have been characterised, from enzymatic reactions to gene regulation. However, some processes, such as the assembly of iron-sulphur clusters, still contain reactions about which we know little (for reviews, see [37, 38]). Thus, any model developed to analyse the relationship between iron homeostasis and oxidative stress damage must be able to deal with several different levels of knowledge about a large number of elements (i.e. biological species: proteins, genes, ions, metabolites...) and reactions. It should also include different levels of regulation in different cellular compartments. Furthermore, we wanted to model cellular processes as closely to the molecular mechanisms as possible. For example, if a protein is absent in some specific condition we could have modelled this fact using a negative statement. Instead, we would rather model the underlying molecular mechanism that are easily transcribed using only positive statements. In our example, it means describing why this protein is absent: is the protein degraded in this specific condition? Is the promoter of the corresponding gene occupied by some transcription factor? Is the mRNA degraded?

Our modelling approach is derived from the Biocham language [39], a rule-based formalism. Each biological species (elements) in the model is either "ON" (present) or "OFF" (absent). The elements interact through biological reactions, formalised as rules. We have extended the Biocham approach by adding a weighting to each reaction, making it possible to define a probability for each reaction to simulate the large differences encountered between certain reaction rates (see Methods). This Boolean probabilist approach is designed to mimic the biological reality as closely as possible, given the qualitative nature of the model.

The list of elements and reactions

The model was built from a core of reactions involved in iron homeostasis. We included in the model the main biochemical pathways, carbon, nitrogen, sulphur and phosphate metabolism, taking into account the cellular distribution of the components and the reactions (see Methods and Figure 1). The connectivity between reactions is depicted in Figure 2. The extracellular space contains the nutrients imported into the cell when the appropriate permeases or uptake systems are synthesized. The nutrients include glucose as a source of carbon. Its uptake may be modified by assigning different weights to the import reaction. When transcribed, the nuclear genes produce proteins that are exported to the cytoplasm and then targeted to their final subcellular compartment (transcription and translation are modelled as one reaction because no alteration of the mRNAs, for example active degradations, was introduced). Mitochondria perform the main energy-related metabolic reactions. Substrates and products shuttle between cytoplasm and mitochondria through appropriate translocators. A limited number of functions were attributed to the vacuole, mostly relating to the homeostasis of metals and phosphate. Our model currently includes 642 elements and 1007 reactions (see Table 1 and Additional file 1 for the detailed list of reactions). The weights of the reactions were set to the default weight of 1 unless changing the weight of a reaction was necessary to model a phenomenon realistically (see Methods for more details and Table 2).
Figure 1

Overview of the content of the model. Main pathways included in the model and their cellular localisation.
Figure 2

Graph of the model elements. Each element is a node, two elements are connected if they are involved in the same reaction. Blue nodes are constant elements, red nodes are non-constant elements. Green lines are the connections involving oxygen. The node size increases with the number of elements connected to it.

Table 1

Number of elements and reactions of the model

642 Elements:





Proteins and protein complexes


Other chemical components and metabolites


Unlimited elements (always ON)

1007 Reactions:



Gene transcription/translation


Translocation of elements between compartments


Protein complex synthesis/cofactor binding


Sulphur metabolism/glutathione biosynthesis


Haem/siroheme/FeS/PLP synthesis


ATP/GTP synthesis


Oxidative stress


Element degradation


NADP/NADPH synthesis


Allantoin (NH3)

Table 2

List of the reactions with weights other than the default value of one.



Source and/or comment

Basal transcription/translation of condition or transcription factor dependent genes


Should be significantly lower than regulated expression

Active Degradation


Should be significantly higher than the production

Aft1/2p::nucleus = [Glt1/Aco1p-FeS::c, Grx3p::nucleus] Aft1/2p::c



Low affinity Fet4p metal transport



Glutathione/Cysteine mitochondrial import


Balance of flux between compartments in wt

Neoglucogenesis enzymatic reactions



Reaction catalysed by Zwf1p


Major source of NADPH. [76]

TCA cycle (reverse)



Atp production by Atp-synthase



O2m::c = [Sod1p-Cu-Zn::c] H2O2::c



O2m::m = [Sod2p::m] H2O2::m



O2m::(m or c) + H2O2::(m or c) Ot::(m or c)



O2m::(m or c) + H::(m or c) H2O2::(m or c)



Ot::(m or c) + Aco1p-FeS::(m or c) Aco1p::(m or c)


Less probable than aconitase reaction but more probable than spontaneous degradation.

O2m::(m or c) + Aco1p-FeS::(m or c Aco1p::(m or c)


Same as "Ot::(m or c) + Aco1p-FeS::(m or c)".

Ot::(m or c) + any iron-sulphur cluster containing protein protein without the iron-sulphur cluster


Less probable than the same reaction involving Aco1p-FeS

O2m::(m or c) + any iron-sulphur cluster containing protein protein without the iron-sulphur cluster


Less probable than the same reaction involving Aco1p-FeS

Ot::c + Met5p-siroheme-FeS::c Met5p-siroheme::c


Less probable than other FeS proteinsdue to presence of siroheme

O2m::c + Met5p-siroheme-FeS::c Met5p-siroheme::c


Less probable than other FeS proteins due to presence of siroheme

H2O2::c = [Ctt1p-Heme::c] O2


Higher turnover number than most enzymes but less than SODs [80]

H2O2::m = [Cta1p-Heme::m] O2


Higher turnover number than most enzymes but less than SODs [80]

#O2::ext #O2::ext + O2


Provides saturating oxygen.

#glucose::ext #glucose::ext + glucose::c


Limited glucose uptake as described in [36]

Degradation of protein or peptide elements


Must be significantly less than the default.

Degradation of other elements


Less than proteins degradation.

m or c = mitochondria or cytoplasm.


The output of one simulation is "ON" (present) or "OFF" (absent) for each 642 element at each iteration (simulations of 20 millions of iterations or until a steady state was reached, See Simulations in the Methods section for more details). In order to describe the average evolution of each element, each model was simulated 100 times. We then defined, for each element, a parameter called PoP (Percentage of Presence). The PoP is the number of iterations where the element is "ON" during the last million of iterations, averaged over 100 simulations. For example, a PoP for the element NADPH::cytoplasm of 24.3% means that, on average, during the last millions of iterations, this element is "ON" for 24.3% of the iterations. Recall that the model is a Boolean one, and one should not try to interpret or relate these values to quantities such as concentration or number of molecules. However, because an element can be involved into several reactions (and consequently consumed when one of these reactions is triggered), this parameter reveal the "availability" and the "usability" of the element during the time evolution of the system. Hence, it is a simple parameter that convey important biological information.

We also defined a reference model, referred to as the wild type model (WT model) which correspond to a model where all genes are "ON" during the simulations. In this WT model, our list of reactions and their weights allow the PoPs to be: 1/biologically meaningful according the data from the literature and, 2/high enough so that no reaction is always impossible during WT simulations. However some elements, such as the ROS, should remain low in the WT model (PoP less than 1%) as observed in wild type cells. A higher PoP for these elements would be interpreted as related to a situation of oxidative stress. To reach this goal, in some specific cases, the default weights should be altered. As an example, when every weights in the WT model are set to 1, the PoP of the hydroxyl radical (the most deleterious ROS) is higher than 70%. It means that this element would be too often encountered in the system which is obviously biologically irrelevant. In our WT model, this element has a PoP lower than 1% due to appropriate weights on the reactions involved in handling ROS (see Table 2 for a list).

In the following, we used the PoP to compare the outputs of one model (for example when simulating an in silico mutation) versus the reference model (WT model).

The distribution of the PoP for the various elements (Figure 3) indicated that the vast majority of elements remained present at steady-state levels from a few percent (e.g. toxic forms of oxygen) to 100% (e.g. all the genes). The PoP of most proteins was relatively high (higher than 50%), given their low probability of degradation. Eighty proteins are active in a specific cellular compartment other than the cytoplasm. Since they are produced in the cytoplasm and then translocated to their compartment, the overall PoP for their residence in the cytoplasm is about 66%. This causes the peak around 60% in Figure 3.
Figure 3

Histogram of the PoPs of the elements in the WT model at steady state.

We performed a sensitivity analysis of the outputs of the model (PoP at steady state) when the weights of the reactions were modified by multiplications or divisions up to a factor of 10 (see Methods and Additional file 2 for more details). Varying the weight of a given reaction had no impact on the model when the variation was small (weight +/- 10%). For some reactions, when the variations were larger, the PoP of few elements (up to 20 out of 642) was modified. Each of these variations has only a local impact on the graph (as defined in Figure 2). This analysis shows that the model is robust to modifications of the weights of the reactions within a dynamic range of 100, that reflects most biological situations.

Validation of the model

Our validation strategy involved evaluating the consistency of simulations of the model at different levels. We first considered the effects of removing every gene, thus generating an in silico mutant. The simulations of the mutants were compared to experimental phenotypes of the corresponding biological mutants. A second validation involved analysing the anaerobiosis to aerobiosis transition. Therefore, we compared the simulations of our model with and without oxygen to experimental metabolic flux data.

Evolution of representative elements in selected mutants

For validation of the model, we compared the results of the simulations for the WT model (all elements present at the beginning of the simulations) with those for models in which specific changes were introduced.

We analyzed the effects of removing selected genes, thus generating in silico mutants. The resulting mutated genes have the same designations as the corresponding biological mutants, but with an "is-" prefix (e.g. is-hem15 represents the model with the HEM15 gene switched "OFF" during the simulations - see Methods for more details).

We show the evolution of the PoP during the simulations for selected elements in representative mutants in Figure 4. Figure 4-A describes a simple situation in which deletion of the gene encoding ferrochelatase (is-hem15) leads to an increase in the PoP of protoporphyrin IX, the substrate of the enzyme, paralleling the situation occurring in vivo in the corresponding biological mutant [40]. Similarly, in a model in which superoxide dismutase was eliminated (is-sod1, Figure 4-B), the PoP of its substrate, the superoxide anion [41] increased strongly. As a consequence, the PoP of the hydroxyl radical also increased, due to the Fenton and Haber-Weiss reactions.
Figure 4

Mean PoP of a number of elements in selected mutants. The gene was turned o after the first million of steps. A: PoP of protoporphyrin-IX when the gene HEM15 is turned o, B: superoxide anion and hydroxyl radical when the gene SOD1 is set to "OFF", C-F: PoP of glutamate, cysteine and glutathion C) when ACO1 is set to "OFF". Cysteine reached zero at steady-state, D: ACO1 is set to "OFF" and a source of glutamate is added, E: HEM1 is set to "OFF", F: HEM1 is set to "OFF" and a source of glutamate is added.

Figure 4-C-F illustrates two more complex biological situations. Yeast mutants lacking the tricarboxylic acid cycle enzyme aconitase are known to be auxotrophic for glutamate [42], due to a lack of alpha-ketoglutarate and impaired nitrogen assimilation. The is-aco1 mutant (Figure 4-C) displayed a decrease in the PoP of glutamate and, as a consequence, a decrease in the PoP of glutathione (GSH), the synthesis of which requires glutamate (GSH is a γ-glutamyl-cysteinyl-glycine tripeptide). However, the PoP of cysteine was not affected in the is-aco1 mutant. Indeed, it may even have increased slightly, as cysteine was no longer used for GSH synthesis. If we "supplemented" our model with glutamate (Figure 4-D), fixing this element as "always present", then the is-aco1 mutant no longer displayed defective GSH synthesis. Similarly, mutants affected in the early steps of the haem biosynthesis pathway are known to be auxotrophic for cysteine/methionine, due to a defect in the synthesis of siroheme, the cofactor of sulphite reductase [43]. Therefore they display impaired sulphur assimilation. An analysis of simulations for the is-hem1 mutant (Figure 4-E) showed a decrease in the PoP of cysteine when hem1 was removed from the model. Accordingly, the PoP of GSH also decreased if siroheme synthesis was impaired. More surprisingly, the model predicted a decreased in the PoP of glutamate for the is-hem1 mutant. The provision of unlimited amounts of glutamate did not restore cysteine and glutathione synthesis in the is-hem1 model (Figure 4-F). Furthermore, glutamate synthesis was restored by providing cysteine as an unlimited element in the is-hem1 model. The decrease in the PoP of glutamate in the is-hem1 model is therefore a consequence of the decrease in the amount of cysteine. Indeed, Figure 4-E shows that the PoP of cysteine decreases before the PoP of glutamate. (for more analysis of phenotypes, see "Systematic analysis of predicted mutant phenotypes" below). The "phenotypes" of such mutants were found to be consistent with published data for the phenotypes of the corresponding biological mutants, providing the first demonstration that our model was of high quality.

Systematic analysis of predicted mutant phenotypes

We simulated the model in which each element (gene or element supplied as unlimited) was removed individually from the model (100 simulations for each of the 198 elements removed). Whenever possible, we compared the is-phenotypes to the manually curated phenotypes reported in the Saccharomyces Genome Database (SGD). We extracted from the simulations and from SGD, data corresponding to documented phenotypes related to "auxotrophies", "chemical compound accumulation", "oxidative stress resistance", "respiratory growth" and "nitrogen source utilization" (see Additional file 3). The simulations of the different is-mutants gave rise to a remarkably accurate qualitative description of most of the experimental phenotype. There are some unavoidable discrepancies that may be traced back to different problems.

The "auxotrophies" phenotypes were found consistent in 18/21 simulated mutants. We did not evidence a methionine auxotrophy in the is-zwf1 mutant, because in the corresponding model, the decrease in PoP of NADPH is not important enough (PoP of NADPH::cytoplasm: WT = 24.3+/-0.3%, is-zwf1 = 17.7+/-0.3%). NADPH remains present due to the NADH kinase activity and a constant supply in NAD in the model.

We analysed the sensitivity of the is-mutants to oxidative stress (SGD "oxidative stress resistance"). Under our simulation conditions, the mutants were not challenged with large amounts of exogenous hydrogen peroxide, paraquat or diamide as oxidizing agents (all these conditions are usually required to evidence in vivo sensitivity to oxidative stress). Therefore, the simulations are not expected to show any variation in the PoP of the ROS. Indeed, several is-mutant do not accumulate ROS, while the corresponding biological mutants were reported as sensitive to oxidative stress. However, under our simulation conditions many mutants (13/26 in silico mutants) showed an increase in the PoP of the ROS. This means that our model allows to point out the mutations the more prone to induce indogenous oxidative stress.

The "respiratory growth" phenotype reported in SGD is most often strictly related to the capability of the cells to grow on non-fermentative sources of carbon (glycerol, ethanol or lactate). All our simulations were run using glucose as the sole source of carbon. We should therefore not be able to evaluate the "respiratory growth" phenotype from our simulations. However, we show here that using the variations of the PoP of protons in the intermembrane space of mitochondria is indeed a very good is-marker of respiration defect. As much as 63/92 mutants descriptions could be validated this way, while 27/92 descriptions remained similar in is-mutants and is-WT and only 1/92 discrepancy was found between model and biology.

The "chemical compound accumulation" phenotype was very discriminating, with 33/45 descriptions in full agreement with the biology, 3/45 descriptions showing no change between is-mutant and is-WT. Our model did not allow to evidence accumulation of zinc in a is-zrt1 mutant, but this phenotype reported in SGD was conditional to the presence of aluminum ions in the growth medium of the biological mutant [44]. Some of the is-adk1 mutant phenotypes did not fit to the published ones, mostly because in this is-mutant, it is diffcult to quantitatively describe the imbalance between fermentation and respiration. The biological mutant exhibits apparently contradictory phenotypes, with increased ethanol accumulation and increased respiratory growth. Interestingly, our model predicts that the is-yfh1 mutant is depleted in oxidized glutathione, while SGD reports increased GSSG. However, recent results show indeed that the prediction from the model is correct [45] and this impacts on NADPH content due to higher G6PDH (Zwf1p) activity. All the mutants defective in "nitrogen sources utilization" were found in the simulations (5/5).

On a total number of 145 situations where the is-phenotypes were comparable to the SGD data, 91.7% were consistent (see Additional file 3). These results indicated that through a large scale analysis of the phenotypes expected from the biology, we were able to recover the corresponding is-phenotypes.

Analysis of the anaerobiosis to aerobiosis transition, or "how to deal with dynamic processes?"

We then analyzed the changes in several independent elements in the simulations of physiologically relevant transitions. One of our goals was to analyse the oxidative stress response. We therefore modeled the transition from anaerobiosis to aerobiosis. We analyzed simulations in which a key element, oxygen, was removed from the model. Oxygen is required to re-oxidise the reduced equivalents produced by cell metabolism. However, S. cerevisiae can grow in the absence of oxygen, due to its ability to shift from a respiratory to a fermentative metabolism. NADH is then oxidised by alcohol dehydrogenases, leading to ethanol production. We therefore compared the simulations of our model with different levels of oxygen to the experimental data of Jouhten et al.[36], who measured metabolic fluxes in yeast cells grown in a chemostat with a limited glucose supply in the presence of 20.9% oxygen (aerobic conditions) and under anaerobic conditions (0% oxygen). These experimental conditions were taken into account by weighting "glucose import" to a low value of 0.1. We compared the experimental fluxes with the reaction frequencies at steady state of key carbon metabolism reactions, in both the presence and absence of oxygen (Figure 5). All the in silico reactions displayed patterns of variation highly similar to those observed in vivo, with the exception of the reaction catalysed by Oac1p, a mitochondrial bidirectional oxaloacetate transporter with broad specificity for various anions. There is a simple explanation for the discrepancy observed for this specific reaction (of 10 considered). In anaerobiosis, cytoplasmic oxaloacetate is less produced than in aerobiosis. However, it is also less consumed than produced. Because the ratio of the production over the consumption is higher in anaerobiosis than in aerobiosis, oxaloacetate PoP is higher. Therefore, since one of the few remaining possible reactions involving oxaloacetate is transports by OAC1, the model predicts an increase in this transport in the absence of oxygen. Another important consideration when analysing results for the N2-O2 transition is that the model includes a number of regulatory mechanisms to describe the biological response of the yeast cell to oxygen deprivation accurately. These regulations were of two kinds. First, we modelled the induction, under anaerobic conditions, of the Rox1 responsive regulon [46, 47]. Second, we modelled the effects of glucose repression on the alcohol dehydrogenase system, by adjusting the rate of degradation of the corresponding enzymes as a function of the presence or absence of oxygen [48]. This reflects the post-transcriptional regulation mechanisms involved in these processes.
Figure 5

Variations of the fluxes in selected reactions with and without oxygen. Occurrence of 10 reactions (WT model) in the simulations, compared with experimental fluxes [36] (red). A-I: frequency of the reaction at steady state (green), J: PoP of ethanol at steady state (green).

This global validation of the model provides important evidence that the frequencies of utilization of the reactions during the simulations are qualitatively consistent with the metabolic fluxes observed experimentally.

We show here that our probabilistic Boolean modelling strategy may provide a useful description of the dynamics of biological systems. Further analysis is now required to determine the extent to which this key notion of fluxes, estimated as reaction frequencies at steady-state, can be generalised.

Global analysis of the output of the simulations

We evaluated the potential of the model to identify emerging properties through a global analysis of the simulations of the in silico mutants. In silico mutations were performed on 187 genes of the model. Indeed, it is impossible to predict the outputs of such a large number of simulations. We also included in this meta-analysis the simulations of models in which each of the 11 fixed elements was removed, independently, representing changes in the composition of the growth medium We evaluated whether it was possible to uncover from the simulations of the different models some specific properties of this complex system. A Bayesian classification algorithm (implemented at Autoclass@IJM, [49]) was then used to cluster the outputs of all simulations, defining classes of is-mutants with similar phenotypes and classes of elements varying in a similar manner in the simulations. Genes and elements were clustered into discrete sets (Figure 6 and Table 3). These mutations were grouped into seven classes that were consistent to some extent with the known biology of the system used (hem mutants, met mutants, cys mutants) and with additional mutations in genes with no straightforward relationship to the other members of the classes. This suggests that these mutations may define previously masked metabolic states, some of which it may be impossible to test experimentally, as several elements may be simultaneously missing, leading to cell death in vivo. A detailed analysis of each class showed that most simulations were correlated with biological data, but it was not always possible to predict the class to which a given gene would be assigned. Many is-mutations had only a limited impact on the steady-state PoP of most of the elements of the model when simulated with the predefined set of unlimited elements (Additional file 4). However, several classes of is-mutations had very pronounced phenotypes.
Figure 6

Variations of the PoP, at steady state, of each element in the model (columns), in each in silico mutant (rows), with respect to WT simulations. For each mutant and for each element, was calculated. Positive values (red) indicate that the PoP was higher in the mutant simulations and negative values (green) imply that the PoP was higher in the WT simulations. For a complete content of all clusters see additional le 2 (cdt file of this clustering). See methods for details.

Table 3

Composition of the clusters of mutants from Fig. 6 (rows).

Cluster 1


Cluster 2


Cluster 3


Cluster 4


Cluster 5


Cluster 6


Cluster 7


Cluster 8


Cluster 9


Cluster 10


Cluster 11


Cluster 12


Cluster 13


Each mutant is referred to by its gene name or the constant element that was turned "OFF". COX1A13 = genes COX1 to COX13; ATP1A5 = genes ATP1 to ATP5; ATP710A23 = genes ATP7 and ATP10 to ATP23; ATP689 = genes ATP6, ATP8 and ATP; LCCA = long chain carboxylic acid; ZN2 = Zn2+; FE3 = extracellular Fe3+; COA = coenzyme A.

One of these phenotypes was an increase in the PoP of toxic oxygen species. Two clusters of is-mutants displayed a large increase in the PoP of the hydroxyl radical. The first cluster (class 2) had defective sulphate assimilation, sulphur amino-acid synthesis and iron-sulphur cluster biogenesis or assembly. Among them were the sulphate permease (is-sul1), enzymes of the methionine/cysteine biosynthesis pathways (is-met1, 2, 3, 4, is-met8, 14, 16, 25) and enzymes involved in the mitochondrial assembly of iron-sulphur clusters (is-isu1, is-atm1 and is-yfh1, which is known to be extremely sensitive to oxidative stress). The cluster was particularly interesting, because it also grouped together is-mutations in genes and the depletion of elements involved in controlling major metabolic pathways (is-sam1, homoserine, CH3): CH3, representing the entry point to one-carbon metabolism, S-adenosyl methionine production (SAM1), and alpha-ketoglutarate dehydrogenase (KGD1 and KGD2), which connects carbon and nitrogen metabolism in yeast.

The second cluster (class 6) included is-mutants in which the early steps of haem synthesis (HEM1 to 4) were impaired, with phenotypes similar to wild-type models run in the absence of oxygen (except, of course, for the accumulation of ROS). A similar phenotype was observed for is-mutants with defective pyridoxal phosphate cofactor synthesis (PLP is required for ALA synthase activity). Interestingly, this class of is-mutants also included the WT minus Zn conditions, consistent with the requirement of Zn as a cofactor for ALA dehydratase (HEM2). The is-adk1 mutant, lacking adenylate kinase, also belonged to this class. All these is-mutants had a high PoP of toxic forms of oxygen. Most sulphur-containing compounds had a low PoP, due to the inability of these is-mutants to assimilate sulphate.

These two classes of mutations were strongly correlated with an iron-related phenotype, the accumulation in the nucleus of the iron-responsive transcription factor Aft1p, thus mimicking the "AFT1-up" phenotype described for biological Δyfh1 [50, 51], Δatm1 [50] and Δcbf1 [52] mutants (Cbf1 is a global transcriptional regulator of the cys/met regulon [53]).

Is-mutants affected in the late steps of haem synthesis (HEM12 to 15, class 3) had very different phenotypes. They lacked cytochromes and displayed an increase in the PoP of hydrogen peroxide due to the lack of catalase, but their sulphur metabolism was unaffected.

These results show that our model is robust to in silico mutations, in that it predicts coherent phenotypes. To our knowledge, this is the first classification of such a large number of in silico mutations by a single model generating coherent and biologically relevant results.

Focusing on changes in the PoP of a given element makes it possible to explore multiple alterations to gene expression, such as loss of function (gene always "OFF") or overexpression (setting a higher weight). In this context, we might speculate as to which mutations might restore the wild-type level of a missing element in a given mutant, for example. Approaches of this type should help to target experimental studies trying to identify functional interactions between genes, but should also facilitate the exploration of new hypotheses relating to certain poorly documented aspects of the control of cellular homeostasis.

Exploration of alternative hypothesis: an example of a use of the model

Although iron homeostasis in yeast is fairly well described, some aspects of its metabolism remain unclear. As an example, it has been reported that some mutant strains affected in the iron-sulphur biosynthesis pathway accumulate iron within mitochondria as iron-phosphate amorphous precipitates or nanoparticules ([54] (Δyfh1), [55] (Δatm1), [56] (Δyah1)). Iron-phosphate complexes (FePi) are virtually insoluble in aqueous media at physiological pH (Ksp 9.91 × 10-16) [57] and iron is therefore not available for biosynthesis. There should be some mechanism in wild-type cells that either prevent the formation of FePi complexes or allows the re-mobilization of iron for haem or iron-sulphur biosynthesis.

We used our model to explore alternative hypothesis related to FePi accumulation within wild-type and iron suphur deficient cells (the Δyfh1 mutant as an example). We first postulated the presence within cells of a highly efficient reaction that reverts FePi aggregates to free Fe and Pi. This reaction is catalyzed in the model by an unknown element "X" and the weight of its reaction is set to 1000. As shown in Figure 7-C, removing this reaction from the model led to a dramatic increase in the PoP of mitochondrial FePi, a decrease of both mitochondrial Fe and Pi and a steadily increasing oxidative stress concomitant to the loss of iron-sulphur clusters.
Figure 7

Variation of selected elements PoP under different hyptothesis regarding the accumulation of iron-phosphate aggregates. A: Initial Model - gene YFH1 set to "OFF", B: Hypothesis 1 - gene YFH1 set to "OFF", C: Hypothesis 1 - Element X set to "OFF". Values at steady-state: Oxydized Mir1p = 90%, non oxidized Mir1p = 20%, hydroxyl radical = 99%, FeS = 1%, D: Hypothesis 2 - gene YFH1 set to "OFF".

We wanted however to evaluate the evolution of the model in situations where the unknown element "X" was present, using different hypothesis of physiological relevance. The first hypothesis takes into account the well described alteration of the function of the major Pi transporter, Mir1p, when a single cysteine residue is modified (see Additional file 5). Mir1p shifts from an anti-port (Pi/Pi or Pi/OH-) activity to a symport (Pi/H+) activity [58]. We therefore simulated this change of activity when Mir1p is the target of an oxidative modification under stress conditions. Figure 7-B shows that while the wild-type cells are unaffected by this modification, a yfh1 mutant, known to be very sensitive to oxidative stress and to accumulate FePi, indeed exhibits a slight but significative increase in Fe, Pi and FePi PoP within mitochondria, while ROS PoP strongly increase and FeS dramatically drops.

However, in the model, Mir1p is not the only transporter involved in the exchange of Pi between the cytoplasm and the mitochondria. Sal1p (ATP-ADP/Pi antiporter) and Dic1p (dicarboxylic acids/Pi antiporter) also contribute to the Pi homeostasis. We therefore simulated a model where all the three transporters may have their activities altered by an oxidative stress. The phenotypes were similar to those observed under the previous hypothesis but the PoP of FePi was only slightly higher (data not shown), indicating a modest cumulative effect of the oxidative alterations of the three transporters.

FePi nanoparticles are dissociated in vitro when treated with hydrosul fite, a strong reductant [56]. We evaluated the evolution of the model when "X" was identified as a thiol as a potential reducing agent. In our model, this thiol may be either cysteine or reduced glutathione. Removing cysteine impacts dramatically on the model and precludes the evaluation of its role in FePi formation. Therefore, we tested the hypothesis that "X" may be glutathione. Figure 7-D shows that under this hypothesis, most of the phenotypes of the biological yfh1 mutant are appropriatly described in the model since in this mutant, the glutathione PoP is 0% at steady state: there is a strong increase of the PoP of FePi in addition to the previously observed high PoP of ROS, loss of iron-sulphur clusters and drop in mitochondrial glutathione. This example shows how the model may be used to investigate poorly understood aspects of the cell biology. Differents hypothesis were explored to evaluate the mechanism of the accumulation of iron-phosphate within mitochondria. The simulations indicated that, although perturbation of phosphate transport accross the membrane could be considered, the hypothesis of a reductive dissociation of FePi by thiols is more likely. Indeed, many of the phenotypes of the mutant simulated are adequately rendered in this latter case. We are currently trying to test this hypothesis by using combination of biological mutants such as the Δgsh1. However, this study is made difficult since a Δgsh1 strain was shown to have limited growth without exogeneous glutathione. The lack of glutathione triggers an apoptotic response after few cell divisions [59]. The Δgsh1 strain accumulated mitochondrial iron but the activity of some FeS enzymes within mitochondria remained normal (the cytosolic FeS pool was deeply depleted) [60]. Interestingly, the double mutant Δgsh1 Δatm1 was described as non-viable.

The example of the investigation of the mechanisms of accumulation of FePi within mitochondria in pathological situations illustrates the power of this approach. Several hypotheses are easily explored by introducing the appropriate reactions and elements into the initial model. Testing different hypotheses paves the way to different type of experiments such as genetic/biochemical identification of the "X" element that makes iron soluble even in the presence of phosphate, with special emphasis on glutathione or other thiol containing compounds status within the cells. It will be also necessary to measure the activity of the various phosphate exchangers in mitochondria purified from different mutant strains, or from mitochondria subjected to strong oxidative stress. We believe these applications of the model will contribute to a better understanding on what contribute most to the control of iron homeostasis. This approach is to be extended to the analysis of other pathways of important physiological relevance, in normal or pathological cells.

Application of this methodology to other systems

The new extension of the Biocham language introduced and validated here could potentially be extended to many other studies. Systems biology papers frequently contain "interaction maps" for cellular or biochemical processes, as defined in a number of studies (see for example [16, 61]) based on the formalism of Molecular Interaction Maps http://​discover.​nci.​nih.​gov/​mim. Several authors have suggested that their interaction maps could serve as a first step towards the development of a computational model and have highlighted the need for appropriate tools to translate these maps into models [62]. Our approach appears highly suitable for this purpose. Each edge of an interaction map could be translated into a rule, as used in this paper. The weights of rules may be defined as a function of biological knowledge, as shown here. This last step seems to be manageable, because these maps are generated by qualified biologists and are based on extensive mining of published texts. The definition of weights simply provides a means of translating a biologist's intuition and knowledge into a numerical value; this value does not need to be obtained through experimentation (although it could be) and requires definition only with respect to the other rules. Our sensitivity analysis shows that only the orders of magnitude for the weights are relevant. Biologists should then simulate and analyse their models. We believe that the analysis of simulations will prove to be a powerful tool for the analysis of large biological systems [63]).

Such a qualitative approach may be seen as a first step towards more detailed mathematical models. Indeed, no quantitative value is required and only a rough description for each interaction is needed. Despite providing only an approximation of the interaction maps, we believe that this method can provide biologists with meaningful information.


Iron homeostasis, oxidative stress and metabolism are tightly linked within cells, with far-reaching implications at the level of the whole organisms. Indeed, disturbances in these processes are increasingly frequently being linked with disease. However, it remains unclear how the disruption of part of this system affects the rest of the network. Previous studies have modelled iron-related processes at the subcellular level (iron-sulphur cluster assembly in yeast mitochondria [64, 65]), at the cellular level in E. coli[33] (modelling of the genetic regulation of iron-related genes) or at the organism level in humans (regulation of iron levels; [30, 34, 35, 62]. Other models dealing with oxidative stress-related systems have been published [66, 67] but, to our knowledge, this is the first large model to include iron homeostasis as well as oxidative stress and some aspects of metabolism.

Using a Boolean probabilistic methodology based on the rule-based approach proposed by Biocham, we have constructed a detailed model of this complex system in yeast. Our methodology uses:

  • Only positive statements (no OR or NOT) to reflect the molecular nature of biological interactions,

  • Weighted reactions to reproduce large differences in certain reaction rates,

  • Stochastic simulations, to reproduce the variability of biological systems.

This new methodology combines the adaptability of Boolean models (no need for quantitative data for every element or reaction, possibility of including "fuzzy knowledge" etc.) with the need to reflect the reality of biological processes as closely as possible.

The resulting model consists of 642 elements and 1007 reactions, including iron homeostasis, oxygen-related reactions and the main carbon, nitrogen, sulphur and phosphate pathways. It includes several levels of regulation: gene expression ("housekeeping" or conditionally expressed genes), proteins production, protein targeting to specific subcellular locations and degradation, nutrient supply and transformation through the main metabolic pathways and the uptake and use of ions. Combined with our simulation algorithm, this model can be used to investigate the effects of specific perturbations on simulations, facilitating direct comparison with the "wild-type" situation. These perturbations may be changes in the "growth conditions" or gene deletions.

The analysis of simulations for selected in silico mutants showed the consistency of the model to specific experimental phenotypes. Comparison of the simulations of the wild-type model run with and without oxygen to experimental metabolic fluxes measured under similar conditions showed that most aspects of the adaptive responses of the system were accurately addressed. In a global analysis, in silico mutations generated for each gene of the model were compared with the "wild-type" model. A clustering analysis for all the simulations led to the identification of clear phenotypic profiles, thus providing new insights into some metabolic response to stress conditions. Finally, the model was also used to explore several new hypotheses in order to better understand some unexpected phenotypes in given mutants. All these results show that this model, and the underlying modelling strategy, are powerful tools for improving our understanding of iron homeostasis in relation to the main features of metabolism, energy production and oxidative stress.


Building the model

The list of reactions

As a starting point, the S. cerevisiae bibliome was searched for references related to a list of genes known to be involved in iron homeostasis. The bibliome consist of all articles concerning S. cerevisiae found on PubMed.

From these articles, a list of reactions was manually inferred and selected on the basis of our knowledge.

The new genes cited in these articles that we felt to be critical for model accuracy were then used to direct a new PubMed search. This process was repeated until we were unable to identify any new reactions strictly related to iron homeostasis.

The same process was used for inorganic phosphate homeostasis, for some oxidative stress reactions and for any other processes that we felt to be important for our model that were not described in databases.

We searched the Swissprot database for a list of proteins requiring iron-sulphur clusters or haem as a cofactor. We then selected metabolic pathways involving these proteins that could be directly or indirectly linked to iron homeostasis or oxidative stress. The reactions describing these pathways were expressed according to SGD pathways [68] on the Saccharomyces Genome Database website http://​www.​yeastgenome.​org.

However, when including a given pathway, we did not systematically describe all its steps. If a pathway included several steps producing intermediate metabolites not required by any other pathway included in our model, we wrote the whole pathway as one reaction, unless we had reasons to include the intermediate reactions. For example, siroheme biosynthesis from uroporphyrinogen-III involves three reactions. The two intermediate metabolites produced, precorrin-2 and sirohydrochlorin, are not required by any other reaction already included in the model, so we could have expressed the whole process as one reaction. However, the first step involves S-adenosyl-methionine and S-adenosyl-homocysteine, which are already included in several reactions in our model, and the last step involves iron. Thus, if siroheme cannot be produced - in a mutant for instance - we want to be able to determine whether this deficiency is related to 2-S-adenosyl-methionine synthesis or iron availability. We therefore included siroheme biosynthesis as two reactions, the first producing precorrin-2 from uroporphyrinogen-III and 2-S-adenosyl-methionine and the second producing siroheme from precorrin-2 and NAD.

Finally, we searched the yeast metabolome (described in SGD pathways) for reactions that might link several metabolites already included in our model. This is the reason for which we included alanine degradation, for example, in the model.

The weights of the reactions

The default value for the weight of the reactions was one. However some reactions were given a weight lower than one (for most degradation reactions, the weight is 0.01), or higher than one. For example, the reaction catalysed by the Sod1 superoxide dismutase was given a weight of 100, to take into account the extremely high abundance of this protein in yeast cells (519,000 molecules per cell; [69]), its very high catalytic efficiency and turnover number. A list of reactions given weights other than one is provided in Table 2. We simulated our model with all weights set to 1, and it was unable to produce realistic outputs (data not shown). For example, the hydroxyl radical elements have PoP lower than 1% in the complete model, whereas the PoP is higher than 70% in the model without weights, which is not realistic (the parameter PoP is defined below, in the Simulations subsection). The WT model corresponds to a model where the outputs are biologically meaningfull according the data of the literature and our own experience. We performed a sensitivity analysis of the outputs of the model (PoP at steady state) when the weights of the reactions are modified. Each weight was multiplied by a coefficient k. The differences between each PoP of the initial model and the PoP of the model with the modified weight were computed. Our results show that the model is robust to modifications of the weights when the coefficient k ranges from 0.1 to 10. See Additional file 2 for more details.



We used a modified version of the Biocham asynchronous Boolean simulation algorithm, which can be summarised as follows:
  1. 1.

    Initial state: the list of elements that are "ON" at the beginning of the simulation.

  2. 2.

    Based on the list of elements that are "ON", the list of possible reactions is inferred: A reaction is possible if its reactants and modifiers are "ON".

  3. 3.

    A reaction is randomly selected from the list of possible reactions

  4. 4.

    The products of the selected reaction are set to "ON".

  5. 5.

    The reactants are randomly either set to "OFF" or left "ON".

  6. 6.

    A new list of elements that are "ON" is computed.

  7. 7.

    Steps 2 to 6 are repeated for each simulation step.


In the Biocham algorithm, the reactants are not always set to "OFF", so it is possible to reselect the reaction in subsequent steps (see above step 5). This possibility reflects the presence of more than one molecule of each sort in biological systems. If the same reaction is selected over and over again, all the reactants will eventually be turned "OFF" and the reaction will cease to be possible unless other reactions set the reactants back to "ON". Indeed, in biological systems, all elements may be considered to exist in limited quantities.

We also had to overcome a limitation of Biocham to obtain simulations as close as possible to real biological systems: this algorithm does not mimic the large differences in reactions rates observed in real biological systems. Although we do not precisely know the rate of all the reactions in the model, we can reasonably state that the degradation of an enzyme is far less likely to occur than the reaction this enzyme catalyzes. Another example is related to the functioning of the TCA cycle: we know that some reactions of the cycle are technically reversible but the reaction always runs in one direction in practice, because this direction is thermodynamically more favorable. We needed to mimic this situation, to make our model as realistic as possible. We therefore decided to extend the algorithm, by weighting the reactions. From the weights of the possible reaction, at step 2, in the previous algorithm, we calculated the probability of a reaction being selected as the weight of the reaction divided by the sum of the weights of all possible reactions. A reaction is therefore more likely to occur if it has a high weighting.

If the experimental rates of all reactions in the model were known, we could set their weights accordingly. As these rates were not all known, we defined relative weights, using a default rate of one. The weight of a given reaction was modified if we had quantitative or qualitative experimental knowledge relating to its reaction rate (see Table 2 for a list of reactions with weights other than one).

In this paper, all the elements were set to "ON" at the beginning of the simulations, defining thus the "initial conditions". We tested different initial conditions (for the non constant elements only) which always resulted in the same steady-state (data not shown).

Mean profiles

The simulation algorithm was developed in the C (for computation) and Python (for file/data management) languages on a Linux workstation. The simulations were performed on a cluster of 40 nodes (DualCore AMD 64 bits Opteron bi-Processor, 2Go RAM, PBS/Maui scheduler). Each type of model was simulated 100 times, until a steady state was reached (see below). From these 100 simulations and for each element, we defined the PoP as the percentage of simulations in which the element was "ON" (present), for each simulation step. This calculation process is referred to as the profile of an element. Figure 4 shows examples of such profiles in which the mean PoP was also calculated every 100,000 steps.

Steady states

We checked whether steady state had been reached, by performing an ANOVA on the last 2 million steps, for every element (the null hypothesis is that the PoPs are independent from the iterations). If a significant result was obtained (p-value below 0.001), then the simulation was rerun for a larger number of steps. As background noise in the simulation can generate false positives, if the means of the last million steps of the new simulation of each element was equal to the means of the previous simulations, then we considered steady state to have been reached.

For some mutants, some elements have PoPs that increase or decrease slowly. For example, in some mutants the reactive oxygen species have increasing PoP. Most of them will have a PoP of 100 at steady-state but some will reached this maximum value faster than others. We advocate that comparing the simulations of the mutants after the same number of steps may be related to some biological properties of the mutants. To be able to take into consideration these differences in our systematic in silico mutants analysis, we used the PoP of the elements after 20 millions of steps even though the simulations were not at steady-state yet.

74% of the simulations reached steady-state before 20 million steps and among the remaining 26% only some very slowly increasing or decreasing elements have not reached steady-state yet.

In silico mutations and mutants

Model simulations began with all elements "ON" and continued until steady state was reached: this situation corresponds to "wild-type" simulations. To simulate a mutant for a given gene, a wild-type simulation was carried out for one million of steps. The gene in question was then turned "OFF" and the simulation was continued for 19 more millions of steps for the systematic mutants analysis or until steady-state was reached for the other analysis. This method for simulating mutants mimics experimental Tet-OFF mutants, in which the transcription of a given gene is controlled by a tetracycline- responsive promoter and can be turned off by adding tetracycline to the growth medium [70].

The model was simulated according to this procedure with, for each set of simulations, one of the unlimited elements deleted. Each type of model in which an unlimited element was turned "OFF" was referred to as a "mutant". Note that none of the reactions were modified.


For each "mutant" model, the mean of the last million of steps of the simulation was calculated for each element. Then, for each element and for each "mutant", a distance to the "wild-type" was calculated as follows: (formula based on A.Ultsch RelDi [71]). In the resulting matrix, each column is an element of the model and each row a "mutant". This matrix was then clustered, using Autoclass@IJM (real location data type, relative error = 0.01), which identifies classes of "mutants" causing similar changes in the simulations. The matrix was then transposed and clustered again: this generated classes of elements changing in a similar fashion in the different "mutants". Figure 6 shows the matrix clustered for both elements and "mutants".

Figure 6 is annotated: the rows are annotated with the most significantly enriched gene ontology process (calculated with GoTermFinder [72]), the columns were annotated manually (only the columns containing many elements which PoP changes significantly were annotated - see Additional file 6 for the complete clustering).

Phenotypes analysis

The file containing the manually curated phenotypes of the yeast mutants was retrieved from SGD (file edited the 03/13/10). From this file, we extracted the phenotypes associated with mutants of genes present in the model. From this selection, we further extracted the phenotypes that could be compared with our simulation results (e.g. the auxotrophies related to molecules present in the model). Six types of phenotypes were extracted: "auxotrophies" (cysteine, methionine, heme, glutamate), "chemical compound accumulation" (of elements present in the model), "oxidative stress resistance", "respiratory growth" and "nitrogen source". Then, we compared these phenotypes with the PoP of the corresponding elements in the WT model and in the mutants. For the "oxidative stress resistance", we did not simulate the mutants with an additional source of stress, as where observed the experimental phenotypes. Instead, we looked for the production of ROS in the mutant as compared to the WT in standard simulations. Therefore, only the constitutively stressed mutants showed similar phenotypes in our simulations and in vivo. As for the "respiratory growth" phenotypes, we compared the PoP of the protons in the intermembrane space (noted "Hinter::mitochondria" in the model), because in the model an increase or a decrease in this element PoP can be directly linked to a defect in the respiratory metabolism. For the "nitrogen source utilization" phenotypes, we compared the PoP of the products of the utilization of the nitrogen source. See Additional file 3 for the full results.



Fiona Achcar holds a grant-in-aid from the French Ministére de l'enseigneur supérieur et de la recherche. This project was partly funded by grants from the Digiteo Research Park (Saint-Aubin, France) and the Microsoft High Performance Computing for Education program. The English text was edited by Alex Edelman and Associates.

Authors’ Affiliations

Modelling in Integrative Biology, Institut Jacques Monod - UMR7592 - CNRS - Univ. Paris-Diderot
Mitochondria, Metals and Oxidative Stress, Institut Jacques Monod - UMR7592 - CNRS - Univ. Paris-Diderot


  1. Millen KJ, Gleeson JG: Cerebellar development and disease. Curr Opin Neurobiol. 2008, 18: 12-9. 10.1016/j.conb.2008.05.010PubMed CentralView ArticlePubMed
  2. Orth M, Schapira AH: Mitochondria and degenerative disorders. Am J Med Genet. 2001, 106: 27-36. 10.1002/ajmg.1425View ArticlePubMed
  3. Schapira AH: Mitochondrial dysfunction in neurodegenerative diseases. Neurochem Res. 2008, 33 (12): 2502-9. 10.1007/s11064-008-9855-xView ArticlePubMed
  4. Tabrizi SJ, Schapira AH: Secondary abnormalities of mitochondrial DNA associated with neurodegeneration. Biochem Soc Symp. 1999, 66: 99-110.View ArticlePubMed
  5. Xu JX: Radical metabolism is partner to energy metabolism in mitochondria. Ann N Y Acad Sci. 2004, 1011: 57-60. 10.1196/annals.1293.006View ArticlePubMed
  6. Murphy MP: How mitochondria produce reactive oxygen species. Biochem J. 2009, 417: 1-13. 10.1042/BJ20081386PubMed CentralView ArticlePubMed
  7. Evans MD, Dizdaroglu M, Cooke MS: Oxidative DNA damage and disease: induction, repair and significance. Mutat Res. 2004, 567: 1-61. 10.1016/j.mrrev.2003.11.001View ArticlePubMed
  8. Grimsrud PA, Xie H, Griffin TJ, Bernlohr DA: Oxidative stress and covalent modification of protein with bioactive aldehydes. J Biol Chem. 2008, 283 (32): 21837-41. 10.1074/jbc.R700019200PubMed CentralView ArticlePubMed
  9. Niki E: Lipid peroxidation: physiological levels and dual biological effects. Free Radic Biol Med. 2009, 47 (5): 469-84. 10.1016/j.freeradbiomed.2009.05.032View ArticlePubMed
  10. Kwong JQ, Beal MF, Manfredi G: The role of mitochondria in inherited neurodegenerative diseases. J Neurochem. 2006, 97 (6): 1659-75. 10.1111/j.1471-4159.2006.03990.xView ArticlePubMed
  11. Ohta S, Ohsawa I: Dysfunction of mitochondria and oxidative stress in the pathogenesis of Alzheimer's disease: on defects in the cytochrome c oxidase complex and aldehyde detoxification. J Alzheimers Dis. 2006, 9 (2): 155-66.PubMed
  12. Valko M, Morris H, Cronin MT: Metals, toxicity and oxidative stress. Curr Med Chem. 2005, 12 (10): 1161-208. 10.2174/0929867053764635View ArticlePubMed
  13. Heimpel H, Anselstetter V, Chrobak L, Denecke J, Einsiedler B, Gallmeier K, Griesshammer A, Marquardt T, Janka-Schaub G, Kron M, Kohne E: Congenital dyserythropoietic anemia type II: epidemiology, clinical appearance, and prognosis based on long-term observation. Blood. 2003, 102 (13): 4576-81. 10.1182/blood-2003-02-0613View ArticlePubMed
  14. Lee PL, Beutler E: Regulation of hepcidin and iron-overload disease. Annu Rev Pathol. 2009, 4: 489-515. 10.1146/annurev.pathol.4.110807.092205View ArticlePubMed
  15. Nicolas G, Bennoun M, Porteu A, Mativet S, Beaumont C, Grandchamp B, Sirito M, Sawadogo M, Kahn A, Vaulont S: Severe iron deficiency anemia in transgenic mice expressing liver hepcidin. Proc Natl Acad Sci USA. 2002, 99 (7): 4596-601. 10.1073/pnas.072632499PubMed CentralView ArticlePubMed
  16. Calzone L, Gelay A, Zinovyev A, Radvanyi F, Barillot E: A comprehensive modular map of molecular interactions in RB/E2F pathway. Mol Syst Biol. 2008, 4: 173- 10.1038/msb.2008.7PubMed CentralView ArticlePubMed
  17. Kohn KW, Aladjem MI, Kim S, Weinstein JN, Pommier Y: Depicting combinatorial complexity with the molecular interaction map notation. Mol Syst Biol. 2006, 2: 51-PubMed CentralPubMed
  18. Gerard C, Goldbeter A: Temporal self-organization of the cyclin/Cdk network driving the mammalian cell cycle. Proc Natl Acad Sci USA. 2009, 106 (51): 21643-8. 10.1073/pnas.0903827106PubMed CentralView ArticlePubMed
  19. Noble D: Biophysics and systems biology. Philos Transact A Math Phys Eng Sci. 2010, 368 (1914): 1125-39. 10.1098/rsta.2009.0245View Article
  20. Tyson JJ, Albert R, Goldbeter A, Sible J: Biological switches and clocks. J R Soc Interface. 2008, 5 (Suppl 1): S1-8. 10.1098/rsif.2008.0179.focusPubMed CentralView ArticlePubMed
  21. Alvarez-Buylla ER, Chaos A, Aldana M, Benitez M, Cortes-Poza Y, Espinosa-Soto C, Hartasanchez DA, Lotto RB, Malkin D, Santos GJE, Padilla-Longoria P: Floral morphogenesis: stochastic explorations of a gene network epigenetic landscape. PLoS One. 2008, 3 (11): e3626- 10.1371/journal.pone.0003626PubMed CentralView ArticlePubMed
  22. Frazier JM, Chushak Y, Foy B: Stochastic simulation and analysis of biomolecular reaction networks. BMC Syst Biol. 2009, 3: 64- 10.1186/1752-0509-3-64PubMed CentralView ArticlePubMed
  23. Shmulevich I, Aitchison JD: Deterministic and stochastic models of genetic regulatory networks. Methods Enzymol. 2009, 467: 335-56. full_text full_textPubMed CentralView ArticlePubMed
  24. Blinov ML, Faeder JR, Goldstein B, Hlavacek WS: A network model of early events in epidermal growth factor receptor signaling that accounts for combinatorial complexity. Biosystems. 2006, 83 (2-3): 136-51. 10.1016/j.biosystems.2005.06.014View ArticlePubMed
  25. Faeder JR, Blinov ML, Hlavacek WS: Rule-based modeling of biochemical systems with BioNetGen. Methods Mol Biol. 2009, 500: 113-67. full_text full_textView ArticlePubMed
  26. Davidich MI, Bornholdt S: Boolean network model predicts cell cycle sequence of fission yeast. PLoS One. 2008, 3 (2): e1672- 10.1371/journal.pone.0001672PubMed CentralView ArticlePubMed
  27. Mai Z, Liu H: Boolean network-based analysis of the apoptosis network: irreversible apoptosis and stable surviving. J Theor Biol. 2009, 259 (4): 760-9. 10.1016/j.jtbi.2009.04.024View ArticlePubMed
  28. Mangla K, Dill DL, Horowitz MA: Timing robustness in the budding and fission yeast cell cycles. PloS One. 2010, 5 (2): e8906-[PMID: 20126540], 10.1371/journal.pone.0008906PubMed CentralView ArticlePubMed
  29. Schlatter R, Schmich K, Vizcarra IA, Scheurich P, Sauter T, Borner C, Ederer M, Merfort I, Sawodny O: ON/OFF and beyond-a boolean model of apoptosis. PLoS Comput Biol. 2009, 5 (12): e1000595- 10.1371/journal.pcbi.1000595PubMed CentralView ArticlePubMed
  30. Formanowicz D, Sackmann A, Formanowicz P, Blazewicz J: Petri net based model of the body iron homeostasis. J Biomed Inform. 2007, 40 (5): 476-85. 10.1016/j.jbi.2006.12.001View ArticlePubMed
  31. Hawari AH, Hussein ZA: Simulation of a Petri net-based Model of the Terpenoid Biosynthetic Pathway. BMC Bioinformatics. 2010, 11: 83- 10.1186/1471-2105-11-83PubMed CentralView ArticlePubMed
  32. Li C, Nagasaki M, Ueno K, Miyano S: Simulation-based model checking approach to cell fate specification during Caenorhabditis elegans vulval development by hybrid functional Petri net with extension. BMC Syst Biol. 2009, 3: 42- 10.1186/1752-0509-3-42PubMed CentralView ArticlePubMed
  33. Semsey S, Andersson AMC, Krishna S, Jensen MH, Massé E, Sneppen K: Genetic regulation of fluxes: iron homeostasis of Escherichia coli. Nucleic Acids Research. 2006, 34 (17): 4960-4967. 10.1093/nar/gkl627PubMed CentralView ArticlePubMed
  34. Sackmann A, Formanowicz D, Formanowicz P, Koch I, Blazewicz J: An analysis of the Petri net based model of the human body iron homeostasis process. Computational Biology and Chemistry. 2007, 31: 1-10. 10.1016/j.compbiolchem.2006.09.005View ArticlePubMed
  35. Sarria B, Dainty JR, Fox TE, Fairweather-Tait SJ: Estimation of iron absorption in humans using compartmental modelling. European Journal of Clinical Nutrition. 2005, 59: 142-144.http://​www.​ncbi.​nlm.​nih.​gov.​gate1.​inist.​fr/​pubmed/​15252423 10.1038/sj.ejcn.1602030View ArticlePubMed
  36. Jouhten P, Rintala E, Huuskonen A, Tamminen A, Toivari M, Wiebe M, Ruohonen L, Penttila M, Maaheimo H: Oxygen dependence of metabolic fluxes and energy generation of Saccharomyces cerevisiae CEN.PK113-1A. BMC Syst Biol. 2008, 2: 60- 10.1186/1752-0509-2-60PubMed CentralView ArticlePubMed
  37. Lill R, Muhlenhoff U: Iron-sulfur protein biogenesis in eukaryotes: components and mechanisms. Annu Rev Cell Dev Biol. 2006, 22: 457-86. 10.1146/annurev.cellbio.22.010305.104538View ArticlePubMed
  38. Rouault TA, Tong WH: Iron-sulphur cluster biogenesis and mitochondrial iron homeostasis. Nat Rev Mol Cell Biol. 2005, 6 (4): 345-51. 10.1038/nrm1620View ArticlePubMed
  39. Chabrier-Rivier N, Chiaverini M, Danos V, Fages F, Schächter V: Modeling and querying biomolecular interaction networks. Theoretical Computer Science. 2004, 325: 25-44. 10.1016/j.tcs.2004.03.063http://​www.​sciencedirect.​com/​science/​article/​B6V1G-4C74FM9-1/​2/​978d950405c118e8​82fd1008fba7ab23​ 10.1016/j.tcs.2004.03.063View Article
  40. Urban-Grimal D, Labbe-Bois R: Genetic and biochemical characterization of mutants of Saccharomyces cerevisiae blocked in six different steps of heme biosynthesis. Mol Gen Genet. 1981, 183: 85-92. 10.1007/BF00270144View ArticlePubMed
  41. Liu XF, Elashvili I, Gralla EB, Valentine JS, Lapinskas P, Culotta VC: Yeast lacking superoxide dismutase. Isolation of genetic suppressors. J Biol Chem. 1992, 267 (26): 18298-302.PubMed
  42. Ogur M, Roshanmanesh A, Ogur S: Tricarboxylic Acid Cycle Mutants in Saccharomyces: Comparison of Independently Derived Mutants. Science. 1965, 147: 1590- 10.1126/science.147.3665.1590View ArticlePubMed
  43. Gollub EG, Liu KP, Dayan J, Adlersberg M, Sprinson DB: Yeast mutants deficient in heme biosynthesis and a heme mutant additionally blocked in cyclization of 2, 3-oxidosqualene. J Biol Chem. 1977, 252 (9): 2846-54.PubMed
  44. Tamura S, Yoshimura E: Promotion of Zn(2+) uptake by Al (3+) in a Saccharomyces Cerevisiae mutant that lacks the ZRT1 gene encoding a high-affinity Zn transporter. Biological Trace Element Research. 2008, 124 (3): 262-268. 10.1007/s12011-008-8145-4View ArticlePubMed
  45. Auchère F, Santos R, Planamente S, Lesuisse E, Camadro J: Glutathione-dependent redox status of frataxin-deficient cells in a yeast model of Friedreich's ataxia. Human Molecular Genetics. 2008, 17 (18): 2790-2802. [PMID: 18562474],http://​www.​ncbi.​nlm.​nih.​gov.​gate1.​inist.​fr/​pubmed/​18562474View ArticlePubMed
  46. Kastaniotis AJ, Zitomer RS: Rox1 mediated repression. Oxygen dependent repression in yeast. Adv Exp Med Biol. 2000, 475: 185-95. full_text full_textView ArticlePubMed
  47. Zitomer RS, Lowry CV: Regulation of gene expression by oxygen in Saccharomyces cerevisiae. Microbiol Rev. 1992, 56: 1-11.PubMed CentralPubMed
  48. Lutstorf U, Megnet R: Multiple forms of alcohol dehydrogenase in Saccharomyces cerevisiae. I. Physiological control of ADH-2 and properties of ADH-2 and ADH-4. Arch Biochem Biophys. 1968, 126 (3): 933-44. 10.1016/0003-9861(68)90487-6View ArticlePubMed
  49. Achcar F, Camadro J, Mestivier D: AutoClass@IJM: a powerful tool for Bayesian classification of heterogeneous data in biology. Nucleic Acids Research. 2009, W63-67. 37 Web Server,
  50. Babcock M, de Silva D, Oaks R, Davis-Kaplan S, Jiralerspong S, Montermini L, Pandolfo M, Kaplan J: Regulation of mitochondrial iron accumulation by Yfh1p, a putative homolog of frataxin. Science. 1997, 276 (5319): 1709-12. 10.1126/science.276.5319.1709View ArticlePubMed
  51. Foury F, Talibi D: Mitochondrial control of iron homeostasis. A genome wide analysis of gene expression in a yeast frataxin-deficient strain. J Biol Chem. 2001, 276 (11): 7762-8. 10.1074/jbc.M005804200View ArticlePubMed
  52. Measday V, Baetz K, Guzzo J, Yuen K, Kwok T, Sheikh B, Ding H, Ueta R, Hoac T, Cheng B, Pot I, Tong A, Yamaguchi-Iwai Y, Boone C, Hieter P, Andrews B: Systematic yeast synthetic lethal and synthetic dosage lethal screens identify genes required for chromosome segregation. Proc Natl Acad Sci USA. 2005, 102 (39): 13956-61. 10.1073/pnas.0503504102PubMed CentralView ArticlePubMed
  53. Thomas D, Jacquemin I, Surdin-Kerjan Y: MET4, a leucine zipper protein, and centromere-binding factor 1 are both required for transcriptional activation of sulfur metabolism in Saccharomyces cerevisiae. Mol Cell Biol. 1992, 12 (4): 1719-27.PubMed CentralView ArticlePubMed
  54. Lesuisse E, Santos R, Matzanke BF, Knight SAB, Camadro J, Dancis A: Iron use for haeme synthesis is under control of the yeast frataxin homologue (Yfh1). Human Molecular Genetics. 2003, 12 (8): 879-889. 10.1093/hmg/ddg096View ArticlePubMed
  55. Miao R, Kim H, Koppolu UMK, Ellis EA, Scott RA, Lindahl PA: Biophysical characterization of the iron in mitochondria from Atm1p-depleted Saccharomyces cerevisiae. Biochemistry. 2009, 48 (40): 9556-9568. 10.1021/bi901110nPubMed CentralView ArticlePubMed
  56. Miao R, Martinho M, Morales JG, Kim H, Ellis EA, Lill R, Hendrich MP, Münck E, Lindahl PA: EPR and Mössbauer spectroscopy of intact mitochondria isolated from Yah1p-depleted Saccharomyces cerevisiae. Biochemistry. 2008, 47 (37): 9888-9899. 10.1021/bi801047qView ArticlePubMed
  57. Lide DR: Solubility Product Constants. CRC Handbook of Chemistry and Physics. 2009, CRC Press, 90,http://​www.​hbcpnetbase.​com
  58. Ligeti E, Fonyó A: Mitochondrial phosphate carrier. Functional role of its SH groups and interrelations within the carrier unit. European Journal of Biochemistry/FEBS. 1987, 167: 167-173. 10.1111/j.1432-1033.1987.tb13319.xView ArticlePubMed
  59. Madeo F, Fröhlich E, Ligr M, Grey M, Sigrist SJ, Wolf DH, Fröhlich KU: Oxygen stress: a regulator of apoptosis in yeast. The Journal of Cell Biology. 1999, 145 (4): 757-767. 10.1083/jcb.145.4.757PubMed CentralView ArticlePubMed
  60. Sipos K, Lange H, Fekete Z, Ullmann P, Lill R, Kispal G: Maturation of cytosolic iron-sulfur proteins requires glutathione. The Journal of Biological Chemistry. 2002, 277 (30): 26944-26949. 10.1074/jbc.M200677200View ArticlePubMed
  61. Kohn KW, Aladjem MI, Weinstein JN, Pommier Y: Molecular interaction maps of bioregulatory networks: a general rubric for systems biology. Mol Biol Cell. 2006, 17: 1-13. 10.1091/mbc.E05-09-0824PubMed CentralView ArticlePubMed
  62. Hower V, Mendes P, Torti FM, Laubenbacher R, Akman S, Shulaev V, Torti SV: A general map of iron metabolism and tissue-speci c subnetworks. Mol Biosyst. 2009, 5 (5): 422-43. 10.1039/b816714cPubMed CentralView ArticlePubMed
  63. Hlavacek WS: How to deal with large models?. Molecular Systems Biology. 2009, 5: 240- 10.1038/msb.2008.80PubMed CentralView ArticlePubMed
  64. Alves R, Herrero E, Sorribas A: Predictive reconstruction of the mitochondrial iron-sulfur cluster assembly metabolism. II. Role of glutaredoxin Grx5. Proteins. 2004, 57 (3): 481-492. 10.1002/prot.20228View ArticlePubMed
  65. Alves R, Herrero E, Sorribas A: Predictive reconstruction of the mitochondrial iron-sulfur cluster assembly metabolism: I. The role of the protein pair ferredoxin-ferredoxin reductase (Yah1-Arh1). Proteins. 2004, 56 (2): 354-366. 10.1002/prot.20110View ArticlePubMed
  66. Qi Z, Miller GW, Voit EO: Computational Systems Analysis of Dopamine Metabolism. PLoS ONE. 2008, 3 (6): e2444.- 10.1371/journal.pone.0002444PubMed CentralView ArticlePubMed
  67. Zhou L, Cortassa S, Wei A, Aon MA, Winslow RL, O'Rourke B: Modeling Cardiac Action Potential Shortening Driven by Oxidative Stress-Induced Mitochondrial Oscillations in Guinea Pig Cardiomyocytes. Biophysical Journal. 2009, 97 (7): 1843-1852. 10.1016/j.bpj.2009.07.029PubMed CentralView ArticlePubMed
  68. Caspi R, Foerster H, Fulcher CA, Hopkinson R, Ingraham J, Kaipa P, Krummenacker M, Paley S, Pick J, Rhee SY, Tissier C, Zhang P, Karp PD: MetaCyc: a multiorganism database of metabolic pathways and enzymes. Nucleic Acids Research. 2006, D511-516. 34 Database,
  69. Ghaemmaghami S, Huh W, Bower K, Howson RW, Belle A, Dephoure N, O'Shea EK, Weissman JS: Global analysis of protein expression in yeast. Nature. 2003, 425 (6959): 737-741. 10.1038/nature02046View ArticlePubMed
  70. Mnaimneh S, Davierwala AP, Haynes J, Moffat J, Peng WT, Zhang W, Yang X, Pootoolal J, Chua G, Lopez A, Trochesset M, Morse D, Krogan NJ, Hiley SL, Li Z, Morris Q, Grigull J, Mitsakakis N, Roberts CJ, Greenblatt JF, Boone C, Kaiser CA, Andrews BJ, Hughes TR: Exploration of essential gene functions via titratable promoter alleles. Cell. 2004, 118: 31-44. 10.1016/j.cell.2004.06.013View ArticlePubMed
  71. Ultsch A: Improving the identification of differentially expressed genes in cDNA microarray experiments. Classification - the Ubiquitous Challenge. 2005, 378-385. full_text.http://​[http://​apps.​isiknowledge.​com.​gate1.​inist.​fr/​fullrecord.​do?​product=​WOS&​search m ode=​GeneralSearch&​qid=​1&​SID=​U2pB75G28PK3AI@kBaf&​page=​1&​doc=​1]View Article
  72. Boyle EI, Weng S, Gollub J, Jin H, Botstein D, Cherry JM, Sherlock G: GO::TermFinder - open source software for accessing Gene Ontology information and finding significantly enriched Gene Ontology terms associated with a list of genes. Bioinformatics. 2004, 20 (18): 3710-3715. 10.1093/bioinformatics/bth456PubMed CentralView ArticlePubMed
  73. Chen OS, Crisp RJ, Valachovic M, Bard M, Winge DR, Kaplan J: Transcription of the yeast iron regulon does not respond directly to iron but rather to iron-sulfur cluster biosynthesis. The Journal of Biological Chemistry. 2004, 279 (28): 29513-8. 10.1074/jbc.M403209200View ArticlePubMed
  74. Dix DR, Bridgham JT, Broderius MA, Byersdorfer CA, Eide DJ: The FET4 gene encodes the low affinity Fe(II) transport protein of Saccharomyces cerevisiae. The Journal of Biological Chemistry. 1994, 269 (42): 26092-26099.PubMed
  75. Lin SS, Manchester JK, Gordon JI: Enhanced gluconeogenesis and increased energy storage as hallmarks of aging in Saccharomyces cerevisiae. The Journal of Biological Chemistry. 2001, 276 (38): 36000-7. 10.1074/jbc.M103509200View ArticlePubMed
  76. Minard KI, McAlister-Henn L: Sources of NADPH in yeast vary with carbon source. The Journal of Biological Chemistry. 2005, 280 (48): 39890-39896. 10.1074/jbc.M509461200View ArticlePubMed
  77. Voet D, Voet JG: Regulation of the Citric Acid Cycle. In Biochemistry. 1995, 557-559. Wiley, 2,
  78. Hinkle PC, McCarty RE: How cells make ATP. Scientific American. 1978, 238 (3): 104-117. 10.1038/scientificamerican0378-104View ArticlePubMed
  79. Auchère F, Rusnak F: What is the ultimate fate of superoxide anion in vivo?. Journal of Biological Inorganic Chemistry: JBIC: A Publication of the Society of Biological Inorganic Chemistry. 2002, 7 (6): 664-7.View Article
  80. Aebi H: Catalase in vitro. Methods Enzymol. 1984, 105: 121-6. full_text full_textView ArticlePubMed


© Achcar et al; licensee BioMed Central Ltd. 2011

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://​creativecommons.​org/​licenses/​by/​2.​0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.