Extension of the yeast metabolic model to include iron metabolism and its use to estimate global levels of iron-recruiting enzyme abundance from cofactor requirements.

Metabolic networks adapt to changes in their environment by modulating the activity of their enzymes and transporters; often by changing their abundance. Understanding such quantitative changes can shed light onto how metabolic adaptation works, or how it can fail and lead to a metabolically dysfunctional state. We propose a strategy to quantify metabolic protein requirements for cofactor-utilising enzymes and transporters through constraint-based modelling. The first eukaryotic genome-scale metabolic model to comprehensively represent iron metabolism was constructed, extending the most recent community model of the Saccharomyces cerevisiae metabolic network. Partial functional impairment of the genes involved in the maturation of iron-sulphur (Fe-S) proteins was investigated employing the model and the in silico analysis revealed extensive rewiring of the fluxes in response to this functional impairment, despite its marginal phenotypic effect. The optimal turnover rate of enzymes bearing ion cofactors can be determined via this novel approach; yeast metabolism, at steady state, was determined to employ a constant turnover of its iron-recruiting enzyme at a rate of 3.02 × 10 -11  mmol·(g biomass) -1 ·h  -1 .

Metabolic models with high predictive ability are important tools for the investigation and engineering of metabolism (Aung, Henry, & Walker, 2013). Metabolic models can provide reasonable predictions when direct measurements of network components are infeasible, as is demonstrated by predictions on flux distributions using genome-scale models (Orth, Thiele, & Palsson, 2010). Analyses with stoichiometric models can be used to predict flux, but only a limited number of studies exist on their use for estimating absolute enzyme abundances, for example in yeast (Nilsson & Nielsen, 2016;Sánchez et al., 2017); an alternative approach would be to make use of the cofactors to determine the relative abundances of those enzymes that use those cofactors.
To exploit this relationship between metabolic enzymes and their cofactors, we propose a network-based strategy to determine optimal enzyme abundance by model predictions based on the metabolic requirements of their cognate cofactors. We have studied a family of cofactors that are not, themselves, metabolic intermediates and which do not have donor functional groups. We chose to study the iron cofactor family, including iron in its ionic form (Fe 3+ , Fe 2+ ) and complex ion entities such as the haem family (sirohaem, haems A, B, C and O), and the iron-sulphur clusters (2Fe-2S, 4Fe-4S), because these iron-containing cofactors bind nearly 10% of the documented enzymes and transporters in the metabolic network. Baker's yeast, the first eukaryote to have its genome sequenced (Goffeau et al., 1996), has long been a favourite model organism (Botstein, Chervitz, & Cherry, 1997). Cellular activities including DNA replication, recombination and repair, RNA transcription and translation, intracellular trafficking, as well as the enzymatic activities of general metabolism, and mitochondrial biogenesis are conserved from yeast to human (Barrientos, 2003). The availability of comprehensive and powerful genome-scale models of the yeast metabolic network (Aung et al., 2013) for almost two decades (Famili, Forster, Nielsen, & Palsson, 2003) made yeast an ideal model for our study. Although there is a wealth of knowledge about iron utilisation and homoeostasis in yeast (De Freitas et al., 2003;Lill & Mühlenhoff, 2008;Miethke & Marahiel, 2007), this information has not been integrated into the curated genome-scale metabolic model, thus limiting the model's usefulness for in silico studies. This issue has recently been highlighted by a comprehensive study on the comparative analysis of yeast metabolic models (Heavner & Price, 2015).
In this study, we have extended the genome-scale metabolic model of Saccharomyces cerevisiae to include iron metabolism, and made the first comprehensive mathematical representation of any inorganic ion in a model of a network in a eukaryote; this model is publicly available (BioModels [Chelliah et al., 2015] access no: MODEL1709260000). The new model has been benchmarked against known environmental responses (viz., changes in the availability of iron and copper) and genetic modifications relating to iron metabolism. The model has allowed the identification of the extensive rewiring of metabolic fluxes to cope with the hemizygosity of essential genes involved Fe-S cluster maturation. The model permitted the establishment of the requirements for iron-family cofactors, based on how the fluxes were distributed through the metabolic network. These requirements were then used as proxies for calculating the metabolic requirement of those metabolic enzymes and transporters that employ iron species as cofactors.

| Primary metabolic model, simulation environment and model annotation
The primary model for the incorporation of iron metabolism was selected as the most recent stoichiometric model of the S. cerevisiae metabolic network (v7.6; Aung et al., 2013). The extended model (Yeast7.Fe) is provided as Supporting Information S1 and the details on modification and extension of the existing model are provided in Supporting Information S2 in supplemental material. The extended genome-scale model of S. cerevisiae is available in a COBRA compatible-SBML format (v.4). A specific growth rate of 0.1 hr −1 was used unless otherwise specified. The iron and copper supplementations of the system was determined from a standard defined medium for S. cerevisiae (Baganz, Hayes, Marren, Gardner, & Oliver, 1997). Iron and copper limitations were imposed on the system by scaling the extracellular availability of the respective ions in the medium down to 10% of the original upper bound.
Details regarding the extension annotations and the simulation environment and approaches are provided in Supporting Information S3, and the simulation code is provided in Supporting Information S4.

| Cofactor representation
For any metabolic reaction catalysed by a cofactor requiring enzyme (N): where the stoichiometric coefficients of the substrates (A and B) and products (C and D) are denoted by their cognate lower case letters, a cofactor (X) of enzyme N enters the metabolic reaction as a substrate and leaves as a catalytically unreactive product: Details regarding the representation of cofactors are provided in Supporting Information S3. The distribution of fluxes throughout the metabolic network would be correctly influenced by the requirements for, and the availability of, iron in this representation. Iron-containing cofactors could be required for some of the enzymes necessary for catalysing a reaction, and a dedicated amount of cofactors should be reserved for this use, without being actually involved in the reaction itself. The effect that this modification had on the distribution of fluxes, including the growth rate, will be discussed in Section 3.

| Iron metabolism in yeast metabolic models
Although the predictive power of the yeast metabolic network model has improved substantially over the years, it is still limited by the F I G U R E 1 Schematic representation of iron metabolism in the yeast model. Pathways and metabolites are represented by uppercase and lowercase letters, respectively. Directionality of the fluxes through the pathway is specified by arrows (for single reaction steps) and block arrows (for lumped consecutive reaction steps). Metabolic enzymes are shown in teal colour. The cell and organelle boundaries are represented as double dashed lines; the mitochondrion, nucleus, vacuole and ER have cartoon representations. The minimal representation of iron metabolism in the existing Yeast 7.6 model is provided in (a). The details on the reductive, non-reductive and xenosiderophore-bound iron uptake, intracellular transport and storage of iron, haem and sirohaem biosynthetic and degradation pathways are provided in (b). Details regarding the biogenesis of Fe-S clusters in the mitochondrion (ISC machinery), and the maturation of apoenzymes (A) into Fe-S cluster-bound holoenzymes (H) in the mitochondrion (ISC machinery), in the cytosol and in the nucleus (CIA machinery) are provided in (c). An empty scaffold (ES) and its sulphonylated form (SS) were introduced as pseudometabolites in the Fe-S cluster formation mechanism. The regulation of iron uptake via the iron regulon, employing the negative feedback from Fe-S cluster biogenesis, is demonstrated in (d). The shuttling of the signals representing the availability of mitochondrial iron (PS) and its depletion (AS) were introduced as pseudo-metabolites to modulate and activate the reductive iron uptake routes. For simplifications of the function and activity of the iron regulon, see text [Color figure can be viewed at wileyonlinelibrary.com] omission or incomplete representation of some pathways in the model. Our previous work on the Yeast 7 metabolic network model highlighted the fact that the high metabolic burden (characterised by high reaction fluxes, and often indicative of the low efficiency of their cognate enzymes (Bonarius, Hatzimanikatis, Meesters, Schmid, & Tramper, 1996) carried by pathways involved in energy generation and processes imposed limitations on the model's predictive ability (Dikicioglu, Kırdar, & Oliver, 2015). A recent analysis that we carried out by constraining this model by fluxes calculated using the intracellular concentrations of intermediates in the purine nucleotide biosynthetic pathway, as determined by HPLC analysis (Hesketh, Vergnano, Wan, & Oliver, 2017), demonstrated that the prediction of growth rate was at least half or twice the experimentally determined value. The distribution of the fluxes indicated that the iron uptake and utilisation pathways were inactive because they were inadequately represented in the model and completely disconnected from the rest of the metabolic network ( Figure 1a). We addressed each of the deficiencies through extensive literature curation and will explain how iron metabolism was incorporated into the network model in the following sections.

| Uptake, intracellular transport and storage of iron
Both reductive (Fe 3+ -transporting, high affinity) and non-reductive (Fe 2+ -transporting, low affinity) iron uptake mechanisms were incorporated into the metabolic network. High-affinity iron uptake, as both free and xenosiderophore-bound iron, was represented by the reductive pathway. The intracellular transport of iron and complex iron entities, such as haem and sirohaem, across intracellular boundaries and the storage of residual components and excess iron were also considered ( Figure 1b). Iron export is not known to be exhibited by S. cerevisiae (Haas, Eisendle, & Turgeon, 2008) and so was excluded.
The concentration limit to determine the lower bound of the flux through the low-affinity iron uptake reactions was set at 1 µM (Lesuisse, Blaiseau, Dancis, & Camadro, 2001). The high-and lowaffinity uptake systems were modelled as working independently of one another. However, in the absence of a functional low-affinity system, the high-affinity system will be used, even in the presence of abundant iron (and vice versa); these situations are allowed for in the model (Haas et al., 2008;Lesuisse et al., 2001). Specific ARN family transporters involved in the uptake of each iron-bound xenosiderophore and the fate of each xenosiderophore in the yeast cell was The mechanism of copper recruitment by the reductive pathway of high-affinity iron uptake (Philpott, 2006) necessitated the incorporation of copper uptake into the metabolic network. Since copper metabolism was not represented in the primary metabolic network of yeast to any degree, the network was extended to incorporate the uptake of copper and its function in high-affinity iron uptake. It is important to note that the representation of copper was not exhaustive; only activities that are relevant to iron metabolism were considered in this reconstruction. The threshold for switching between high-and low-affinity copper transport was set as 20 µM (Hassett, Dix, Eide, & Kosman, 2000). The threshold concentrations were used to calculate a threshold flux boundary to be used in the model. Two assumptions were made in modelling the uptake of copper across the cell envelope: (a) Although included as a unique species in the model, the high-affinity copper transporter Ctr3p was not associated with the copper uptake reaction in the model along with Ctr1p. This exclusion was necessary since the CTR3 gene, found in strains of the S288C lineage that were used in this study, has been inactivated by insertion of a Ty2 transposon (Knight, Labbé, Kwon, Kosman, & Thiele, 1996). (b) Fre1p copper reductase was similarly excluded from gene-reaction associations since it was reported to be active only during the first 3-4 hr post-inoculation (Georgatsou & Alexandraki, 1994). The in silico analyses conducted in this study were all carried out using pseudo-steady-state assumption; by that time, Fre2p was thought to determine copper reductase activity.

| Biosynthesis and recycling of complex iron entities
The biosynthesis and recycling of haem, sirohaem and Fe/S clusters were fully implemented in the genome-scale model of the yeast metabolic network. The synthesis of 5-aminolevulinate from glycine and succinyl-CoA via the Shemin pathway (Ferreira & Gong, 1995) was introduced and 5-aminolevulinate uptake was excluded from the model. Haem catabolism in the endoplasmic reticulum, and the storage of bilirubin in the vacuole were also introduced to the model de novo (Figure 1b).

| The iron regulon
Two different regulatory mechanisms for iron uptake via the iron regulon were implemented in this model: (a) the yeast-specific haem-regulated positive feedback route (Lill & Mühlenhoff, 2008), and (b) the mitochondrial Fe-S cluster biogenesis-associated negative feedback route.
Haem deficiency, indicated by the deactivation of the enzymes encoded by the essential genes (HEM1-4, HEM12, HEM13, HEM15) of the pathway, was reported to be an indicator of low-iron uptake in yeast (Lill & Mühlenhoff, 2008). We coupled the essential enzymes of haem biosynthesis to the reaction representing low-affinity iron uptake catalysed by Fet4p to account for this positive feedback mechanism by associating these essential enzymes with the lowaffinity iron uptake reaction. This coupling with low-affinity iron uptake ensured that haem biosynthesis would not be over-activated by the network unless extracellular iron was abundant.
The negative feedback on iron uptake by Fe/S cluster biogenesis is shut down upon iron depletion. This was modelled by introducing two Aft1p, which is responsible for relaying the iron depletion signal to the cell boundary, was reported to have an additional role in creating iron resources for the cell by binding Cth2p to facilitate the degradation of the mRNAs for iron-containing enzymes (Cherry et al., 2012). This route was excluded from the network since the model does not treat enzymes as either the substrates or products of reactions. Aft1p was also reported to mediate saving iron by activating Vth1p-mediated biotin uptake since biotin synthesis was reported to be iron consuming (Shakoury-Elizeh et al., 2004). This route was also excluded since biotin co-enzyme metabolism was not considered to be of direct interest to our system.

| Iron family cofactor considerations of metabolic enzymes
Iron family cofactors, copper ions and pyridoxine (the last is used in only one reaction) were used as substrates and untransformed reactants in those reactions that were catalysed by enzymes activated by these cofactors, as described in "Cofactor representation." This allowed us to establish the connectivity between iron metabolism, described above, and the primary metabolic network of yeast. For this purpose, the copper, iron, haem, sirohaem, pyridoxine and Fe-S requirements of the metabolic network were Reactions, which were reported to take place at the membrane but did not have a detailed mechanism explained, were represented to occur across the membrane. This device simplified the model without compromising on the details of iron metabolism. The enzymes associated with these reactions were localised to the membrane to highlight the transmembrane nature of the reaction.
Pseudo-metabolites were introduced such that they did not interfere with the material, energy, or redox balances of the network. They were introduced in coupled reactions to avoid accumulation, and these reactions were unbounded so that the system was not constrained by the flux limitations through these reactions when the underdetermined system was optimised for a given objective.
Only 14  Fe, showed that they were associated with genes that were significantly enriched for the nucleoside phosphate metabolism process term (p < 10 -38 ) along with other metabolic processes in line with our observations on the analysis involving the purine intermediates used as flux constraints. We carried out a sensitivity analysis to investigate robustness by selecting the fluxes that displayed a change more than 15-35% (in 5% increments), and observed that the genes encoding the enzymes that catalysed these reactions were significantly enriched for the same, or very similar processes (see Supporting Information S5 in supplemental material).
The variability of the flux distributions was taken into consideration in conducting this analysis (for details of the method, see Supporting Information S3 in supplemental material).
The analysis of the fluxes in the biomass-modified model, Y7.
FeBM, indicated that the magnitude of only 8% of the fluxes were altered as a response to this change and that the magnitude of the change was minimal (Supporting Information S10). The modifications were proportional enrichment in the fluxes for the production of the iron entities that were incorporated into the biomass equation.
Despite being far from providing a complete picture due to the problems regarding data availability discussed above, this exercise of

| Iron-recruiting enzyme requirements of the metabolic network model
The stoichiometry of the cofactors introduced into the reactions was observed to be closely related to the predictions of growth rate, which meant that such predictions were very sensitive to the metabolic requirement for iron in the network. The total iron requirement of the cell was calculated from the iron composition of a standard defined medium for yeast (Chiu & Segrè, 2008), and the recruitment of iron or complex iron entities as cofactors by the enzymes was simulated based on the constraint imposed by the iron transport flux. Fine-tuning for the sustainable iron recruitment indicated that the stoichiometric coefficients for these terms need to be in the order of magnitude of 10 -14 if the supply of iron is not to reduce the rate of growth (Figure 2). The metabolic requirement of total turnover for iron-recruiting enzymes was determined as 3.02 × 10 -11 mmol·(g biomass) −1 ·h −1 based on this model. The

| Benchmarking the model against environmental and genetic challenges
Having established the iron-recruiting enzyme requirement of the metabolic network, we extended the investigation to understand how the in silico system adapted to environmental challenges that were tailored specifically to exploit this model. The response of the network to different levels of extracellular iron and copper availability was investigated through flux distributions. For this purpose, the upper and lower bounds of the fluxes through several reactions were set to zero to mimic reports available in the literature (Haas et al., 2008;Lesuisse et al., 2001). The high-affinity iron uptake system was activated by blocking the routes through low affinity Fe 2+ uptake routes. The absence of extracellular xenosiderophores was also taken into consideration for the simulation of an S. cerevisiae monoculture. The predictions on flux distributions were in agreement with the expected physiological outcome on how the yeast metabolism responded to low or high abundance of copper (Cankorur-Cetinkaya et al., 2013;Vest et al., 2016) or iron (Holmes-Hampton, Jhurry, McCormick, & Lindahl, 2013;Vest et al., 2016; Table 3).
We then investigated how Y7.Fe could be used to study the impact of the deletion of the genes ARH1, ATM1 or YFH1 or a reduction in the functionality of their protein products, which are all essential components of the ISC machinery. All three genes are essential, although only arh1 deletants could be identified as inviable employing the model. This implies that YFH1 (the Fe-S cluster scaffold protein) and ATM1 (the Fe-S apocluster transporter) have essential functions outside of the metabolic network, perhaps for the supply of Fe-S clusters for essential non-enzyme proteins.
Despite their essentiality, the hemizygous mutants of these genes in the diploid BY4743 genetic background did not produce any significant growth defects (p < 0.01) with the specific growth rate for all mutants remaining within ± − 0.43 0.01 h 1 as also indicated by the model predictions. Our experiments demonstrated that the hemizygous mutants did not display significant differences (p < 0.01) with respect to their utilisation of glucose or ammonium as carbon and nitrogen sources, respectively, nor in their production of ethanol or glycerol. On average, the hemizygous mutants were observed to consume more iron per unit optical density equivalent number of cells although this difference was not observed to be significant due to the high variance between replicates. However, a significant difference in the intracellular distribution of copper was observed between these mutants and the wild type was observed. That said, measurements of intracellular haem, copper, reduced or total iron content of the mutants gave no further insight into the essential roles of ARH1 , YFH1 and ATM1 inside or outside metabolism (Supporting Information S5).
The reactions with non-zero fluxes were more often catalysed by enzymes encoded by essential genes (2.32-fold over-enrichment; p < 10 -40 ) in Y7.Fe than in Y7.6. In line with this observation, ARH1 was widely associated with reactions having non-zero fluxes in the metabolic network. ARH1 being an essential gene also for the model, the reorganisation of the fluxes in the network could only be investigated by introducing an "in silico reduction of function." We used flux as a proxy for the functional capacity of the enzyme catalysing that reaction and lowered the flux value by 50% to mimic the impact of heterozygosity. The flux bounds that were set to 50% of their wild-type flux value did not relate to nutrient uptake fluxes; F I G U R E 2 Indispensable role of iron for yeast. This plot demonstrates how growth rate predictions of the metabolic network model are affected by the amount of iron recruited by the metabolic enzymes. The stoichiometric coefficient of iron ions and complex iron entities to be recruited by iron requiring enzymes without impairing growth substantially can be determined based on the constraints imposed by the experimentally permissible limits of iron uptake therefore, they did not change the nutrient limitation of the system and this has been confirmed by the glucose uptake flux being maintained at its upper bound limit under both conditions. Consequently, the predicted growth rate was reduced by 50% in response to limiting the flux through the reaction catalysed by ARH1 at 50% of its original value, the expected outcome for the remaining non-zerofluxes was to be reduced by 50%, or to remain unchanged. However, 14% of all fluxes were observed to be rewired or had changed by a factor other than a 50% reduction.
Some 22% of the enzymes and transporters represented in Y7.Fe The extended model of the yeast metabolic network now serves as a functional platform to study rare disorders associated with the iron metabolism, as well as those of other ions to the extent of their involvement in iron-associated pathways. ARH1, studied extensively above, encodes the yeast homolog of the adrenodoxin reductase, and mutations in the human gene are reported to be responsible for auditory neuropathy and optic atrophy (Paul et al., 2017). The model also allowed us to simulate the impairment of CCC2, the coppertransporting ATPase, mutations in whose human homolog, ATP7B, is responsible for Wilson's disease (Cankorur-Cetinkaya et al., 2013;Júlvez, Dikicioglu, & Oliver, 2018). Using the model, we were able to demonstrate a critical impairment associated with this disordernamely, that under conditions of high oxygen availability, the deletion of CCC2 impaired growth in the absence of copper supplementation (Table 3).

| DISCUSSION
The complex network of interactions between enzymes, metabolites and their regulators defines transitions between metabolic states.
Metabolic state shifts mediate adaptive responses that allow cells to respond to environmental or genetic perturbations. Such shifts in the state of metabolism may be desirable in some cases, such as biotechnological applications, or most unwelcome in other situations such as switching to a diseased state that may induce cancer or other conditions. Genome-scale models of the metabolic network are of central importance since they allow us to achieve a better understanding of how the fluxes can be rewired in response to an internal or external input, and thus have been frequently used in biotechnological applications as a prediction tool, since the aim is often to obtain a preferred flux distribution (Kerkhoven, Lahtvee, & Nielsen, 2014). Iron has an important role in a range of biotechnological applications from improving human nutrition (Puig, Andres-Ccolas, Garcia-Molina, & Penarrubia, 2007) to the production of recombinant proteins (Eck et al., 2018), and healthcare applications concerning iron storage abnormalities (Valerio, 2007). The incorporation of an extensive, well-curated iron metabolism in metabolic networks may greatly enhance the opportunities metabolic models offer for applications in red, green, or white biotechnology.
This paper has presented a comprehensive extension of the genome-scale stoichiometric model of the yeast metabolic network by incorporating the involvement of ion cofactors, which are not, themselves, consumed in the metabolic reactions in which they are involved. We selected the iron ion as our test case as it is involved in nearly 10% of all enzymatic steps in the yeast metabolic network. A substantial number of non-metabolic proteins also recruit iron cofactors. Although iron species to be used by non-metabolic proteins can be produced by the metabolic network model, they remain as dead-end metabolites since these models do not include non-metabolic biological processes. In its current form, introducing the empirical iron requirement of yeast as a constraint for the metabolic network would yield an inevitable overestimation of the iron requirement for metabolism, since all of the iron transported across the cell boundary can theoretically be utilised metabolically.
Nevertheless, the quantitation of distinct iron-containing cofactor species contributing to the biomass would allow an improved definition of biomass composition to be established, which would indirectly account for the iron species requirements of non-metabolic proteins, resolving this problem of overestimation.
Although this study was specific to iron, we contend that our work has provided a new approach for handling co-enzymes and co-factors in stoichiometric models of metabolism, and that this approach can be generalised for other such entities. We made use of pseudo-metabolites, which were not produced or consumed by the cell, in modelling the system. We used this formalism since these intermediates, functioning in cyclic interconversion pathways, allowed a more explicit description of the molecular mechanisms represented in the model. The Y7.Fe model allowed us to define the optimal turnover rate of iron cofactors in the metabolic network, something which is not possible to achieve experimentally using current analytical technologies. The method provided a means to estimate enzyme fluxes from reaction fluxes and cofactor availability. Proteomics data available at the global scale can, at present, only provide us with the relative quantitation of some enzymes; therefore, empirical protein abundance data do not afford us the means to evaluate our predictions. However, we note once again that ion cofactor availability, rather than empirically determined enzyme turnover rates, was more prominent as a limiting factor in determining the predictions that can be made using the Y7.Fe model. Therefore, we believe our method provides a very reasonable estimate of the ironrecruiting enzyme requirements of the metabolic network under defined circumstances, which is currently not possible with any other method.
Our model represents an unconventional way of extending metabolic network models to incorporate non-metabolic information.
The ion cofactors are "non-metabolic" in the sense that they are conserved moieties that are not, themselves, transformed by metabolic reactions. Although they can be considered as "pseudo- g. glucose) availability, but also by the availability of an essential micronutrient, that is iron. The ability to incorporate the roles of such ion cofactors could prove essential for constructing models of biological systems at the level of the whole cell or organism.