A Robust Model for Circadian Redox Oscillations

Because of their potential to oxidize and damage cellular protein and lipids, ROS levels must be under tight regulation. Prx3 is one of the major antioxidant proteins involved in H₂O₂ removal in mitochondria. A conserved cysteine (Cys) residue in its active site is oxidized by H₂O₂ to Cys-sulfenic acid (Prx3-SOH, reaction k1 in Figure 1B). This intermediate can react with a second conserved Cys from another Prx3 subunit to produce an intermolecular disulfide bond (reaction k2), which can be reduced by the thioredoxin (Trx) system of the cell (reaction k3). Alternatively, Prx3-SOH can undergo further oxidation, termed hyperoxidation, in an S-sulfinylation reaction, to form a sulfinic acid (Prx3-SO₂H, reaction k4) [9,13]. Prx3-SO₂H is catalytically inactive [9,13], and in adrenal gland cells, brown adipocytes or cardiomyocytes, it can be reduced back to Prx3-SOH by action of sulfiredoxin (Srx) (reaction k5) [11]. In other cell types and other model organisms, the Srx-mediated reduction of Prx-SO₂H is not required for rhythm generation [12].

In mitochondria from brown adipose, adrenal gland and heart tissue, Prx3-SO₂H and Srx levels oscillate with a circadian period [11]. H₂O₂ increase results in the oxidation and inactivation of Prx3 (reactions k1 and k4 from Figure 1B). The consequence is the accumulation of H₂O₂ inside mitochondria and its overflow to the cytosol (k6 in Figure 1B). Cytosolic H₂O₂ activates pathways to control its own production. Among others, it stimulates Srx oxidation and intermolecular disulfide bridge formation with Hsp90 (reaction k7) to promote translocation of Srx to the mitochondria (reaction k8). Mitochondrial Srx can reduce and thus reactivate Prx3-SO₂H (reaction k5) [11]. Srx levels peak in mitochondria ∼8 h after Prx3-SO₂H, and Srx becomes sensitive to degradation (reaction k9) when Prx3-SO₂H levels decrease [11]. The mitochondrial import of Srx thus constitutes a negative feedback that enables a new cycle of H₂O₂ removal and Prx3-SO₂H accumulation (Figure 1B).

In order to understand how fast biochemical redox reactions result in slow 24 h rhythms, we develop a deterministic model containing the biochemical species from the Prx3/Srx system shown in Figure 1B. The complete set of equations of this large model, which we refer to as detailed model, are found in Appendix A. Unfortunately, quantitative details of the kinetic processes are not known and thus estimating parameters in such a large model represents a challenge. For this reason, and in order to gain insight into the design principles of redox rhythms in the Prx3/Srx system, we simplify the detailed model to its core motif (Figure 2A). Details of the model reduction are found in Appendix B.

We identify the core motif that is required for generation of Prx3-SO₂H and mitochondrial Srx circadian oscillations by systematically changing parameters and clamping variables to their mean value [14] (Appendix B). The scheme and equations of this minimal model are shown in Figure 2. It only contains the variables, kinetic parameters, and reactions that are necessary and sufficient for self-sustained rhythm generation.

The core oscillator model consists of only five variables and just one oxidation event (Figure 2A). A represents the active but partially oxidized Prx3 (Prx3-SOH) that gets further hyperoxidized by the danger 1 (D1, mitochondrial H₂O₂) to produce the inactive Prx3-SO₂H, I. As the A→I reaction occurs, A levels decrease. The result is that the ability of mitochondria to remove D1 also decreases. When mitochondria can no longer eliminate the oxidant, D1 accumulates and overflows to the cytosol, where it is referred to as danger 2 (D2, cytosolic H₂O₂). D2 activates the translocation of the rescuer R (mitochondrial Srx), that reduces I back to A, closing the negative feedback loop. The Prx3-SH oxidation to Prx3-SOH and the D2-induced oxidation of cytosolic Srx (reactions k1 and k7 in Figure 1B) were found to be dispensable for rhythm generation. Thus, these biochemical species, as well as the reactions they take part in, are omitted from the minimal model.

We formulate ordinary differential equations (ODEs) to describe the behavior of the system over time. Production and removal terms are modeled with mass action kinetics (Figure 2B). The nonlinear term in the model arises from the quasi-steady state approximation [15,16] performed on the equation for variable A, whose production and removal reactions are modeled as bilinear terms (details provided in Appendix B). We assume constant D1 production (parameter p in the model) and constant Prx3 pool A+I over time. Moreover, we only consider the negative feedback performed by R-induced reactivation of A, which is stimulated by the D2 increase. However, other additional feedback loops might exist, as seen by studies that have shown that cytosolic H₂O₂D2 can feed back and decrease D1 production by activation of the stress-activated protein kinase p38 MAPK, among other pathways [17].

The dynamics of a system depend greatly on parameter values. Unfortunately, most kinetic rates in the Prx3/Srx sytem are not known in quantitative detail and therefore finding reasonable parameters represents a major challenge. Nevertheless, some physiological constraints can be taken into account to narrow down the plausible range of parameter values. First of all, the system needs to oscillate with a circadian period and the phase relationship between I and R should be approximately 8 h, according to previously published data [11]. Other constraints come from enzymatic assays performed with Prx3 and Srx, as well as from studies that have focused on H₂O₂ signaling and transport [11,18,19,20,21,22,23,24]. They have estimated rates of Prx3 hyperoxidation, Srx-mediated Prx3-SO₂H reduction, and Srx degradation and have allowed the derivation of a reasonable rate for H₂O₂ translocation across membranes (Table 1).

In order to obtain a reasonable set of parameters that satisfies the physiological constraints, we systematically vary all model parameters and perform bifurcation and sensitivity analyses. Bifurcation diagrams show under which parameter values the system oscillates and sensitivity analyses address how the oscillation period varies over the oscillatory range as a parameter is changed.

Bifurcation diagrams are presented for a selected choice of parameters in Figure 3A–C. The bifurcation plots show that the system oscillates for values of d,e,q<1. The limit cycle emerges at d=0.02 and disappears at d=0.5 (Figure 3A, bifurcation plots for parameters e and q are shown in Appendix C). We also find that, in contrast to small d,e,q parameter values, the system requires a 1000–10,000-fold higher a oxidation rate to enter into the oscillatory regime. The Hopf bifurcation occurs at a=120 (Figure 3B) and once this value is reached, oscillations persist for a wide range of kinetic parameters, with little effects in oscillation amplitude or period, indicating robustness of the model. The period sensitivity analyses reveal that the oscillation period depends strongly on the translocation rates d and e and on the R degradation parameter q, being most sensitive to changes in parameter d (Figure 3C and Appendix D).

These findings show that a fast A inactivation reaction (high a value) in combination with a slower negative feedback loop (1000–10,000-fold lower d,e,q rates) constitute the backbone of the redox oscillator model (Figure 3D). We thus find a plausible set of model parameter values that produce robust circadian oscillations with the expected characteristics: a period of 24.2 h and a I/R phase difference of 8.7 h, as determined by maxima estimation of the R and I dynamics (Figure 4). Moreover, the parameter choice is in agreement with the high Prx3 oxidation rate measured experimentally [18] and with the ∼10,000-fold slower physiological H₂O₂ translocation (Table 1).

According to our results, only a small fraction of active Prx3 A is needed to keep D1 levels close to 0, as seen by the higher I levels compared to A during the whole simulated time (Figure 4). The results also predict that D2 levels peak 3.9 h after D1 (or 4.5 h after I), an observation that could be relevant in the context of the canonical TTFL and will be further addressed in the Discussion section.

Interestingly, insights from the redox oscillator model apply as well to larger and more complex models. Firstly, negative feedback loops (D1−D2−R in the minimal model) as well as nonlinear terms (given by the A equation, Figure 2B) are required to achieve self-sustained oscillations [25]. Secondly, overcritical delays about quarter to half of a period are necessary for oscillations to arise [25,26,27] (in the minimal model, slow d,e,q rates are required to yield the 8.7 h delay). Lastly, translocation and degradation rates have profound effects on the period [28,29,30,31]. It should be noted, nevertheless, that the choice of parameters is not tailored to specific kinetic data. Kinetic parameters might differ among tissues and might also depend on physiological conditions.

We have shown that, in order to oscillate, the mitochondrial Prx3-SO₂H/Srx (I/R) system requires a fast D1-induced inactivation of A followed by a lengthy and slower negative feedback loop D1−D2−R that reactivates A (Figure 3D). As seen by the equations from Figure 2B, the dynamics of A,I, and D1 depend on the fast a oxidation rate compared to the slower dynamics of D2 or R. This is reflected in the relaxation- or triangular-like waveforms of A,I, and D1 (Figure 4). Such relaxation oscillations are generally characterized by two consecutive processes that occur on different timescales. There are intervals of time, during which little happens (A=0 or I=1 in the system’s dynamics), followed by time intervals with considerable changes. The result is the triangular-like dynamics observed for A, I, and D1 (Figure 4).

It is known that oscillations depend on negative feedbacks, but that in the absence of appropriate delays or nonlinear terms, negative feedback circuits often settle into a stable steady state termed homeostasis [25,32,33]. Several computational studies have demonstrated that sufficiently strong nonlinearities are required to generate self-sustained oscillations [31,34,35,36]. Such nonlinear terms often exhibit sigmoidal or switch-like characteristics. Surprisingly, no switch-like curve is evident from the only nonlinear term of the model abD1RbR+aD1 (Figure 2B and Appendix B), nor from the response of any variable to upstream components (data not shown). Nonetheless, we find that this term is at the core of the oscillations, as its linearization results in the stabilization of the steady state. The large value of a allows the production of self-sustained rhythms, even in the absence of an explicit kinetic switch.

In order to emphasize the importance of the nonlinear term in the context of oscillations, and given the switch-like dynamics of A,I, and D1, we still use the term switch. Instead of using it to describe how the steady state of a first component responds to a second component which is upstream of the former (often done in the literature [36,37,38,39]), we focus here on temporal switches that mark transitions between phases. We see that the model dynamics can be split into three phases separated by two temporal switches (Figure 5).

The first phase is characterized by the inactivation of A as a result of D1 removal. As long as there is active Prx3 (A) in the system, mitochondrial H₂O₂ (D1) levels are kept in check and D1 is removed. But the cost of the D1 clearance is the inactivation of A to I, and thus A levels decrease as I increases (first box in Figure 5B). When most of the active Prx3 pool A has been hyperoxidized and inactivated to I, the first switch occurs and the system progresses to the next phase. As a consequence of the A inactivation, its ability to remove D1 molecules decreases and thus D1 accumulates and leaks to the cytosol. In other words, the flux of D1 changes from inactivating A to leaking, what gives this first temporal switch its name: inactivation/leakage (I/L) switch. As a result of the I/L switch, D1 accumulates, leaks and activates the D1−D2−R negative feedback. Consequently, D2 and R levels increase (phase 2, second box in Figure 5B). At the end of the second phase, the increase in R is accompanied by a decrease in D1. Once D1 and R reach critical values in the nonlinear term, the negative feedback becomes effective and the dynamics of A, that are governed by the nonlinearity (Figure 2B), change suddenly. This triggers the activation (A) switch and the sharp increase in A levels (phase 3, third box in Figure 5B). The system is thus switched to its active state and a new round of D1 removal can start. Furthermore, R becomes sensitive to degradation. The series of two temporal switches separate the two timescales in the relaxation oscillations. The I/L switch sets the start of the slow timescale (the negative feedback loop), whereas the A switch marks the beginning of the fast clearance of D1.

We have designed the first model for the complex biochemical system of Prx3-SO₂H/Srx redox oscillations and we find that the loop Prx3-SOH – Prx3-SO₂H, together with the negative feedback mediated by Srx (Figure 1B), is necessary and sufficient for the generation of oscillations. The design principles of this oscillator are (i) a fast Prx3-SOH inactivation followed by (ii) a slow and delayed negative feedback loop, where mitochondrial H₂O₂ leaks to the cytosol to promote the translocation of cytosolic Srx to mitochondria, approximately 9 h after the inactivation of Prx3-SOH. This simple backbone reproduces the previously described circadian oscillations of mitochondrial Prx3-SO₂H and Srx as well as the delay in mitochondrial Srx import.

Although quantitative details of the kinetic parameters in the Prx3/Srx system remain largely unknown, we have found a reasonable set of parameters that is consistent with the biology they describe. Previous biochemical studies have estimated that Prxs are 1000–10,000 times more reactive to H₂O₂ than other reduced cellular thiols [18], and that thus they can confine the diffusion of the oxidant [23]. Under the assumption that 0.01% of the membrane area contains transport proteins (Aqp) that can passively transport H₂O₂ out of the compartment [22], the translocation rate of H₂O₂ is 10,000-fold smaller than its reduction rate (Table 1). This is in agreement with the parameter values, a=1000 and d=0.2. Since the import of Srx to mitochondria also requires a protein transporter, we have assumed this rate to be in the same order of magnitude as the translocation of H₂O₂ to the cytosol. The experimental determination of the rate at which a protein is transported through a membrane is challenging, because the transport depends on a number of factors including pH, membrane voltage, protein length, protein 3D conformation, membrane fatty acid composition, and ATP levels (if the transport is active), among others. Hence the estimation of translocation kinetics is only preliminary. Furthermore, the choices for Srx degradation rate q and hyperoxidized Prx3 reduction rate by Srx b are consistent with previous experimental studies, which have estimated a half-life of sulfiredoxin of 4–5 h [11,24], and a rate of reduction of hyperoxidized Prx by Srx of approximately 2 M⁻¹s⁻¹, 1000–10,000 times smaller than its oxidation rate [19] (Table 1).

It is known that relaxation oscillations typically depend on positive feedback loops [40,41]. However, even in the absence of explicit positive feedback loops, the redox oscillator model still produces oscillations of Prx3-SOH, Prx3-SO₂H and mitochondrial H₂O₂ that are of relaxation type. Two temporal switches characterize their triangular-like waveform and divide the system’s dynamics into three phases. A summary is shown in Figure 6. This scheme with three phases and two switches is reminiscent to that of the cell cycle. A series of biochemical switches controls transitions between the various phases of the cell cycle. They maintain its orderly progression and act as checkpoints to ensure that each phase has been properly completed before progression to the next phase. Such switches have been shown to generate decisive, robust (and potentially irreversible) transitions and trigger stable oscillations [33,42,43]. It should be noted, however, that the switches that regulate cell cycle transitions are the consequence of a complicated network of positive and negative feedback loops and result in bistability, as opposed to our simple model that (strikingly) contains only one negative feedback loop. The redox model highlights the importance of negative feedbacks and the diversity of switches in the context of oscillations.

The detailed model (Appendix A) also oscillates along the lines of the design principles identified for the reduced model. It shows the same characteristics: relaxation oscillations, temporal switches, and expected period and I/R phase delay, supporting the robustness of the minimal core model. According to biochemical assays, the rate constant for Prx3-SH oxidation to sulfenic acid Prx3-SOH is in the order of 10⁷ M⁻¹ s⁻¹ [18]. From the simulations of the detailed model we see that when this rate (k1 in the model) increases, the amplitude of Prx3-SO₂H and Srx oscillations also increases (Appendix A). This suggests that although the loop Prx3-SH – Prx3-SOH – Prx3-SS – Prx3-SH is not required for oscillations, it fine-tunes the dynamics of the Prx3/Srx oscillatory system. Large Prx3-SO₂H and Srx amplitudes are in agreement with experimental data. Western blots show intense bands at the peak of Prx3-SO₂H or Srx oscillations and no bands at the trough [11].

The oscillator model presented here applies to Prx oscillations in mammalian adrenal gland, heart, and brown adipose tissue. Prx daily rhythms have been found not only in mammals, but also in a variety of eukaryotes, cyanobacteria, and archaea, making these oscillation an evolutionarily conserved circadian rhythm marker [10]. A key unanswered question is what determines Prx oscillations in other cell types and phyla. Srx indeed accounts for the oscillations in adrenal gland, brown adipose tissue, and heart [11], but other organisms that display oscillations in Prx do not express Srx homologs (e.g., C. elegans and N. crassa), suggesting that other mechanisms might be in place. The same is true for mouse red blood cells. Although they express Srx, Prx-SO₂H rhythms are largely unaltered in erythrocytes derived from Srx knockout mice and, by contrast, application of proteasome inhibitors eliminated the decay phase of the Prx-SO₂H rhythm in wild-type mouse erithrocytes [12]. These findings suggest that proteasomal degradation might account for the decay phase and hence for the negative feedback of the Prx-SO₂H oscillations in other cell types and phyla.

Given that Prx proteins and their oxidation rhythms are highly conserved among species, Prx rhythms seem to be evolutionarily ancient. It is also likely that they emerged before the circadian TTFL, which is believed to have evolved independently in separate species [10]. It is an important goal to now understand the relationship between the canonical TTFL and Prx oscillations, as it is becoming evident that there is a reciprocal interplay between the two. In one direction, the period length of Prx-SO₂H rhythms is extended in embryonic fibroblasts derived from Cry1/2^−/− mice [8]. In behaviorally arhythmic Drosophila mutants and N. crassa mutants displaying a lengthened period, Prx oscillations are perturbed in phase [10]. In the other direction, Arabidopsis mutants lacking Prx genes display TTFL rhythms that are altered in phase or amplitude [10], and knockdown of specific PRX isoforms in human U2OS cells affects the period length and amplitude of clock gene rhythms [44].

The modeling approach described in this study assumes constant H₂O₂ production. Nevertheless, the sources of mitochondrial H₂O₂ in tissues where the redox clockwork is present have been shown to be, at least in part, under the control of the circadian TTFL. In adrenal gland, for example, the major source of H₂O₂ is steroidogenesis, and this process is regulated by the TTFL [45,46]. In heart or brown adipose tissue, in which oxidative metabolism is very high, the respiratory component becomes an important source of H₂O₂, and respiration rate has also been shown to be rhythmically controlled by the canonical TTFL [47].

Moreover, the Prx3/Srx redox oscillator model contains only one negative feedback, namely the reduction of Prx3-SO₂H (I) by Srx (R), which is stimulated by the cytosolic H₂O₂ (D2) increase. Cytosolic H₂O₂, however, can signal as a second messenger and activate the p38 MAPK pathway to decrease mitochondrial H₂O₂ production [11,17]. The periodic H₂O₂ release to the cytosol is expected to act on other cytosolic or nuclear targets that are not yet identified. Canonical clock proteins might be directly or indirectly affected by the cytosolic H₂O₂ increase. This idea is supported by a study from 2014, which showed that the interaction of the two mammalian core clock proteins PER2 and CRY1 is redox-sensitive [48], as well as by other studies that have shown that the cellular redox poise can regulate the TTFL oscillator [49,50]. The H₂O₂ leakage might thus represent a potential coupling node between both the redox and the TTFL clock. The reciprocal interplay between the Prx system and the local TTFL clock might allow synchronization between local metabolic activity and systemic circadian regulation. The presence of Srx as part of the Prx clockwork in adrenal gland, brown adipose tissue or heart, but not in the liver, for instance, might indicate differences in how metabolic signals are transmitted to the TTFL.

Mathematical models are most useful when they make predictions that can be tested (and hopefully validated) in vivo. Because the Prx3/Srx oscillator model is not quantitative in detail due to the limited knowledge of rate constants, it only allows us to make some basic predictions. Interfering with the expression of Srx or specific aquoporin (Aqp) isoforms involved in the H₂O₂ translocation across mitochondrial membranes would, intuitively, disrupt Prx3/Srx oscillations. Moreover, the simulations show that H₂O₂ peaks in the cytoplasm ∼4.5 h after the peak of Prx3-SO₂H (or ∼4 h after mitochondrial H₂O₂). This observation might represent an attractive target for redox sensors that could measure when H₂O₂ peaks in the cytosol in relation to its peak in mitochondria. The peak of cytoplasmatic H₂O₂ at a proper time might be of importance to transmit appropriately timed signals about the redox state of the cell to the canonical TTFL oscillator. It is likely that a stable phase relationship is required between both oscillators for their optimal function. Our model also predicts that the system should be able to oscillate even under conditions of constant mitochondrial H₂O₂ production. Such a self-sustained nature of the Prx3 cycle could be easily tested in isolated mitochondria led to respire in a constant environment.

It should be noted that the redox model reproduces the circadian rhythms at the level of mitochondrial membranes. Mitochondria vary in number depending on the cell type. For example, a human liver cell contains 1000–2000 mitochondria, whereas a cardiac myocyte contains 7000–10,000 mitochondria [51]. Although single mitochondria might be competent redox oscillators, it seems plausible that they synchronize their redox circadian cycles to each other, as seen by the clear bands in published western blots of Prx3-SO₂H or Srx [11], in order to transmit the information about the redox state of the cell to the TTFL clock as an ensemble.

Finally, we close with a short remark about the beauty of simple and generic models. We have elaborated on the basic ingredients that constitute the Prx3/Srx redox oscillator in mammalian mitochondria from adrenal gland, heart, and brown adipose tissue. We have found the core motif that is necessary and sufficient to generate self-sustained rhythms that reproduce previously published experimental data. Our model has been implemented using plausible and simple, yet generic mathematical, representations. This five-variable network might be relevant for other simplified oscillatory systems that follow the same logic of fast (in)activation followed by slow negative feedback, such as the mitotic oscillator involving cyclin and Cdc2 kinase interactions or the Ca²⁺ oscillations based on Ca²⁺-induced Ca²⁺ release [52]. The fact that nature has converged on common mathematical structures underlines the links between similar dynamic phenomena occurring in widely different biological settings. It even constitutes a step forward in the synthetic design of reliable biological clocks [53,54].