Reactive oxygen species (ROS) play a complex role in cellular biology. Initially viewed merely as harmful byproducts of metabolism, ROS are now known to serve additional functions as intracellular regulators of various signaling pathways. At low levels, ROS function as a reactive second messenger mainly by reversible oxidation of key amino acids of target proteins. High ROS levels, by contrast, cause damage to various biomolecules (lipids, proteins, and DNA) and have been linked to pathologies including neurodegeneration (Smith et al., 1997; Giasson et al., 2000), atherosclerosis (Patetsios et al., 2001; Guzik et al., 2006), and renal disease (Nishikawa et al., 2000; Dounousi et al., 2006). ROS play an especially complex role in cancer, with various oncogenes and tumor suppressors influencing, and influenced by, the redox environment of the cell (Meng et al., 2002; Leslie et al., 2003).

It is challenging to study experimentally how endogenous and exogenous sources of ROS are handled by the cell. In addition, the context- and dose-dependent intracellular consequences of ROS can result in confounding observations. ROS has been found to stimulate proliferation in some cell types under certain experimental conditions (Ruiz-Ginés et al., 2000; Nieto et al., 2002) and inhibit proliferation (Qu et al., 2011) or induce apoptosis (Pierce et al., 1991) in others. Such contextual and experimental complexities make it difficult to understand ROS holistically by experimentation alone.

Computational modeling approaches can tackle this problem by simulating the concurrent dynamics of many variables, including those that are difficult to access experimentally (Janes and Lauffenburger, 2013). Our review here covers the handful of models described thus far for ROS production and ROS regulation of signaling pathways (Figure 1). We start with a brief introduction of ROS, the biological processes that generate them, and the signaling pathways that sense ROS and detoxify the cell. We next focus on models that simulate ROS production by the mitochondria, the predominant intracellular source of ROS, and by membrane-bound enzymes, which link extracellular signaling to intracellular ROS production. We also discuss models that simulate ROS regulation of various signaling pathways, giving a broader view of the influence ROS have on intracellular signaling. Finally, we discuss the need for a systems-level analysis of ROS signaling to provide a generalizable framework that accounts for the many downstream cellular effects of ROS.

Reactive oxygen species are chemically reactive species derived from the incomplete reduction of oxygen. Common examples are hydrogen peroxide (H₂O₂), superoxide (O2−), and hydroxyl radicals (HO⋅). For most cell types, ROS production mainly occurs through mitochondrial oxidative phosphorylation. Inherent leakiness of the electron transport chain (ETC), specifically from complex I and complex III, causes some electrons to flow out of the pathway and partially reduce oxygen to the superoxide anion. Other endogenous sources of ROS include the membrane-bound enzyme NADPH oxidase, which produces ROS in response to various ligands, along with other enzymes such as xanthine oxidase, cyclooxygenases, and nitric oxide synthase. To detoxify excess ROS, cells are equipped with antioxidants, including the thiols glutathione and thioredoxin, which are regulated by the transcription factor NRF2 (Itoh et al., 1997; Chorley et al., 2012). Tight regulation of ROS by antioxidant systems is necessary to balance the generation of ROS and to curtail oxidative damage to biomolecules.

Reactive oxygen species are difficult to measure reliably within cells. For example, early measurements of cellular oxidative state used dyes such as 2′,7′-dichlorofluorescin (DCFH) which were later shown to create additional radical species (Rota et al., 1999). The secondary reactions and instability of the dye made long-term imaging with DCFH impossible. To address these deficiencies, researchers engineered fluorescent reporter proteins such as HyPer, RoGFP, and RxYFP (Hanson et al., 2004; Belousov et al., 2006; Meyer and Dick, 2010), which undergo redox-sensitive conformational changes that elicit a change in fluorescence. Genetically encoded fluorescent proteins enable live-cell imaging and are further capable of localizing ROS production to specific sub-cellular compartments (Malinouski et al., 2011; Swain et al., 2016). Nevertheless, these sensors might miss low concentrations of H₂O₂ due to the endogenous enzyme peroxiredoxin, which is ∼100-fold more active toward H₂O₂ compared to introduced probes (Ezeriņa et al., 2014).

Given these challenges, some groups have taken a more computational approach to calculating the kinetics associated with H₂O₂ (Brito and Antunes, 2014). Lim et al. built a reaction-diffusion model to study localization of H₂O₂, which is important for control and specificity of redox signaling. They incorporated cytoplasmic diffusion into their reduced kinetic model of H₂O₂ clearance, in which peroxiredoxin is the dominant scavenging molecule (Lim et al., 2015). Using modeled concentration profiles obtained after bolus addition of H₂O₂ to the extracellular medium, the authors determined order-of-magnitude estimates for intracellular H₂O₂ diffusion through the cytosol, with a length scale of 4 μm and a time scale of 1 ms (Lim et al., 2016). The short length scale and rapid time scale indicate that H₂O₂ degradation and signaling are localized to the area where H₂O₂ is produced, contradicting the common modeling assumption of a well-mixed cytoplasm. This finding could explain discrepancies observed between bolus addition versus steady intracellular generation of H₂O₂ (Sobotta et al., 2013; Cheong et al., 2015). Rapid H₂O₂ scavenging also has implications for intracellular signaling, as H₂O₂ reactivity is limited to molecules in the subcellular vicinity.

Physiologically, ROS production increases under hypoxic conditions. Hypoxia causes a decrease in the maximum reaction rate of complex IV, which is thought to cause excess electron leakage from other components of the ETC, such as complex III (Chandel et al., 2000). After a return to normoxic conditions, ROS remain at higher hypoxic levels. The stable switch in ROS production is relevant during organ transplantation and other surgeries requiring an ischemic period. To better understand the mechanism behind this bistability, Selivanov et al. (2009) modeled the Q cycle mechanism of Complex III in the mitochondrial respiratory chain as the primary mechanism of ROS generation. Complex III can take on as many as 400 redox states due to its binding to quinones. The authors elaborated a system of differential equations describing the evolution of all of the redox states of Complex III. Model simulations predicted that Complex III can exist in two different steady states, a low ROS-producing state and a high ROS-producing state. This bistability is dependent upon the initial conditions of the system, specifically the redox state as predicted by levels of semiquinone and free ubiquinol. If starting in a highly reduced state, the overall system remains reduced, whereas if it starts in a less reduced state, Complex III progresses to a steady state with low semiquinone concentration and thus low ROS production. The overall system evolves to the high ROS producing state either by an increase in succinate concentration, causing the reduction of ubiquinone to ubiquinol, or a decrease in oxygen content. Once switched to this high ROS-producing state, Complex III persists in that state even after a return to lower succinate concentration or normoxic conditions. The sustained increase in ROS production provides a mechanism that may contribute to reperfusion injury after ischemia.

The model predictions were experimentally validated in isolated rat brain mitochondria incubated with succinate with or without the addition of ADP. The addition of ADP, and subsequent synthesis of ATP, switches mitochondria to a low ROS-producing state and thereby lowers mitochondrial membrane potential. Once all the ADP is consumed, membrane potential increases to pre-ADP levels, but ROS production remains at the lower initial level. These results agree with model predictions that two levels of ROS production could coexist under the same set of parameters and give rise to metabolic heterogeneity in an isogenic population of cells.

Reactive oxygen species are also produced by Complex I of the ETC. Gauthier et al. (2013a) built a computational model of the ETC focusing on ROS production by both complex I and III (Figure 1). Simulations were used to study the control of ROS production in cardiac myocytes under different metabolic conditions. The model is composed of non-linear ordinary differential equations describing the oxidation states of the various forms of ubiquinone, produced by complex I electron transfer, and the three subunits of complex III: cytochrome b, cytochrome c1, and the iron sulfur protein. The authors investigated how mitochondrial membrane potential, matrix pH, and ROS scavenging affect ROS production and control. When membrane potential increased 20 mV higher than unstressed cells to above ∼150 mV (Padmaraj et al., 2014), the model predicted that complex III ROS production as a function of membrane potential switches from zeroth order (constant production) to first order (exponential production). Increased membrane potential leads to a reduction in the Q cycle reaction rate, or conversion of ubiquinol to ubiquinone, causing a substantial increase in ROS production rate.

The model predictions agree with experimental results reporting a threshold membrane potential of 153 mV, above which ROS production increases dramatically (Korshunov et al., 1997). When simulating an increase in mitochondrial matrix pH, the model predicted that ROS production from complexes I and III increases during forward electron transport. This pH-dependent mechanism of ROS generation was experimentally observed by Selivanov et al. (2008) who found that an increase in pH from 6 to 7 caused a threefold increase in ROS production rate. During reverse electron transport, where electrons flow toward complex I in the presence of a weak reducing agent, complex I ROS production also increased with matrix alkalinization. The model therefore correctly predicted the dependence of ROS production on both mitochondrial membrane potential and matrix pH.

ROS levels are governed not only by production, but also by clearance through scavenging mechanisms. To gain a more complete picture of ROS dynamics, Gauthier et al. added glutathione and thioredoxin-mediated ROS scavenging mechanisms to the model (Aon et al., 2012). ROS production decreased to a minimum level as the mitochondrial environment became more oxidized and then rose again as the scavenging systems became depleted. This result agrees with the redox-optimized ROS balance hypothesis (Aon et al., 2010), which states that ROS levels are lowest at an intermediate mitochondrial redox potential. Together, the authors’ minimal model of ROS regulation produced results that matched many independent experiments describing different regimes of ROS production, providing support for the hypothesis that the cellular redox state influences the rate of ROS production.

Applying their model, Gauthier went on to investigate ROS production and scavenging in the context of heart failure (Gauthier et al., 2013b). Integrating the ETC-ROS model discussed above into a mitochondrial energetic-redox model (Kembro et al., 2013) allowed the authors to test the hypothesis that mitochondrial Ca²⁺ mismanagement leads to high levels of ROS during heart failure. In agreement with this hypothesis, their model showed that under conditions of mismanaged mitochondrial Ca²⁺, NADH levels decrease drastically under simulated cardiac pacing, highlighting the link between compromised Ca²⁺ and NADH regulation. Lower amounts of NADH lead to lower NADPH levels and a decreased ability to reduce scavenging enzymes for reuse, causing ROS to accumulate. ROS levels in mitochondria increase under increased load (Satapati et al., 2015), but surprisingly the model predicted that ROS production actually decreases under these conditions. The net increase in ROS abundance stems from an even-further reduction in ROS scavenging, which causes ROS accumulation in the cell. Therefore, in the setting of heart failure, preservation or restoration of scavenging enzymes may prove more effective than efforts to block ROS production (Linke et al., 2005).

In another model of ROS production by ETC complexes I and III, Bazil et al. (2016) described ROS generation by oxidative phosphorylation coupled to ATP demand. They updated an existing kinetic model of oxidative phosphorylation (Beard, 2005) to include ROS generation by complexes I and III and first-order scavenging by superoxide dismutase and peroxidase. Model simulations agreed with previous findings that free radical production by complex III is higher than complex I production under physiological conditions (Gauthier et al., 2013a). As ATP demand increases, the steady state production of ROS also increases, in line with experimental observations (Liu, 2010). The authors further applied their model to study reverse electron transport that is seen during reperfusion. Simulating ischemia/reperfusion led to bistability in ROS production (Selivanov et al., 2009, 2011) only when the activity of complex II was increased. Complex II activity requires electrons to be supplied to the quinone pool by the dehydrogenation of succinate to fumarate. The predicted importance of complex II agrees with work by Chouchani et al. (2014) showing that succinate is a main driver of mitochondrial ROS production upon reperfusion (Chouchani et al., 2014). Therefore, complex II inhibition during reperfusion could prove useful to decrease ROS production and reperfusion injury (Valls-Lacalle et al., 2016).

In a phenomenon known as ROS-Induced ROS Release (RIRR), damaged mitochondria produce an increased amount of ROS, which causes surrounding mitochondria to increase ROS production through a positive-feedback loop. Park et al. (2011) used an agent-based model describing inter-mitochondrial signaling to study the role of mitochondrial network dynamics in mitochondria-driven ROS production (Figure 1). Simulations were performed with three different mitochondrial networks: uniformly distributed mitochondria, as seen in cardiomyocytes; irregularly distributed mitochondria, as in neurons; and sparsely distributed mitochondria, as found in white blood cells (Maianski et al., 2004). The simulations introduced hydrogen peroxide as an initial oxidative stress, causing mitochondria in the surrounding area to produce more ROS. Mitochondrial ROS diffuse stochastically by random walk in 2D space, amplifying the local ROS response. Depending on the initial oxidative stress insult and the mitochondrial network dynamics, ROS production is blocked by antioxidant enzyme systems or becomes amplified by RIRR, which propagates ROS through the entire cell. The goal of the model was to predict the percent reactive mitochondria resulting from an initial oxidative stress input and the initial hydrogen peroxide concentration that causes RIRR.

The model indicated that ROS propagation is faster in the cardiomyocyte model than in the irregular distribution model, as shown by a higher dose dependence of reactive mitochondria as a function of initial oxidative stress. In addition to mitochondrial distribution, the model predicted that the density of mitochondria affects the response to oxidative stress inputs. Cells with a low density of mitochondria have considerable ROS propagation after low levels of oxidative stress, while cells with a high density of mitochondria only show strong ROS propagation after high levels of oxidative stress. The authors hypothesized that these differing responses to oxidative stress are due to differences in ROS signal transduction between mitochondrial networks. They further simulated the addition of different antioxidants to find that superoxide-scavenging antioxidants block ROS propagation more effectively in the cardiomyocyte model, while antioxidants that detoxify hydrogen peroxide are more effective in the irregular-distribution and low-density models of mitochondria. These results suggest that mitochondrial network configuration influences which molecular species is used to propagate ROS in the cell.

Cyclosporin A (CsA) is an immunosuppressant, which indirectly causes oxidative stress (Parra Cid et al., 2003) and adaptively activates the NRF2 pathway in the kidney. Hamon et al. (2014) fused an in vitro pharmacokinetic model (Wilmes et al., 2013) of CsA distribution in cultured renal epithelial cells with a dynamical model of NRF2 signaling originally designed to capture the cellular response to xenobiotics (Zhang et al., 2009). The authors adapted the NRF2 model to accommodate ROS as a state variable generated in proportion to cytosolic CsA (Hamon et al., 2014). In the revised model, ROS are detoxified by glutathione peroxidase and further act as an oxidant and inhibitor of KEAP1, which degrades NRF2. Last, CsA was forbidden from interacting with the aryl hydrocarbon receptor as xenobiotics do in the original model, because experimental data was lacking for such an interaction. To parameterize the fused model, Bayesian inference was used together with transcriptomic, proteomic, and metabolomic data collected from cells treated with different concentrations of CsA dosed daily for 2 weeks. The model predicted that low doses of CsA yielded widespread oscillations throughout the network as cells metabolized the administered CsA and detoxified ROS before the next administration. At high doses, however, the cell is overwhelmed and the modeled network locks into an elevated state of ROS adaptation. These predictions were not followed up experimentally, but the work of Hamon et al. (2014) nonetheless illustrates how toxicologic models can be repurposed for ROS specifically.

Aside from the ETC, ROS also play a key role in pathogen clearance. Neutrophils utilize ROS to attack bacteria engulfed within a phagosome. The source of this ROS burst is not from mitochondrial respiration but from the NADPH oxidase complex at the plasma membrane (Henderson and Chappel, 1996) (Figure 1). Levels of ROS oscillate in the neutrophil (Wymann et al., 1989); however, the mechanism behind these oscillations was unclear. Olsen et al. (2003) proposed that the oscillations arose from interactions among myeloperoxidase, melatonin, NADPH, and NADPH oxidase. To explore the oscillatory behavior, they built a two-compartment, differential equation model of the phagosome and the cytoplasm (Olsen et al., 2003). Without NADPH oxidase activity, model simulations exclusively produced damped oscillations that converged to a steady-state; by contrast, addition of NADPH oxidase elicited sustained oscillations similar to those reported experimentally (Petty, 2001). The authors triggered ROS oscillations in neutrophils with the activating chemotactic peptide FMLP and showed that pre-incubation with an inhibitor of NADPH oxidase blocked oscillations. Their model further predicted that melatonin would change the amplitude of the ROS oscillations measured. Pre-incubation of FMLP-activated neutrophils with melatonin confirmed the predicted increases in ROS amplitude. Computational and experimental modeling of NADPH oxidase in this setting allowed the authors to understand the basis of melatonin “priming” previously observed in neutrophils (Recchioni et al., 1998), underlining the power of pairing in silico and in vitro experiments.

Reactive oxygen species are also generated as a secondary byproduct of multiple signal transduction cascades (Meier et al., 1989; Thannickal and Fanburg, 1995). Haack et al. (2015) built a model of the WNT/β-catenin signaling pathway that included membrane-related processes as well as ROS signaling. The authors sought to explain experimental results showing that disruption of membrane lipid rafts inhibits WNT/β-catenin signaling, and also that ROS activate WNT signaling in the context of differentiation of human neural progenitor cells (Funato et al., 2006; Rharass et al., 2014). The three-compartment model is based on mass-action kinetics and includes: a membrane model, in which WNT binds to receptor LRP6 causing its phosphorylation within lipid rafts; an intracellular model, in which AXIN binds phosphorylated LRP6 to prevent degradation of β-catenin; and a redox model, in which ROS increase the concentration of DVL-bound AXIN, making AXIN unable to degrade β-catenin. Simulations were initiated with a burst of ROS that was shown experimentally to coincide with the beginning of neural progenitor differentiation induced by growth-factor withdrawal. The model predicted an immediate, transient β-catenin stabilization resulting from redox-dependent DVL/AXIN binding, followed by a sustained β-catenin response arising from lipid raft-dependent canonical WNT signaling. The immediate β-catenin accumulation was observed experimentally by the authors in lipid raft-deficient cells that maintained a transient β-catenin response. By including ROS signaling, the extended WNT/β-catenin signaling model correctly captured experimental β-catenin nuclear dynamics during early neuronal differentiation.

Dwivedi et al. (2015) looked at modulation of cell signaling by ROS in another setting, using the IL-4 signaling pathway as a redox-regulated case study. The authors sought to identify the most important mechanisms of redox regulation in the IL-4 pathway, which is important for regulating the effector T-cell response. The activated IL-4 receptor complex upregulates ROS through NADPH oxidase (Sharma et al., 2008), which influences signal transduction that proceeds through JAKs and culminates in the phosphorylation of STAT6. To identify the combination of regulatory mechanisms that best recapitulated the dynamics of IL-4 induced STAT6 phosphorylation, the authors turned to Monte Carlo analysis of an IL-4 ordinary differential equation model. Phosphorylated STAT6 dynamics were best captured by a model that incorporated a protein tyrosine phosphatase whose activity and nucleocytoplasmic shuttling were ROS sensitive. ROS regulation of phosphatase activity and localization, along with other ROS-independent mechanisms, were included in the systems-level model of IL-4 signaling, with parameters fit to experimental data in IL-4-stimulated Jurkat cells. The model predicted diminished STAT6 phosphorylation following IL-4 stimulation and ROS inhibition, which was confirmed experimentally by NADPH oxidase inhibition of IL-4-stimulated cells. Transient oxidation of protein tyrosine phosphatases was also observed experimentally by oxidized protein tyrosine phosphatase immunoprecipitation of extracts from IL-4-stimulated Jurkat cells. The authors’ systems-level model provides a framework for investigating additional modes of receptor-initiated oxidation not previously explored.

Smith and Shanley (2013) adapted an existing differential equation model of insulin signaling (Sedaghat et al., 2002) to incorporate ROS and study the interplay between insulin signaling and oxidative stress (Figure 1). Insulin-stimulated ROS production was assumed to occur through activation of NADPH oxidase and is about fivefold higher than the background level of mitochondrially produced ROS (Mahadev et al., 2001). ROS deactivate the phosphatases PTEN and PTP, activate the kinases JNK and IKK, and are detoxified by cytoplasmic SOD2. This model was used to make predictions about ROS, FOXO, SOD2, and insulin receptor abundances over long time scales.

When hydrogen peroxide was added as an oxidant to the system with or without insulin stimulation, the model predicted surprisingly different responses. Hydrogen peroxide alone caused modest glucose uptake and insulin alone caused strong glucose uptake, but hydrogen peroxide and insulin stimulation together were antagonistic, causing only moderate glucose uptake. In the model, this dampening effect of oxidative stress occurs because hydrogen peroxide and insulin together activate protein kinases (e.g., JNK, IKK), which cause hyperphosphorylation of IRS1 and decrease its ability to form the IRS1-PI3K complex that stimulates glucose uptake. The model also predicted that the FOXO-mediated antioxidant response depends critically on the extent of oxidative stress. With low oxidative stress, the antioxidant enzyme SOD2 is upregulated by FOXO through a JNK-mediated mechanism, but SOD2 is downregulated at higher levels of stress through an IKK-mediated mechanism. Although simplified in its handling of ROS and antioxidant pathways, this integrated systems model captures some of the complexities of oxidative stress for an important metabolic pathway.

The previously discussed models incorporated ROS into a single canonical signaling pathway, but the generation and handling of ROS pervades multiple pathways and can lead to counterintuitive cell outcomes (Wang et al., 2011). Zeigler et al. (2016) incorporated ROS as part of a much larger signaling network to identify regulators of cardiac fibrosis. A cardiac fibroblast signaling network was designed to study drivers of fibrosis, which is implicated in many cardiac pathologies (Moreo et al., 2009). The network was compiled from experimental data on 10 pathways that are known to be important in cardiac injury, such as the IL-1 and TGFβ pathways. The model was formed using a logic-based differential equation approach (Kraeutler et al., 2010), whereby species are represented as differential equations with rates of change dictated by Hill functions and truth tables comprised of interacting biomolecules. ROS production is controlled by the activity of NADPH oxidase and feeds into the truth tables of JNK and ERK activation, linking ROS generation at the plasma membrane to downstream intracellular responses.

In the model, reducing ROS levels had far-reaching and context-dependent effects on the network. Under baseline conditions, reductions in ROS caused a decrease in matrix metalloproteinase 9 (MMP9), which is important for the breakdown of extracellular matrix. By contrast, in an environment with high TGFβ signaling, like in a myocardial infarction, reducing ROS led to an increase in MMP9 activity. Therefore, a therapeutic prediction of the model is that antioxidant treatment for fibrosis would be more beneficial in the high TGFβ environment of myocardial infarct.

The computational models of ROS biology covered in this review largely focus on ROS handling within the cell or on ROS modulation of canonical signaling pathways. In the future, we anticipate more sophisticated models that combine ROS handling and signaling concurrently. A prime test bed for such an integrated approach would be the NF-κB pathway, which is activated by ROS (Schreck et al., 1991) and is responsible for inducing scavenging enzymes such as SOD2 (Wong and Goeddel, 1988). Finn and Kemp (2012) have assembled a provisional model of IKKβ S-glutathionylation in the setting of antioxidants and chemotherapy-induced ROS. The coupling of signaling, production, and scavenging could give rise to feedback networks that explain the variable oxidative stress observed in some settings among single cells in very similar microenvironments (Wang et al., 2011).

There is also a need to build multiscale models that place ROS in the broader context of developing tissues, tumors, and infections (Schwarz, 1996; Bajikar and Janes, 2012; Wang et al., 2012). The dynamics of proliferation and death impinge on metabolism and signal transduction, which culminate to impact the redox state of cells in the population. Crosstalk between these cellular pathways may require different classes of modeling than those implemented so far (Anderson et al., 2006; Chitforoushzadeh et al., 2016). Advances in measurement will likewise expand the scope of targets modified by ROS (Kim et al., 2015) and reveal the extent to which molecular crosstalk is underappreciated.

Integration of ROS signaling into larger networks may allow researchers to predict outcomes of drug treatments that affect ROS generation, which causes drug resistance in some cancer contexts (Nieborowska-Skorska et al., 2014; Okon et al., 2015). A deeper understanding of ROS network dynamics could generate combinatorial treatments that avoid neutralizing drug efficacy. In the broader human population, there are many polymorphisms that affect ROS generation and scavenging, such as p22^phox C242T and SOD2 A16V (Bresciani et al., 2013; Meijles et al., 2016). These variants may tune how ROS interacts with other signaling networks, contributing to heterogeneous patient responses during therapy.

ROS are a fact of life that cannot be ignored. Like a living cell, investigators must find ways to deal with ROS holistically and achieve our goals despite their presence. The tools for pharmacologic modulation of ROS are predominantly limited to antioxidants. Systems modeling of ROS may one day provide a venue for exploring more-precise interventions that account for the complex biological processes involved.

EP and CS conceived of the review with input from KJ. EP and CS drafted the review with revisions provided by KJ.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.