- Research article
- Open Access
- Published:

# A mathematical model of iron import and trafficking in wild-type and Mrs3/4ΔΔ yeast cells

*BMC Systems Biology*
**volume 13**, Article number: 23 (2019)

## Abstract

### Background

Iron plays crucial roles in the metabolism of eukaryotic cells. Much iron is trafficked into mitochondria where it is used for iron-sulfur cluster assembly and heme biosynthesis. A yeast strain in which Mrs3/4, the high-affinity iron importers on the mitochondrial inner membrane, are deleted exhibits a slow-growth phenotype when grown under iron-deficient conditions. However, these cells grow at WT rates under iron-sufficient conditions. The object of this study was to develop a mathematical model that could explain this recovery on the molecular level.

### Results

A multi-tiered strategy was used to solve an ordinary-differential-equations-based mathematical model of iron import, trafficking, and regulation in growing *Saccharomyces cerevisiae* cells. At the simplest level of modeling, all iron in the cell was presumed to be a single species and the cell was considered to be a single homogeneous volume. Optimized parameters associated with the rate of iron import and the rate of dilution due to cell growth were determined. At the next level of complexity, the cell was divided into three regions, including cytosol, mitochondria, and vacuoles, each of which was presumed to contain a single form of iron. Optimized parameters associated with import into these regions were determined. At the final level of complexity, nine components were assumed within the same three cellular regions. Parameters obtained at simpler levels of complexity were used to help solve the more complex versions of the model; this was advantageous because the data used for solving the simpler model variants were more reliable and complete relative to those required for the more complex variants. The optimized full-complexity model simulated the observed phenotype of WT and Mrs3/4ΔΔ cells with acceptable fidelity, and the model exhibited some predictive power.

### Conclusions

The developed model highlights the importance of an Fe^{II} mitochondrial pool and the necessary exclusion of O_{2} in the mitochondrial matrix for eukaryotic iron-sulfur cluster metabolism. Similar multi-tiered strategies could be used for any micronutrient in which concentrations and metabolic forms have been determined in different organelles within a growing eukaryotic cell.

## Background

The complexity of biochemical processes in growing eukaryotic cells is enormous, often rendering the corresponding genetic phenotypes difficult to understand at the chemical level. One means of analyzing such systems is to develop ordinary-differential-equation (ODE^{Footnote 1})-based kinetic models [1,2,3]. In principle, such models can reveal on a quantitative basis whether observed phenotypic behavior could emerge from a proposed system of reacting chemical players using a particular set of kinetic and thermodynamic parameters. This is a huge advantage relative to the common practice of describing complex biochemical processes as a cartoon or scheme. Another advantage of math-based kinetic models is that all assumptions are explicit and available for public inspection whereas cartoons and schemes generally include hidden assumptions. The major *disadvantage* of such kinetic models is that a complete and accurate dataset, including rate-law expressions, rate-constants, and reactant concentrations, are required to solve them and to endow them with predictive power. Rarely is all such information available, and available information is often less quantitative than desired.

A common approach to circumventing this problem is to employ *simple* models (in terms of numbers of components and reactions) that nevertheless remain capable of generating observed cellular behavior and of explaining genetic phenotypes. Designing such models involves deciding which species and reactions to include, which to leave out, and which to combine into groups. Such decisions often boil-down to whether including an additional component or reaction is “worth” (in terms of generating the desired behavior) an additional adjustable parameter. Simple models with few adjustable parameters simplify reality but they can also provide fundamental insights into reality - by penetrating through the entangled and bewildering complexity of a highly complex system.

Iron is critical for all eukaryotic cells [4, 5]. It is present in many forms including heme centers, iron-sulfur clusters (ISCs), nonheme mononuclear species, and iron-oxo dimeric centers. Such centers are commonly found in the active-sites of metalloenzymes. Iron plays a major role in energy metabolism; e.g. there are iron-rich respiratory complexes located on the inner membrane of mitochondria. Mitochondria are the primary site in the cell where ISCs are assembled, and the only site where iron is installed into porphyrins during heme biosynthesis. For these reasons, mitochondria are a major ‘hub’ for iron trafficking.

The cytosol also plays an important role in iron trafficking, in that nutrient iron enters this region prior to being distributed to the organelles. Most of the iron that enters the cytosol is probably in the Fe^{II} state, but neither the oxidation state nor the concentration of cytosolic Fe has been established [6]. The vacuoles are another trafficking ‘hub’ in yeast, as much of the iron imported into these cells (when grown on iron-sufficient media) is stored in these acidic organelles [7, 8]. Vacuolar iron is predominately found as a mononuclear nonheme high spin (NHHS) Fe^{III} species, probably coordinated to polyphosphate ions [9].

Iron is tightly regulated in cells, and some insightful mathematical models involving iron metabolism, trafficking and regulation have been developed. Twenty years ago, Omholt et al. designed and analyzed a model of the IRP/IRE iron regulatory system in mammalian cells [10]. More recently, Mobilia et al. developed a similar model that assumed scarce or unavailable data; they also developed new methods to represent data by constrained inequalities [11, 12]. Chifman and coworkers developed an ODE-based model for iron dysregulation in cancer cells in which the roles of the IRP-based regulation, the iron storage protein ferritin, the iron export protein ferroportin, the labile iron pool, reactive oxygen species, and the cancer-associated Ras protein were emphasized [13], as well as a logical-rule-based mathematical model of iron homeostasis in healthy mammalian cells [14]. Mitchell and Mendes used ODE’s to model iron metabolism and regulation in a liver cell and its interaction with blood plasma [15]. They emphasized the role of iron-regulating hormone hepcidin and the regulatory and storage systems mentioned above, and they simulated the effects of iron-overload disease. Their model was complex - involving 66 adjustable parameters many of which were not experimentally determined. None of the above models included iron-sulfur cluster (ISC) synthesis, the role of mitochondria (or other organelles), and none modeled growing cells. In terms of biological emphasis, the model of Achcar et al. [16] is most relevant to the current study. They developed a model of iron metabolism and oxidative stress in yeast cells using a Boolean approach using weighted reactions. Their model included ISC assembly, as well as organelles such as mitochondria, vacuoles, cytosol, and nucleus. However, their model was exceedingly complicated (642 components and 1007 reactions) and was not ODE based [16]. They modeled the development of Fe^{III} (phosphate) oxyhydroxide nanoparticles in mitochondria of mutant cells lacking ISC assembly proteins (e.g. Yfh1, the yeast frataxin homolog), similar to the emphasis of our previous model [17]. They included a reaction in which an unidentified species X converted nanoparticles into free iron, and hypothesized that X might be glutathione. In contrast, our model emphasized the role of oxygen in controlling nanoparticle formation.

The iron content of yeast cells and the major organelles involved in iron trafficking have been analyzed using Mössbauer (MB) spectroscopy, the most powerful spectroscopic tool for interrogating the iron content of biological samples [18]. If the absolute iron concentration of ^{57}Fe-enriched cells and organelles are known, the absolute concentrations of major groups of iron-containing species in such cells can be calculated using percentages obtained by MB. Such data is used here to develop an advanced mathematical model of iron import and trafficking in eukaryotic cells.

In WT cells, much iron enters mitochondria through Mrs3 and Mrs4, paralogous inner-membrane proteins [19, 20]. These “high affinity” iron-importers contain a small tunnel that allows a low-molecular-mass cytosolic iron species to enter the matrix. We have recently discovered a low-molecular-mass species in mitochondria, designated Fe_{580}, which might serve as feedstock for ISC assembly [21]. Iron can enter mitochondria through alternative pathways, including one that involves Rim2 [22].

Iron import in yeast is regulated according to the ISC activity occurring in mitochondria [23]. When this activity is attenuated, for example by mutations in the ISC assembly machinery, the rate of nutrient iron imported increases. In yeast, iron regulation involves the *Iron Regulon*, a group of 20–30 genes whose expression is controlled by transcription factors Aft1/2 [24, 25]. This includes the Fet3/Ftr1 complex on the plasma membrane through which much iron enters the cell.

Yfh1 helps catalyze ISC assembly in mitochondria [26]. This and other ISC mutant cells accumulate large quantities of iron in the form of Fe^{III} nanoparticles [27, 28]. These cells import excessive iron because the iron regulon is activated in response to insufficient mitochondrial ISCs. Excess iron (in the form of nanoparticles) accumulates in mitochondria because the rate of iron import into the organelle increases due to activation of the iron regulon. The *net* rate of iron import into vacuoles is reduced such that these organelles contain little iron in ISC mutants. Actually, the iron *export* rate is probably *increased* in these mutants. The vacuolar membrane contains an iron-export complex (Fet5/Fth1) that is homologous to the Fet3/Ftr1 iron import complex on the plasma membrane; both are controlled by the iron regulon [29].

We have developed a simple model (Fig. 1, bottom panel) to illustrate the changes in iron import and trafficking that occur in ISC mutants relative to in WT cells [17]. The core assumption of the model is that the matrix of healthy WT mitochondria is *micro-aerobic*, containing less O_{2} than in an aerobic state due to the ability of the respiratory complexes on the inner membrane to quickly reduce much of the O_{2} that would otherwise diffuse into the matrix. Although dissolved [O_{2}] concentrations in the mitochondrial matrix have not been measured directly, three lines of evidence indicate that this space is micro-aerobic. Firstly, in vitro ISC assembly assays must be performed anaerobically because some of the proteins involved are O_{2}-sensitive [30]. Secondly, numerous other enzymes in the matrix, including aconitase, biotin synthase, and lipoic acid synthase are O_{2}-sensitive [31, 32]. Thirdly, the nitrogenase iron protein which is exquisitely O_{2}-labile remains active when installed in the mitochondrial matrix [33]. According to our model, in ISC mutant mitochondria, the lack of ISCs and heme centers cause a deficiency of respiratory complexes; this condition allows O_{2} to diffuse into the matrix and react with a pool of Fe^{II}, forming nanoparticles.

Relative to WT cells, Mrs3/4ΔΔ cells (to be called **ΔΔ** hereafter) grow slowly under iron-deficient conditions but at WT rates in iron-sufficient media [5, 19, 20]. The iron concentration of ΔΔ cells is higher than in comparable WT cells, indicating that the iron regulon is activated. We recently found that mitochondria from iron-deficient ΔΔ cells are dominated by nanoparticles whereas the iron content of mitochondria from iron-sufficient ΔΔ cells are similar to WT mitochondria – i.e. dominated by the ISC and heme centers that are found in respiratory complexes. WT mitochondria also contain a substantial amount of a NHHS Fe^{II} that probably arise from Fe_{580} [34]. Fe_{580} is present in mitochondria from iron-replete ΔΔ and both iron-deficient and iron-replete WT cells. However, our previous model [17] was unable to reproduce the ΔΔ phenotype.

In this paper, we present an improved ODE-based model of iron trafficking and regulation in yeast, and use a multi-tiered strategy to solve it at an expanding steady-state. This model was able to explain both the ΔΔ *and* ΔYfh1 phenotypes while requiring fewer adjustable parameters relative to our previous model.

## Methods and results

As is typical of modeling biochemical processes within cells, the challenge was to generate a useful and insightful model despite sparse and imperfect data [11, 12]. Our strategy for doing this was to optimize the model at different levels of complexity. Model variants ranged from one that consisted of a single iron species and no cellular compartments to one that involved nine species in three cellular compartments. The parameters used to optimize the simpler variants were transferrable to the more complex models. This was an important insight because the data needed to solve simpler systems tend to be more reliable and complete relative to those required to solve more complex variants. A similar strategy could be applied for models involving the trafficking of other micronutrients. The only requirements are that the concentrations and metabolic forms of the micronutrient in the cell and in major organelles be known (at some reasonable level of accuracy) for different growth conditions and/or genetic strains.

The complete chemical model is shown in Fig. 1, bottom panel. We initially solved it (to be referred to as *C*_{9}, the “nine-component” model, including components *C*, *CIA*, *F2*, *F3*, *VP*, *FM*, *FS*, *MP*, and *O*_{2}) at three simpler levels of complexity called *C*_{1} (the “one-component” model, with component *Fe*_{cell}), *C*_{3} (the “three-component” model, including *Fe*_{cyt}, *Fe*_{mit}, and *Fe*_{vac}), and *C*_{4} (the “four-component” model, including *C*, *CIA*, *Fe*_{vac}, and *Fe*_{mit}). These model variants are illustrated in Fig. 1, top and middle panels. We solved C_{9} in this way because the data required to solve the simpler versions were more reliable and complete than those required to solve the C_{9} variant. Importantly, *the parameters that were optimized using the simpler versions could be transferred to the more complex variants*. This minimized the number of adjustable parameters that had to be assigned using less reliable or incomplete data. As far as we are aware, this multi-tiered modeling strategy has not been employed previously within the context of ODE-based models involving the trafficking of iron or any micronutrient within a growing eukaryotic cell. Code for all model variants was written using *Mathematica 10* software (wolfram.com). Initial concentrations for each iron component was 10 μM, and initial [O_{2}] was 0 μM. ODEs were solved to steady-state using the NDSolve routine.

### Development of the C_{1} model

Consider a population of cells growing exponentially on a nutrient form of iron called *N* which enters the cell through a transporter on the plasma membrane (Fig. 1, top panel, red circle). In the experimental results used here in fitting [34], N consisted of 0, 1, 10 and 40 μM ferric citrate plus 1 μM endogenous iron as found in minimal medium. Let *V*_{cell} represent the collective cell volume (within a culture) at time *t*. When cells are growing exponentially, *V*_{cell} will increase according to the relationship.

where *α*_{cell} is the growth rate. During exponential growth *α*_{cell} is constant in time. The optical density at 600 nm of an exponentially growing culture is proportional to *V*_{cell} such that the slope of the {ln(OD_{600}) vs. time} plot affords *α*_{cell}. This parameter has been determined for WT and ΔΔ cells grown in medium containing [N] = 1, 2, 11, and 41 μM ([34] and Table 1). The 8 “data-based” determinations of *α*_{cell} will be called *α*_{cell − dat}.

For simulations, a continuous α function between *N* = 1–41 μM was required. Plots of *α*_{cell} vs. [N] exhibited saturation behavior, suggesting the Michaelis-Menten type function

The desired continuous α function was obtained by fitting (2) against the *α*_{cell − dat}values using the error function.

This error function normalized absolute differences between simulations and data to the average of simulation and data values. This formulation weighed the error associated with each datapoint evenly without regard to the magnitude of the point. Best-fit *α*_{max} and *K*_{α}values are given in Table 2, and plots of *α*_{cell} are shown in Fig. 2a. The simulated growth rates of WT and ΔΔ cells increased as the concentration of iron in the medium [N] increased, mirroring the experimental growth rates with acceptable fidelity (apart from the point associated with ΔΔ cells at [N] = 41 μM). *Acceptable fidelity* is a qualitative term which means that simulations “trended” with the data – i.e. increasing when the data increased, decreasing when they decreased, and remaining flat when they remained flat. Due to the limited amount of data and our trial-and-error method of optimizing, this term is more appropriate than other more quantitative descriptors. The high iron concentration of ΔΔ cells grown with 40 μM ferric citrate is consistent with Mössbauer spectra which are very high-intensity. ΔΔ cells are iron-dysregulated and they import large amounts of iron, and the model developed here can only account for a portion of that accumulation. A more complex model could fit the point better but we opted to keep the model simple.

In the C_{1} model, all iron in the cell is considered to be a single component called *Fe*_{cell}. The concentration [Fe_{cell}] is a function of moles (*n*_{Fecell}) and volume *V*_{cell}, namely [Fe_{cell}] = *n*_{Fecell}/*V*_{cell}. Since the cell is growing as chemistry is occurring, the time-dependent change of [Fe_{cell}] is given by the partial derivative

The first term on the right-hand-side of the last equation of (4) describes the rate of iron import at constant volume (\( N\overset{R_{cell}}{\to }{Fe}_{cell} \)) – i.e. for chemistry occurring in a no-growth cell. The second term reflects dilution due to the growth of cells at constant moles of *Fe*_{cell} - .i.e. for a growing cell devoid of chemistry. Yeast cells lack iron exporters (unlike mammalian cells that contain ferroportin) and no export process is known. Under the *expanding-steady-state* condition, as would exist for a population of exponentially growing cells, the left-hand-side of (4) equals zero, [Fe_{cell}] is constant and the import rate *R*_{cell} equals the dilution rate,

[Fe_{cell}] was measured in WT and ΔΔ cells grown under the four concentrations of [N] [34], and the product of this and corresponding *α*_{cell − dat} values afforded the “data-based” *R*_{cell-dat} values listed in Table 1. These values are shown as the circles in Fig. 2c. Plots of [Fe_{cell}] vs. log_{2}[N] are given in Fig. 2b.

We next assigned a rate-law expression to *R*_{cell} that depended solely on [N], such that a continuous *R*_{cell-sim} function could be generated at all [N]. The iron-importer on the plasma membrane of yeast cells is saturatable by nutrient iron [35], and so we assigned the rate-law for *R*_{cell-sim} to the Hill function

where *sens* is a Hill coefficient allowing for cooperative iron import. *R*_{cell-sim} was optimized by minimizing an ERR function similar to Eq. (3). The resulting optimized *R*_{cell-sim} simulation parameters are given in Table 2. The *R*_{cell-sim} equation was used to generate an ODE (based on the last equation in (4)) that could be used in kinetic modeling (see Additional file 1: Equation S1). However, the current study focuses on the expanding steady-state condition, and so ODEs S1 and S2 were solved at infinitely long times for [N] ranging from 1 to 41. Plots of steady-state *R*_{cell-sim} vs. log_{2}[N] are shown in Fig. 2 bottom panel. As expected, the simulated rate of iron import increased in both WT and ΔΔ cells as the concentration of iron in the medium increased, with higher rates for ΔΔ cells since they accumulate more iron. The [N]-dependent increase in iron import rate is counterbalanced by the [N]-dependent increase in cell growth rate.

### Development of the C_{3} model

In the C_{3} model, cell volume was divided into mitochondria, vacuoles and all remaining compartments, such that

Here, “*cyt*” refers to cytosol *plus* all organelles besides mitochondria and vacuoles; there is insufficient published information to justify subdividing cyt into additional cellular compartments. This collective compartment includes the iron content of the nucleus which contains a significant number of [Fe_{4}S_{4}] containing proteins [36]. Topologically, cyt was treated as though it was exclusively cytosol i.e. surrounding mitochondria and vacuoles and being surrounded by the plasma membrane.

Each cellular compartment in the C_{3} model was presumed to contain a single iron species, called *Fe*_{cyt}, *Fe*_{mit}, and *Fe*_{vac}. The conservation of matter requires that

where *f*_{cyt}, *f*_{mit}, and *f*_{vac} are fractional volumes e.g. *f*_{mit} = *V*_{mit}/*V*_{cell}. In an expanding steady-state, these fractional volumes will be constant such that

All of the derivative terms in (9) are zero in an expanding steady state. For the C_{3} model, *N* is imported into the cytosol forming *Fe*_{cyt} (\( N\overset{R_{cyt}}{\to }{Fe}_{cyt} \)). Some *Fe*_{cyt} is imported into mitochondria (\( {Fe}_{cyt}\overset{R_{mit}}{\to }{Fe}_{mit} \)) and some into vacuoles (\( {Fe}_{cyt}\overset{R_{vac}}{\to }{Fe}_{vac} \)). The rest remains in cyt. Based on this scheme, the time-dependent changes of the concentrations of the Fe species in each region are

Equation (10) follows from the proposed mechanism in which iron first flows into the cytosol and then cytosolic iron *Fe*_{cyt} flows into mitochondria and vacuoles. Volume ratios in the second and third equations of (10) are required to conserve mass as *Fe*_{cyt} moves from one region to another. Under an expanding steady-state, the left-hand-sides of (10) equal zero and

The growth rate of each cellular region will equal the growth rate of the cell multiplied by the fractional volume of that compartment,

Substituting (12) into (11), and using (1) and (8) affords

Published fractional volumes were used to help solve these equations. The cellular content of fermenting exponentially growing, nonbudding *S. cerevisiae* was reconstructed in 3D, and volume fractions were determined [37]. Mitochondria and vacuoles occupied 1.7% and 5.8% of cell volume, respectively. Another study reported that the same two organelles occupied 1.6 and 7.8%, respectively [38]. In a third study, vacuoles in yeast strain W303 (the same as used in our studies) accounted for 10% of cell volume [39]. And in *respiring* yeast cells, mitochondrial volume was 10%–12% of cell volume [40]. Since the model developed here is of iron trafficking in respiring W303 yeast cells, we assumed *f*_{mit} = 0.1, *f*_{vac} = 0.1, and *f*_{cyt} = 0.8.

The relationships given in (13) are connected to (5). Substituting the last two equations of (13) into the first, and then simplifying and comparing to (5) affords the relationship

This equation connects C_{1} and C_{3} models. The rate of iron import into cyt (*R*_{cyt}) equals the data-based rate of Fe import into the cell (*R*_{cell}) divided by the volume fraction *f*_{cyt}. These rates describe the change of iron *concentrations* within the cell or cytosol, not the change in the number of moles of *N* imported. Since *V*_{cyt} < *V*_{cell}, [Fe_{cyt}] will increase faster than [Fe_{cell}], in proportion to the ratio *V*_{cell}*/V*_{cyt}. This is true even though the same number of moles of iron per unit time is imported. The rate-law expression for *R*_{cyt-sim} should also involve a Hill expression, with the same *K*_{N} and *sens* as in (6) but with a maximal velocity that is 1.25-times (1/*f*_{cyt}) faster.

The C_{3} model could not be solved fully until [Fe_{cell}] was separated into [Fe_{cyt}], [Fe_{mit}], and [Fe_{vac}] components for each of the 8 growth/strain conditions investigated. To do this, we relied on the conservation-of-matter Eq. (8), published MB spectra, and on iron concentrations for WT and ΔΔ cells and organelles [34, 41]. The spectra were separated into contributions from the eight iron-containing components specified by the C_{9} model. Then we combined particular components into cytosol, mitochondria, or vacuoles locations (as defined by the model of Fig. 1). Finally, we summed the iron concentrations for all of the components assigned to each compartment to afford our best data-based estimates of [Fe_{cyt}], [Fe_{mit}], and [Fe_{vac}]. Results are given in Table 1.

### Development of the C_{9} model

Before explaining how MB spectra were decomposed, we introduce the components of the C_{9} model. Component *C* represents a NHHS Fe^{II} complex presumed to be present in the cytosol. This component can move into vacuoles and mitochondria, but it can also stay in cyt and react to form component *CIA* (the Cytosolic Iron Sulfur Assembly machine), a second cytosolic iron species. *CIA* represents the sum of the ISCs and low-spin Fe^{II} heme groups in this collective compartment. Numerous ISCs are found in the cytosol and nucleus [37, 42], justifying the inclusion of *CIA* in the model. *FM* represents the pool of NHHS Fe^{II} ions in mitochondria, *FS* represents ISCs and heme centers in the organelle, and *MP* refers to mitochondrial nanoparticles. Components *FM*, *FS*, and *MP* have all been characterized experimentally. *F2* and *F3* are NHHS Fe^{II} and Fe^{III} species in vacuoles, and *VP* represent vacuolar nanoparticles; they have also been characterized experimentally [9, 43]. When *C* enters the vacuoles, this component becomes *F2* some of which oxidizes to *F3*. Some *F3* converts into *VP*. When *C* enters mitochondria, it converts into *FM*, which serves as feedstock for *FS*. *FS* metal centers are viewed as being installed into the respiratory complexes on the inner membrane, which then catalyze the reduction of *O*_{2} to water. *FM* can also react with model component *O*_{2} in the matrix to generate *MP* and ROS. ROS exhibits the exact behavior of *MP* so is not formally included in the model.

### Decomposing Mössbauer features into modeling components

MB spectroscopy detects all of the ^{57}Fe in samples. However, resolution is limited so the spectra under consideration were subdivided into just four groups of iron centers, including NHHS Fe^{III}, NHHS Fe^{II}, the central doublet (CD), and Fe^{III} oxyhydroxide nanoparticles. The CD represents [Fe_{4}S_{4}]^{2+} clusters and low-spin Fe^{II} heme centers; the two cannot be resolved. Other minor spectral features (HS Fe^{II} heme groups and [Fe_{2}S_{2}] clusters) can be resolved and quantified, but we decided to bundle them with the CD since they are not individually represented in the model. The absolute concentrations associated with each group were obtained by multiplying the associated percentages by [Fe_{cell}]. The conservation of mass requires that

These MB features were assigned to the following combinations of modeling components.

Then these species were organized into the three cellular compartments by summing contributions as described by (17).

The one component of the C_{9} model that could not be determined in this way with reasonable accuracy was component *C*. Thus, we relied on published reports to estimate the concentration of cytosolic NHHS Fe^{II}. Petrat et al. [44] used a fluorescent chelator to quantify the concentration of labile iron in hepatocytes and liver endothelial cells at 5–7 μM, and we assumed similar values for iron-sufficient WT yeast cells. We further assumed that the concentration of cytosolic Fe^{II} increases with increasing nutrient iron concentrations, and that [C] in iron-sufficient ΔΔ cells is higher than in WT cells (because the absence of Mrs3/4 should block import of *C* into mitochondria). Within these constraints, we assigned the concentrations of *C* to those listed in Table 3 so as to minimize the ERR function when the other iron concentrations in the table were used. MB spectral decompositions, along with these relationships and assumptions, were sufficient to generate concentrations for all other modeling components (Table 3).

#### Solving the C_{3} model

Once [Fe_{cyt}], [Fe_{mit}] and [Fe_{vac}] were determined, we determined rates of import into each compartment, *R*_{cyt}, *R*_{mit}, and *R*_{vac} as defined by (13). Data-based import rates *R*_{cyt-dat}, *R*_{mit-dat}, and *R*_{vac-dat} for the 8 conditions are shown as circles in Fig. 3 and are tabulated in Table 1. The rate of iron import into “cyt” but not exported into mitochondria or vacuoles equals *R*_{cyt} – *R*_{mit} - *R*_{vac}. According to these rates, iron flows faster into the cyt of ΔΔ cells, and slower into their mitochondria, relative to in WT cells. This makes sense because the absence of Mrs3/4 in ΔΔ cells should hinder *Fe*_{cyt} from entering mitochondria.

We next assigned rate-law expressions to *R*_{mit-sim} and *R*_{vac-sim}. We considered two forms for rate-laws, namely a mass-action form *R*_{i} = *k*_{i}[C]^{n} and a Hill form *V*_{i}[C]^{n}/{*K*_{M}+[C]^{n}} where [C] indicates the concentration of *C* as defined by the C_{9} model components. The latter form was used only when the simpler mass-action form was unable to generate reasonable simulations of the relevant data-based rates. The simple mass-action form was acceptable for *R*_{mit-sim} whereas *R*_{vac-sim} required a Hill term. The terms were optimized using an ERR function. The following rate-law expressions were ultimately selected.

One potentially confusing aspect of solving this system was that we used *C* (a component of the C_{4} and C_{9} models) rather than *Fe*_{cyt} (a component of the C_{3} model) as substrate for these processes. This was done so that the resulting rate-constant k_{mit} and R_{vac-max} would not change when solving the C_{9} model. Had we used *Fe*_{cyt} instead of *C* as substrate in the C_{3} model, *k*_{mit} and *R*_{vac-max} would have been too small for the C_{4} and C_{9} models in which [Fe_{cyt}] = [C] + [CIA]. The rate of iron flowing into mitochondria depends only on [C], not on [C] + [CIA]. Optimized *R*_{mit-sim} and *R*_{vac-sim} values were used along with *R*_{cyt-sim} (obtained from the C_{1} model), to construct a full set of ODEs (Additional file 1: Equations S3–S5) describing the C_{3} model. Once combined in this way, all of the parameters associated with the three rates *R*_{cyt-sim}, *R*_{mit-sim}, and *R*_{vac-sim} were re-optimized against data-based rates using an ERR function. To do this, each parameter was increased and decreased by 10% as all other parameters were fixed; candidate values that lowered ERR were then fixed as the next parameter on the list was varied. The process was repeated for a second round except that each parameter was adjusted ±5%. In the third and final round, each parameter was adjusted ±1%. The final plots are shown in Fig. 3.

### The C_{4} model

We next solved the C_{4} model which was identical to C_{3} except that [Fe_{cyt}] was separated into [CIA] and [C] components. To obtain [CIA], we subtracted the values of [C] given in Table 3 from [Fe_{cyt}], resulting in the *CIA* concentrations listed in Table 3. These values were multiplied by *α*_{cell} to generate *R*_{cia-dat}. We assumed a Hill expression to generate an *R*_{cia-sim} function that minimized differences with *R*_{cia-dat} with acceptable fidelity.

#### Solving the C_{9} model

At this point, the C_{9} model could be solved. The derivative of (17) is

According to the mechanism of Fig. 1, bottom panel, the rates of change in the concentrations of the two cyt iron species are

Adding the two equations of (20) affords the first equations of (19) and (10). The rate of change of the concentrations of the iron-containing species in the mitochondria is given by (21),

Summing the equations of (21) affords the second equations of (19) and (10). Similarly for vacuoles,

Summing the equations of (22) affords the third equations of (19) and (10). Thus, the ODE system for the iron-components of the C_{9} model “collapses” down to that of the C_{3} model when the components of the three regions are summed appropriately. In an expanding steady state, the left-hand-sides of (20), (21), and (22) equal zero such that

Data-based and simulation-based values of *R*_{cyt}, *R*_{mit}, *R*_{vac}, and *R*_{cia} have already been obtained. Using the experimental values of *α*_{cell} and the values of model-component concentrations listed in Table 3, we constructed data-based rates for the formation of each C_{9} component using data from the 4 nutrient conditions in WT and ΔΔ cells (Table 4). *R*_{vp-dat} was then used along with *α*_{cell} and [F3] to generate *R*_{23-dat} as defined in (23). The next step was to assign a rate-law expression to each of the remaining rates associated with the C_{9} model as listed in (23) – expressions that depended solely on other C_{9} components. Once assigned, a system of ODEs could be defined in these terms (Additional file 1: Equations S6–S14) and integrated numerically to afford our final simulations.

We first assigned rate-law expressions for the remaining C_{9} components that did not involve O_{2}, namely *R*_{vp} and *R*_{isu}. The expressions *k*_{vp}⋅[F3] and *k*_{isu⋅}[FM] were sufficient to simulate *R*_{vp-dat} and *R*_{isu-dat} with acceptable fidelity. The simple rate-law *R*_{23-sim} = *k*_{23⋅}[F2] was unable to simulate the data. The problem was that cells grown under low-iron conditions have an unusually high concentration of NHHS Fe^{II}, only a small percentage of which can be assigned to *FM* in mitochondria. Under these conditions, it seemed unlikely that this Fe^{II} could be cytosolic, as there should be low concentrations of [C] (as given in Table 3). The only remaining option (in our model) that could account for the extra Fe^{II} was component *F2* in vacuoles. As cells become iron-sufficient, this effect disappears as [F3] increases. We presumed that the extra *F2* converted into *F3* under these conditions. To coordinate this behavior with increasing cellular iron-sufficiency, we incorporated a Reg_{+FS} function into the *R*_{23-sim} rate-law expression, as we have done previously [17]. In summary, the following rate-law expressions were used in solving the C_{9} model.

The Reg_{+FS} function is the parenthetical term associated with *R*_{23} in Eq. (24). It may be viewed as a valve that regulates the rate by which *F2* converts to *F3*. The valve is almost fully open (value near to 1) when [FS] is much greater than the set-point concentration [FS]_{sp}, and is nearly closed (value near to 0) when [FS] is much less than [FS]_{sp}.

### Effect of O_{2}

Oxygen plays a critical role in the C_{9} model as it reacts with *FM* to generate *MP*. O_{2} is constantly diffusing into the matrix (in accordance with rate *R*_{O2}) and is reduced to H_{2}O by cytochrome c oxidase on the inner membrane. We used [FS] as a proxy for oxidase activity such that the rate of respiration (*R*_{res}) was assumed to be proportional to both [FS] and [O_{2}]. Collectively, these processes determine the dissolved O_{2} concentration in the matrix, as described by

Under an expanded steady-state condition

R_{O2} was presumed to be proportional to the *difference* in the O_{2} concentration in the cytosol (called [O_{2}]_{cyt} – assumed to be fixed at 100 μM) and the concentration of O_{2} in the matrix ([O_{2}]). With rate-law expressions included, (26) becomes

Rearrangement yields

Since all numbers in (28) are positive, the term in the numerator serves to increase [O_{2}] while those of the denominator serve to decrease it. [FM], [FS], and α_{cell} are controlled by other aspects of the model, and so those parameters were not altered in order to generate the behavior desired for [O_{2}] vs [N] in WT vs. ΔΔ cells. This behavior was essentially controlled by the three unassigned parameters, *k*_{O2}, *k*_{res}, and *k*_{mp} contained in Eq. (28).

One major objective for this study was to explain how ΔΔ mitochondria transition from a diseased state (dominated by nanoparticles, *MP* in the model) when cells are grown in low-iron media, to a healthy state (dominated by ISCs and heme centers, *FS* in the model) when they are grown in high-iron media. We also wanted WT mitochondria to be healthy regardless of the iron concentration in the growth medium. The major molecular-level differences between WT and ΔΔ cells are the rates at which iron enters mitochondria (*R*_{mit}) and cells (*R*_{cyt}), and the growth rates – all of which have been set by solving the simpler versions of the model. The key to achieving the desired behavior, according to our model, was to vary the concentration of O_{2} in the matrix. In WT mitochondria, [O_{2}] should be low at all [N] whereas in ΔΔ mitochondria, [O_{2}] should be high at low [N] and low at high [N]. We needed to generate an abrupt decline of [O_{2}] in ΔΔ mitochondria as [N] increases while keeping [O_{2}] low in WT mitochondria at all [N]. *And* we needed to make this happen only by adjusting *k*_{O2}, *k*_{res}, and *k*_{mp}.

The [O_{2}] concentration in the matrix has not been measured directly. We estimated [O_{2}] to be in the ballpark of 1–10 μM for iron-sufficient WT mitochondria as this value is similar to the *K*_{M} for O_{2} reduction by cytochrome c oxidase [45]. We had [MP] vs. [N] data that could be used to help optimize these parameters (especially *k*_{mp}), but they were insufficient.

We also considered the known behavior of Yfh1Δ cells, which we have explained using a similar model [17]. Yfh1Δ mitochondria contain excessive levels of nanoparticles. The previous model explained the excessive nanoparticles as being due to a lack of *FS* (respiratory complexes), which allows for O_{2} to diffuse into the matrix, react with *FM*, and generate *MP*. This behavior (obtained by setting *R*_{isu} = 0) provided another constraint on possible solutions for the current problem. Another consideration was that respiring WT cells grown at all [N] do not accumulate *MP* in their mitochondria.

After extensive trials, we obtained values of *k*_{O2}, *k*_{res}, and *k*_{mp} (listed in Table 2) that generated the best overall behavior. However, despite our efforts to satisfy all of these constraints, we could not completely eliminate the formation of *MP* in WT mitochondria while also having *MP* accumulate at high levels in Yfh1Δ mitochondria. Two additional changes were required to do this, namely increasing *k*_{mit} of WT cells 2-fold and decreasing *k*_{mit} of ΔΔ cells 1.3-fold, both relative to the values obtained by solving the simpler C_{3} model. The adjustment of *k*_{mit} in ΔΔ cells was minor whereas the adjustment for WT cells implies that the concentration of iron in WT mitochondria is actually 2-fold higher than given by the data used for simulations.

In summary, solving the C_{9} model while achieving the desired behavior with O_{2} required that we increase *k*_{mit} of WT cells 2-fold and decrease k_{mit} of ΔΔ cells 1.3-fold relative to values obtained in the C_{3}/C_{4} models. This explains the different values of *k*_{mit} in Table 2. The faster import rate from cytosol into the mitochondria for the C_{9} model also caused a slight decline of cytosolic iron ([C] and [CIA]).

### Final optimization and sensitivity analysis

Once each parameter was optimized individually as described above, we re-optimized the entire system by changing one component at a time while holding the others fixed, as described above. For the C_{1}, C_{3}/C_{4}, and C_{9} model variants, the best-fit ERR values were 0.32, 0.39, and 0.72, respectively. A sensitivity analysis was performed for each parameter by taking the average of the ±1% ERR values, and normalizing the average to the optimal ERR for that parameter [17]. This procedure is calculated using the equation

Highest sensitivity values (Table 2) indicate parameters with the greatest impact on the overall fit of the model; *nisu* (Hill coefficient for ISC assembly), α_{max} (growth rate), *K*_{N} (the Michaelis-Menten constant for nutrient iron import), and *f*_{cyt} (fractional cytosol volume) were the most sensitive.

Simulation plots showing the concentrations of each iron-containing component of the C_{9} model (except for nanoparticles) is shown in Fig. 4. As expected, simulated concentrations of most components increased as the nutrient iron concentration increased. Simulated concentrations of cytosolic and vacuolar components in ΔΔ cells were higher than in WT cells, whereas the simulated concentrations of mitochondrial components *FS* and *FM* in ΔΔ cells were lower than in WT cells. Vacuolar iron is dominated by *F2* under iron-deficient conditions and by *F3* under iron-sufficient ones.

### Effect on O_{2} on mitochondrial nanoparticles

Simulations of mitochondrial O_{2} and nanoparticle concentrations are shown in Fig. 5. There is a nonlinear effect that simulates the observed behavior of ΔΔ mitochondria. Mitochondria from iron-deficient ΔΔ cells contain mostly nanoparticles and are responsible for the slow-growth defect. However, mitochondria from these cells recover when ΔΔ cells are grown in iron-sufficient medium. The plot simulates this recovery. As [N] increases, [O_{2}] levels decline because increasing concentrations of respiratory complexes *FS* prevent O_{2} from diffusing into the matrix and reacting with *FM*. This allows more *FS* to be made which allows even less O_{2} into the matrix. This vicious cycle leads to the observed nonlinear behavior. The same behavior is observed for the formation of nanoparticles (Fig. 5, panel b). Other traces to either side of the best-fit [O_{2}] trace represent the effect of increasing/decreasing one parameter while keeping all others fixed. Since the percentage change for each parameter was the same, the parameters that influence the shape of the plot more dramatically are located on the extremes; those that have little influence on the plot are located near the central optimized curve. A similar effect is obtained by lowering *R*_{isu} (Fig. 5, panel c) which simulates the effect of lowering the Yfh1 concentration in yeast mitochondria (or the frataxin concentration in human mitochondria). WT mitochondria do not exhibit this nonlinear effect because they can exclude O_{2} from the matrix at all [N] considered.

## Conclusions

### Comparison to previous model

The model developed here represents a major advance relative to our previous model [17]. Both simulate iron import and trafficking in a growing yeast cell, both include three regions (cytosol, mitochondria, and vacuoles), and both import a single nutrient iron form *N*. The major difference between the two is their level of complexity, method of optimization, and predictive power. The previous model included ~ 38 adjustable parameters (Table S2 of [17]) whereas the current C_{9} model includes only 25 (Table 2). The previous model was optimized by guessing an initial set of values and minimizing an error function. The current model was optimized using the multi-tiered approach detailed above. Perhaps the most important difference is that the current model predicts the nonlinear O_{2}-dependent behavior described above while the previous model does not. Our current model was solved at different levels of complexity. We solved the simpler variants first, and discovered that the parameters obtained could be transferred to the more complex variants. This multi-tiered strategy was helpful because the parameters obtained by fitting the simpler models used more reliable data. Although the heuristics used to find optimal parameter sets do not mathematically guarantee that they globally minimize the error function, our approach yielded an acceptable solution that is compatible with known biology.

Another strategic difference in modeling approaches was that we excluded all but one Reg function in the current model. This made the current model more responsive to changing parameters and allowed better comprehension of inherent behavior.

In the end, only four parameters differed between ΔΔ and WT simulations, namely *R*_{cyt-max}, *k*_{mit}, *k*_{vp}, and *K*_{α}. *All other assigned parameter values were identical between the two genetic strains.* The ability of the model to reproduce ΔΔ and WT behavior with such few differences in terms of modeling parameters is remarkable. Moreover, we can easily rationalize why at least half of these parameters should be different. A 4.6-fold reduction of *k*_{mit} for ΔΔ cells makes sense because Mrs3/4, the high-affinity importers into mitochondria, are deleted in this strain. *R*_{cyt-max} is 2-times higher for ΔΔ cells because iron is dysregulated in these cells so expression of the Ftr1/Fet3 complex on the plasma membrane should be higher. Explaining why *K*_{α} should be 30-times higher in ΔΔ cells is more difficult. *K*_{α} is a *K*_{M}-like parameter which reflects the sensitivity of the growth rate to changes in the nutrient iron concentration [N]. For some reason, the growth of iron-deficient ΔΔ cells is 30 times less sensitive to increases in [N] than are comparable WT cells. Perhaps this reflects difficulties in flowing sufficient iron into iron-deficient ΔΔ mitochondria to support robust respiratory cell growth. Why *k*_{vp} is 2-fold higher in ΔΔ cells is even more difficult to explain; it implies that the rate of vacuolar nanoparticle formation is faster in ΔΔ cells than in WT cells. The actual mechanism of vacuolar nanoparticle formation is undoubtedly more complicated than is represented in our current model. However, it is a tribute to the model that it has the ability to highlight this effect.

### Predictive power of the model

Mathematical models *might* have predictive power, but this is not guaranteed. This ability is related to how closely the assumed mechanism and assigned kinetic parameters correspond to reality. We have attempted to make our model predictive by keeping it simple and well-grounded experimentally. This was a challenge given the complexity of the process under investigation and the limited amount of relevant data available.

Our model can be used to predict the effect of O_{2} on iron metabolism in yeast cells. It predicts that the iron in mitochondria of Yfh1-deficient cells that have been grown under micro-aerobic (low O_{2}) conditions should predominantly be *FM* (i.e. NHHS Fe^{II}). We are currently examining a Yfh1-deficient strain of yeast, and found that NHHS Fe^{II} rather than nanoparticles are indeed observed in these cells (data not shown). Our model also predicts that O_{2} should not affect vacuolar iron (it should still be present mainly as *F3* (Fe^{III}) under micro-aerobic conditions). However, this prediction is not realized by our current experiments (data not shown), highlighting a deficiency in this particular aspect of the model. We believe that this iterative approach of prediction→testing→remodeling will gradually yield major new insights in understanding iron import, trafficking, and regulation in eukaryotic cells. We are currently using this approach in our studies of the Yfh1-deficient strain.

The same strategy could be applied to model the import and trafficking of any micronutrient within any growing eukaryotic cell, from yeast to human cells. The concentration of the nutrient (or its derivatives) in whole cells and in various organelles and cytosol should be known as should exponential growth rates. Obvious candidates include other metals such as Cu, Mn, Zn, Mo, Co. The same approach could be used to examine the import and trafficking of Pt anticancer drugs into growing human cells. A better understanding of how such drugs are trafficked intracellularly might provide new insights for treating cancer.

## Notes

The abbreviations used are: CD, central doublet; CIA, cytosolic iron-sulfur assembly complex; ISC, iron-sulfur cluster; MB, Mössbauer; NHHS, nonheme high-spin; ODE, ordinary differential equation; ROS, reactive oxygen species; WT, wild-type. C

_{1}, C_{3}, C_{4}, and C_{9}refer to model variants; C_{1}includes one component called*Fe*_{Cell}; C_{3}contains three components called*Fe*_{cyt},*Fe*_{mit}, and*Fe*_{vac}; C_{4}contains four components called*C*,*CIA*,*Fe*_{mit}and*Fe*_{vac}; C_{9}contains nine components called*C*,*CIA*,*F2*,*F3*,*VP*,*FM*,*FS*,*MP*, and*O*_{2}. Parameters*f*_{cyt},*f*_{mit}, and*f*_{vac}are volume fractions representing the ratios V_{cyt}/V_{cell}, V_{mit}/V_{cell}, and V_{vac}/V_{cell}, respectively. ΔΔ is the Mrs3Δ/Mrs4Δ double-deletion strain.*α*_{cell}is the growth rate of the cell, and*N*is the nutrient iron species that is imported into the cell.

## References

Chen WW, Niepel M, Sorger PK. Classic and contemporary approaches to modeling biochemical reactions. Genes Dev. 2010;24:1861–75.

Anbumathi P. Chemical kinetics based modelling approaches and tools for investigating complex biological systems. Asian J Chem. 2014;26:2511–7.

Mizera A, Czeizler E, Petre I. Methods for biochemical model decomposition and quantitative submodel comparison. Israel J Chem. 2011;51:151–64.

Outten CE, Albetel AN. Iron sensing and regulation in

*Saccharomyces cerevisiae*: ironing out the mechanistic details. Curr Opin Microbiol. 2013;16:662–8.Lill R, Dutkiewicz R, Freibert SA, Heidenreich T, Mascarenhas J, Netz DJ, Paul VD, Pierik AJ, Richter N, Stuempfig M, Srinivasan V, Stehling O, Mühlenhoff U. The role of mitochondria and the CIA machinery in the maturation of cytosolic and nuclear iron-sulfur proteins. Eur J Cell Biol. 2015;94:280–91.

Philpott CC, Ryu MS. Special delivery: distributing iron in the cytosol of mammalian cells. Front Pharmacol. 2014;5:173. https://doi.org/10.3389/fphar.2014.00173.

Li LT, Kaplan J. A mitochondrial-vacuolar signaling pathway in yeast that affects iron and copper metabolism. J Biol Chem. 2004;279:33653–61.

Raguzzi F, Lesuisse E, Crichton RR. Iron storage in

*Saccharomyces cerevisiae*. FEBS Lett. 1988;231:253–8.Cockrell AL, Holmes-Hampton GP, McCormick SP, Chakrabarti M, Lindahl PA. Mössbauer and EPR study of iron in vacuoles from fermenting

*Saccharomyces cerevisiae*. Biochemistry. 2011;50:10275–83.Omholt SW, Kefang X, Andersen O, Plahte E. Description and analysis of Switchlike regulatory networks exemplified by a model of cellular Iron homeostasis. J Theor Biol. 1998;195:339–50. https://doi.org/10.1006/jtbi.1998.0800.

Mobilia N, Donze A, Moulis JM, Fanchon E. A model of the cellular iron homeostasis network using semi-formal methods for parameter space exploration. EPTICS. 2012;92:42–57. https://doi.org/10.4204/EPTCS.92.4.

Mobilia N, Donze A, Moulis JM, Fanchon E. Producing a set of models for the iron homeostasis network. EPTICS. 2013;125:92–8. https://doi.org/10.4204/EPTCS.92.4.

Chifman J, Kniss A, Neupane P, Williams I, Leung B, Deng Z, Mendes P, Hower V, Torti FM, Akman SA, Torti SV, Laubenbacher R. The core control system of intracellular iron homeostasis: a mathematical model. J Theor Biol. 2012;300:91–9.

Chifman J, Arat S, Deng ZY, Lemler E, Pino JC, Harris LA, Kochen MA, Lopez CF, Akman SA, Torti FM, Torti SV, Laubenbacher R. Activated oncogenic pathway modifies Iron network in breast epithelial cells: a dynamic modeling perspective. PLoS Comput Biol. 2017;13:e1005352. https://doi.org/10.1371/journal.pcbi.1005352.

Mitchell S, Mendes P. A computational model of liver iron metabolism. PLoS Comput Biol. 2013;9:e1003299.

Achcar F, Camadro JM, Mestivier D. A Boolean probabilistic model of metabolic adaptation to oxygen in relation to iron homeostasis and oxidative stress. BMC Syst Biol. 2011;5:51 http://www.biomedcentral.com/1752-0509/5/51.

Wofford JD, Lindahl PA. Mitochondrial iron-sulfur-cluster activity and cytosolic iron regulate iron traffic in

*Saccharomyces cerevisiae*. J Biol Chem. 2015;290:26968–77.Chakrabarti M, Lindahl PA. The utility of Mössbauer spectroscopy in eukaryotic cell biology and animal physiology. In: Rouault TA, editor. Iron-sulfur clusters, vol. 4. Berlin: Walter de Gruyter Publisher; 2014. p. 49–75.

Foury F, Roganti T. Deletion of the mitochondrial carrier genes MRS3 and MRS4 supresses mitochondrial iron accumulation in a yeast frataxin-deficient strain. J Biol Chem. 2002;277:24475–83.

Mühlenhoff U, Hoffmann B, Richter N, Rietzschel N, Spantgar F, Stehling O, Uzarska MA, Lill R. Compartmentalization of iron between mitochondria and the cytosol and its regulation. Eur J Cell Biol. 2015;94:292–308.

McComick SP, Moore MJ, Lindahl PA. Detection of labile low-molecular-mass transition metal complexes in mitochondria. Biochemistry. 2015;54:3442–53.

Yoon H, Zhang Y, Pain J, Lyver ER, Lesuisse E, Pain D, Dancis A. Rim2, a pyrimidine nucleotide exchanger, is needed for iron utilization in mitochondria. Biochem J. 2011;440:137–46.

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. J Biol Chem. 2004;279:29512–8.

Poor CB, Wegner SV, Li HR, Dlouhy AC, Schuermann JP, Sanishvili R, Hinshaw JR, Riggs-Gelasco PJ, Outten CE, He C. Molecular mechanism and structure of the Saccharomyces cerevisiae iron regulator Aft2. Proc Natl Acad Sci U S A. 2014;111:4043–8.

Courel M, Lallet S, Camadro JM, Blaiseau PL. Direct activation of genes involved in intracellular iron use by the yeast iron responsive transcription factor Aft2 without its paralog Aft1. Mol Cell Biol. 2005;25:6760–71.

Martelli A, Puccio H. Dysregulation of cellular iron metabolism in Friedreich ataxia; from primary iron-sulfur cluster deficit to mitochondrial iron accumulation. Front Pharmacol. 2014;5:130. https://doi.org/10.3389/fphar.2014.00130.

Lesuisse E, Santos R, Matzanke BF, Knight SAB, Camadro JM, Dancis A. Iron use for haeme synthesis is under control of the yeast frataxin homologue (Yfh1) human Mol. G E N. 2003;12:879–89.

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:9556–68.Urbanowski JL, Piper RC. The iron transporter fth1 forms a complex with the Fet5 iron oxidase and resides on the vacuolar membrane. J Biol Chem. 1999;274:38061–70.

Muhlenhoff U, Richhardt N, Gerber J, Lill R. Characterization of iron-sulfur protein assembly in isolated mitochondria - a requirement for ATP, NADH, and reduced iron. J Biol Chem. 2002;277:29810–6.

Gardner PR, Nguyen DD, White CW. Aconitase is a sensitive and critical target of oxygen poisoning in cultured mammalian cells and in rat lungs. Proc Natl Acad Sci U S A. 1994;91:12248–52.

Ollagnier-de Choudens S, Sanakis Y, Hewitson KS, Roach P, Baldwin JE, Münck E, Fontecave M. Iron-sulfur Center of Biotin Synthase and Lipoate Synthase. Biochemistry. 2000;39:4165–73.

López-Torrejón G, Jiménez-Vicente E, Buesa JM, Hernandez JA, Verma HK, Rubio LM. Expression of a functional oxygen-labile nitrogenase component in the mitochondrial matrix of aerobically grown yeast. Nat Commun. 2016;7:1–6. https://doi.org/10.1038/ncomms11426.

Moore MJ, Wofford JD, Dancis A, Lindahl PA. Recovery of mrs3Δmrs4Δ

*Saccharomyces cerevisiae*cells under iron-sufficient conditions, and the role of Fe_{580}. Biochemistry. 2018;57:672–83.Kwok EY, Severance S, Kosman DJ. Evidence for iron channeling in the Fet3p-Ftr1p high-affinity iron uptake complex in the yeast plasma membrane. Biochemistry. 2006;45:6317–27.

Andreini C, Banci L, Rosato A. Exploiting bacterial operons to illuminate human iron-sulfur proteins. J Proteome Res. 2016;15:1308–22.

Yamaguchi M, Namiki Y, Okada H, Mori Y, Furukawa H, Wang J, Ohkusu M, Kawamoto S. Structome of

*Saccharomyces cerevisiae*determined by freeze-substitution and serial ultrathin-sectioning electron microscopy. J Electron Microscopy. 2011;60:337–51.Wei DG, Jacobs S, Modla S, Zhang S, Young CL, Cirino R, Caplan J, Czymmek K. High-resolution three-dimensional reconstruction of a whole yeast cell using focussed-ion beam scanning electron microscopy. Biotechniques. 2012;53:41–8.

Chan YHM, Marshall WF. Organelle size scaling of the budding yeast vacuole is tuned by membrane trafficking rates. Biophys J. 2014;106:1986–96.

Stevens BJ. Variation in number and volume of mitochondria in yeast according to growth conditions: study based on serial sectioning and computer graphic reconstruction. Biol Cell. 1977;28:37–56.

Holmes-Hampton GP, Jhurry ND, McCormick SP, Lindahl PA. Iron content of Saccharomyces cerevisiae cells grown under iron-deficient and iron-overload conditions. Biochemistry. 2013;52:105–14.

Aras S, Bai M, Lee I, Springett R, Hüttemann M, Grossman LI. MNRR1 (formerly CHCHD2) is a bi-organellar regulator of mitochondrial metabolism. Mitochondrion. 2015;20:43–51.

Park J, McCormick SP, Cockrell AL, Chakrabarti M, Lindahl PA. High-spin ferric ions in

*Saccharomyces cerevisiae*vacuoles are reduced to the ferrous state during adenine-precursor detoxification. Biochemistry. 2014;53:3940–51.Petrat F, deGroot H, Rauen U. Subcellular distribution of chelatable iron: a laser scanning microscopic study in isolated hepatocytes and liver endothelial cells. Biochem J. 2001;356:61–9.

Krab K, Kempe H, Wikstrom M. Explaining the enigmatic K(M) for oxygen in cytochrome c oxidase: a kinetic model. Biochim Biophys Acta. 2011;1807:348–58.

## Acknowledgements

We thank Ivan V. Surovtsev (Yale University) for suggesting Eq. (12).

### Funding

This work was supported by the National Institutes of Health (GM127021), the National Science Foundation (MCB-1817389) and the Robert A. Welch Foundation (A1170).

### Availability of data and materials

All data generated or analysed during this study are included in this published article and its supplementary information file. The *Mathematica* code used in simulations will be provided to researchers upon request 1 year after publication.

## Author information

### Authors and Affiliations

### Contributions

JDW helped design the model. He wrote the computer code, optimized the model, and prepared figures. PAL derived the math equations, helped design and optimize the model, offered advice, and wrote much of the manuscript. Both authors approved the final version of the manuscript.

### Corresponding author

## Ethics declarations

### Ethics approval and consent to participate

Not applicable.

### Consent for publication

Not applicable.

### Competing interests

The authors declare that they have no competing interests.

### Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Additional file

### Additional file 1:

This file includes ODEs for the different model variants. (DOCX 61 kb)

## Rights and permissions

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

## About this article

### Cite this article

Wofford, J.D., Lindahl, P.A. A mathematical model of iron import and trafficking in wild-type and Mrs3/4ΔΔ yeast cells.
*BMC Syst Biol* **13**, 23 (2019). https://doi.org/10.1186/s12918-019-0702-2

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/s12918-019-0702-2

### Keywords

- Ordinary differential equations
- Nanoparticles
- Iron-sulfur clusters
- Mitochondria
- Vacuoles
- Cytosol
- Kinetics
- Oxygen