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

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 FeII mitochondrial pool and the necessary exclusion of O2 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. Electronic supplementary material The online version of this article (10.1186/s12918-019-0702-2) contains supplementary material, which is available to authorized users.


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 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 Fig. 1 Strategy for optimizing a model of nutrient iron import, trafficking and regulation in growing eukaryotic cells. Top panel: C 1 model in which all iron in the cell is treated as a single species and the cell is considered to be homogeneous. Middle panel: C 3 model in which the cell is divided into three regions and each region is assumed to contain a single type of iron species. C 4 model includes the reaction forming CIA. Bottom panel: C 9 model in which the cell remains divided into 3 regions but the number of iron-containing species is expanded to 8 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 mitochondriai.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 ironreplete ΔΔ 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 datai.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 → R cell 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, Table 1 Growth rates, iron concentrations, and import rates in ΔΔ and WT cells grown under different nutrient conditions. [N] refers to the μM concentration of iron in the respiring medium, as described [34]. The untreated medium was assumed to contain 1 μM of endogenous iron. For each entry, the top number is datum or data-based estimates (R …-dat ) and the bottom number is the corresponding simulated value (R …-sim [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 ironimporter 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.   Table 1. For data points in this figure and in Fig. 3, the errors estimated in [34] from two determinations are within the marks. Overall errors may be higher 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 → R cyt Fe cyt ). Some Fe cyt is imported into mitochondria ( Fe cyt → R mit Fe mit ) and some into vacuoles ( Fe cyt → R vac 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) [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 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  . 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 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 ironsufficient, 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

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 ] cytassumed 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 molecularlevel differences between WT and ΔΔ cells are the rates at which iron enters mitochondria (R mit ) and cells (R cyt ), and the growth ratesall of which have been set by solving the simpler versions of the model. . 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)

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   N]). This is the central plot. Other plots in Panel a, on either side of the optimized plot, show that only certain parameters affect curve shape and position. For the other traces, the color-coded parameters shown on the side were altered ±10% of their optimized values while holding all other parameters fixed. Changing other parameters yields traces (e.g. k 23 in the white dashed line) that had no effect on the plots and so the trace overlaid the optimal trace in the center. Panel c is a plot of [MP] vs. R isu-max , the maximum rate of FS formation. Low R isu-max values simulate the slow rate of ISC assembly in yfh1Δ cells, while higher values reflect WT conditions 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.

Endnotes
1 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.