Evolution of the repression mechanisms in circadian clocks

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. Supplementary Information The online version contains supplementary material available at (10.1186/s13059-021-02571-0).

The reactions in the detailed model are: (S1) The corresponding mass actions equations are: The KaiBC complex sequesters KaiA when the S site is phosphorylated on KaiC (BS and BST ), where n = 2 indicates the strength of sequestration. To study the stability of our model under ATP variations, we adopt the competitive inhibition Figure S1: Simulation of the detailed model. Oscillations of KaiA (Yellow), KaiBC complexes in different stages (Blue and Magenta) along with the overall phosphorylation level (Black). The proportion of phosphorylated KaiC oscillates roughly between 73% and 25% on a circadian period. KaiA is mostly sequestered by KaiBC complexes (mainly by B-S) except for a short phase, which correspond to a rapid rise in the phosphorylation level. See Table S1 for detailed description of the parameters.
Here, the parameters k 0 ps and k 0 cat are the constant rates under 100% ATP, and K I is the strength of the inhibition from ADP. Simulations of the model under various [ATP]/([ATP]+[ADP]) ratios show that the cyanobacterial clock is robust to ATP variations ( Figure 2E-F, main text), which is consistent with previous experimental results [12].

Bifurcation analysis and stochastic simulations
Next, we generate a bifurcation diagram for the detailed model using XPP-AUTO. Our bifurcation diagram is similar to that in [13]. In particular, oscillations occur when the ratio between the concentration of KaiA and KaiC is within a bounded range while the relative concentration of KaiB has little effect on whether there are oscillations ( Figure S2).
We also simulated our model with the Gillespie algorithm (kinetic Monte Carlo). By keeping the KaiC total concentration a constant C T = 1µM , the number of KaiC molecules is directly related to the total volume in our simulation. We simulate the discrete model for various C T and compare the corresponding phosphorylation levels of KaiC ( Figure S3). When the number of molecules is low (C T = 10, 20, 50), the trajectory exhibits randomness without any sustained oscillations ( Figure  S3). As the number of KaiC molecules increases, however, the corresponding oscillations become more stable. When the number of molecules C T becomes large enough, the general shape of the profile stays almost the same even if the number is doubled from C T = 5000 to C T = 10000. The amplitude of the oscillations from stochastic simulations approaches that of the deterministic model as the number of molecules increases. Moreover, the phase difference between the stochastic and deterministic simulations decreases significantly as the number of molecules increases. Therefore, we predict that to observe stable and synchronized oscillations, the total amount of KaiC proteins must be above a certain threshold.
3 Linear stability analysis of the core model Here, we present a linear stability analysis of the simplified model. We continue to use [T ], [ST ], and [S] as the state variables for KaiC concentrations in various phosphorylation states, and the variable [A] for the concentration of KaiA. Constants C T and A T indicate the total amount of the KaiC and KaiA proteins, respectively, and K d is the dissociation constant in the sequestration of KaiA through KaiC-S. We assume that sequestration happens on a faster timescale, thus reaching equilibrium quickly. Due to the two time scales and the binding of KaiA having no impact on other reactions of S, we solve for the concentration of free KaiA. Figure S2: Two parameter bifurcation diagram of the KaiC system.. The parameter atoc is the fraction of KaiA concentration to KaiC concentration. The parameter btoc is the fraction of KaiB concentration to KaiC concentration. Points within the region enclosed by the blue curve give initial distributions that produce oscillations while we expect no oscillations outside this region.
First, we reduce the number of variables by applying the following conservation law (S14) After applying Equation (S14), the core model reduces to the following set of equations.
K d must be small enough for the system to generate sustainable oscillations. Therefore, since K d << A T , we take the limit as [S], 0}. Then, Equations (S15)-(S18) reduce to: Solving the above system for the equilibrium solutions, we obtain from Equation (S20) and Equation (S21), respectively. Solving for the equilibrium solution of Equation (S19), we obtain a quadratic in terms of [S]: Additionally, we assume k 4 << k 1 since KaiA-enhanced phosphorylation is much faster than (auto)phosphorylation. Hence, we drop the last term in Equation (S22). After dividing k 1 through all terms, we find two steady state solutions: Next, we compute the Jacobian at the equilibrium [S] * = A T : The characteristic function of the Jacobian, J 1 , is then: From the secant condition [40], we conclude that oscillations occur when In the main text, we write condition (S28) as C T − (1 + r)A T ≥ , where r = k 4 k 2 + k 4 k 3 and = 8k 4 k 1 . To verify our analysis, we simulate the model for k 2 = k 3 = k 4 = 0.1 and compare our results with Next, we analyze the Jacobian matrix at the other equilibrium [S] * = C T 3 . Without introducing too many computational details, we assume that k 2 = k 3 = k 4 , which means that autodephosphorylation and autophosphorylation of KaiC are on the same time scale.
The Jacobian of [S] * = C T 3 is: The characteristic polynomial of the Jacobian, J 2 , is a cubic function: The Routh-Hurwitz Stability Criterion [41] states that a necessary condition for a cubic function h(x) = x 3 + a 2 x 2 + a 1 x + a 0 to be stable (i.e., all roots have negative real part) is a 0 > 0, a 2 > 0, and a 2 a 1 > a 0 .
Applying the criterion to Equation (S29), we have a 0 > 0 and a 2 > 0. We check the third condition accordingly: Thus, we know that all three roots of the characteristic function always have negative real parts, implying that the steady state is stable. Therefore, there are no oscillations around the steady state [S] * = C T 3 .
To conclude, we found a necessary condition for the system to generate oscillations, which we also verified by simulations.

Derivation of the repression function in the novel phospholock model
From Figure 4A in the main text, we have the following reaction network.
Network (S31) generates the following system of ODEs.
System (S32) generates the following steady-state equations.
Additionally, we have the following two conservation laws in total activator and total repressor concentrations.
A T = A + AR + ARp (S36) where A T denotes the total concentration of A and R T denotes the total concentration of R.
Using Eqn. (S35), we get Combining the conservation law (S36) with Eqn. (S38), we get an updated conservation laws Moreover, substituting ARp based on Eqn. (S35) into the steady-state equation (S33) gives For simplification, we setK Then, we can rewrite the conservation laws as Substituting the conservation laws (S44) and (S45) into Eqn. (S41) gives the quadratic equatioñ Expanding out Eqn. (S46) gives The solution of Eqn. (S47) is We take the negative solution above to get that Finally, we replace the solution for AR in Eqn. (S48) into the conservation law (S44) to get the concentration of free activator.
In the main text, we replace R T with R as done previously [14,16].

Derivation of the repression function in the adjusted phospholock model with activator phosphorylation
From Figure 4D in the main text, we have the following reaction network.
Network (S49) generates the following system of ODEs.
System (S50) generates the following steady-state equations.
Additionally, we have the following two conservation laws in total activator and total repressor concentrations.
A T = A + AR + ApRp + Ap (S55) where A T denotes the total concentration of A and R T denotes the total concentration of R.
Moreover, substituting Ap based on Eqn. (S58) into the steady-state equation (S51) gives whereK 2 is defined as in (S43). For simplification, we set Then, we can rewrite the conservation laws as Substituting the conservation laws (S63) and (S64) into Eqn. (S61) gives the quadratic equatioñ Expanding out Eqn. (S65) gives The solution of Eqn. (S66) is We take the negative solution above to get that Finally, we replace the solution for AR in Eqn. (S67) into the conservation law (S63) to get the concentration of free activator. ( 6 Sensitivity of transcription-rate function decreases with increasing k 3 values Additionally, we performed a comprehensive analysis of how the peak in sensitivity varies with respect to k 3 . Previous studies illustrated that a higher sensitivity increases the likelihood of oscillations [14]. We saw for individual parameter sets, for which System (N) exhibited oscillations, that the magnitude of the peak sensitivity decreases as k 3 increases ( Figure S4A). Thus, we further investigated whether this is true in general. In particular, for 1000 randomly generated parameter sets for which System (N) exhibited oscillations, we calculate the magnitude of the peak sensitivity first at the k 3 value of the parameter set. Then, we vary k 3 values (x-axis, Figure S4B) and report that variation as percent change from the original k 3 value of the parameter set. Next, we compute the magnitudes of the peak sensitivity for the varied k 3 values and report those as percent variation from the original magnitude (y-axis, Figure S4B). Thus, at 100%, which is the original k 3 value, the magnitude should be 100%, i.e., equal to the original amount. As k 3 increases relative to the Figure S4: The robustness of oscillations decreases as the rate of dissocation increases.(A) The sensitivity of the transcription-rate function of System (N) with one parameter set for which the system exhibited oscillations. We vary the k 3 value to represent an increase in the dissociation rate of phosphorylated repressors. As k 3 increases (more green), the magnitude of the peak sensitivity decreases, i.e., the robustness of oscillations decreases. (B) We randomly generated 1000 parameter sets for which the Neurospora system exhibits oscillations. For each parameter set, we vary the k 3 value relative to the original k 3 value of the parameter set (x-axis records the percent variation). Then, we plot the variation in the magnitude of the peak sensitivity. We compute this variation for each parameter set and plot the mean change. As k 3 increases relative to the original value, the magnitude of the peak sensitivity decreases, indicating that faster dissociation of phosphorylated repressor and activator decreases the robustness of the oscillations.
original k 3 value, the magnitude in the peak sensitivity decreases relative to the original magnitude ( Figure S4B). Similarly, when the k 3 value is higher than the original value, the magnitude is also higher than the original amount ( Figure S4B). Thus, in general, as the rate of dissociation increases, the magnitude of the peak sensitivity decreases. Therefore, increasing the dissociation rate of the phosphorylated repressor and activator decreases the robustness of the oscillations. This phenomenon may explain why CRY1 binds with BMAL1/CLOCK on DNA independently of PER in the early morning [41].