Interaction between FRQ and CK1 proteins controls the timing of daily rhythms in Neurospora

To determine the mechanism of period phenotypes in these mutants, we first compared the interaction between FRQ and WC-2, a critical step in the circadian negative feedback process. Immunopreciptation of WC-2 showed that the amount of FRQ precipitated with WC-2 was similar after normalized with the Input FRQ level in the mutants and in the wild-type strain (Fig. 1c). Because CK1a forms a stable complex with FRQ and it not only phosphorylates FRQ but is also recruited by FRQ to mediate WC phosphorylation to close the circadian negative feedback loop^13,44,48, we examined the FRQ-CK1a interaction in these strains. Because FRQ phosphorylation is rhythmic in constant darkness (DD), we performed immunoprecipitation experiments using cultures grown in constant light (LL) which has constant FRQ phosphorylation profiles. Using a CK1a antibody that specifically recognizes the endogenous CK1a protein (Supplementary Fig. 1a), we demonstrated that the amount of FRQ associated with CK1a (normalized with the Input FRQ) was less in the mutants than the wild-type strain even though total FRQ levels in the mutants were higher (Fig. 1d). This indicated that the FRQ-CK1a interaction was significantly reduced in both mutant strains compared to the wild-type strain. In addition, the FRQ-CK1a interaction was lower in the M9 mutant (period 32.1 h) than in the M10 mutant (period 29.8 h). Thus, the reduced FRQ-CK1a interaction but not FRQ stability can potentially explain the long period phenotypes of both mutants. These results also suggest that the FRQ phosphorylation events impaired in these two mutants promote the FRQ-CK1a interaction.

We next examined two classic frq period mutants: frq¹ and frq⁷⁴⁹. frq¹ is a short period mutant and frq⁷ is a long period mutant, each due to a missense mutation in the FRQ open reading frame (Supplementary Fig. 1b). Their period phenotype was previously attributed to differences in FRQ stability: FRQ is more stable in frq⁷ and is less stable in frq¹⁵⁰. To evaluate the FRQ-CK1a interaction, we introduced a Myc-tagged CK1a expression construct into the mutants and wild-type strain. By analysis of total cell extracts prepared from cultures grown in LL, we found that levels of FRQ were reduced in the frq¹ strain and elevated in the frq⁷ strain compared to the wild-type strain (Fig. 1d, e and Supplementary Fig. 1c). In addition, FRQ-WC interaction was not affected in these mutants (Supplementary Fig. 1d). Immunoprecipitation of Myc-CK1a showed that compared to the wild-type strain, the FRQ-CK1a interaction was significantly increased in the frq¹ mutant but was decreased in the frq⁷ mutant (Fig. 1e).

To confirm this conclusion, the recombinant FLAG-tagged CK1a was incubated with Neurospora extracts from either the wild-type or mutant strains. Precipitation with the FLAG antibody showed that there was more FRQ associated with FLAG-CK1a in the frq¹ extract than in the frq⁷ extract (Fig. 1f). These results indicate that the FRQ-CK1a interaction negatively correlates with period length in these classic clock mutants. Due to the role of FRQ-CK1a interaction in mediating FRQ phosphorylation and its subsequent degradation, the altered FRQ-CK1a interaction explains the different FRQ stability in frq¹ and frq⁷ mutants. These results suggest that the altered FRQ-CK1a interaction, but not FRQ stability, is the primary cause of period phenotypes in these clock mutants.

Immunoprecipitation of CK1a from extracts of mutants showed that the period length negatively correlated with the FRQ-CK1a interaction: The longer the period length of the mutant, the weaker the interaction between FRQ and CK1a (Fig. 2d). In contrast, the FRQ-WC-2 interaction was not affected in these mutants (Supplementary Fig. 2a). These results suggest that the FRQ-CK1a interaction is a major determinant of circadian period length in Neurospora.

Comparison of FRQ phosphorylation profiles for cultures grown in constant light (LL) showed that FRQ became hypophosphorylated in strains carrying mutations in the FRQ FCDs (Fig. 2e). In addition, the degree of impact on FRQ phosphorylation correlated with the FRQ-CK1a interaction: The longer the period length, the more hypophosphorylated FRQ became. The 4A12 mutant was the most hypophosphorylated, and no FRQ-CK1a complex was detected in this mutant. FRQ in the L488V mutant (48 h period) was more hypophosphorylated than those in other mutants, consistent with its dramatic decrease of FRQ-CK1a interaction. Thus, unlike typical kinase-substrate reactions, CK1a must first bind tightly to FRQ prior to its phosphorylation of FRQ.

We compared the FRQ rhythmic phosphorylation profile between the WT and one of the FCD mutants (Q494N) in DD (Supplementary Fig. 2b). Consistent with the long period phenotype of the Q494N strain, the phase of its FRQ phosphorylation rhythm was markedly delayed compared to the WT strain. In addition, we found that period length correlated with FRQ stability in the above mutant strains (Fig. 2f). Furthermore, we also inducibly over-expressed CK1a in the WT and FCD mutants (Q494N and V320I) to compare their period length changes due to an increase of CK1a level. Although the induction of CK1a resulted in period shortening in all strains, the FCD mutants were more sensitive to CK1a over expression than WT (Supplementary Fig. 2c). These results suggest that the binding equilibrium of the FRQ-CK1a complex formation is important for period determination. Together, these results indicate that the FRQ-CK1a interaction is a critical step that determine the efficiency of FRQ phosphorylation, its stability and the period length.

Depending on the site, mutations of FRQ phosphorylation sites lengthen or shorten the circadian period length^9,14,15. The importance of the FRQ-CK1a interaction on period determination suggested that FRQ phosphorylation can also regulate FRQ-CK1 interaction. Supporting this conclusion, the FRQ-CK1a interaction was impaired in both M9 and M10 mutants (Fig. 1d), indicating that FRQ phosphorylation at these sites promotes the interaction between FRQ and CK1a, resulting in their long period phenotypes.

FRQ phosphorylation profiles exhibit a robust circadian rhythm: Newly synthesized FRQ is hypophosphorylated during early subjective morning, is progressively more phosphorylated during the day, and becomes hyperphosphorylated in the middle of subjective night⁵². Thus we compared the relative amount of FRQ associated with CK1a at different time points in constant darkness (DD). Consistent with a previously reported Mass spectrometry data¹⁴, our IP experiments showed that the relative amount of FRQ associated with CK1a exhibited a circadian rhythm (Fig. 3b): It was lowest when FRQ was hyperphosphorylated (DD28) and highest when FRQ phosphorylation was intermediate (DD22). When FRQ was hypophosphorylated (DD16 and DD34), the amount was intermediate. This result suggests that the effect of FRQ phosphorylation can have opposing effects on complex formation: progressive FRQ phosphorylation during the early stage promotes the interaction between FRQ and CK1a, but extensive phosphorylation inhibits the interaction. In contrast, the FRQ-WC-2 interaction was not significantly affected by FRQ phosphorylation levels (Supplementary Fig. 3a).

S513, S519, S683, S685, S885, and S887 are some of the previously identified FRQ phosphorylation sites^9,14,15. The strain with the mutation S513R and the previously described M13 strain (mutations at S683 and S685) resulted in long period phenotypes (Supplementary Fig. 3b). In contrast, the strain with the mutation S519A and the previously described M17 strain (mutations at S885 and S887) led to period shortening (Supplementary Fig. 3b). Compared to the wild-type strain, the relative level of the FRQ-CK1a complex was reduced in the S513R and M13 mutants but was increased in the S519A and M17 mutants (Fig. 3c, d), consistent with their period phenotypes. Moreover, the FRQ-WC-2 interaction was not affected in these mutants (Supplementary Fig. 3c–d).

FRQ is phosphorylated by CK2, but CK2 does not form a stable complex with FRQ^35,36,53. The disruption of the ckb gene, which encodes the regulatory subunit of CK2, led to hypophosphorylation of FRQ and a lengthening of the period by about 6 h (Supplementary Fig. 3b)³⁵. As expected, the FRQ-CK1a interaction was also significantly reduced in the ckb^RIP mutant, explaining its long period phenotype (Fig. 3e). Thus, FRQ phosphorylation at different sites can either positively or negatively regulate the FRQ-CK1a interaction. In addition, the period phenotypes of the FRQ phosphorylation site mutants correlate with the FRQ-CK1a interaction in different clock mutants.

We hypothesized that phosphorylation influences the conformation of FRQ, which in turn affects the binding of FRQ with CK1a. To test this, we compared limited trypsin digestion sensitivity of FRQ in Neurospora extracts from different time points. Hyperphosphorylated FRQ (DD28) was more sensitive to trypsin digestion than was hypophosphorylated FRQ (DD16) (Fig. 3f). This result suggests that phosphorylation induces changes of FRQ structure, which likely impacts FRQ-CK1a complex formation.

To determine whether the FRQ-CK1a interaction affects the circadian negative feedback process, we performed chromatin immunoprecipitation (ChIP) assays using a WC-2 antibody. The relative WC-2 enrichment at the frq C-box was constantly elevated in the 4A12 mutant compared to the wild-type strain at different time points in DD (Fig. 4c). As expected, intermediate levels of WC-2 enrichment at the frq C-box were observed in the Q494N mutant. In addition, the WC-2 enrichment was also higher in the frq⁷ strain and lower in the frq¹ strain compared to the wild-type strain (Fig. 4d), consistent with previously reported FRQ levels in these mutant strains⁵⁴. These results suggest that the FRQ-CK1a interaction determines the efficiency of WC phosphorylation by CK1a and, therefore, the DNA binding activity of WCC.

To confirm the role of the FRQ-CK1a interaction in determining FRQ activity in negative feedback process, we compared the ability of wild-type FRQ and FCD mutants in suppressing frq transcription. To allow side-by-side comparison and to avoid the interference from the endogenous FRQ protein, we inducibly expressed FRQ (wild-type, 4A12 or Q494N) in a frq⁹ strain from a construct induced by quinic acid (QA); this allowed control of FRQ expression in a QA concentration-dependent manner^33,55. The frq⁹ strain has a frame-shift mutation in the frq open reading frame that results in a truncated non-functional FRQ protein and constant high frq RNA level due to the loss of circadian negative feedback process⁴⁹. FRQ protein was induced to higher levels as QA concentration increased, and the FRQ levels were similar at each QA concentration in the three qa-FRQ strains (Fig. 4e). Northern blot analysis was performed to detect the endogenous frq⁹ mRNA. As expected for a normal circadian negative feedback process, the induction of wild-type FRQ expression resulted in a marked decrease of frq⁹ mRNA level at 1 × 10⁻⁴ M QA, and the frq mRNA level was very low at 1 × 10⁻³ M QA (Fig. 4e). In contrast, induction of FRQ.4A12 expression, even at 1 × 10⁻³ M QA, failed to repress frq⁹ mRNA levels (Fig. 4e), indicating the complete loss of circadian negative feedback loop due to lack of FRQ-CK1a complex formation. On the other hand, the induction of FRQ.Q494N resulted in a modest decrease of frq⁹ mRNA levels (Fig. 4e). Together, these results demonstrate that the FRQ-CK1a interaction determines FRQ activity in the circadian negative feedback loop.

Thus, the FRQ-CK1a interaction should influence the period length independent of FRQ stability by determining the time needed to close the negative feedback loop. Supporting this conclusion, we show that the mutants with stronger FRQ-CK1a interactions exhibit short period phenotypes while the mutants with weaker interactions are long period even when FRQ stability was not affected. Consistent with our findings, it was recently shown that mutations of CRY, a negative element in the mammalian circadian negative feedback loop, can result in period length changes independent of its stability in mouse cells¹⁹.

Our model contains 10 variables and 35 parameters. Because most of the parameter values are unknown, we performed an unbiased parameter search by using a Bayesian analysis with the Markov chain Monte Carlo estimation (Supplementary Fig. 4a–c), which generated multiple sets of default parameters that produce robust circadian oscillations of about 22 h (please see detailed description in “Methods”). Distribution of parameters from these sets of parameters are summarized as boxed plots in Supplementary Fig. 5, which indicated that each parameter has different spread with distinct mean values. We randomly chose three sets of independent parameters for in-depth sensitivity analysis to determine the sensitivity of period length to each parameter (Supplementary Tables 1, 2). To our surprise, the same top five parameters that determine the period in this system were identified with each set of parameters (Fig. 5b). These parameters can be divided into two categories: the rates that determine dynamics of the positive feedback loop that promote the association between FRQ and CK1 (α₁, α₂, α_c) as the top three parameters and the rates of frq mRNA and FRQ_p-CK1a degradations (δ₁, δ₅). Our representative data using the parameter set 1 shows that the rate of association between FRQ and CK1a, α₁, was the most sensitive parameter; reductions in α₁ caused periods of more than 60 h (Fig. 5c). Reduction of the rate of CK1 synthesis (α_c) and the rate of phosphorylation FRQ-CK1a complex (α₂) also increased the period to more than 40 h (Fig. 5d, e).

In our model, the rate of degradation of frq mRNA (δ₁) and FRQ degradation after phosphorylation (δ₅) moderately influenced the period, up to 36 h in case of δ₁ and only 27 h for δ₅ (Fig. 5f, g). In contrast, most of the other degradation rates, including the rates of degradation of hyperphosphorylated FRQ (δ₈ and δ₉), showed a large region of oscillatory domain with minimal changes in period (Fig. 5h and Supplementary Fig. 6). This is consistent with the fact that the hyperphosphorylated FRQ inhibits the FRQ-CK1a interaction, and thus it is not as functional as hypophosphorylated FRQ in the negative feedback loop. Thus, these in silico experiments suggest that FRQ degradation only has a modest role in period determination.

The dramatic period effect due to reduction of the rate of association between FRQ and CK1a (α₁) occurs because the system is approaching to a bifurcation point (i.e., a point in parameter space where the stability of a system changes altering the behavior of the system) at which an unstable steady state with a stable limit cycle collides with a saddle point (called a Saddle Node on Invariant Circle)^58,59, generating an oscillation with an infinite period (Fig. 5i). Biologically, the observed dramatic increase of period occurs when any of the rates that determine the accumulation of CK1-FRQ_p (α₁, α₂, α_c), which triggers positive feedback loop, are very low, resulting in a critical slowdown of the system. In other words, a certain amount of CK1-FRQ_p has to accumulate to trigger the positive feedback loop promoting the association of FRQ and CK1a. Therefore, if α₁, α₂, or α_c are low, then it takes much longer time to accumulate a critical concentration of CK1-FRQ_p that will induce the positive feedback loop. In contrast, this unique dynamical behavior is not observed with other parameters in the model. Bifurcation analysis of δ₁ showed that the oscillatory domain is bounded by two Hopf bifurcation points (Fig. 5j, solid green circles) with modest monotonical decrease of period with increasing δ₁^58,59. Representative one-parameter bifurcation diagrams for δ₁ and δ₈ are shown for comparison with that of α₁ in Fig. 5j, k.

The level of CK1 is known to be relatively constant, and the phosphorylation of FRQ by CK1a (α₂) depends on the prior interaction between these two proteins (α₁). Therefore, our computational analyses revealed that the regulated complex formation between FRQ and CK1a is the major parameter in determining the period of the Neurospora circadian clock. In addition, our in silico analyses suggest that the molecular wiring diagram is one of the critical constraints that determine appropriate parameter space for robust oscillations.

The 87–3 (bd, a) strain or a frq complementation strain (KAJ120) was used as the wild-type strain in this study. S513R, S519A, M9, M10, M13, M17, and ckb^RIP strains were created previously^9,15,35. The frq^KO (bd, his-3) strains were used as host strains for constructs with his-3-targeting vectors. The control strain was a frq^KO strain transformed with KAJ120, a vector containing the wild-type frq gene. All constructs bearing FCD1 and FCD2 mutations were transformed into the frq^KO strain. Constructs with the qa-2 promoter driving expression of Myc.His.CK1a and pBT6 were co-introduced into the KAJ120, S513R, S519A, and 4A12 strains by homologous recombination. Positive transformants were verified by western blot analyses. Constructs with the qa-2 promoter driving expression of Myc.His.CK1a.hph were introduced into the 87–3, frq¹, frq⁷, and ckb^RIP strains at the csr locus by homologous recombination. Constructs with the qa-2 promoter driving expression of Myc.His.CK1a.bar were introduced into the KAJ120, V320I and Q494N strains at the csr locus by homologous recombination. pBA50 constructs with the qa-2 promoter driving expression of WT, Q494N, or 4A12 was introduced into the frq⁹ (bd, his-3) strain at the his-3 locus by homologous recombination. Positive transformants were identified by western blot analyses, and homokaryon strains were isolated by microconidia purification using 5-µm filters (Millipore, Cat#SLSV025LS).

Liquid cultures were grown in minimal media (1× Vogel’s, 2% glucose). When QA (Sigma-Aldrich, Cat#138622-100 G) was used to activate the qa-2 promoter, liquid cultures were grown in 1 × 10⁻² M QA (pH 5.8), 1× Vogel’s, 0.1% glucose, and 0.17% arginine. For rhythmic experiments, Neurospora was cultured in petri dishes in liquid medium for 2 days. The Neurospora mats were cut into discs and transferred into medium-containing flasks and were harvested at the indicated time points.

The medium for race tube assays contained 1× Vogel’s salts, 0.1 or 0.03% glucose, 0.17% arginine, 50 ng/mL biotin (Sigma-Aldrich, Cat#B4639-500MG), and 1.5% agar with or without QA. After entrainment of 24 h in the LL condition, race tubes were transferred to the DD conditions and marked every 24 h.

Protein extraction, quantification, and western blot analyses were performed as previously described^51,52,66. Briefly, tissue was ground in liquid nitrogen with a mortar and pestle and suspended in ice-cold extraction buffer (50 mM HEPES (pH7.4), 137 mM NaCl, 10% Glycerol) with protease inhibitors at final concentrations: Pepstatin A (1 μg/mL), Leupeptin (1 μg/ml), PMSF (1 mM). For certain experiments with PPase inhibitors (made fresh): 25 mM NaF, 10 mM Na₄P₂O₇.10H₂O, 2 mM Na₃VO₄, 1 mM EDTA. After centrifugation, protein concentration was measured using protein assay dye (Bio-Rad, Cat#500-0006). For western blot analyses, equal amounts of total protein (40 μg) were loaded in each protein lane. After electrophoresis, proteins were transferred onto PVDF membrane (Immobilon-P, EMD Millipore, Cat#IPV1100010), and western blot analyses were performed. To analyze the phosphorylation profile of WC-2, 7.5% SDS-PAGE gels containing a ratio of 149:1 acrylamide/bisacrylamide were used. Otherwise, 7.5% SDS-PAGE gels containing a ratio of 37.5:1 acrylamide/bisacrylamide were used.

RNA was extracted with Trizol (Ambion, Cat#15596018) in accordance with the manufacturer’s protocol, and then further purified with 2.5 M LiCl (Sigma-Aldrich, L4408-500G)⁶⁷. For northern blot analyses, equal amounts of total RNA (20 μg) were loaded onto agarose gels. After electrophoresis, the RNA was transferred onto nitrocellulose membrane (GVS North America, Cat#1212590). The membrane was probed with [³²P] UTP (PerkinElmer, Cat#BLU007H250UC) -labeled RNA probes specific for frq 5′UTR. RNA probes were transcribed in vitro from PCR products by T7 RNA polymerase (Ambion) as per the manufacturer’s protocol. The frq primer sequences used for the template amplification were frq Northern F, 5′-TAATACGACTCACTATAGGGGAAGCCTGCGATGCTACTCAC-3′, and frq Northern R, 5′-CATCGCATCCACTTCGCACAC-3′. The list of primers used in this study is shown in Supplementary Table 3.

IP analyses were performed as previously described⁵¹. Briefly, Neurospora proteins were extracted as described above. For each immunoprecipitation reaction, 2 mg protein and 1.5 μL c-Myc antibody (Roche, Cat#11667203001), 1 μL WC-2 antibody⁴³, or 1 μL CK1a antibody were used. After incubation with antibody for 3 h, 40 μL GammaBind G Sepharose beads (GE Healthcare) were added, and samples were incubated for 1 h. For the FLAG pull-down assay, 4 mg purified FLAG-CK1a fusion protein and 2 μL FLAG antibody (Sigma-Aldrich, Cat#F1804) were incubated with 2 mg Neurospora extracts. After incubation with antibody for 3 h, 40 protein GammaBind G Sepharose beads were added, and samples were incubated for 1 h. Immunoprecipitated proteins were washed three times using extraction buffer before Western blot analysis. IP experiments were performed using cultures harvested in constant light (LL) except for the rhythmic experiments. Quantification of relative FRQ-CK1a or FRQ-WC-2 interaction levels are based on the ratio of IP to Input and normalized with CK1a or WC-2 levels, respectively.

For CHX (Sigma-Aldrich, Cat#C7698-5G) treatments, discs made from 2-day-old mycelial mats of individual strains were grown in liquid media for 1 day under constant light conditions. CHX was then added to a final concentration of 10 μg/mL, and tissues were harvested at the indicated time points.

For the trypsin-sensitivity assay, protein extracts were diluted to a protein concentration of 2.5 μg/μL. A 100-μL aliquot was treated with trypsin (final concentration 1 μg/mL) at 25 °C. A 20 μL sample was taken from the reaction at each time point (0, 5, 15, and 30 min) after addition of trypsin. Protein samples were mixed with protein loading buffer and resolved by SDS-PAGE. To compare trypsin sensitivity of FRQ from different strains, experiments were performed side by side, and the protein samples were transferred to the same membrane for western blot analysis.

ChIP assays were performed as previously described^48,68. Briefly, Neurospora tissues were fixed with 1% formaldehyde for 15 min at room temperature with shaking. Glycine was then added at a final concentration of 125 mM. The cross-linked tissues were ground and resuspended in lysis buffer (50 mM HEPES, pH 7.5, 137 mM NaCl, 1 mM EDTA, 1% Triton X-100, 0.1% deoxycholate, 0.1% SDS). Chromatin was sheared using sonication to 500–1000 base-pair fragments. For each immunoprecipitation reaction, 1 mg protein was used. WC-2 antibody was generated previously⁴³. The ChIP reaction was carried out with 2 μL WC-2 antibody. Immunoprecipitated DNA was enriched using GammaBind G Sepharose beads and washed sequentially with low-salt washing buffer (0.1% SDS, 1% Triton X-100, 2 mM EDTA, 150 mM NaCl, 20 mM Tris–HCl, pH 8), high-salt washing buffer (0.1% SDS, 1% Triton X-100, 2 mM EDTA, 500 mM NaCl, 20 mM Tris–HCl, pH 8), LNDET buffer (0.25 M LiCl, 1% NP40, 0.1% deoxycholate, 1 mM EDTA,10 mM Tris–HCl, pH 8), and TE buffer (1 mM EDTA, 10 mM Tris–HCl, pH 8). Elution buffer (1% SDS, 0.1 M NaHCO3) was used to retrieve protein-DNA complex. Finally, eluted DNA was purified and quantified by real-time qPCR. The primer sets used were frq C-box F, 5′-GTCAAGCTCGTACCCACATC-3′, and frq C-box R, 5′-CCGAAAGTATCTTGAGCCTCC-3′. Occupancies were normalized by the ratio of ChIP to Input DNA.

The GST-CK1a fusion protein (containing CK1a amino acids Q299-E384) was expressed in BL21 cells, and the recombinant proteins were purified and used as the antigens to generate rabbit polyclonal antisera⁶⁹.

Quantification of western blot and northern blot data were performed using Image J software. All studies were performed at least three independent experiments (different days and different preparations). Error bars are standard deviations for IP assays and standard error of means for race tube assays. Statistical significance was determined by Student’s t-test.

Variables used are:

[F_m] = concentration of frq mRNA

[F] = concentration of FRQ

[C] = concentration of CK1

[F:C] = concentration of complex FRQ:CK1

[F_p:C] = concentration of complex FRQ_p:CK1

[W_a] = concentration of WCC_a

[W_i] = concentration of WCC_i

[F_p:C:W_a] = concentration of complex FRQ_p:CK1:WCC_a

[F_pp:C:W_i] = concentration of complex FRQ_pp:CK1:WCC_i

[F_pp] = concentration of double-phosphorylated FRQ

Parameters used are listed below:

α₁, α₂, α₃, α₄: forward reaction rates

α₅: dissociation rate

β₁, β₂, β₃, β₄: backward reaction rates

α_f: synthesis rate of FRQ protein

α_c: synthesis rate of CK1

α_w: synthesis rate of WWCa

δ₁, δ₂, δ₃, δ₄, δ₅, δ₆, δ₇, δ₈, δ₉, δ₁₀: degradation rates

σ_1,σ₂: transcription rates

σ₃: activation rate of autocatalysis

σ₄: F_pp:C:W_i dependent dissociation of F:C

σ₅: conversion of W_i to W_a

K₁, K₂, K₃, K₄: thresholds

: Hill coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m,n,\ell ,s$$\end{document} m n ℓ s , , ,

System of ordinary differential equations are listed below: 1 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{{\mathrm{d}}[F_m]}}{{{\mathrm{d}}t}} = \sigma _1\frac{{[W_a]^m}}{{K_1^m + [W_a]^m}} + \sigma _2\frac{{[F_p:C:W_a]^n}}{{K_2^n + [F_p:C:W_a]^n}} - \delta _1[F_m]$$\end{document} d [ ] F m d t [ ] W a m K 1 m + [ ] W a m [ ] F p W a : : C n K 2 n + [ ] F p W a : : C n = + − σ 1 σ 2 δ 1 [ ] F m 2 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{{\mathrm{d}}[F]}}{{{\mathrm{d}}t}} \, = \, \alpha _f[F_m] - \alpha _1[F][C] + \beta _1[F:C] - \sigma _3\frac{{[F][C][F_p:C]^\ell }}{{K_3^\ell + [F_p:C]^\ell }} \\ \,+ \sigma _4[F:C]\frac{{[F_{\mathrm{pp}}:C:W_i]^s}}{{K_4^s + [F_{\mathrm{pp}}:C:W_i]^s}} - \delta _2[F]$$\end{document} d [ ] F d t [ ] F [ ] C [ ] F p : C ℓ K 3 ℓ + [ ] F p : C ℓ [ ] F pp W i : : C s K 4 s + [ ] F pp W i : : C s = − + − + − α f α 1 β 1 σ 3 σ 4 δ 2 [ ] F m [ ] F [ ] C [ ] F C : [ ] F C : [ ] F 3 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{{\mathrm{d}}[C]}}{{{\mathrm{d}}t}} \, = \, \alpha _c - \alpha _1[F][C] + \beta _1[F:C] - \sigma _3\frac{{[F][C][F_p:C]^\ell }}{{K_3^\ell + [F_p:C]^\ell }} \\ +\, \alpha _5[F_{pp}:C:W_i] \!+ \sigma _4[F:C]\frac{{[F_{\mathrm{pp}}:C:W_i]^s}}{{K_4^s \!+ [F_{\mathrm{pp}}:C:W_i]^s}} - \delta _3[C]$$\end{document} d [ ] C d t [ ] F [ ] C [ ] F p : C ℓ K 3 ℓ + [ ] F p : C ℓ [ ] F pp W i : : C s K 4 s + [ ] F pp W i : : C s = − + − + + − α c α 1 β 1 σ 3 α 5 σ 4 δ 3 [ ] F [ ] C [ ] F C : [ ] F p p W i : : C [ ] F C : [ ] C 4 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{\mathrm{d}[F:C]}}{{\mathrm{d}t}} \, = \, \alpha _1[F][C] - \beta _1[F:C] + \sigma _3\frac{{[F][C][F_p:C]^\ell }}{{K_3^\ell + [F_p:C]^\ell }} - \sigma _4[F:C]\frac{{[F_{\mathrm{pp}}:C:W_i]^s}}{{K_4^s + [F_{\mathrm{pp}}:C:W_i]^s}} - \alpha _2[F:C] + \beta _2[F_p:C] - \delta _4[F:C]\\ \\ = \, \alpha _1[F][C] + \sigma _3\frac{{[F][C][F_p:C]^\ell }}{{K_3^\ell + [F_p:C]^\ell }} - \sigma _4[F:C]\frac{{[F_{\mathrm{pp}}:C:W_i]^s}}{{K_4^s + [F_{\mathrm{pp}}:C:W_i]^s}} + \beta _2[F_p:C] - (\alpha _2 + \beta _1 + \delta _4)[F:C]$$\end{document} d [ ] F C : d t [ ] F [ ] C [ ] F p : C ℓ K 3 ℓ + [ ] F p : C ℓ [ ] F pp W i : : C s K 4 s + [ ] F pp W i : : C s [ ] F [ ] C [ ] F p : C ℓ K 3 ℓ + [ ] F p : C ℓ [ ] F pp W i : : C s K 4 s + [ ] F pp W i : : C s = − + − − + − = + − + − α 1 β 1 σ 3 σ 4 α 2 β 2 δ 4 α 1 σ 3 σ 4 β 2 [ ] F [ ] C [ ] F C : [ ] F C : [ ] F C : [ ] F p : C [ ] F C : [ ] F [ ] C [ ] F C : [ ] F p : C ( ) α 2 β 1 δ 4 + + [ ] F C : 5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{{\mathrm{d}}[F_p:C]}}{{{\mathrm{d}}t}} = \alpha _2[F:C] - \beta _2[F_p:C] - \alpha _3[F_p:C][W_a] + \beta _3[F_p:C:W_a] - \delta _5[F_p:C]\\ = \alpha _2[F:C] - \alpha _3[F_p:C][W_a] + \beta _3[F_p:C:W_a] - (\beta _2 + \delta _5)[F_p:C]$$\end{document} d [ ] F p : C d t = − − + − = − + − α 2 β 2 α 3 β 3 δ 5 α 2 α 3 β 3 [ ] F C : [ ] F p : C [ ] F p : C [ ] W a [ ] F p W a : : C [ ] F p : C [ ] F C : [ ] F p : C [ ] W a [ ] F p W a : : C ( ) β 2 δ 5 + [ ] F p : C 6 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{{\mathrm{d}}[W_a]}}{{{\mathrm{d}}t}} = \alpha _w - \alpha _3[F_p:C][W_a] + \beta _3[F_p:C:W_a] + \sigma _5[W_i] - \delta _6[W_a]$$\end{document} d [ ] W a d t = − + + − α w α 3 β 3 σ 5 δ 6 [ ] F p : C [ ] W a [ ] F p W a : : C [ ] W i [ ] W a 7 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{{\mathrm{d}}[F_p:C:W_a]}}{{{\mathrm{d}}t}} \, = \alpha _3[F_p:C][W_a] - \beta _3[F_p:C:W_a] - \alpha _4[F_p:C:W_a] + \beta _4[F_{\mathrm{pp}}:C:W_i] - \delta _7[F_p:C:W_a]\\ = \alpha _3[F_p:C][W_a] + \beta _4[F_{\mathrm{pp}}:C:W_i] - (\beta _3 + \alpha _4 + \delta _7)[F_p:C:W_a]$$\end{document} d [ ] F p W a : : C d t = − − + − = + − α 3 β 3 α 4 β 4 δ 7 α 3 β 4 [ ] F p : C [ ] W a [ ] F p W a : : C [ ] F p W a : : C [ ] F pp W i : : C [ ] F p W a : : C [ ] F p : C [ ] W a [ ] F pp W i : : C ( ) β 3 α 4 δ 7 + + [ ] F p W a : : C 8 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{{\mathrm{d}}[F_{\mathrm{pp}}:C:W_i]}}{{{\mathrm{d}}t}} \, = \alpha _4[F_p:C:W_a] - \beta _4[F_{\mathrm{pp}}:C:W_i] - \alpha _5[F_{\mathrm{pp}}:C:W_i] - \delta _8[F_{\mathrm{pp}}:C:W_i]\\ = \alpha _4[F_p:C:W_a] - (\beta _4 + \alpha _5 + \delta _8)[F_{\mathrm{pp}}:C:W_i]$$\end{document} d [ ] F pp W i : : C d t = − − − = − α 4 β 4 α 5 δ 8 α 4 [ ] F p W a : : C [ ] F pp W i : : C [ ] F pp W i : : C [ ] F pp W i : : C [ ] F p W a : : C ( ) β 4 α 5 δ 8 + + [ ] F pp W i : : C 9 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{{\mathrm{d}}[F_{\mathrm{pp}}]}}{{{\mathrm{d}}t}} = \alpha _5[F_{\mathrm{pp}}:C:W_i] - \delta _9[F_{\mathrm{pp}}]$$\end{document} d [ ] F pp d t = − α 5 δ 9 [ ] F pp W i : : C [ ] F pp 10 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{{\mathrm{d}}[W_i]}}{{{\mathrm{d}}t}} = \alpha _5[F_{\mathrm{pp}}:C:W_i] - \sigma _5[W_i] - \delta _{10}[W_i] = \alpha _5[F_{\mathrm{pp}}:C:W_i] - (\sigma _5 + \delta _{10})[W_i]$$\end{document} d [ ] W i d t = − − = − α 5 σ 5 δ 10 α 5 [ ] F pp W i : : C [ ] W i [ ] W i [ ] F pp W i : : C ( ) σ 5 δ 10 + [ ] W i

As a positive control, we measure the frq gene expression using a luciferase reporter fused with the frq promoter in constant dark (DD) conditions in three independent experiments (Supplementary Fig. 4a). The independent experiments result in different frq expression levels, each with slightly increasing trend. To fit our mathematical model assuming a periodic curve of frq, we removed the estimated trend from the original expression level by using a second-order polynomial regression, then standardized the detrended series (Supplementary Fig. 4b).

To estimate the model parameters in the above mathematical model, we adopt a Bayesian analysis with the Markov chain Monte Carlo (MCMC) estimation. In the statistical inference, we view the experimental measurement of frq as a realization of the underlying biological process. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde y_{it}$$\end{document}ỹit denote the detrended and standardized expression level measured at time t in the i-th experiment, i = 1, 2, 3 and t = 56, …, 121 (Supplementary Fig. 4b). The simulation model output for time t is denoted by y_t and its standardized value as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde y_t$$\end{document}ỹt, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde y_t = (y_t - \bar y)/s_y$$\end{document}ỹt=(yt−ȳ)∕sy where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar y$$\end{document}ȳ and s_y are the sample mean and standard deviation of the model outputs {y_t}, respectively.

We put a simple assumption that the deviation of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde y_{it}$$\end{document}ỹit from the simulation output follows the Gaussian random error, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde y_{it} = \tilde y_t + \varepsilon _t$$\end{document}ỹit=ỹt+εt where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon _t \sim N(0,\sigma ^2)$$\end{document}εt~N(0,σ2). To restrict the solution to be a periodic curve, we put the informative prior distribution: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde y_t - \tilde y_{t + 22} \sim N\left( {0,\tau ^2} \right)$$\end{document}ỹt−ỹt+22~N0,τ2 for 56\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\le t \le$$\end{document}≤t≤ 99. We set the tuning parameter (hyperparameter) τ as 0.1, so that the fitted values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{ \tilde y_t\}$$\end{document}{ỹt} do not result in a nonperiodic curve. For other parameters, we put weakly informative priors, so the parameter estimation is mainly driven by the experimental data rather than the prior distribution. Then, we run MCMC to explore the region of the candidate parameter sets whose simulation output \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{ \tilde y_t\}$$\end{document}{ỹt} is similar to the experimental observation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{ \tilde y_{it}\}$$\end{document}{ỹit}, i = 1, 2, 3. After 1000 iterations of burn-in period, we store the MCMC draws up to 4000 iterations, resulting in numerous sets of model parameter estimates. Distribution of parameters is summarized in Supplementary Fig. 5 as boxed plots with a horizontal black line indicating the mean of each parameter. We observe that each parameter has different spread with distinct mean values. For example, α₅ shows a large spread of parameter space ranging from close to 0 to about 11 with a mean value close to 0. This indicates that α₅ is not a critical parameter in the system, because oscillations can be generated even when its value is close to 0 (with a wide range of values of α₅). In contrast, δ₁ and δ₅ show tighter region of parameter space, which suggest that they may play critical roles in generating autonomous oscillations. However, distribution of these parameters does not indicate how the circadian period is altered as a function of each parameter.

Supplementary Fig. 4c shows that the model-fitted frq gene expression levels (solid curves) determined by the randomly selected three different sets of estimated parameters are similar to each other. Supplementary Table 1 shows that the parameter values are distinct which show uncertainty in the parameter estimation. In other words, there exist multiple sets of parameters whose simulation outputs (the frq curves) are identical.

For each parameter set, we performed one-parameter bifurcation analyses to investigate the dynamical changes in the Neurospora circadian clock system such as changes in period and molecular levels. Bifurcation diagrams were obtained by varying a single parameter, while the rest of the parameters were held fixed as their default values; see Fig. 5c–k and Supplementary Fig. 6.

We also evaluated the period sensitivity of each parameter for three parameter sets as follows:\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\it{\phi }} = \frac{{\Delta P}}{{P_0}}$$\end{document}ϕ=ΔPP0where P₀ is the period of default parameter values and ΔP is the difference between the maximum period and minimum period as each parameter varies from 0 to 100% more than the default parameter value. For each parameter set, the ranking of period sensitivity is listed in Supplementary Table 2.

Further information on research design is available in thelinked to this article. Nature Research Reporting Summary

Supplementary Information Peer Review File Reporting Summary