Dysregulation of iron metabolism in cancer is well documented and it has been suggested that there is interdependence between excess iron and increased cancer incidence and progression. In an effort to better understand the linkages between iron metabolism and breast cancer, a predictive mathematical model of an expanded iron homeostasis pathway was constructed that includes species involved in iron utilization, oxidative stress response and oncogenic pathways. The model leads to three predictions. The first is that overexpression of iron regulatory protein 2 (IRP2) recapitulates many aspects of the alterations in free iron and iron-related proteins in cancer cells without affecting the oxidative stress response or the oncogenic pathways included in the model. This prediction was validated by experimentation. The second prediction is that iron-related proteins are dramatically affected by mitochondrial ferritin overexpression. This prediction was validated by results in the pertinent literature not used for model construction. The third prediction is that oncogenic Ras pathways contribute to altered iron homeostasis in cancer cells. This prediction was validated by a combination of simulation experiments of Ras overexpression and catalase knockout in conjunction with the literature. The model successfully captures key aspects of iron metabolism in breast cancer cells and provides a framework upon which more detailed models can be built.
This section provides biological background about iron metabolism and its connection to some oncogenic pathways.
We describe a discrete dynamic model of the network in Fig 1, based on an encoding of the regulatory logic for each node through a “logical” update rule. This type of model is qualitative, in the sense that each species can assume a finite set of states rather than quantitative concentrations of molecular species. For this study, we adopted a ternary logic, an extension of Boolean logic. Our choice of ternary logic was motivated by the fact that iron levels cannot be viewed as either ON = 1 or OFF = 0. The iron homeostasis pathway is the major focus of our study and both low and high levels of iron are detrimental, so that it is tightly controlled. Additionally, IRP2 at both low and high activity levels does alter the iron pathway [68]. With only two states it would not be clear when IRP2 operates at low activity levels, as it would be represented the same way as normally active protein. For our model to be able to differentiate between iron homeostasis (normal levels of iron) and low/high iron levels as well as activities/concentration levels of various proteins, we chose to represent each species by three levels: low, normal and high. In the language of logical models the state of a particular species is described by 0 if the species is low/inactive, by 1 if at normal/intermediate activity, and by 2 if high/active. In analogy to the Boolean formalism, we can compute the future state of a species at time step t + 1 using the states of other species at time step t. Fundamental OR and AND gates for two species X and Y are defined as max{X, Y} and min{X, Y}, respectively, where X, Y ∈ {0,1, 2}. To differentiate from the Boolean OR and AND gates, we denote these gates by Max and Min, respectively. The NOT gate (denoted here by X ¯) is defined by inverting the input, i. e. , leaving 1 unchanged and inverting 0 and 2. For a concrete example, consider heme in Fig 1. It is produced through ALA synthase (ALAS1) but inhibited by HO-1. Then the logical function (update rule) that predicts how much heme is present at time t + 1 can be computed as follows: heme (t+1) =Min (ALAS1 (t), HO−1¯ (t) ). This means that, if HO-1 was 0 (low) and ALAS1 was 2 (high) at time t, then heme will be 2 (high) at time t + 1. Based on the biological knowledge described in the previous section, we translated the interactions of the normal cell network into logical functions (see Table 1). One caveat about logical models that is not present for Boolean models is that species can change for example from a low state to a high one in one time step, skipping intermediate concentrations. This is biologically unrealistic. Thus, to address the continuity issue we have also implemented a methodology commonly used for logical models that takes into account the previous state of the regulated species (see [69] for details). For purposes of simulation, we converted the logical rules into polynomial functions to obtain a so-called polynomial dynamical system (PDS). A description of the construction of the PDS and the entire system can be found in the Materials and Methods section and in the supplemental file S1 PDS, respectively. To analyze the dynamic properties of the model we simulated the entire state space and computed the basins of attraction of the system. For this purpose, we used an encoding of the model as a polynomial dynamical system, as described above, and customized scripts written in Perl and Python (see Materials and Methods section). The size of the model’s state space is 324 = 282,429,536,481, where 24 is the number of species in the network and 3 is the number of states (low, medium, high) per species. We employed a synchronous update schedule for the species in the network; all species were updated simultaneously based on the states of their input species at the previous time step. Each state leads to another state, eventually converging to a steady state or a limit-cycle (a set of recurring states), which are called attractors. A collection of initial states that lead to a particular attractor is termed the basin of attraction. Under this scheme, each state belongs to the basin of attraction of only one attractor: a point attractor (steady state) or a cycle attractor (limit-cycle). These attractors correspond to different phenotypes in the biological context and can describe various behaviors of the system such as homeostasis. We simulated the normal cell model and also investigated the long-term behavior of this model under different conditions, namely, the effects of knockout (k/o) or overexpression (o/e) of one or more species. To simulate these experimental conditions, we set the update rule for a particular species to a constant equal to 0 or 2, respectively. In other words, regardless of the input (regulators), the species of interest will always stay at the chosen level. Our results are summarized in Fig 2, which shows all the species and their long-term behavior. Simulations were performed by exhaustively enumerating the transitions of the model on all possible 324 states. The normal cell network has no cycle attractors and reaches a unique stable steady state (point attractor) indicating that all species are at their respective normal levels regardless of the initial starting state (Fig 2 top line of the heat map labeled Normal).
It is well-known that iron metabolism in breast epithelial cells is differentially regulated as cells transition to malignancy. Determining the causes for this altered phenotype is complicated by the complexity of iron regulation and its connection to several other processes, such as response to oxidative stress and changes in iron consumption [80], as well as crosstalk with oncogenic pathways. Integrating these different influences on the iron phenotype in normal and malignant cells can benefit greatly from a systematic approach through dynamic mathematical modeling, beyond the network approach taken in [80]. The model presented here is a first step toward a comprehensive understanding of the iron phenotype of cells as it changes in breast cancer. We have chosen to construct a qualitative model of an intracellular iron network (Fig 1) to capture its fundamental dynamic features (attractors). The main reason for our choice of modeling platform is that our current knowledge of the kinetics involved in these different processes as well as mechanisms underlying these complex reactions is very limited, so that a quantitative model, such as a system of ordinary differential equations is more challenging to construct. We have validated our model using both experimental data and information from the literature not used in model construction. In particular, we have experimentally validated the model prediction that IRP2 overexpression in the normal cell network only alters the iron homeostasis pathway, leaving the other model components unchanged. Also, our model agrees with the current literature that overexpression of mitochondrial ferritin (Ftmt) increases both IRPs and TfR1, decreases cytosolic Ft and reduces cytosolic and mitochondrial iron pools [71]. In addition, we have shown that shutting down trafficking of iron into the mitochondria, together with Ras overexpression and Cat reduced bioactivity, does lead to the observed cancer phenotype of the iron homeostasis pathway. However, it might be possible that further refinements of the model can lead to the required phenotype by altering only the oncogenic pathway. Not all known information about the normal and cancer phenotypes can be captured by the model, however. This is likely due to the fact that some key features of this system are not represented completely, such as an iron-sensing regulator in the mitochondria and iron-sulfur cluster (ISC) synthesis. It has been suggested that frataxin, a nuclear-encoded mitochondrial protein, may act as an iron-sensing regulator and even function as a switch between heme and ISC synthesis [81–83]. At this stage, one cannot determine whether it is frataxin or some other iron sensor/regulator, but we have suggested the possibility of a mitochondrial iron-sensing node in our current model (depicted as a question mark in Fig 4 (a) ). This adjusted normal cell model also reaches a unique stable steady state agreeing with our model discussed in the Results section (Fig 4 (b) ). Additionally, we simulated two more models using the following perturbations: (i) knockout of the sensor node and (ii) overexpression of Ras and of the sensor node, and low bioactivity of CAT (2nd and 3rd rows in Fig 4 (b) ). Interestingly, the sensor node k/o model agrees with experimental data that in frataxin k/o mice heme is decreased, TfR1 is upregulated and iron uptake via Mfrn is increased, leading to cytosolic iron-deficiency and mitochondrial iron overload [84,85]. This strongly indicates that there is a sensor/regulator, and thus further refinements of the model can provide insight into mitochondrial iron regulation and utilization, and potentially suggest new experiments that can validate new connections. The latter model produced the same cancer phenotype of the iron homeostasis pathway (see Eq (2) ) and also implied that cancer cells have reduced heme biosynthesis. Furthermore, we note that the latter model allows Ftmt and ALAS1 from the iron utilization pathway to have high expression levels (compare 9th row in Fig 2 to row 3 in in Fig 4 (b) ). While we do not have much evidence about Ftmt in cancer, there are some studies about ALAS1 in lung cancer. It was found that ALAS1 protein levels were substantially increased in non-small-cell lung cancer cells compared to normal cells [86]. This suggests the possibility to expand the cancer phenotype of the iron homeostasis pathway to the iron utilization pathway. Of course, one can simulate a model by setting various proteins to their respective observed levels, but then we gain no information about the drivers that change iron metabolism in cancer. Ideally, we would like to include other pathways implicated in breast cancer to capture different molecular subtypes of breast cancer and iron cancer phenotypes associated with them.
We begin by defining a set of rules that describe various relations between molecular species, from which we then build the entire model. If species X is, inducing species Y (X → Y) or species X is inhibiting species Y (X ⊣ Y) then we represent these relationships via a transition table as depicted in Table 2. Notice that inhibition in Table 2 is just a logical NOT gate, denoted here by X ¯. The other two fundamental gates, OR and AND, for two species X and Y regulating species Z (X → Z ← Y), are defined as max{X, Y} and min{X, Y} respectively, for X, Y ∈ {0,1, 2}, and denoted here by Max and Min. We can express the above gates as polynomials over a finite field on three elements, F 3. If we limit the exponent of each variable in a polynomial to be less than or equal to 2, then one can show that any logical rule constructed from these three operations has a unique polynomial representation, using x¯=2+2xMax (x, y) =x2y2+x2y+xy2+2xy+x+yMin (x, y) =2x2y2+2x2y+2xy2+xy. (3) One can check that polynomials given by Eq (3) agree with definitions of fundamental gates as described in the paragraph above, e. g. , max{1,2} = 2 and Max (1,2) = (1) 2 (2) 2 + (1) 2 (2) + (1) (2) 2 + 2 (1) (2) + 1 + 2 = 2, where the right-hand side is computed modulo 3. Various adjustments to the strength of a particular regulation can be made by altering entries in the Table 2. For example, it has been suggested that IRP1, when active, contributes less to the regulation of ferritin (Ft) than IRP2 [68] (see Table 3). These tables mean that when IRP2 = 2 (active) it will inhibit Ft, whereas when IRP1 = 2 (active) it will have a lesser effect on Ft. Thus, we can represent regulation of Ft by IRP2 in Table 3 using Eq 3: IRP2¯=2+2⋅IRP2. Now, for IRP1 regulating Ft according to this new adjustment, one can also find a polynomial representing Table 3 (left table). For convenience, whenever we use an adjusted regulation we will place an asterisk (*) in front of the variable inside the logic gate. *IRP2¯=2+2⋅ (IRP1) 2. To match current biological knowledge we have adjusted regulation of IRP1 for TfR1 and Fpn as well [68]. The transition table for Fpn is similar to Ft. For IRP1 regulating TfR1: when IRP1 = 2, then Tfr1 = 1, while, when IRP1 is 0 or 1 then TfR1 is also 0 or 1, respectively. Additionally, we modified regulation of Keap1 by Nrf2 to reflect current literature [44]. For Nrf2 regulating Keap1 we have that when Nrf2 = 0 then Keap1 = 1, while when Nrf2 is 1 or 2 then Keap1 is also 1 or 2, respectively. To make sure that we preserve continuity (i. e. , each species changes at most one unit in one time step), we are going to employ methodology as described in [69]. The underlying reasoning is that this can be accomplished by taking into account the previous state (e. g. , concentration or activity) of the regulated species, in effect adding a self-regulation loop to each network node. The future value of the regulated species under continuity is computed as follows. Let fxi be the update function for xi. To ensure that each variable changes at most 1 unit, define a function h (xi, fxi) for the future value of the variable xi: h (x i, f x i) = x i + 1 if f x i > x i x i if f x i = x i x i - 1 if f x i < x i (4) LIP, heme, and ROS do not undergo self-degradation/self-regulation and hence we do not apply continuity to these species. In order to compute final polynomials, we are going to make use of the following property of finite fields: Remark 0. 1 If h: F p n → F p is any function then there is a polynomial g: F p n → F p so that h (x) = g (x) for all x ∈ F p n. One can find g by using the following formula, g (x) = ∑ c ∈ F p n h (c) ∏ j n (1 - (x j - c j) p - 1), (5) where h (c) is the update function as defined by Eq (4), c is a vector of input variables, and the right-hand side is computed modulo p. All of these logic gates, transition tables describing different strength of regulation and continuity, are then appropriately translated into final polynomial functions over a finite field with three elements. These polynomial functions then form what is called a polynomial dynamical system (PDS) over a finite field. Below, we fully describe a construction of the update function for ferritin (Ft) in our network (see Fig 1). The entire PDS system can be found in the supplemental file S1 PDS. The attractors of the models were found using 2 algorithms: the attractor finder by random sampling (Algo. 1) that is written in Perl and the attractor finder by iterating over all possible states (Algo. 2) that uses a custom written Python package. The codes can be found at https: //github. com/LoLab-VU/LogicalModel. Models that were used for simulations are located in the same directory under NewModels_2015_12_18 and NewModels_2015_8_17 folders. Supporting file S1 Simulations provides additional o/e and k/o simulation results using attractor finder by random sampling. The index for the order of variables is available from row 21 to row 45. After 3,000 random sampling, the basin size of the attractor is specified in the table. The first program requires a model file, the number of states and a sampling size, which is 100,000 here. We randomly selected 100,000 states and stored the attractor states to have a broad perspective on the possible attractors of a model. It was utilized to test which overexpression and knockout models could be potential cancer models. To ensure that we know all attractors of the models of interest, we ran the second program, which requires a model file and number of states. Optional arguments include start and end states and an option to create images of attractor states. The model file is parsed and compiled into an executable function with Cython [87]. We iterated through all possible states of each model (3N), storing only the attractor states. Simulations were performed in parallel using mpi4py [88] running on large cluster computers. Algorithm 1 Pseudo code for attractor finder by random sampling 1: procedure For i in 100,000 ▹Iterate over 100,000 randomly selected states 2: sampled = empty set 3: state = changebase (random (3N) ) 4: while state ∉ sampled do 5: sampled. add (state) 6: state = update (state) ▹Update function is the compiled model 7: state. pop () ▹Returns the last state added to sampled Algorithm 2 Pseudo code for attractor finder by iterating over all possible states 1: procedure For i in 3N ▹Iterate over all possible states 2: sampled = empty set 3: state = changebase (i) 4: while state ∉ sampled do 5: sampled. add (state) 6: state = update (state) ▹Update function is the compiled model 7: state. pop () ▹Returns the last state added to sampled MCF10A, non-tumorigenic immortalized human mammary epithelial cells were obtained from the Wake Forest University Comprehensive Cancer Center Tissue Culture Core facility. The cells were maintained in a suggested condition by ATCC. To overexpress IRP2 in MCF10A cells, the lentiviral vector pSL2-IRP2 [68] was applied. Briefly, MCF10A cells were infected with the concentrated viral particles from pSL2-IRP2 and pLS2 empty vector (as a control). The infection efficiencies for both infections were over 90% based on GFP fluorescence in cells. The cell lysates were harvested for subsequent analysis seven days after infection. Western blotting was performed as previously described [68]. Antibodies: GAPDH (Fitzgerald), TfR1 and c-Myc (Invitrogen), IRP2 and EGFR (Santa Cruz Biotechnology), Keap1 (Cell Signaling Technology), HO-1 and IL-6 (Abcam), ferritin H ([89]).
Iron is required for cellular metabolism and growth, but can be toxic due to its ability to cause high oxidative stress and consequently DNA damage. To prevent damage, all organisms that require iron have developed mechanisms to tightly control iron levels. Dysregulation of iron metabolism is detrimental and can contribute to a wide range of diseases, including cancer. This paper presents a predictive mathematical model of iron regulation linked to iron utilization, oxidative stress, and the oncogenic response specific to normal breast epithelial cells. The model uses a discrete modeling framework to generate novel biological hypotheses for an investigation of how normal breast cells become malignant cells, capturing a breast cancer phenotype of iron homeostasis through overexpression and knockout simulations. The new biology discovered is (1) IRP2 overexpression alters the iron homeostasis pathway in breast cells, without affecting the oxidative stress response or oncogenic pathways, (2) an activated oncogenic pathway disrupts iron regulation in breast cancer cells.
heme
medicine and health sciences
breast tumors
oxidative stress
cancers and neoplasms
oncology
physiological processes
mathematics
homeostasis
algebra
mitochondria
bioenergetics
cellular structures and organelles
polynomials
proteins
breast cancer
chemistry
biochemistry
carcinogenesis
iron
cell biology
post-translational modification
physiology
biology and life sciences
physical sciences
energy-producing organelles
chemical elements
Activated Oncogenic Pathway Modifies Iron Network in Breast Epithelial Cells: A Dynamic Modeling Perspective