- Research
- Open access
- Published:

# Evolution of the repression mechanisms in circadian clocks

*Genome Biology*
**volume 23**, Article number: 17 (2022)

## Abstract

### Background

Circadian (daily) timekeeping is essential to the survival of many organisms. An integral part of all circadian timekeeping systems is negative feedback between an activator and repressor. However, the role of this feedback varies widely between lower and higher organisms.

### Results

Here, we study repression mechanisms in the cyanobacterial and eukaryotic clocks through mathematical modeling and systems analysis. We find a common mathematical model that describes the mechanism by which organisms generate rhythms; however, transcription’s role in this has diverged. In cyanobacteria, protein sequestration and phosphorylation generate and regulate rhythms while transcription regulation keeps proteins in proper stoichiometric balance. Based on recent experimental work, we propose a repressor phospholock mechanism that models the negative feedback through transcription in clocks of higher organisms. Interestingly, this model, when coupled with activator phosphorylation, allows for oscillations over a wide range of protein stoichiometries, thereby reconciling the negative feedback mechanism in *Neurospora* with that in mammals and cyanobacteria.

### Conclusions

Taken together, these results paint a picture of how circadian timekeeping may have evolved.

## Background

Circadian clocks time most behavioral and physiological events in a 24-h day to optimize fitness in many organisms [1]. Molecular mechanisms are responsible for driving these rhythms at the cellular level [2]. A common theme for the generation of oscillations at the molecular level is the presence of negative feedback [3]. However, the mechanisms by which various organisms implement negative feedback differ. For example, in cyanobacteria, three proteins, KaiA, KaiB, and KaiC, interact to generate near 24-h rhythms in the phosphorylation status of KaiC, giving rise to a purely post-translational clock. On the other hand, eukaryotic organisms produce the required negative feedback through transcriptional activation of a repressor that, after a sufficient amount of repressor is present, inhibits the activator from further promoting transcription.

Despite our increased understanding of the core negative feedback architecture across many organisms, the evolutionary paradigm of circadian clocks remains unclear. Here, we develop novel mathematical models of the clocks in cyanobacteria, *Neurospora*, and mammals that reveal convergent principles by which diverse organisms generate oscillations. Our cyanobacterial clock models reflect protein sequestration between KaiA and KaiC [4, 5] while transcription regulates proper stoichiometric balance. Additionally, we present a mathematical model of the eukaryotic clock incorporating a novel “phospholock” mechanism that couples phosphorylation with protein sequestration. We find that the addition of phosphorylation increases the robustness of oscillations at biologically reasonable protein affinities. Furthermore, we show that, when adding activator phosphorylation, our model allows for oscillations at lower stoichiometric ratios of repressor to activator, reconciling experimental results from *Neurospora* with sequestration models.

Altogether, our analysis reveals that the circadian clock mechanisms among prokaryotes and eukaryotes may have evolved towards similar ends. That is, as the clock’s core feedback moved from post-translational regulation in cyanobacteria to a transcription-translation feedback loop in *Neurospora* and higher organisms, the regulation of stoichiometry, which is essential for robust oscillations, moved from a transcription-translation feedback loop to a post-translational mechanism.

## Results

### A detailed model of the cyanobacterial clock

Cyanobacteria are among the simplest organisms to exhibit circadian oscillations in cellular components. Proper functioning of the cyanobacterial circadian clock relies on interactions among three proteins: KaiA, KaiB, and KaiC [6]. Central to the protein interactions are two indispensable phosphorylation sites on KaiC-Ser-431 (S431) and KaiC-Thr-432 (T432). Recently, crystal structural analysis of the protein complexes in the KaiABC system revealed a detailed sequestration mechanism of KaiA by KaiBC protein complexes (Fig. 1A) and confirmed in silico [7, 8]. First, KaiA acts as an enzyme to enhance the phosphorylation of KaiC on the T432 site and, subsequently, the S431 site. Then, once site S431 is phosphorylated, KaiC undergoes a conformational change from a pre-hydrolysis state to a post-hydrolysis state creating a hub [9]. Simultaneously, KaiB undergoes a fold-change transition into an active state that is captured by the post-hydrolysis KaiC on the N-terminal (CI) domain. The KaiBC complex recruits KaiA proteins and prevents KaiC activation on the C-terminal (CII) domain. In other words, KaiBC complexes sequester free KaiA proteins.

Here, we model the sequestration mechanism of KaiA by the KaiBC complex with a system of ordinary differential equations (ODEs) assuming mass-action kinetics (see Eqs. –(8) in Additional file 1: Supplementary Information). Our model considers KaiC proteins in four states: unphosphorylated KaiC (denoted by *U*), S431 phosphorylated KaiC (denoted by *S*), T432 phosphorylated KaiC (denoted by *T*), and S431 and T432 phosphorylated KaiC (denoted by *ST*). A similar notation is used in the phosphorylation status of KaiBC. Our model assumes that KaiB binding only happens when the S431 site is phosphorylated [10, 11].

Simulations of this detailed model suggest that the components oscillate over a wide range of reaction rates. Analyzing the temporal profiles of the KaiABC oscillator, we see that the fraction of phosphorylated KaiC oscillates with a roughly 24-h period (Fig. S1, Additional file 1: Supplementary Information). Also, an asymmetry between the KaiC phosphorylation and dephosphorylation phases are consistent with previous experimental results [10, 12]. Analysis of the phases of these profiles suggests that the dynamics are indeed what we expect: first, free KaiA activates unphosphorylated KaiC (*U*); next, activated KaiC phosphorylates on the T432 site (*T*), which then slowly phosphorylates on the S431 site (*ST*). The total phosphorylation level of KaiC increases accordingly during the activation phase. Meanwhile, KaiB begins to bind with the S431 and T432 phosphorylated KaiC. The resulting KaiBC complexes in turn sequester KaiA through tight binding. Once sequestration depletes KaiA sufficiently, phosphorylation stops and dephosphorylation dominates, leading the system back to a highly unphosphorylated state.

Our detailed model reproduces important qualitative and quantitative results, which we illustrate by plotting circadian data from [12] alongside our model simulations (Fig. 2). In particular, we observe an asymmetric circadian rhythm where KaiC proteins spend less time in phosphorylation than dephosphorylation (Fig. 2A, B). Also, our model shows a robust period of around 24 h with less than 5% fluctuations for several ATP/ADP ratios (Fig. 2C–F). The CI domain of KaiC binds with KaiB, and thus sequesters KaiA, at night as shown experimentally in [12] (Fig. 2G–H). Taken together, our model recapitulates experimental results revealing the importance of the KaiC CI domain for the proper functioning of the cyanobacterial circadian clock.

Finally, we perform a bifurcation analysis (Fig. S2, Additional file 1: Supplementary Information) and stochastic simulations (Fig. S3, Additional file 1: Supplementary Information) of our detailed model. This bifurcation analysis reveals a result similar to that in [13] where the relative ratio of KaiA and KaiB to KaiC plays an important role in the circadian clock (Fig. S2, Additional file 1: Supplementary Information). Stochastic simulations reveal that, as the total number of KaiC molecules in the system increases, the stability of the oscillations increases. When the total number of KaiC molecules is low (*C*_{T}≤50), the trajectories exhibit randomness without sustained oscillations (Fig. S3A, Additional file 1: Supplementary Information). As *C*_{T} increases, the stochastic simulations approach that of the deterministic system (Fig. S3, Additional file 1: Supplementary Information).

### A simple model of the cyanobacterial clock reveals tight binding is necessary for oscillations, similar to the mammalian system

Here, we propose a simple model of the cyanobacterial clock where only the KaiA and KaiC proteins are explicitly included (see Fig. 1B for the schematic). Our model is based on the detailed cyanobacterial model from above which accounts for a wide range of experimental observations. In our simplified model, given below in System (C), the KaiA-enhanced phosphorylation of the *T* site is much faster than autophosphorylation of the *S* site. The tight binding between KaiA and the unphosphorylated KaiC ensures that phosphorylation proceeds rapidly even in the presence of low amounts of free KaiA (the first equation in System (C)). As the phosphorylation continues, the number of KaiC proteins phosphorylated on the *S* site increases (the second equation in System (C)). In turn, sequestration of KaiA increases through tight-binding (the last equation in System (C)).

The variable *C*_{T} is the total amount of the KaiC protein in the system. In the last equation, the free amount of KaiA, denoted by [*A*], is represented as a function of S431 phosphorylated KaiC that is derived under equilibrium assumptions similar to those in [14] (see Supplementary Information, Section 3 for a detailed derivation). As in the detailed model, the variable [*T*] gives the amount of KaiC protein phosphorylated on the *T* site. Similarly, the variables [*S*] and [*S**T*] reflect the amount of KaiC protein phosphorylated on the *S* site and both the *S* and *T* sites, respectively.

Our simulations demonstrate that KaiA sequestration through tight binding with S431 phosphorylated KaiC is indispensable for generating oscillations for various values of the binding affinity of KaiA and S431 phosphorylated KaiC (*K*_{d}). In particular, as the *K*_{d} value decreases (i.e., the binding affinity increases), the cusp of the KaiA concentration sharpens, and the magnitude of the sensitivity of the KaiA concentration to the amount of *S* increases (Fig. 3A, B). Increasing sensitivity, or order of reaction [15], of the activator KaiA with respect to S431 phosphorylated KaiC corresponds to an increasing likelihood of oscillations [15, 16]. Moreover, for increasing values of *K*_{d}, the fraction of parameter sets that generate oscillations decreases (Fig. 3C). This result shows that the system is more likely to generate stable oscillations when there is stronger KaiA sequestration (i.e., smaller *K*_{d}). A similar dependence on smaller *K*_{d} values is present in the interaction of the activator and the repressor in the mammalian system [14].

### A necessary stoichiometric ratio for oscillations in the simple cyanobacterial model

Additional theoretical analysis (see Supplementary Information, Section 3) and simulations confirm that a necessary condition for oscillations to occur in our model is

where *r* is a constant related to the rates of phosphorylation and dephosphorylation, and *ε* is a small positive constant. Consistent with existing experimental and modeling results [13, 17, 18], our simulations suggest that a balanced molar ratio between KaiA and KaiC abundance is crucial for generating sustainable oscillations. The observation that KaiB abundance does not affect the circadian rhythm in the same way is also consistent with the model design. In our simulations, we choose parameter values such that *r*=2 and *ε*<<1. The condition effectively turns into the linear criterion:

Simulation results verify that parameters violating Condition (1) rarely generate oscillations (Fig. 3D). A similar mechanism has been proposed for the mammalian circadian clock where repressors and activators are in a balanced molar ratio [14].

### Simple cyanobacterial model is functionally equivalent to the previous Kim-Forger model

System (C) reflects the interplay between phosphorylation and sequestration in the KaiABC oscillator in cyanobacteria. Correctly timed and ordered phosphorylation is necessary for robust circadian rhythm generation. Surprisingly, even though the mechanisms that generate rhythms are different among cyanobacteria and higher-order eukaryotes, the form of the repression function is similar. In particular, recall the simple Kim-Forger model of the mammalian circadian clock from [14, 19]:

Here, *M* models the mRNA of the repressor protein, *P*_{C} is the cytosolic repressor protein, and *P* is the nuclear repressor protein. The parameter *K*_{d} is the binding affinity of the activator *A* and the nuclear repressor, *P*. The function *f*(*P*,*A*,*K*_{d}) models the amount of free activator as a function of total activator, estimated by *A*, and repressor (*P*) present [14]. Note that the repression functions in Systems (C) and (M) are equivalent, indicating that the action of repression is conserved across the organisms even though the specific mechanisms of repression are different.

### Additional transcription-translation feedback loop adds robustness to oscillations in cyanobacteria

We also find that an additional transcription-translation feedback loop (TTFL) plays a significant role in sustaining robust oscillations. Besides the KaiABC system, the histidine kinases SasA and CikA regulate the output signaling from the post-translational oscillator to the transcriptional activity in the in vivo oscillator. The integrated roles of SasA, CikA, and RpaA together with KaiABC in the cyanobacterial circadian clock were summarized recently [20]. When KaiC phosphorylation reaches its peak, SasA binds to the CI domain of the ST-phosphorylated KaiC [9, 21], autophosphorylates, and transfers the phosphate group to RpaA, turning on the transcription factor [22]. During the dephosphorylation phase, when KaiB binds to the CI domain to form the KaiBC complex, SasA is released from KaiC. CikA is then recruited by the KaiBC complex, dephosphorylates RpaA, and thus, inhibits the transcription. In other words, the transcription of the *kaiBC* gene is activated when most of the KaiC proteins are highly phosphorylated and is inhibited when KaiC proteins are mostly in the *S* state. Therefore, we model the TTFL as an inhibition scheme where S431 phosphorylated KaiC acts like an inhibitor of the *kaiBC* gene (Fig. 1C, System (TTFL) below).

The variable *M* is the amount of mRNA transcribed from the *kaiBC* gene.

We investigate the role of the transcription-translation feedback loop (System (TTFL)) as an additional negative feedback loop in the KaiC system by comparing with a modified system with a constitutive source of transcription (Fig. 1D, System (PTR) below). See Table S1 (Additional file 1: Supplementary Information) for a detailed description of the parameters and the ranges of parameters used.

We simulate both Systems (TTFL) and (PTR) with 22,500 randomly generated parameter sets varying the transcription rate, *K*_{s}, and the total KaiA concentration, *A*_{T}. The TTFL model generates circadian oscillations under a wide range of parameters (98.84%, Fig. 3E), while the PTR model, without the additional feedback loop, only exhibits oscillations when satisfying Condition (1) (60.33%, Fig. 3F). Altogether, our analysis and simulations suggest that the TTFL, as an additional negative feedback loop in the cyanobacterial clock, can help sustain the required molar ratio balance through a homeostatic mechanism.

### A novel “Phospholock” model of the eukaryotic clock

In light of recent results that reveal a complex interplay between binding of the activator and repressor as well as the timing and ordering of phosphorylation of the repressors [23–25], we introduce an extension of System (M) for higher-order eukaryotes that incorporates additional phosphorylation of the repressor after binding to the activator. In particular, the repressor complex, after binding with the activator, is phosphorylated, subsequently leading to differential binding affinities between the activator and repressor complexes at various phosphorylation states and dissociation of the two complexes occurring after sufficient phosphorylation. We call this mechanism the “phospholock” as the additional phosphorylation of the repressor in turn keeps the activator and repressor complex sequestered together more so than pure protein sequestration. Notably, after phosphorylation is incorporated, the mathematical model of the repression mechanism is equivalent to those Systems (C) and (M) up to additional coefficients as shown by the model below (System (E)).

In this phospholock model, the activator complex, A, promotes the transcription of the repressor gene. The repressor mRNA, M, is then translated to the repressor protein r, which forms a final complex R with other proteins, e.g., CK1. Next, the repressor complex binds to A, leading to the subsequent phosphorylation of the repressor complex. After sufficient phosphorylation, the activator and repressor complexes dissociate, leaving the activator free to promote transcription of the repressor gene again. In this way, we incorporate phosphorylation, which has been shown to be essential to oscillations, with the previously established protein sequestration repression mechanism.

We describe the phospholock mechanism with a simple model below (see Supplementary Information, Section 4, for the derivation of the repression function, *f*(*R*)).

In Fig. 4B, we plot the repression function *f*(*R*) for a specific parameter set where *K*_{d}=1 while varying the ratio of phosphorylation to dephosphorylation (i.e., \(k_{2} = {\frac {k_{2f}}{k_{2r}}}\)). As *k*_{2} increases, the cusp of the repression function sharpens (Fig. 4B, top). Moreover, the magnitude of the peak sensitivity increases as *k*_{2} increases (Fig. 4B, bottom).

Next, we investigate how the interplay between the *K*_{d} value and the ratio of the strength of phosphorylation to dephosphorylation (*k*_{2}) affects the robustness of oscillations. As *K*_{d} values decrease, the likelihood that the system generates oscillations increases in general (Fig. 4C). Moreover, at small *K*_{d} values, the system is more likely to generate oscillations when dephosphorylation is stronger than phosphorylation, i.e., *k*_{2}<1. However, as *K*_{d} increases, e.g., *K*_{d}≥0.1, the system is more likely to generate oscillations when phosphorylation is stronger than dephosphorylation, i.e., *k*_{2}>1 (Fig. 4C).

### Additional activator phosphorylation in the phospholock model allows for wider stoichiometric ratios

While System (E) represents the mechanisms of the circadian clock in Drosophila and mammals, we now investigate the mechanism in *Neurospora*. The repression mechanism in the *Neurospora* circadian clock historically has been regarded as a Hill-type repression mechanism over a pure protein sequestration mechanism [19]. In this Hill-type repression, the repressor inactivates the activator by catalyzing several successive phosphorylations [26]. Moreover, the higher number of phosphorylations results in a higher chance of oscillations occurring [19]. Even though the assumptions used to derive the Hill function are very restrictive [16], mathematical models of the *Neurospora* clock [27–29] continue to use a Hill-type repression term to model the mechanism by which FRQ represses WCC over protein sequestration for two key reasons: 1) phosphorylation of the activator is a crucial element in the generation of oscillations [30–32], and 2) the stoichiometric ratio of repressor to activator in *Neurospora* is not consistent with required stoichiometric ratios in protein sequestration models [33].

However, when we add phosphorylation of the activator to the phospholock model (Fig. 4D), we find that oscillations are possible at a wide range of stoichiometric ratios. System (E), now with added phosphorylation of the activator, becomes

Note that, even with the additional phosphorylation of the activator, the repression functions in Systems (C), (E), and (N) are equivalent up to changes in parameter coefficients. See Supplementary Information, Section 5, for the derivation of the repression function.

Next, we simulate thousands of parameter sets that generate oscillations for increasing *k*_{3} values and then compute the stoichiometric ratio given by each parameter set. As *k*_{3} increases above zero, the range of stoichiometric ratios that sustain oscillations increases (Fig. 4E). In fact, the stoichiometric ratio can attain values very close to zero (Fig. 4E), values well within the range of measured levels of nuclear FRQ and WCC in *Neurospora* [34]. Thus, adding in the phosphorylation of the activator can significantly decrease the number of repressors required for the system to exhibit oscillations. In this case, phosphorylation acts to inactivate WCC without requiring the presence of FRQ. Moreover, increasing *k*_{3} actually decreases the robustness of the oscillations in this case (see Supplementary Information, Section 6).

## Discussion

Here, we present mathematical models of the cyanobacterial, mammalian, and Neurospora circadian clocks that reflect the crucial interplay between phosphorylation and sequestration of activators and repressors. Our models ultimately reveal a common repression function that generates oscillations at the molecular level independently in the three systems. This consistent repression function is particularly surprising considering the divergent mechanisms by which oscillations are generated, e.g., purely post-translational (cyanobacteria) versus a transcription-translation feedback loop (eukaryotes).

Moreover, our results reveal stoichiometric balancing requirements for oscillations and how these requirements may be satisfied in the three systems. For example, an additional transcription-translation feedback loop (TTFL) in cyanobacteria plays a significant role in sustaining robust oscillations by maintaining the required molar ratio (Condition (1)). This additional feedback loop matches the proposed role of the Rev-erb *α* TTFL in the mammalian clock [14, 19] that acts synergistically with the core feedback loop to maintain the molar ratio between activators and repressors. For example, when the molar ratio is perturbed to greater than 1 (i.e., more repressor), the additional feedback loop promotes activator production, restoring the 1:1 molar ratio required for robust oscillations [19]. Similarly, when the molar ratio is perturbed in the cyanobacterial system to greater than its required range (i.e., more free KaiC), then RpaA is dephosphorylated leading to low transcription of the *kaiBC* genes and a decrease in free KaiC, restoring the molar ratio required for robust oscillations.

To compare oscillations across species, we do not specifically account for the structure of some clock complexes. Accounting for this in molar ratios (e.g., the KaiC hexamer) shows that the requirements we see are consistent with other previous work (e.g., [18]). That being said, the mechanisms we discuss are far from a complete picture of circadian timekeeping. In fact, factors like co-operative binding or ultrasensitivity also play key roles [35, 36].

In *Neurospora*, it has been shown experimentally that phosphorylation of the activator is crucial to the generation of rhythms, and stoichiometric ratios of the repressor to activator do not agree with the ranges required from the protein sequestration paradigm. However, when additional phosphorylation of the activator is added to the phospholock model, the reconciliation among the possible models emerges. Additional phosphorylation of the activator allows for another method by which the repressor may inactivate the activator, namely successive phosphorylations. Once the activator is phosphorylated, it must be dephosphorylated to promote transcription again. Thus, fewer repressors are required to attain the required amount of inactivation of the activator than when there is no phosphorylation of the activator. That is, with no phosphorylation, once the repressor is degraded, free activators may promote transcription, so another repressor is required to sequester the activator. However, if the activator is phosphorylated, then another repressor does not need to bind to the activator. Thus, fewer repressors are required leading to a wider range in the stoichiometric ratios that generate oscillations.

Recent work using mice showed biochemically that PER promotes the nuclear entry of CRY and the phosphorylation of CLOCK [37]. This additional phosphorylation of CLOCK prevents it from binding to DNA, resulting in what the authors call “displacement-type repression” [37]. This is mechanistically identical to previous work showing that FRQ-mediated phosphorylation of specific sites on WC-1 and WC-2 is required for repression of WCC [32]. The two mechanisms may be identical, i.e., phosphorylation of the activator complex effected by the repressor complex brings about repression of the activator by preventing it from binding to DNA.

In previous protein sequestration models, the robustness of oscillations relied heavily on a tight binding affinity between the activator and repressor. In this way, the repressor sufficiently sequesters the activator, rendering it unable to promote transcription. However, these previous models required unrealistic binding affinities. In the phospholock model, the inactivation of the activator by the repressor is balanced between binding of the two components and phosphorylation/dephosphorylation of the repressor. That is, at low binding affinities of the activator and repressor, the phospholock model reveals that stronger phosphorylation adds robustness to the oscillations by keeping the activator inactivated. In contrast, at high binding affinities, stronger dephosphorylation adds robustness by freeing some activators to promote transcription.

The phospholock model also revealed that rapid degradation of the repressor is not required for sustained oscillations, confirming previous experimental results [23, 24, 38]. That is, simulations revealed that decreasing the *k*_{3} value, a surrogate marker of the degradation strength of the phosphorylated repressor, increases the sensitivity of the repression function, thereby increasing the robustness of oscillations. Therefore, our results are consistent with the important feature, now shown experimentally in both *Neurospora* [38] and mice [23], that rapid degradation of the phosphorylated repressors is not required for rhythm generation.

Altogether, our novel models of the cyanobacteria and eukaryotic circadian clocks corroborate the emerging theory that molecular clock mechanisms in divergent organisms may have evolved convergently [39]. In particular, a simple mechanism of the cyanobacteria clock is as follows: KaiA is required for the phosphorylation of KaiC on the *T* site. When the subsequent *S* site of KaiC is phosphorylated, KaiC binds to KaiB and gets “activated”. The “activated” KaiC can then bind and inactivate KaiA with the help of KaiB. We assume that KaiC, when bound to KaiB, is very efficient in inhibiting KaiA, such that effectively all KaiA is sequestered even in the presence of a small amount of KaiC with the *S* site phosphorylated. After the phosphorylation is completed and KaiC gets dephosphorylated at the *S* site, there is eventually not enough KaiC to sequester all KaiA, and KaiA is released. The free KaiA then efficiently phosphorylates the unphosphorylated KaiC, starting the cycle again. The eukaryotic systems work via a similar mechanism. An activator complex, A, activates the transcription of repressor proteins, which is functionally identical to KaiA activating KaiC on the *T* site. Next, the eukaryotic repressor proteins sequester and inactivate the activating complex. Similarly, the phosphorylation of KaiC on the *S* site allows KaiC to sequester and inactivate KaiA. Finally, the repressor proteins in eukaryotes undergo several phosphorylations and then dissociate, freeing the activating complex, completing the cycle. Similarly, as KaiC is dephosphorylated below a certain threshold, it releases KaiA, which then facilitates KaiC phosphorylation on the *T* site, completing the cycle. Despite differences among the three clocks, it is surprising that the systems’ core mechanisms share a similar architecture for generating robust oscillations.

Follow-up work could consider simulated evolution in silico experiments that may reveal efficient and robust mechanisms for generating oscillations across independent species. Such mechanisms could help further delineate how clocks may have evolved. Additionally, the plant molecular circadian clock is much more complex than that found in other organisms, at least in its TTFL structure. Further work should explore how sequestration and the mechanisms to balance molar ratios we discuss may play a role in the plant circadian clock. Further experiments could measure or change the molar ratios of activators and repressors to determine how the phospholock mechanism we propose extends the range where sustained oscillations are found.

## Conclusion

Our mathematical models reveal circadian clocks in cyanobacteria and some eukaryotes employ both protein sequestration and phosphorylation to generate oscillations. An additional transcription-translation feedback loop in cyanobacteria adds robustness to the oscillations over constitutive transcription of the *kaibc* genes. Interestingly, the Vrille or Rev-erb transcription-translation feedback loops in Drosophila and mammals, respectively, may have rediscovered this cyanobacterial mechanism to add robustness oscillations. In fungi and mammals, phosphorylation acts as a lock to keep repressor and activator complexes coupled and thus adds robustness to the oscillations, especially at higher *K*_{d} values (lower protein affinity). Furthermore, the phospholock in *Neurospora* includes phosphorylation of the activator, which drives wider stoichiometric ratios of repressor to activator. Thus, our model reveals how *Neurospora* fits into the previous protein sequestration model for mammalian circadian clocks. Taken together, these results paint a picture of how circadian timekeeping may have evolved [40–42].

## Methods

### Computing the fraction of parameter sets for which System (C) generates oscillations

For 14 *K*_{d} values increasing from 10^{−4} to 10^{−1}, we randomly generate 10^{5} parameter sets based on the parameter values in Table S1 (Additional file 1: Supplementary Information). In particular, we fix the *k*_{2}, *k*_{3}, *k*_{4}, *V*_{trsp}, *V*_{d}, and *V*_{m} parameters and select all other parameters from a random uniform distribution based on the ranges in Table S1 (Additional file 1: Supplementary Information). Then, we use the fast Fourier transform (fft) function in MATLAB to assess whether the parameter sets are sufficient for System (C) to exhibit oscillations.

Next, we simulate the TTFL and PTR models each 22,500 times after randomly generating parameter sets when varying *K*_{s} and *A*_{T} concentrations based on the ranges found in Table S1 (Additional file 1: Supplementary Information). We calculate the period by first applying the FFT on the time course and then identifying the strongest frequency in the spectrum.

### Plotting sensitivity of repression functions

We plot the repression function *f*([*S*],*A*_{T},*K*_{d}) with respect to [*S*] from System (C) for decreasing values of *K*_{d} (10^{−2} to 10^{−4}) using the parameter values in Table S1 (Additional file 1: Supplementary Information) (Fig. 3A). Additionally, in Fig.3B, we plot the sensitivity of *f*([*S*],*A*_{T},*K*_{d}) (equal to \({\frac {\mathrm {d}\log f([S], A_{T}, K_{d})}{\mathrm {d}\log [A]}}\) [15, 16]) computed in MATLAB R2020b.

Similarly, we plot the repression function *f*(*R*) from System (E) using the parameter values found in Table 1 for \(k_{2} = \frac {k_{2f}}{k_{2r}}\) values ranging from 1 to 10^{4} (Fig. 4B, top) for repressor amounts from 75 to 150. We also plot the sensitivity of *f*(*R*) (equal to \({\frac {\mathrm {d}\log f(R)}{\mathrm {d}\log R}}\) [15, 16]) for the same range of repressor values and *k*_{2} values (Fig. 4B, bottom) computed in MATLAB R2020b.

### Counting parameter sets for which the phospholock model generates oscillations

We randomly generated 10^{4} parameter sets where *A*_{T}, *α*_{1}, *β*_{1}, …, *α*_{3}, and *β*_{3} were sampled from a uniform distribution between 0 and 100. Then, we sampled the parameters *k*_{1f} and *k*_{2f} from a uniform distribution between 0 and 1. To guarantee *K*_{d}=1, we set *k*_{1r}=*k*_{1f}. Additionally, we sampled *k*_{3} from a uniform distribution between 0 and 0.1.

Next, we simulated the phospholock model for each randomly generated parameter set using the function ode23tb in MATLAB. Specifically, we randomly initialized the initial conditions and ran the model for 300 time units. Then, we took the mean of the M values and simulated the model again. The second time, however, we ran an event location procedure where the ode solver saves all times that satisfy: (1) *M* is equal to the mean calculated from the initial run and (2) *M* is increasing. Then, we concluded that the model oscillated if

where *t**e*_{1} is the last time that the event took place, *t**e*_{2} is the second to last, and *t**e*_{3} is the third to last. We checked several plots of the oscillating parameter sets to make sure that our procedure selected parameter sets that generate oscillations.

Finally, we simulated the 10^{4} randomly generated parameter sets varying the phosphorylation strengths by setting *k*_{2r}=100·*k*_{2f} (strongest dephosphorylation, Fig. 2C; left), *k*_{2r}=10·*k*_{2f}, *k*_{2r}=*k*_{2f}, …, *k*_{2r}=0.0001·*k*_{2f} (strongest phosphorylation, Fig. 2C; right). Each time, we counted how many parameter sets generated oscillations.

### Stoichiometric ratio analysis

Beginning with *k*_{3}=0, we randomly generate parameter sets until we recovered 100 such that the phospholock model with added activator phosphorylation generates oscillations (see Eqn. (56) in Supplementary Information for the updated repression function in this case). In particular, we randomly sample *A*_{T}, *α*_{1}, *β*_{1}, …, *α*_{3}, and *β*_{3} from a uniform distribution between 0 and 100. Next, we sample *k*_{1f}, *k*_{2r}, and *k*_{4f} from a uniform distribution between 0 and 1. We set *K*_{d}=10^{−5} in this case, so *k*_{1r}=10^{−5}·*k*_{1f}. Additionally, we assume phosphorylation is stronger than dephosphorylation in the phospholock, so we set *k*_{2f}=10·*k*_{2r}. Finally, we set autophosphorylation and autodephosphorylation as equal, so *k*_{4r}=*k*_{4f}. Parameter sets were assessed for whether they generated oscillations based on the same method as above.

For each of the 100 parameter sets that exhibit oscillations, we compute the stoichiometric ratio in the following way. First, we solve the following equation for R.

where

The solution to Eqn. (2) gives the repressor concentration at the steady state. Then, the stoichiometric ratio is

where R solves Eqn. (2).

## Availability of data and materials

The codes for the phospholock analyses are provided in the GitHub repository (https://github.com/jptyler/Phospholock) under MIT license [43].

## References

Ouyang Y, Andersson CR, Kondo T, Golden SS, Johnson CH. Resonating circadian clocks enhance fitness in cyanobacteria. P Natl Acad Sci USA. 1998; 95(15):8660–64. https://doi.org/10.1073/pnas.95.15.8660. http://arxiv.org/abs/https://www.pnas.org/content/95/15/8660.full.pdf.

Dunlap JC. Molecular bases for circadian clocks. Cell. 1999; 96(2):271–90. https://doi.org/10.1016/S0092-8674(00)80566-8.

Alon U. Network motifs: theory and experimental approaches. Nat Rev Genet. 2007; 8(6):450–61. https://doi.org/10.1038/nrg2102.

Tseng R, Chang Y-G, Bravo I, Latham R, Chaudhary A, Kuo N-W, LiWang A. Cooperative KaiA-KaiB-KaiC interactions affect KaiB/SasA competition in the circadian clock of cyanobacteria. J Mol Biol. 2014; 426(2):389–402. https://doi.org/10.1016/j.jmb.2013.09.040.

Chang Y-G, Cohen SE, Phong C, Myers WK, Kim Y-I, Tseng R, Lin J, Zhang L, Boyd JS, Lee Y, Kang S, Lee D, Li S, Britt RD, Rust MJ, Golden SS, LiWang A. A protein fold switch joins the circadian oscillator to clock output in cyanobacteria. Science. 2015; 349(6245):324–28. https://doi.org/10.1126/science.1260031 http://arxiv.org/abs/https://science.sciencemag.org/content/349/6245/324.full.pdf.

Nishiwaki T, Satomi Y, Nakajima M, Lee C, Kiyohara R, Kageyama H, Kitayama Y, Temamoto M, Yamaguchi A, Hijikata A, Go M, Iwasaki H, Takao T, Kondo T. Role of KaiC phosphorylation in the circadian clock system of Synechococcus elongatus PCC 7942. P Natl Acad Sci USA. 2004; 101(38):13927–32. https://doi.org/10.1073/pnas.0403906101. http://arxiv.org/abs/https://www.pnas.org/content/101/38/13927.full.pdf.

Clodong S, Dühring U, Kronk L, Wilde A, Axmann I, Herzel H, Kollmann M. Functioning and robustness of a bacterial circadian clock. Mol Syst Biol. 2007; 3:90. https://doi.org/10.1038/msb4100128.

AXMANN IM, LEGEWIE S, HERZEL H, MINIMAL CIRCADIAN CLOCK MODEL A:54–64. https://doi.org/10.1142/9781860949920_0006. http://arxiv.org/abs/https://www.worldscientific.com/doi/pdf/10.1142/9781860949920_0006. https://www.worldscientific.com/doi/abs/10.1142/9781860949920_0006.

Tseng R, Goularte NF, Chavan A, Luu J, Cohen SE, Chang Y-G, Heisler J, Li S, Michael AK, Tripathi S, Golden SS, LiWang A, Partch CL. Structural basis of the day-night transition in a bacterial circadian clock. Science. 2017; 355(6330):1174–80. https://doi.org/10.1126/science.aag2516. http://arxiv.org/abs/https://science.sciencemag.org/content/355/6330/1174.full.pdf.

Rust MJ, Markson JS, Lane WS, Fisher DS, O’Shea EK. Ordered phosphorylation governs oscillation of a three-protein circadian clock. Science. 2007; 318(5851):809–12.

Lin J, Chew J, Chockanathan U, Rust MJ. Mixtures of opposing phosphorylations within hexamers precisely time feedback in the cyanobacterial circadian clock. P Natl Acad Sci USA. 2014; 111(37):3937–45. https://doi.org/10.1073/pnas.1408692111. http://arxiv.org/abs/https://www.pnas.org/content/111/37/E3937.full.pdf.

Phong C, Markson JS, Wilhoite CM, Rust MJ. Robust and tunable circadian rhythms from differentially sensitive catalytic domains. P Natl Acad Sci USA. 2013; 110(3):1124–29. https://doi.org/10.1073/pnas.1212113110. http://arxiv.org/abs/https://www.pnas.org/content/110/3/1124.full.pdf.

van Zon JS, Lubensky DK, Altena PRH, ten Wolde PR. An allosteric model of circadian KaiC phosphorylation. P Natl Acad Sci USA. 2007; 104(18):7420–25. https://doi.org/10.1073/pnas.0608665104. http://arxiv.org/abs/https://www.pnas.org/content/104/18/7420.full.pdf.

Kim JK, Forger DB. A mechanism for robust circadian timekeeping via stoichiometric balance. Mol Syst Biol. 2012; 8(1):630. https://doi.org/10.1038/msb.2012.62. http://arxiv.org/abs/https://www.embopress.org/doi/pdf/10.1038/msb.2012.62.

Thron CD. The secant condition for instability in biochemical feedback control–I, The role of cooperativity and saturability. B Math Biol. 1991; 53(3):383–401. https://doi.org/10.1016/S0092-8240(05)80394-5.

Forger DB. Biological Clocks, Rhythms, and Oscillations. Cambridge: The MIT Press; 2017.

Kageyama H, Nishiwaki T, Nakajima M, Iwasaki H, Oyama T, Kondo T. Cyanobacterial circadian pacemaker: Kai protein complex dynamics in the KaiC phosphorylation cycle in vitro. Mol Cell. 2006; 23(2):161–71. https://doi.org/10.1016/j.molcel.2006.05.039.

Mori T, Sugiyama S, Byrne M, Johnson CH, Uchihashi T, Ando T. Revealing circadian mechanisms of integration and resilience by visualizing clock proteins working in real time. Nat Commun. 2018; 9(1):3245. https://doi.org/10.1038/s41467-018-05438-4.

Kim JK. Protein sequestration versus Hill-type repression in circadian clock models. IET Syst Biol. 2016; 10:125–13510.

Swan JA, Golden SS, LiWang A, Partch CL. Structure, function, and mechanism of the core circadian clock in cyanobacteria. J Biol Chem. 2018; 293(14):5026–34. https://doi.org/10.1074/jbc.TM117.001433.

Gutu A, O’Shea EK. Two antagonistic clock-regulated histidine kinases time the activation of circadian gene expression. Mol Cell. 2013; 50(2):288–94. https://doi.org/10.1016/j.molcel.2013.02.022.

Takai N, Nakajima M, Oyama T, Kito R, Sugita C, Sugita M, Kondo T, Iwasaki H. A KaiC-associating SasA–RpaA two-component regulatory system as a major circadian timing mediator in cyanobacteria. P Natl Acad Sci USA. 2006; 103(32):12109–14. https://doi.org/10.1073/pnas.0602955103. http://arxiv.org/abs/https://www.pnas.org/content/103/32/12109.full.pdf.

Ode KL, Ukai H, Susaki EA, Narumi R, Matsumoto K, Hara J, Koide N, Abe T, Kanemaki MT, Kiyonari H, Ueda HR. Knockout-rescue embryonic stem cell-derived mouse reveals circadian-period control by quality and quantity of cry1. Mol Cell. 2017; 65(1):176–90. https://doi.org/10.1016/j.molcel.2016.11.022.

Dunlap JC, Loros JJ. Just-So Stories and Origin Myths: Phosphorylation and Structural Disorder in Circadian Clock Proteins. Mol Cell. 2018; 69(2):165–68. https://doi.org/10.1016/j.molcel.2017.11.028.

Baker CL, Kettenbach AN, Loros JJ, Gerber SA, Dunlap JC. Quantitative proteomics reveals a dynamic interactome and phase-specific phosphorylation in the

*Neurospora*circadian clock. Mol Cell. 2009; 34(3):354–63. https://doi.org/10.1016/j.molcel.2009.04.023.Gonze D, Abou-Jaoudé W. The Goodwin model: Behind the Hill function. PLoS ONE. 2013; 8(8):1–15. https://doi.org/10.1371/journal.pone.0069573.

Hong CI, Jolma IW, Loros JJ, Dunlap JC, Ruoff P. Simulating dark expressions and interactions of frq and wc-1 in the Neurospora circadian clock. Biophys J. 2008; 94(4):1221–32. https://doi.org/10.1529/biophysj.107.115154.

Bellman J, Kim JK, Lim S, Hong CI. Modeling reveals a key mechanism for light-dependent phase shifts of Neurospora circadian rhythms. Biophys J. 2018; 115(6):1093–102. https://doi.org/10.1016/j.bpj.2018.07.029.

Upadhyay A, Brunner M, Herzel H. An inactivation switch enables rhythms in a Neurospora clock model. Int J Mol Sci. 2019; 20(12):2985. https://doi.org/10.3390/ijms20122985.

He Q, Shu H, Cheng P, Chen S, Wang L, Liu Y. Light-independent phosphorylation of white collar-1 regulates its function in the neurospora circadian negative feedback loop. J Biol Chem. 2005; 280(17):17526–32. https://doi.org/10.1074/jbc.M414010200. http://arxiv.org/abs/http://www.jbc.org/content/280/17/17526.full.pdf+html.

He Q, Cha J, He Q, Lee H-C, Yang Y, Liu Y. CKI and CKII mediate the FREQUENCY-dependent phosphorylation of the WHITE COLLAR complex to close the Neurospora circadian negative feedback loop. Gen Dev. 2006; 20(18):2552–65. https://doi.org/10.1101/gad.1463506. http://arxiv.org/abs/http://genesdev.cshlp.org/content/20/18/2552.full.pdf+html.

Wang B, Kettenbach AN, Zhou X, Loros JJ, Dunlap JC. The phospho-code determining circadian feedback loop closure and output in Neurospora. Mol Cell. 2019; 74(4):771–7843. https://doi.org/10.1016/j.molcel.2019.03.003.

Schafmeier T, Haase A, Káldi K, Scholz J, Fuchs M, Brunner M. Transcriptional feedback of

*Neurospora*circadian clock gene by phosphorylation-dependent inactivation of its transcription factor. Cell. 2005; 122(2):235–46. https://doi.org/10.1016/j.cell.2005.05.032.Hong CI, Ruoff P, Loros JJ, Dunlap JC. Closing the circadian negative feedback loop: Frq-dependent clearance of wc-1 from the nucleus. Gen Dev. 2008; 22(22):3196–204. https://doi.org/10.1101/gad.1706908. http://arxiv.org/abs/http://genesdev.cshlp.org/content/22/22/3196.full.pdf+html.

Hong L, Lavrentovich DO, Chavan A, Leypunskiy E, Li E, Matthews C, LiWang A, Rust MJ, Dinner AR. Bayesian modeling reveals metabolite-dependent ultrasensitivity in the cyanobacterial circadian clock. Mol Syst Biol. 2020; 16(6):9355. https://doi.org/10.15252/msb.20199355. http://arxiv.org/abs/https://www.embopress.org/doi/pdf/10.15252/msb.20199355.

Heisler J, Swan JA, Palacios JG, Sancar C, Ernst DC, Spangler RK, Bagshaw CR, Tripathi S, Crosby P, Golden SS, Partch CL, LiWang A. Structural mimicry confers robustness in the cyanobacterial circadian clock. bioRxiv. 2020. https://doi.org/10.1101/2020.06.17.158394. https://www.biorxiv.org/content/early/2020/06/19/2020.06.17.158394.full.pdf.

Cao X, Yang Y, Selby CP, Liu Z, Sancar A. Molecular mechanism of the repressive phase of the mammalian circadian clock. P Natl Acad Sci USA. 2021; 118(2). https://doi.org/10.1073/pnas.2021174118. http://arxiv.org/abs/https://www.pnas.org/content/118/2/e2021174118.full.pdf.

Larrondo LF, Olivares-Yañez C, Baker CL, Loros JJ, Dunlap JC. Decoupling circadian clock protein turnover from circadian period determination. Science. 2015; 347(6221). https://doi.org/10.1126/science.1257277. http://arxiv.org/abs/https://science.sciencemag.org/content/347/6221/1257277.full.pdf.

Partch CL. Orchestration of circadian timing by macromolecular protein assemblies. J Mol Biol. 2020; 432(12):3426–48. https://doi.org/10.1016/j.jmb.2019.12.046. Circadian Regulation: from Molecules to Physiology.

Tyson JJ, Othmer HG. The dynamics of feedback cellular control circuits in biochemical pathways. Prog Theor Biol. 1978; 5:1–62.

Yang X. Generalized form of Hurwitz-Routh criterion and Hopf bifurcation of higher order. Appl Math Lett. 2002; 15(5):615–21. https://doi.org/10.1016/S0893-9659(02)80014-3.

Parico GCG, Perez I, Fribourgh JL, Hernandez BN, Lee H-W, Partch CL. The human CRY1 tail controls circadian timing by regulating its association with CLOCK:BMAL1. P Natl Acad Sci USA. 2020. https://doi.org/10.1073/pnas.1920653117. http://arxiv.org/abs/https://www.pnas.org/content/early/2020/10/23/1920653117.full.pdf.

Tyler J, Lu Y, Dunlap J, Forger D. Phospholock. GitHub. 2021. https://github.com/jptyler/Phospholock. Accessed 30 Dec 2021.

### Review history

The review history is available as Additional file 2.

## Funding

This work was supported by a National Institutes of Health Training Grant (T32 HL007622), by the National Institutes of Health (R35GM118021), and by the National Science Foundation (DMS 1714059 and DMS 2052499).

## Author information

### Authors and Affiliations

### Contributions

JT, YL, JD, and DF designed the study. JT and YL performed the mathematical analysis. JT, YL, JD, and DF analyzed the results and wrote the manuscript. All authors read and approved the final manuscript.

### Corresponding author

## Ethics declarations

### Ethics approval and consent to participate

Not applicable.

### Competing interests

The authors declare they have no competing interests.

## Additional information

### Publisher’s Note

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

## Supplementary Information

**Additional file 1**

A separate pdf file that includes more detailed mathematical derivations and analysis of the models presented in the main text, Figures S1-S4, and Tables S1 and S2.

**Additional file 2**

Review history.

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

## About this article

### Cite this article

Tyler, J., Lu, Y., Dunlap, J. *et al.* Evolution of the repression mechanisms in circadian clocks.
*Genome Biol* **23**, 17 (2022). https://doi.org/10.1186/s13059-021-02571-0

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/s13059-021-02571-0