A Circadian Clock-Regulated Toggle Switch Explains AtGRP7 and AtGRP8 Oscillations in Arabidopsis thaliana

The circadian clock controls many physiological processes in higher plants and causes a large fraction of the genome to be expressed with a 24h rhythm. The transcripts encoding the RNA-binding proteins AtGRP7 (Arabidopsis thaliana Glycine Rich Protein 7) and AtGRP8 oscillate with evening peaks. The circadian clock components CCA1 and LHY negatively affect AtGRP7 expression at the level of transcription. AtGRP7 and AtGRP8, in turn, negatively auto-regulate and reciprocally cross-regulate post-transcriptionally: high protein levels promote the generation of an alternative splice form that is rapidly degraded. This clock-regulated feedback loop has been proposed to act as a molecular slave oscillator in clock output. While mathematical models describing the circadian core oscillator in Arabidopsis thaliana were introduced recently, we propose here the first model of a circadian slave oscillator. We define the slave oscillator in terms of ordinary differential equations and identify the model's parameters by an optimization procedure based on experimental results. The model successfully reproduces the pertinent experimental findings such as waveforms, phases, and half-lives of the time-dependent concentrations. Furthermore, we obtain insights into possible mechanisms underlying the observed experimental dynamics: the negative auto-regulation and reciprocal cross-regulation via alternative splicing could be responsible for the sharply peaking waveforms of the AtGRP7 and AtGRP8 mRNA. Moreover, our results suggest that the AtGRP8 transcript oscillations are subordinated to those of AtGRP7 due to a higher impact of AtGRP7 protein on alternative splicing of its own and of the AtGRP8 pre-mRNA compared to the impact of AtGRP8 protein. Importantly, a bifurcation analysis provides theoretical evidence that the slave oscillator could be a toggle switch, arising from the reciprocal cross-regulation at the post-transcriptional level. In view of this, transcriptional repression of AtGRP7 and AtGRP8 by LHY and CCA1 induces oscillations of the toggle switch, leading to the observed high-amplitude oscillations of AtGRP7 mRNA.

f T LD & f T LL -Synchronized Oscillations 1: AtGRP7 and AtGRP8 mRNA adopt circadian oscillations in all diurnal data sets studied for wild type plants. We therefore require that they adopt the same period as the core oscillator both under LD and LL conditions.
The corresponding cost function term for the LD condition is defined as contributing a value smaller than unity if the deviation ∆T LD,i := T Core LD − T Slave LD,i of the slave oscillator's averaged period T Slave LD,i from the averaged core oscillator's period T Core LD is smaller than three minutes, i.e. δT LD := 0.05h. The index i = 7 accounts for the period of AtGRP7 mRNA and i = 8 for that of AtGRP8 mRNA. The core oscillator period T Core LD is determined by the time difference between two successive maxima of the LHY/CCA1 protein oscillations P L (t). Likewise we define the periods T Slave LD,i of AtGRP7 and AtGRP8 mRNA oscillations as the time difference between two consecutive maxima of the AtGRP7 and AtGRP8 mRNA concentrations with δT LL = 0.05h.
f LC -Synchronized Oscillations 2: Our circadian core oscillator model gives rise to periodic oscillations of P L (t) under 12h:12h LD as well as LL conditions [7]. We therefore require that the slave oscillator should also periodically oscillate under both circumstances. Similarly as in [1][2][3] we thus define to quantify the entrainment of the slave oscillator by the core oscillator under LD and LL conditions. The first term gives a contribution smaller than unity if the relative variability around its mean value in LD of the concentration maxima K max i of a given oscillating variable K i is smaller than 5%. Likewise, the second term accounts for the LL conditions and the third term for the variability of the phases φ K i (see main text, section "Comparison with Experimental Results and Predictions") among successive LD-cycles.
In summary the term (4) thus penalizes significant variations in phase or height of the oscillation maxima (e.g. chaotic solutions or damped oscillations) while the terms (2) and (3) impose the periodicity of the core oscillator to the slave. f φ -Experimentally observed Phases: The model of our core oscillator was optimized for 12h:12h LD experimental entrainment conditions [7]. We therefore used the data sets from the DIURNAL data bank experimentally obtained under those 12h:12h LD conditions as well. The phase of the AtGRP7 mRNA oscillations was estimated as φ exp M 7 = zt 10h [8]. The AtGRP8 mRNA oscillations precede those of AtGRP7 by two hours [8], similar to what was found in studies for other entrainment conditions [8,9], so we set φ exp M 8 = zt 8h. There is, to the best of our knowledge, no published data on AtGRP7 or AtGRP8 protein concentrations for 12h:12h LD entrainment conditions. In [10] one can find data on protein concentrations under 8h:16h LD entrainment conditions which revealed that AtGRP7 protein peaks approximately four hours after its mRNA. We used this as a gross approximation also under 12h:12h LD conditions and set φ exp , respectively. f φ then contains four summands and reads as where the different reliability of the data is reflected by our error tolerances δφ M 7 = δφ M 8 = 0.5h and δφ P 7 = δφ P 8 = 1.0h, respectively.
f A -Amplitude: In view of [8] we expect that the mRNAs and proteins of AtGRP7 and AtGRP8 exhibit concentrations (number of molecules per volume within the pertinent cell compartment) which are very roughly speaking of the same order of magnitude as the concentrations M L (t) and P L (t) of LHY/CCA1. Therefore we defined the summand f A related to the amplitude A of the oscillations as where g(A) is depicted in Figure S10. This function g(A) has been obtained according to the following considerations: Starting with the Ansatz the parameters c min and c max were chosen so as to generate g(A) values smaller than one for all peak-trough-values 2A within the heuristically predetermined range between 0.3 and 2.0, resulting in c min ≈ 0.30, c max ≈ 2.01.
f ∆ & f ∇ -Waveform: Next we turn to the waveform of the experimentally observed AtGRP7 and AtGRP8 mRNA as well as protein oscillations. Experimental time traces of AtGRP7 and AtGRP8 mRNA show non-sinusoidal behavior [8,10]. After being at its trough value, AtGRP7 and AtGRP8 mRNA quickly recovers, typically within four hours, to its maximal value, and needs a slightly longer time to reach its trough value again, creating nevertheless sharp mRNA peaks. Experimentally observed AtGRP7 protein time traces [10], anti-phasic to its mRNA, are much broader and stay at their trough value much shorter. Similarly as in [1][2][3] we thus define where M max+τ M i denotes the M i concentration τ M hours before or after its concentration maximum and likewise for P max+τ P

7
. Generally speaking f ∆ thus penalizes solutions that do not show the experimentally observed sharp mRNA peak and broad protein peak, respectively. More precisely the first term on the right hand side of equation (8) contributes a value smaller than unity for solutions of the M 7 (t) and M 8 (t) oscillations in LD conditions that drop at least by 2/3 of their maximum value within 4 hours before and 6 hours after the respective peak time, and larger than unity otherwise. The second term works in the same manner for LL conditions. Each summand of the third term similarly contributes a value smaller than unity for solutions of the P 7 (t) oscillations in LL conditions, that drop at least by 2/3 of their maximum value within 12 hours before and 8 hours after the respective peak time, consistent with the broad protein oscillations observed in [10].
The term f ∇ in (1) accounts for the trough-broadness of the AtGRP7 and AtGRP8 mRNA oscillations in LL conditions and is defined as Here ∇ 7 is defined as the time the AtGRP7 mRNA concentration exceeds its minimal value by less than 10 per cent of the peak-trough-difference and similarly for ∇ 8 . The corresponding experimental quantities ∇ exp i are quite noisy, ranging from 8 [8] to 12 [10] hours. We therefore set ∇ exp i = 10h and δ∇ = 2h.
f t 1/2 -Half-Life: In experiments measuring the AtGRP7 mRNA half-life [11,12], plants were grown under 16h:8h LD long-day entrainment conditions and were transferred to a medium supplemented with 150 µg ml −1 cordycepin, two hours before expecting the mRNA peak, in order to suppress transcription. Afterwards, the mRNA abundance was measured with one-hour time intervals, observing that the mRNA is reduced to 50% of its original amount within 3-4 hours. We therefore set t exp 1/2 = 3.5h with an error tolerance of δt 1/2 = 0.5h. We followed the experimental protocol in our simulations exactly. First, we inhibited the transcription of AtGRP7 and AtGRP8 in the model by setting v 7 and v 8 in equations (1) and (4) of the main text to zero. Since the production of AtGRP7 mRNA depends on the normal splicing of its pre-mRNA and since its degradation kinetics are of Michaelis-Menten type, its half-life will depend on the systems initial conditions. This initial concentrations of the slave oscillator were set, in analogy to experiments, to the values two hours before the maximum of the AtGRP7 mRNA oscillation is expected within the last entrainment cycle, i.e. for t ∈ [312h, 336h]. Given these initial conditions, we solved the resulting equations for a maximal time of τ max = 16h and define the half-life t 1/2 as the time a given concentration needs to decline to the half of its initial value (see Figure S3 A/B). The summand f t 1/2 then reads as In case the system did not decay to half of the initial AtGRP7 mRNA concentration within τ max = 16h, we set t 1/2 = τ max , therefore defining a maximal contribution to the cost function f max t 1/2 = 16h−3.5h

0.5h
As there is no experimental data available on the half-life of AtGRP8 mRNA, the cost-function lacks a corresponding term. In [12] the half-life of the alternatively spliced mRNA of AtGRP7 was estimated as 0.5h. Since this alternatively spliced mRNA was not assumed to be able of producing functional protein, it will therefore not feed back to the system described by equations (1)-(6) in the main text and a corresponding cost function-term can therefore be neglected.

C: Search for Self-Sustained Oscillations
We tried to find self-sustained autonomous oscillation in our decoupled slave oscillator system, i.e. P L (t) = 0 in equations (1)-(6) of the main text, by modifying the kinetic parameters of Table 1 without changing our model itself. To this end we searched for autonomous oscillations in the neighborhood of the parameter set in Table 1 by means of a gradient descent method that drives the real part of the complex conjugated pair of eigenvalues of the dynamical fixed point across the imaginary axis, corresponding to a Hopf bifurcation of the fixed point into a periodic solution. Remarkably enough, this algorithm indeed generated self-sustaining oscillations (see Figure 7 B of the main text) by only changing the protein degradation rates m 7,2 and m 8,2 , while leaving all other kinetic parameters unchanged. The period of these oscillations is mainly determined by the protein degradation rate m 7,2 of AtGRP7, exhibiting a decreasing period with increasing m 7,2 , but still remaining below the experimentally observed circadian periodicity. In doing so, we restricted our search to a sufficiently small neighborhood of the reference parameter set from Table 1 in order to maintain a reasonably good agreement with the experimental facts 1-5 from section Parameter Estimation of the main text. However, even within this neighborhood all parameter sets exhibiting oscillations actually disagreed already quite notably from the experimental facts 1-5 (see section Parameter Estimation of the main text).

D: Search for Noise-Induced Oscillations
Following [13,14] we model concentration fluctuations by adding noise in equations (1)-(6) of the main text, i.e. equation (7) of the main text now takes the forṁ x k (t) = g k (t, x(t), p) + σx k (t)ζ k (t), where σ quantifies the noise strengths and ζ k (t) are independent, delta-correlated Gaussian white noises of zero mean. The noisy dynamics (18) were tackled by using XPPAUT (Version 6.11). Figure 7 C demonstrates the possibility of noise-induced self-sustained oscillations for suitable noise strengths σ. For the specific σ value from Figure 7 C a spectral analysis reveals a main peak around T ≈ 21.4 h, a period close to the experimentally observed circadian rhythmicity.