Mechanism for the Generation of Robust Circadian Oscillations through Ultransensitivity and Differential Binding Affinity

Abstract

Biochemical circadian rhythm oscillations play an important role in many signaling mechanisms. In this work, we explore some of the biophysical mechanisms responsible for sustaining robust oscillations by constructing a minimal but analytically tractable model of the circadian oscillations in the KaiABC protein system found in the cyanobacteria S. elongatus. In particular, our minimal model explicitly accounts for two experimentally characterized biophysical features of the KaiABC protein system, namely, a differential binding affinity and an ultrasensitive response. Our analytical work shows how these mechanisms might be crucial for promoting robust oscillations even in suboptimal nutrient conditions. Our analytical and numerical work also identifies mechanisms by which biological clocks can stably maintain a constant time period under a variety of nutrient conditions. Finally, our work also explores the thermodynamic costs associated with the generation of robust sustained oscillations and shows that the net rate of entropy production alone might not be a good figure of merit to asses the quality of oscillations.

Figure 1

KaiC monomer. The following schematic has been inspired from ref (14). The KaiC protein exists as a hexamer, and each monomer consists of 2 domains, CI and CII. The CII domain has two phosphorylation sites, Ser-431 and Thr-432, a KaiA binding site, and a nucleotide binding site (which binds either ATP or ADP). The CI domain binds to KaiB and helps sequester KaiA. Subsequently, the KaiABC complex will be denoted using ^–/A/BCI_TP/DP–^–/ACII_TP/DP^U/T/S/D. Here, TP/DP denotes ATP/ADP attached to the domain: U denotes that none of the sites in CII are phosphorylated, S means that only the serine site is phosphorylated, T means the threonine site is phosphorylated, and D denotes the doubly phosphorylated form. A attached to CI denotes sequestered KaiA; A attached to CII denotes active KaiA acting as an assistant in phosphorylation. B attached to CI implies the inactive form which will start sequestering KaiA.

Figure 2

In penel a, rows are labeled I, II, III, IV, and columns are labeled A, B, C, D. In panel a, the colors in the reaction arrows correspond to those in panel b. Active conformations are denoted using a cyan background, and inactive conformations are denoted using a red background. In our model (panel b), the horizontal axis represents the amount of phosphorylation in the system, with ϕ = 0 and ϕ = 2π corresponding to the completely dephosphorylated state and ϕ = π corresponding to the completely phosphorylated hexamer. The phosphorylation function is a linearly increasing function, 0 at ϕ = 0, 1 at ϕ = π, and then symmetrically decreasing from ϕ = π to 2π. Thus, phosphorylation, . Changes in the phosphorylation levels of the KaiC hexamers give rise to oscillations. KaiA binds to KaiC during the ”day” and promotes phosphorylation, whereas at “night”, KaiB binds to KaiC and sequesters KaiA, thus leading to dephosphorylation. The horizontal rungs in all the states correspond to the phosphotransfer reactions and the hydrolysis of ATP accompanying it, i.e., the red arrows between IA → IIB, and IB → IIC, purple arrows between IIIC → IIIA, and green arrows between IVD → IVA in Figure 2a. The ratio of the forward and backward rates is given by, γ, γ₁, and γ₂ which are all less than 1, because of the fact that these describe reactions coupled to ATP hydrolysis which are inherently irreversible. In the model, α > 1 is responsible for differential affinity, corresponds to % ATP, and k₁ helps in tuning ultrasensitivity. Free KaiA, A_f, provides nonlinearity to the system.

Figure 3

Ultrasensitive response in phosphorylation of KaiC with regard to the total KaiA concentration for K_d0 = 10 and α = 10. The values in the bracket are the Hill coefficients for the response curves (calculated using the method of relative amplification). Values of other parameters are given in Table S2. These kinetics are in the absence of KaiB and P₂ states (ω = ω₁ = 0); i.e., they represent only the active form of KaiC in Figure 2b. Thus, there are no oscillations, and the system always settles into a final steady state.

Figure 4

Comparison between numerical and analytical results for the time-independent solution of P₁ states (eq 3.2) for different α’s. The figure in the inset is a representation of the Markov state network with the P₁ states highlighted. In the main figure, gray corresponds to α = 2, red to α = 4, blue to α = 6, and green to α = 8.

Figure 5

Time period of oscillations for various α and K_d0, i.e., at varying levels of differential affinity and % ATP. k₁ = 0. Other parameters are given in Table S1. Since k₁ = 0, there is no effect of ultrasensitivity. The figure on the left represents time periods calculated by numerically simulating the FPEs. The figure on the right represents the time periods which were calculated from the imaginary part of the maximum positive eigenvalue of the instability matrix W, for small perturbations around the steady-state probability distribution. As can be seen, the analytical solution provides us with a good approximation of the time period as well as the critical α at which oscillations take place for different K_d0 values. The contours in the figure are for the time period of the oscillations.

Figure 6

Comparison between numerical and approximate analytical results for the time-independent solution of P₃ states for the case when k₁ ≠ 0 (eq 3.4). The figure in the inset represents the Markov state network with the P₃ states highlighted. In the main figure, gray corresponds to k₁ = 0, cyan to k₁ = 10^–4, violet to k₁ = 5 × 10^–4, red to k₁ = 10^–3, blue to k₁ = 5 × 10^–3, green to k₁ = 10^–2.

Figure 7

Instability leading to oscillations when changing k₁. The y-axis denotes the maximum eigenvalue of the rate matrix W for the perturbations (refer to Supporting Information). The presence of a positive eigenvalue denotes that the time-independent steady state is unstable. α = 10, and the other parameter values are listed in Table S2.

Figure 8

Value of α required for the onset of oscillations as a function of K_d0. Estimates have been obtained both from our theory and from numerical simulations. We set k₁ = 0 for these calculations.

Figure 9

Value of k₁ required for the onset of oscillations as a function of K_d0. Since k₁ ≠ 0 is only approximately tractable analytically, we have only plotted estimates from numerical simulations.

Figure 10

Time period of oscillations for various K_d0 and k₁ values, i.e., at different levels of % ATP and ultrasensitivity. The white region denotes the parameter space which does not support oscillations. This is also supported by the plot for the amplitude of oscillations, Figure 11. In order to have oscillations at higher values of K_d0, the system requires a higher value of k₁. The contours in the figure are for the time period of oscillations.

Figure 11

Amplitude of oscillations as a function of K_d0 and k₁ at α = 10 and parameters given in Supporting Information, Section S2. The contours in the figure are for the amplitude of oscillations.

Figure 12

Velocity of phosphorylation wavepacket as a function of average angle for k₁ = 0.05, with various K_d0’s and other parameters as given in Table S2. Here, the average angle ⟨ϕ⟩ = ∑_ϕϕP(ϕ), and velocity . The time period of oscillation for the different cycles is denoted along the curves.

Figure 13

Entropy production rate vs α for K_d0 = 5, k₁ = 0, and other parameters given in Table S1. Oscillations start at α = 21. α = 1 corresponds to the absence of differential affinity. In order to have oscillations, an additional 0.113 units of energy are required. This energy goes into building coherence among the KaiABC oscillator population

Figure 14

Entropy production rate vs k₁ for α = 10, K_d0 = 8 and other parameters given in Table S2. Unlike the case with changing α in Figure 13 where the entropy production plateaus very quickly with increasing α, in this case, the entropy production increases almost linearly with increasing k₁. As expected, decreasing K_d0 and increasing k₁ lead to a higher dissipation of energy. Oscillations start at k₁ = 0.03. k₁ = 0 corresponds to the absence of ultrasensitivity in the system. An additional 0.052 units of energy are dissipated in order to have oscillations. This additional energy goes into improving the ultrasensitive response of the system, eventually leading to coherence.