Authors: Haochen Fu, Chenyi Fei, Qi Ouyang, Yuhai Tu
Categories: Biological Sciences, Physical Sciences, biochemical oscillators, circadian clocks, kinetic regulation, temperature compensation, 435
Source: Proceedings of the National Academy of Sciences of the United States of America
Although individual kinetic rates in biochemical reactions are sensitive to temperature, most circadian clocks exhibit a relatively constant period across a wide range of temperatures, a phenomenon called temperature compensation (TC). However, it remains unclear how different biochemical oscillators achieve TC. In this study, using representative biochemical oscillator models with different underlying reaction networks, we demonstrate a general kinetic regulation mechanism for TC regardless of the network structure. We find that by driving the system into a regime far from onset where the period increases strongly with at least one of the kinetic rates in the system to balance its inverse dependence on other rates, robust TC can be achieved for a wide range of parameters in different networks.
Keywords: temperature compensation, biochemical oscillators, circadian clocks, kinetic regulation
Biological systems are subject to temperature changes in their environments. Most biochemical reactions are dependent on temperature T, with their kinetic rate constants k obeying the Arrhenius law (1, 2):
where A is a temperature-independent prefactor, kB is the Boltzmann constant, and E is the activation energy. The temperature sensitivity of a chemical reaction can be characterized by its Q10 factor, which measures the change in the reaction rate when the temperature is increased by 10 K. For typical biochemical reactions, the activation energy is E∼20kBT0 (3), leading to a Q10≈2 under room temperature T0≈300 K—the rate of a biochemical reaction approximately doubles upon a 10 K increase in temperature (around room temperature). This empirical law is also known as Van’t Hoff’s rule from physical chemistry. Such temperature dependence occurs in various biological processes, such as bacterial cell growth (4) and vertebrate muscle contraction (5).
Given the large Q10 values for individual reaction rates, it is remarkable that nearly all circadian clocks can maintain a relatively constant period over a wide range of temperatures, a phenomenon known as temperature compensation (TC) (6–9). Specifically, the Q10 value of the period of a circadian clock is between 0.9 and 1.1 (10); namely, a 10 K increase in temperature only results in a 10% change in the period of the circadian clock, which is much smaller than the typical twofold changes for individual reaction rates. While some previous studies attributed TC to certain genes or proteins whose activities are insensitive to temperature (11–13), biological clocks are not simply insensitive to temperature; for example, it is well known that the clock can be readily entrained by external periodic temperature signals (14–16).
In general, dynamics of biochemical oscillations are governed by biochemical reaction networks with multiple biomolecules interacting through different biochemical reactions. As a result, the period P(k) of the oscillation should be a function of all the reaction rates k={ki} with i(=1,2,...) labeling the individual reactions. Here, we define the “period sensitivities”
to describe the dependence of the period P on the i-th reaction rate ki—reactions are referred to as period-lengthening (Ci>0) or period-shortening reactions (Ci<0), respectively.
From a systems perspective, while each individual reaction rate can be sensitive to temperature with a high Q10, TC could in principle arise from a balance between period-lengthening reactions and period-shortening reactions (17–21). To see this, using Eq. 1, we can approximate Q10 of the period as
where the summation is over all the reactions, Ei is the activation energy of the i-th reaction in the unit of kBT0. From now on, we will refer to the quantity “Q10” as Q10 of the period defined by the first line of Eq. 3. Perfect TC (Q10=1) is achieved by having ∑iCiEi=0, which we will refer to as the “balancing scenario” hereafter. Since the activation energies are positive (Ei>0), perfect TC requires the cancelation of the weighted contributions CiEi from the period-lengthening reactions (Ci>0) and the period-shortening reactions (Ci<0).
Given a reaction network with certain Ci’s, one can fine-tune Ei’s to minimize ∑iCiEi, which may be achieved by natural selection (22). However, it remains an interesting open question whether there are design principles that facilitate the realization of the balancing scenario. Some previous studies tackled the problem from the perspective of the reaction network structure (20, 23–28). For example, a network with positive feedback is likely to have period-lengthening reactions (Ci>0) (23, 29), and a network with robust adaption can result in many zero Ci’s that reduces the effective dimension of the problem (26, 27). Nevertheless, even for a preferred reaction network for TC, the performance of TC still depends on kinetic rates in the network (23), suggesting additional requirements to achieve TC in these systems. On the other hand, it is typically easier for organisms to achieve a biological function, such as TC, by adjusting specific kinetic rate parameters within an existing network structure rather than modifying the entire network, because the former only needs small modifications in molecular structures, whereas the latter typically requires additional molecular components in the system. Therefore, it is important to explore strategies for controlling (regulating) kinetic rate parameters in order to achieve TC.
Indeed, kinetic regulation, by which we mean controlling kinetic rates in a temperature-independent way, in general and in particular, changing the prefactors in Eq. 1 by varying certain molecule concentrations in the system, plays an important role in biochemical oscillators. The nonequilibrium nature of limit-cycle oscillation dynamics dictates that the kinetic rates in the underlying biochemical network break detailed balance (30). Furthermore, the kinetic rates must be increased beyond the onset of oscillation, typically a Hopf bifurcation for limit-cycle oscillations (31), in order to generate the oscillatory behaviors. However, these basic requirements do not necessarily result in temperature-compensated oscillations. Thus, conditions for TC remain elusive.
In this paper, we start our investigation with two simplest nonlinear oscillators—the Van der Pol model and the Brusselator model—representing two fundamental motifs in two-component oscillatory the activator–inhibitor motif and the substrate-depletion motif, respectively (32). We find that kinetic regulation provides a general mechanism wherein increasing certain period-lengthening reaction rate(s), which is possible only far from the onset, can lead to robust TC within a realistic range of activation energies. The energy cost associated with this mechanism is also studied. Finally, we verify the applicability of this general TC mechanism in other oscillatory systems including two realistic biological clocks.
We start by considering one of the simplest and the most generic nonlinear oscillators—the Van der Pol (VdP) oscillator (33), described by
By analyzing its Jacobian matrix near the Hopf bifurcation point, we notice that its network topology belongs to the activator–inhibitor motif (23, 32, 34) (Fig. 1A, see also SI Appendix, section 1.A). To study the performance of TC in the VdP model, we assume that the rate parameters Ω and μ depend on the temperature through the Arrhenius law Eq. 1 with activation energies EΩ and Eμ, respectively. Unless otherwise noted, we study effects of kinetic regulation for TC by varying the temperature-independent prefactors (Aμ, AΩ, etc.) in Eq. 1.
Fig. 1. The Van der Pol (VdP) model can achieve TC by kinetic regulation. (A) The topology of the VdP model is an activator–inhibitor. (B) The period at the room temperature (T0=298 K) and that at 10 K above vary with the “kinetic rate” μ by varying the prefactor of μ; Q10 value is calculated by the ratio of the period at such two temperatures. Ω is fixed at Ω=1 at T=T0. The activation energies are set at Eμ=26.7 and EΩ=13.3 (Eμ≈2EΩ) in the unit of kBT0. (C) The Q10 value in the activation-energy space when μ=1∼Ω and when μ=40≫Ω (both μ and Ω are the values at T=T0). The TC domain 0.9≤Q10≤1.1 is highlighted in red. (D) The TC domain (%), defined by the area of the TC domain divided by the total area of the activation-energy space in (C), varies in terms of μ at T=T0.
The VdP model (Eq. 4) undergoes a Hopf bifurcation at μ=0, with a stable limit cycle for μ>0. For simplicity, we keep Ω(T0)=1 and regulate μ(T0) (unless otherwise noted, all μ values are at T=T0). As shown in Fig. 1B, the period of the limit-cycle oscillation is sensitive to temperature change when μ is small or comparable to Ω. However, if μ becomes much larger than Ω, TC of the period can be achieved, reflected clearly by the Q10 value of the period approaching 1 (Fig. 1B). Such a realization of TC can also be under the constraint of constant period at T=T0 by varying both μ and Ω (SI Appendix, Fig. S1).
The good TC performance for large μ/Ω as shown in Fig. 1B is achieved for a particular choice of activation energies. The natural question is how robust is TC with respect to the choices of the activation energy parameters. To address this question, we calculate Q10 in the activation-energy space of Eμ-EΩ. We limit the range of the activation energies to be between 10 and 30 kBT0, which are relevant for typical biochemical reactions (3) and also sufficiently large so that the reactions are not trivially temperature-insensitive. On each Eμ-EΩ phase diagram, μ and Ω at T=T0 are controlled, and the prefactors of μ and Ω change accordingly with the activation energies. When μ is small or comparable to Ω, Q10 is globally small in the activation-energy space, and there is no regime that satisfies the good TC requirement 0.9≤Q10≤1.1, as shown in Fig. 1C. However, when μ≫Ω, there exists a band domain where 0.9≤Q10≤1.1, named as TC domain, suggesting that TC can be robustly achieved with a variety of combinations of activation energies (Fig. 1C). Indeed, the larger the TC domain, the more robustly TC can be achieved. Thus, the size of the TC domain serves as the most important measure of TC performance for a given kinetic regulation scheme. By this measure, kinetic regulation by increasing μ improves the TC performance in the VdP model, as shown in Fig. 1D.
In biochemical reaction systems that exhibit oscillatory behaviors, if all the reaction rates are changed by the same factor λ, then the period should change by a factor λ−1. From this time scaling invariance, one can obtain an exact sum rule among period sensitivities of all the reactions in the system (21, 35):
which serves as a global constraint for all the period sensitivities. Indeed, the period sensitivities of the VdP model always satisfy Cμ+CΩ=−1 (Fig. 2A). If we define the period sensitivity vector C whose components are the period sensitivities of each reaction, and correspondingly the activation energy vector E, the sum rule Eq. 5 constrains C on a subplane, whereas the condition for perfect TC requires that E is perpendicular to C, i.e.,
Fig. 2. The TC mechanism in the VdP model in the balancing scenario. (A) The period sensitivities vary with μ (Ω=1). (B) The period sensitivity vector C is constrained on the ∑iCi=−1 subplane, whereas the activation energy vector E is constrained in the activation-energy space. Perfect TC requires C·E=0. (B1) When all Ci’s are negative or zero (corresponding to small μ/Ω), the TC domain does not exist in the activation-energy space. (B2) TC domain exists when there exists a large positive Ci (corresponding to large μ/Ω). (C) The behavior of the amplitude and the progression speed on the limit cycle (C1) near the onset (μ/Ω=0.1) and (C2) far from the onset (μ/Ω=100). The color bar shows the inverse of the progression speed, i.e., time per distance, for both panels.
Given that all the activation energies are positive, TC cannot be achieved in a biochemical network where all components of C are negative (Fig. 2B1), which is certainly allowed by the sum rule. Thus, a necessary condition for TC is the existence of period-lengthening reactions with positive period sensitivities, which is a nontrivial requirement in light of the sum rule (Eq. 5). Quantitatively, if the positive components of C are small, the corresponding components of vector E must be extremely large compared to other components, and E is restricted on the margin of the nontrivial activation-energy space. On the other hand, a large positive component of C can move E toward the center of the activation-energy space, allowing TC to be achieved by more combinations of activation energies that are physiologically plausible (red-shaded regime in Fig. 2B2). Therefore, TC can be achieved more robustly in the presence of large positive period sensitivities.
The general analysis based on the period sensitivity sum rule can be used to understand why TC is systematically improved by regulation of “kinetic rates” μ and Ω. As shown in Fig. 2A, when μ≲Ω, Cμ≈0 and CΩ≈−1 (see also SI Appendix, section 1.A for theoretical derivation). Hence, the negative CΩEΩ cannot be balanced by a nearly zero CμEμ, leading to a Q10 substantially smaller than 1 (Fig. 1B). As μ increases to μ≫Ω, Cμ increases to 1 while CΩ decreases to −2 (SI Appendix, section 1.A). In this case, the balance between CΩEΩ and CμEμ becomes possible. Especially, to achieve perfect TC, we only need Eμ=2EΩ, so that CΩEΩ+CμEμ=0, as is the case in Fig. 1B.
How does a large positive period sensitivity emerge from kinetic regulation? In a biochemical oscillator, kinetic rate constants have two 1) controlling the progression speed in the phase space of the oscillation and 2) shaping the limit cycle in particular determining the amplitude of the oscillation. Intuitively, increasing a kinetic rate will increase the progression speed on the limit cycle in the phase space. If the limit-cycle trajectory is unchanged, increasing the kinetic rates should generally lead to a faster period, corresponding to negative period sensitivities, as indicated by Eq. 5. To obtain a positive period sensitivity, such a kinetic rate must change the shape of the limit cycle; more specifically, it must increase the amplitude of the oscillation, so that the increased distance in the phase space will compensate for the increased progression speed (36–38). However, this does not necessarily lead to a large positive period sensitivity. For example, in the VdP model with small μ/Ω, the system is close to onset and the oscillation is highly harmonic with a near-circular limit cycle as shown in Fig. 2C1. According to the normal form of Hopf bifurcation (31), the radius of the circular limit cycle scales as r∼μ1/2, while the near-uniform progression speed of the limit cycle scales as v≈Ωr∼μ1/2 (Fig. 2C1). Therefore, the contributions of μ to the amplitude and progression speed of the limit cycle oscillation cancel out, resulting in an overall small Cμ≈0. In summary, a kinetic rate that positively correlates with the oscillation amplitude is necessary but insufficient to have a large positive period sensitivity.
How does Cμ become large (close to 1) in the μ≫Ω regime for the VdP oscillator? To address this question, we examine the limit-cycle trajectory when μ≫Ω. As shown in Fig. 2C2, the limit cycle is separated into fast and slow phases (the inverse progression speed, i.e., the time per distance traveled, is shown by the color along the trajectory), and the period is mostly determined by the time spent on the slow phases along the two segments of the X-nullcline. The amplitude of the slow phases along the X-nullcline, denoted by ΔYslow, is proportional to μ/Ω (see SI Appendix, section 1.A for derivation). Surprisingly, unlike the situation near the onset, the progression speed in the slow phases is given by dY/dt=ΩX∼Ω, insensitive to the change of μ. Indeed, increasing μ accelerates the progression speed linearly only in the fast phases where the time is negligible (SI Appendix, section 1.A). Consequently, the period of the oscillation P≈2ΔYslow/|dY/dt|∼μ/Ω2 increases linearly with μ, leading to a large positive Cμ.
Overall, the separation of the fast and slow phases along the limit-cycle trajectory in the parameter regime far beyond the oscillation onset, a phenomenon which we call oscillation phase separation (OPS), plays an important role in creating a large positive period sensitivity that is crucial for TC. In the VdP model, the emergence of the slow phase(s) with a μ-insensitive progression speed and a μ-increasing amplitude is the origin of a period-lengthening rate with a large positive period sensitivity. In the rest of the paper, we show that this general mechanism for TC based on OPS applies to other more realistic biochemical oscillatory systems.
While the VdP model studied in previous sections represents the activator–inhibitor motif, in this section, we study TC in another class of oscillators with the substrate-depletion motif (23, 32, 34) as represented by the reversible Brusselator model (39–41). This classical model describes the nonlinear autocatalytic reactions (Fig. 3A)
Fig. 3. TC can be achieved by kinetic regulation in Brusselator. (A) The reversible Brusselator model. The concentrations of A, B, and D are kept constant over time. When the two reverse reaction rates, k−2=k−2,0[D] and k−3, are negligible, the Brusselator becomes its original irreversible version. (B) In the irreversible Brusselator, the period at two temperatures (298 K and 308 K) and correspondingly Q10 value of the period vary with the kinetic rate k2. The activation energies are set as E1=15, E2=25, and E3=E−1=20 (kBT0). (C) In the irreversible Brusselator, the enhancement of TC is caused by the emergence of large positive period sensitivity C2, which increases from nearly 0 to 2 with the increase of k2. (D) The limit cycle of the irreversible Brusselator far from the onset with the color bar showing the inverse of the progression speed. The large positive period sensitivity of k2 is a result of a k2-sensitive amplitude and a k2-insensitive progression speed of the slow phase that dominates the period.
where the forward and reverse kinetic rates depend on temperature via Eq. 1, and A, B, and D have constant concentrations [A], [B], and [D], respectively.
For simplicity, we first neglect the reverse kinetic rates k−2,0 and k−3, and the model becomes the original irreversible Brusselator (18, 39, 40). Sustained oscillation occurs when the kinetic rate k2=k2,0[B] exceeds the Hopf bifurcation critical point k2c, beyond which the amplitude of limit-cycle oscillation further increases with k2. Remarkably, increasing k2 can also enhance the performance of TC (Fig. 3B), and the good TC performance at large k2 is robust as it can be realized in a wide range of activation energies (see SI Appendix, Fig. S3 for details). In the Brusselator model, the period sensitivity C2 for reaction k2 is positive and it increases from nearly 0 to around 2 with increasing k2 (Fig. 3C), which leads to better TC performance at larger values of k2.
To understand the positive period sensitivity of k2 when k2≫k2c, we analyze the limit-cycle trajectory in the X-Z phase space where Z=X+Y denotes the total concentration of X and Y. Similar to the VdP model, the limit-cycle trajectory also exhibits the OPS behavior at large values of k2 with a slow phase and three relatively fast phases as shown in Fig. 3D. During the slow phase, X is kept at a roughly constant but low level X≪k1,0[A]/k−1, resulting in a nearly constant progression speed vslow≈dZ/dt=k1,0[A]−k−1X≈k1,0[A], independent of k2 (see SI Appendix, section 1.B for details). On the other hand, the amplitude of the slow branch ΔZslow, and hence the time on the slow phase τslow≈ΔZslow/vslow, asymptotically increases with k2 as k22 (see SI Appendix, section 1.B for detailed derivations). Indeed, the time in the other phases negatively depends on k2, but this period-shortening effect is negligible as the times spent in these fast phases are much shorter than that in the slowest phase (see SI Appendix, section 1.B for detailed derivations). As a result, the period can be approximated by τslow that scales with k22, leading to the period sensitivity C2≈2 when k2≫k2c.
In a biochemical system, kinetic regulation is often achieved by regulating rate-associated molecular concentrations. For example, in the Brusselator model, regulation of k2 can be easily implemented by changing the concentration [B]. However, it consumes free energy to maintain a molecular concentration at a nonequilibrium level (30). Thus, kinetic regulations of reaction rates incur an energy cost. To study the energy cost of TC enhancement in Brusselator, we consider the full reaction network with finite reverse reaction rates, k−2=k−2,0[D] and k−3. The chemical potential difference between B and D (in units of kBT0) reads
which characterizes the irreversibility of the system (30, 41)—when ΔμDB=0, the system is at equilibrium, and sustained oscillation cannot be achieved. Increasing ΔμDB drives the system away from equilibrium and enables oscillation, which can be achieved by either increasing [B] or decreasing [D]. To quantify the energy cost, we compute the free-energy dissipation per oscillation period (30, 41), ΔW, given by
where Ji+(t) and Ji−(t) denote, respectively, the forward and the reverse fluxes of the i-th reversible reaction pair at time t.
For the reversible Brusselator model, we compute the Q10 of the period, the period sensitivity C2 of k2, and the energy dissipation ΔW for varying reaction rates k2 and k−2, as shown in Fig. 4A–C. Remarkably, either increasing [B] or decreasing [D] can increase the positive period sensitivity C2, and consequently enhance TC (Fig. 4 A and C). A perfect TC is achieved in the limit of extremely large [B] and extremely small [D] where the system is far from equilibrium. On the other hand, the free-energy dissipation per cycle (ΔW) is also increased by [B] and decreased by [D] (Fig. 4B). If we plot Q10 against ΔW (Fig. 4D) for different choices of [B] and [D], we find that the upper bound of Q10 is increased by the energy cost. Our results indicate that kinetic regulation for TC costs energy and biochemical oscillators can exploit free-energy dissipation to enhance TC performance by kinetic regulation.
Fig. 4. In the reversible Brusselator, kinetic regulation consumes free energy to enhance the TC performance. (A–C) Q10, energy dissipation ΔW, and period sensitivity C2 in the parameter space spanned by k2=k2,0[B] and k−2=k−2,0[D] (log scale). The dark blue regime does not have sustained oscillation. Activation energies are set as E1=15, E2=25, and E3=E−1=E−2=E−3=20 (kBT0). (D) Scatter plots of Q10 against the energy cost lnΔW according to A and B, with the color showing the period sensitivity C2 in C. The dashed curve shows the upper bound of the scatter plot.
To study whether the above mechanism of kinetic regulation can be applied to understanding TC in realistic biological systems, we focus on two well-known systems that exhibit circadian the Drosophila circadian clock and the Kai system in cyanobacteria circadian clocks. The former represents the family of transcription–translation circadian clocks, while the latter stands for the posttranscriptional circadian clocks. Both systems have been reported to have TC (6, 42–46).
We first investigate the model proposed by Tyson (32, 47) for the Drosophila circadian clock, which has been reduced to a two-species model consisting of only the per mRNA and total PER proteins (SI Appendix, section 1.C). The architecture of the Tyson model contains a substrate-depletion motif (SI Appendix, section 1.C), the same as the Brusselator model studied in the previous subsections. Similar to the Brusselator model, we show that increasing one of the kinetic rates, specifically the phosphorylation rate of PER monomer k1, drives the system away from the onset and the period sensitivity of k1 has a large positive value at large k1 (see SI Appendix, section 1.C and Fig. S4 for details). As a result, TC can be achieved through the same general OPS mechanism by kinetic regulation of k1.
To confirm that the OPS mechanism of TC also works in more complex high-dimensional systems, we next focus on the van Zon–Hatakeyama (vZH) model for the Kai system (48, 49). In this model, the timekeeper protein KaiC forms a hexamer that can switch between active and inactive conformations (states) (Fig. 5A). Active KaiC has seven phosphorylation states corresponding to the number of phosphorylated monomers from 0 to 6. Each phosphorylation step requires an active KaiC to bind to the enzyme KaiA, with the binding affinity decreasing as the phosphorylation level of KaiC increases, which is crucial for synchronization of individual KaiC hexamers. On the other hand, inactive KaiC can spontaneously dephosphorylate. The transitions between active and inactive states occur only in the fully phosphorylated state (between C6 and C~6) or in the fully dephosphorylated state (between C˜0 and C0).
Fig. 5. TC is achieved by regulation of total KaiA concentration in a model of the Kai system. (A) The van Zon–Hatakeyama (vZH) model. KaiC is a hexamer that can be in the active form or the inactive form. During the active states, KaiC can be phosphorylated when binding to KaiA. The binding affinity of KaiA and KaiC decreases with the phosphorylation level of KaiC (48). Based on the oscillation dynamics, different states of KaiC can be grouped into three X=[C5]+[AC5], Y=∑i=04[Ci]+[ACi], and Z=[KaiC]total−X−Y, such that kp6 of our interest is explicit in the simplified network of X, Y, and Z. (B) Q10 of the period varies with the total KaiA concentration. The dots are experimental data from ref. 50, and the curve is from the simulation of the model with parameters tuned to fit the data (see SI Appendix, section 1.D for details). Activation energies are set to Ep6=26, and the other Ei=6 (kBT0). (C) At an intermediate level of total [KaiA] (ln[KaiA]=0.4), the limit-cycle trajectory projected on the phase subspace of X, Y, and Z shows OPS: S, S^′^, and F. Each phase shows the transition between two variables while the third variable is roughly constant. (D) The relative changes of the amplitude (length of the trajectory), the mean progression speed, and the time spent on each of the three phases are calculated when kp6 is perturbed by ±20%. (E) The period sensitivities of all reaction rates vary with the total [KaiA]. Among them, the period sensitivity of kp6 (the kinetic rate of the sixth phosphorylation step) has a relatively high positive sensitivity in the middle range of total [KaiA].
Previously, Hatakeyama and Kaneko investigated the TC property in this model (48). Their results suggested that achieving TC relies on the competition for limited KaiA among various active KaiC phosphorylation states. However, their results also showed that this limited-KaiA scheme for TC is effective only in specific parameter regime (48). To understand this observation, we hypothesize that the Kai system can achieve TC through kinetic regulation of total KaiA concentrations, since KaiA is involved in all phosphorylation rates. Indeed, recent experiments have shown that TC is achieved at a moderate level of KaiA, but not near the onset of oscillation at a high KaiA level (50). Another in vivo study demonstrates that the TC property is abolished in a mutant lacking active KaiA (51). These findings indicate that TC exists only within an intermediate range of KaiA levels.
To better understand the dependence of TC on the concentration [KaiA], we performed simulations based on the vZH model (see SI Appendix, section 1.D for model details). As we vary the KaiA concentration, there are two onsets of the lower onset is determined by the minimum amount of KaiA needed for phosphorylation reactions; and the higher onset is determined by the maximum KaiA concentration beyond which synchronization fails due to lack of competition among individual KaiC hexamers for KaiA. We computed Q10 for the oscillation period in the range of KaiA concentration between these two onsets. Consistent with the experimental data (50, 51), we found that the vZH model fails to achieve TC (Q10<0.9) when [KaiA] is close to the two onsets; however, the value of Q10 peaks around 1 (perfect TC) at intermediate levels of KaiA away from the two onsets as shown in Fig. 5B. Similar to previous models considered in this paper, the nonmonotonic change in Q10 with increasing KaiA levels relies on the crucial kinetic rate kp6 of the sixth phosphorylation step, whose period sensitivity is negative or close to zero near the onsets of oscillation but becomes the largest positive period sensitivity in the intermediate range of [KaiA] (Fig. 5E).
To investigate how the period sensitivity of kp6 becomes large and positive in the intermediate range of [KaiA], we combine the state variables into three X=[C5]+[AC5], Y=∑i=04[Ci]+[ACi], and Z=[KaiC]total−X−Y (Fig. 5A), and study the limit cycle projected in the three-dimensional phase space (Fig. 5C). The oscillation is evidently separated into three S, S^′^, and F. Each phase represents the conversion between two groups with the third group roughly constant. For example, the S phase is when Z is converted to Y while X remains unchanged. We investigate how the amplitude (i.e., the trajectory length), mean progression speed, and the time of each phase vary with kp6. As shown in Fig. 5D, the amplitudes of all three phases increase with kp6 in a similar behavior; as discussed previously, this amplitude-increasing effect is the necessary condition for a period-lengthening reaction. What distinguishes the three phases is the behavior of the mean progression speeds. While the mean progression speed of the F phase scales with kp6 in a similar behavior to the amplitude, those of the S phase and the S^′^ phase are relatively insensitive to changes in kp6. As a result, the time in the F phase is insensitive to changes in kp6, while the time in S and S^′^ phases increases with kp6. Therefore, OPS is the key to generating a relatively large positive sensitivity for kp6: kp6 only accelerates the F phase, a small part of the period, so that this speeding effect of kp6 cannot fully counterbalance its amplitude-increasing effect. We also observe that this OPS mechanism is absent near the onset of oscillation. In this regime, an increase in kp6 simultaneously elevates both the amplitude and the progression speed of the entire limit cycle, leading to a small period sensitivity (see SI Appendix, Fig. S5 for details).
It has been a long-standing problem whether circadian clocks in various organisms share certain common mechanisms for achieving TC (21, 36, 37, 48, 52–54). In this paper, we provide a new perspective to answer this question. By studying four models of biochemical oscillations, we have demonstrated that kinetic regulation can serve as a general mechanism for achieving TC. We focus on cases where kinetic regulation can be implemented by regulation of molecular concentrations, such as [B] or [D] in the Brusselator model and the total [KaiA] in the vZH model. The underlying mechanism of such kinetic regulation is the emergence of large positive period sensitivities. A large positive period sensitivity of a particular rate k is nontrivial given the general period sensitivity sum rule that constrains the sum of all period sensitivities to be −1. Here, we show that a large positive period sensitivity can be caused by the OPS phenomenon wherein the whole oscillation trajectory is separated into different phases in time and the dominant slow phase(s) has an amplitude that increases with a particular rate k and a progression speed that does not increase with k.
To achieve OPS that is critical for TC, a necessary condition is that the system must be driven far from the onset of limit-cycle oscillation (Hopf bifurcation). The reason is that near Hopf bifurcation the oscillation is sinusoidal (31), leading to a nearly uniform progression speed along the whole limit cycle without OPS. In other words, OPS can only emerge from nonsinusoidal oscillations in the “nonlinear regime” where the progression speed and its sensitivity to various kinetic rates are uneven on the limit cycle. One important class of nonsinusoidal oscillation is known as the relaxation oscillation where the separation of timescales occurs (55, 56). Our results suggest that large positive period sensitivities generally exist in relaxation oscillators, such as the VdP model (Figs. 1 and 2) and the Brusselator (Fig. 3), due to OPS. In another important class of oscillators called the negative-feedback oscillators (28, 57, 58), we only noticed a moderate OPS when the system is driven far away from the onset by kinetic regulation (SI Appendix, section E and Fig. S6). In summary, even though OPS exists in all the oscillators we have studied, the degree of OPS and consequently the degree of TC enhancement by kinetic regulation depend on the structure of the oscillator network.
Naturally, driving a system into its nonlinear regime far from the onset requires extra free-energy dissipation, as we show explicitly in the case of the Brusselator model where higher free-energy dissipation is needed to achieve better TC performance (Fig. 4D). In the Kai system, sustained oscillation consumes adenosine triphosphate (ATP) (45). However, most experiments for TC study in the Kai system maintain a sufficiently high level of ATP (46, 50, 59), which makes the effect of free-energy cost for TC hard to assess. In another study, the Kai system is shown to have metabolic compensation (60); that is, the period is insensitive to the change in ATP level. This raises the question of why the cell maintains a high ATP level and consumes more energy for a circadian clock than the minimum energy required for the onset of oscillation. We speculate that kinetic regulation of the ATP concentration may be used to drive the system away from the onset where the OPS mechanism for TC is at play. If the ATP concentration is lowered to a level close to the onset, the Kai system could lose its TC property. This prediction can be tested in future experiments involving the Kai system and possibly other circadian clocks.
Most biological oscillators contain a small number of molecules and they need to operate in a room-temperature environment with large thermal fluctuations. To suppress noise caused by the stochastic biochemical reactions, biological oscillators deploy various control mechanisms in order to carry out their functions, e.g., accurate oscillation, high entrainability to external signals, and synchronization (34, 41, 61), all of which require kinetic rate regulation and extra free-energy dissipation to drive the system beyond the onset of oscillation. Our present work shows that TC, the robustness of the oscillation period against temperature fluctuation, another important function for a class of biological oscillators, can also be achieved by kinetic regulation with the corresponding energy dissipation. These results suggest that the energy-assisted kinetic regulation scheme may serve as a general mechanism for biological networks in particular biochemical oscillators to control internal and external fluctuations to enhance their functional performance. This not only deepens our understanding of biological clocks but also offers valuable insights into constructing synthetic clocks with relevant biological functions.
Simulations were performed on Matlab using its ordinary differential equation (ODE) solvers ode23s (for the VdP model, Brusselator, and Tyson’s model), ode15s (for the vZH model), and ode78 (for the negative-feedback oscillators for higher computational accuracy).
Additional details of modeling and analysis are described in SI Appendix, SI Text.
The work of Y.T. was partially supported by a NIH Grant (R35GM131734). The work of H.F. was partly supported by NIH Maximizing Investigators’ Research Award (MIRA) (R35GM139622) to Suckjoon Jun. The work of Q.O. was partially supported by NSFC12090054. We thank Michael Rust, Hongli Wang, Fangting Li, and Fangzhou Xiao for their helpful discussions and suggestions. Y.T. would like to thank the Center for Computational Biology at the Flatiron Institute for hospitality while a portion of the work was carried out. H.F. would like to thank Suckjoon Jun for his support on this work.
H.F., C.F., Q.O., and Y.T. designed research; H.F. and C.F. performed research; H.F., C.F., and Y.T. analyzed data; Y.T. initiated project; and H.F., C.F., Q.O., and Y.T. wrote the paper.
The authors declare no competing interest.
Simulation codes used in this paper are available to the readers on GitHub (https://github.com/Eigenboy/TC-kinetic_regulation.git) (62). All other data are included in the manuscript and/or SI Appendix.
Simulation codes used in this paper are available to the readers on GitHub (https://github.com/Eigenboy/TC-kinetic_regulation.git) (62). All other data are included in the manuscript and/or SI Appendix.