Role of the reaction-structure coupling in temperature compensation of the KaiABC circadian rhythm

When the mixture solution of cyanobacterial proteins, KaiA, KaiB, and KaiC, is incubated with ATP in vitro, the phosphorylation level of KaiC shows stable oscillations with the temperature-compensated circadian period. Elucidating this temperature compensation is essential for understanding the KaiABC circadian clock, but its mechanism has remained a mystery. We analyzed the KaiABC temperature compensation by developing a theoretical model describing the feedback relations among reactions and structural transitions in the KaiC molecule. The model showed that the reduced structural cooperativity should weaken the negative feedback coupling among reactions and structural transitions, which enlarges the oscillation amplitude and period, explaining the observed significant period extension upon single amino-acid residue substitution. We propose that an increase in thermal fluctuations similarly attenuates the reaction-structure feedback, explaining the temperature compensation in the KaiABC clock. The model explained the experimentally observed responses of the oscillation phase to the temperature shift or the ADP-concentration change and suggested that the ATPase reactions in the CI domain of KaiC affect the period depending on how the reaction rates are modulated. The KaiABC clock provides a unique opportunity to analyze how the reaction-structure coupling regulates the system-level synchronized oscillations of molecules.


Introduction
The mixture solution of cyanobacterial proteins, KaiA, KaiB, and KaiC, shows the robust structural and chemical oscillations with the period of approximately 24 h when the solution is incubated with ATP in vitro [1][2][3]. This period is insensitive to the temperature change, showing the feature of temperature compensation [1]. An important question is whether this feature has a common molecular mechanism or the same mathematical principle as temperature compensation in generic transcription-translation oscillations (TTO). Circadian clocks in many organisms are driven by the time-delayed negative feedback in the TTO [4,5], whose oscillation period is temperature compensated [6,7]; the ratio of the period in 10˚C difference is 0.9 ≲ Q 10 ≲ 1.1, which is much closer to 1 than the value ≳ 1.5 expected from the temperature dependence of normal biochemical reactions. This temperature compensation has been studied with various theoretical models [8][9][10][11][12][13][14][15][16][17][18][19], but its mechanism remains a challenging problem.
There have been at least three views or approaches to temperature compensation; (i) the balance between opposing reactions, (ii) the correlation between oscillation period and amplitude, and (iii) the role of a critical reaction step. One view is the balance between reactions opposingly working to shorten or lengthen the period upon temperature change [7,9]. Such balancing was mathematically formulated [20] and studied with different models [10][11][12][13][14][15]. In searching for the balancing reactions, various mechanisms were examined, including the balance between negative and positive feedback regulations [9] and the balance between ways of resetting bifurcations [12]. In particular, Lakin-Thomas et al. emphasized that the period length should correlate with the amplitude of oscillations [8], suggesting the possible use of the correlation as a clue to find the reactions responsible for the balance [17]. An alternative approach is to search for the critical reaction step or the molecule that determines temperature compensation. In mammals, phosphorylation of period 2 (PER2) regulated by casein kinase Iε/δ (CKIε/δ) is temperature insensitive [21], and this insensitivity was attributed to the reaction mechanism of the CKIε/δ molecule [22]. Therefore, it is meaningful to analyze the temperature compensation of the KaiABC posttranslational oscillations from the three views discussed for TTO. Previously, two hypotheses on the KaiABC temperature compensation were proposed based on the view (iii) of the critical reaction step and the view (i) on the balance between competing reactions. A critical reaction step is ATP hydrolysis in KaiC. KaiC hexamer is a slow ATPase, and variation of the ATPase activity among KaiC mutants is correlated with the variation of the oscillation frequency of those mutants [23,24]. Because the ATPase activity of KaiC is temperature insensitive [23], this correlation suggested that ATPase activity determines temperature compensation in the KaiABC oscillations. In order to clarify whether such causality exists behind the observed correlation, further experimental and theoretical investigations are necessary. Another hypothesis was based on the assumption of balance among the competing binding reactions of KaiA to KaiC at different phosphorylation levels [25]. With this hypothesis, the population of KaiC at a highly phosphorylated state dominates in low temperature, which increases the free unbound KaiA molecules, and increases the overall binding rate of KaiA to compensate for the decrease of reaction rates at the low temperature. However, the accumulation of KaiC at a highly phosphorylated state in low temperature was not observed experimentally [1]; therefore, exploration for other hypotheses and comparison among them are necessary.
Here, we propose a hypothesis based on the view (ii) of the correlation between oscillation period and amplitude. A recent experimental report revealed that substituting an amino-acid residue near the CI and CII domains interface in KaiC induces a striking change in the period from 15 h to 158 h with a suggested correlation between period and amplitude [26]. The modified period length in the mutant was anti-correlated with the volume of the residue after the substitution [26], indicating that the structural coupling between the CI and CII domains is crucial to determining the period; the period was enlarged when the structural coupling between domains was weaker with the smaller volume of residue at the domain interface. The recent X-ray analysis showed that KaiC undergoes structural transitions depending on the state of two phosphorylation sites in the CII and the nucleotide-binding state in the CI [27]. This observation indicated that the structural transitions of KaiC take place through the allosteric communication between CI and CII domains depending on the phosphorylation reactions in the CII and the ATPase reactions in the CI. On the other hand, the rates of those reactions should depend on the structure. Therefore, it is plausible to assume the feedback relations among reactions and structural transitions. Substituting smaller volume residue at the CI-CII interface should reduce the transition cooperativity and weaken the feedback coupling. The weakening of the negative feedback lengthens the period in general nonlinear oscillators as found in the TTO model [9]; therefore, the observed change in the mutants can be explained if the coupling among reactions and structural transitions of KaiC constitutes the negative feedback relation. Modifying strength of the negative feedback leads to the correlated change in amplitude and period in general nonlinear oscillators [9]. Ito-Miwa et al. compared five examples of the CI-CII interface mutations, showing two of them have smaller amplitude with shorter period [26], consistently with the negative-feedback hypothesis. However, the amplitude change of the rest three mutations was not sufficiently quantified [26], suggesting the need for the further statistical evaluation of the mutant oscillations.
In the present study, we propose that the structural coupling between the CI and CII is weakened through thermal fluctuations, weakening the negative feedback relations. In the higher temperature, the larger thermal fluctuations at the CI-CII interface should obscure the specific atomic interactions at the interface, producing a similar effect to the substitution to the smaller volume residue. This weakening of interactions at the interface reduces the negative feedback strength, enhancing the oscillation amplitude, lengthening the period, and compensating for the thermal acceleration of reactions in the higher temperature. We analyze this hypothesis with a model of the KaiABC oscillator and discuss possible tests of the model prediction. We also analyze the correlation between the ATPase reactions and the oscillation frequency to discuss the role of the ATPase reactions in the Kai system temperature compensation.

Problems at two levels; the molecular and ensemble levels
In modeling the KaiABC oscillator, we need to analyze two mechanisms: how individual KaiC molecules oscillate and how oscillations of many KaiC molecules synchronize to generate the ensemble-level oscillations in solution. Many theoretical works focused on the latter question as the sequential change of phosphorylation level of individual KaiC molecules was assumed in advance; then, the synchronization was explained using various hypotheses [25,[28][29][30][31][32][33][34][35][36][37][38][39].
A plausible hypothesis is based on KaiA sequestration [25,[34][35][36][37][38][39][40][41][42][43][44]; the preferential KaiA binding to particular KaiC states sequestrates KaiA, reducing the KaiA binding rate in the other KaiC states, leading to the accumulation of the population in those states, and producing coherent synchronizaed oscillations. This hypothesis is consistent with the experimental observation that the ensemble oscillations disappear when KaiA is too abundant in the solution [45]. Various states of KaiC were assumed as the KaiA-sequestrating states; some models assumed KaiA is sequestrated into the lowly phosphorylated KaiC in the phosphorylation (P) process [25,[34][35][36]40] or in the dephosphorylation (dP) process [37,38]. The other models assumed that KaiA is sequestrated by forming the KaiC-KaiB-KaiA complexes that appear during the dP process [39,[41][42][43][44]. The present author showed [44] that the assumption of the KaiA sequestration into the KaiC-KaiB-KaiA complexes quantitatively explains the experimental data on how the oscillations are entrained when two solutions oscillating at different phases are mixed [46].
In the present study, we use the hypotheses of the KaiA sequestration into the KaiC-KaiB-KaiA complexes to explain the synchronization. We also model the mechanism of how oscillations are driven in individual KaiC molecules. In this way, we address the questions extending over the two levels, the individual molecular level and the ensemble level, to analyze how the molecular-level feedback coupling determines the ensemble-level oscillations and temperature compensation.

Feedback coupling of reactions and structural transitions at the molecular level
At the molecular level, KaiC forms a hexamer [47][48][49], which is denoted here by C 6 . The KaiC monomer is composed of the N-terminal (CI) and C-terminal (CII) domains [50], which are assembled to the CI and CII rings in C 6 [49] (Fig 1). The NMR [51,52], small-angle X-ray diffraction [53], biochemical analyses [54], and X-ray crystallography [27] showed the cooperative structural transitions of KaiC hexamer between two, the structure in the P phase and the structure in the dP phase. We use the order parameter 0 � X � 1 to describe the transitions between two typical conformations and thermal fluctuations around each conformation. We write X(k, t) � 1 when the kth KaiC hexamer at time t takes the structure in the P phase, and X(k, t) � 0 when it takes the structure in the dP phase.
In the present model, we consider that the structural state is a hub of multifold feedback relations among reactions and structural transitions (Fig 1). We describe individual KaiC molecules with the structural parameter X and coarse-grained variables representing three types of reactions; (1) the binding/unbinding reactions of KaiA and KaiB to/from KaiC, (2) the P/dP reactions in the CII, and (3) the ATPase reactions in the CI. We consider that these three types of reactions directly or indirectly depend on X, and these reactions affect X, constituting the multifold feedback relations.
Binding/unbinding reactions of KaiA and KaiB. KaiA and KaiB bind/unbind to/from KaiC in a coordinated way [55][56][57][58][59][60]. The CII ring of a KaiC hexamer can bind a KaiA dimer during the P process to form C 6 A 2 [36,61]. The cryo-electron microscopy and mass-spectrometry showed that each CI domain can bind a KaiB monomer to form KaiC-KaiB complexes, further binding KaiA dimers to form KaiC-KaiB-KaiA complexes, C 6 B i A 2j with j � i � 6 [60]. The stoichiometry C 6 B i A 2j implies the large capacity of KaiC-KaiB-KaiA complexes to absorb KaiA molecules. We consider the probability of the kth KaiC hexamer forming C 6 A 2 at time t, P C 6 A 2 ðk; tÞ, and the probability forming C 6 B i A 2j , P C 6 B i A 2j ðk; tÞ.
Binding/unbinding of KaiA to/from the CII has a timescale of seconds [36]. Because this kinetics is much faster than the other reactions in KaiC, we describe P C 6 A 2 with the quasiequilibrium approximation, P C 6 A 2 ðk; tÞ ¼ x A ðtÞg C:A ðk; tÞP C 6 B 0 A 0 ðk; tÞ with g C:A (k, t) = h A (k, t)/ f A (k, t), where x A (t) is the concentration of free unbound KaiA dimer at time t. h A (k, t) and f A (k, t) are the binding and unbinding rate constants of a KaiA dimer to and from the CII, respectively. We represent the preferential KaiA binding to the X � 1 structure by assuming that h A (k, t) is an increasing function and f A (k, t) is a decreasing function of X(k, t) (See Methods).
Binding/unbinding of KaiB has a timescale of an hour [36,39]. We describe the slow temporal variation of P C 6 B i A 2j ðk; tÞ by integrating the kinetic equations (Methods), which are represented with the rate constants of KaiB binding and unbinding, h B (k, t) and f B (k, t), respectively, and the rate constants of KaiA binding and unbinding to and from the KaiB, h AB and f AB , respectively. We represent the tendency of preferential binding of KaiB to the X � 0 structure by assuming that h B is a decreasing function and f B is an increasing function of X(k, t). Because KaiA does not directly interact with KaiC in this process, we assume h AB and f AB are independent of X. The unusually slow yet specific binding kinetics of the KaiB should be attributed to the fold transitions of KaiB. KaiB switches between ground state (KaiB gs ) and fold-switched state (KaiB fs ) by changing its secondary structures [62]. When only KaiB fs has a significant binding affinity to the CI and KaiB fs is the excited state with the activation energy of ΔE gs!fs > 0, the small factor exp(−ΔE gs!fs /(k B T)) explains the small h B . P/dP reactions. Each CII domain has two sites, Ser431 and Thr432, to be phosphorylated, which amounts to 12 sites in a KaiC hexamer. For simplicity, we do not distinguish Ser431 and Thr432 in the present expression, describing the phosphorylation level of 12 sites with the parameter 0 � D(k, t) � 1; D(k, t) = 1 when 12 sites in the CII of the kth KaiC hexamer are all The CII has twelve sites to be phosphorylated. The KaiC hexamer undergoes allosteric transitions between two structures; the structure (X � 1) in the phosphorylation (P) phase and the structure (X � 0) in the dephosphorylation (dP) phase. A KaiA dimer can bind on the CII of the X � 1 KaiC with the probability P C 6 A 2 , which promotes the P process to increase the phosphorylation level D. KaiB monomers can bind on the CI of the X � 0 KaiC, and a KaiA dimer can further bind on each KaiB monomer forming KaiC-KaiB-KaiA complexes with the probability P C 6 B i A 2j with j � i � 6. KaiA unbinds from the CII of the X � 0 KaiC, promoting the dP process to decrease D. Based on these observations, the model describes the feedback coupling in the KaiC hexamer by assuming (i) the KaiA binding on the CII stabilizes the X � 1 state, and (ii) the KaiB binding on the CI stabilizes the X � 0 structure. The assumptions (i) and (ii) account for the positive feedback to stabilize the X � 1 and X � 0 states. The model also assumes that (iii) the gradual rise of D destabilizes the X � 1 structure and (iv) the gradual fall of D destabilizes the X � 0 structure. The assumptions (iii) and (iv) support the time-delayed negative feedback to drive the transitions between the X � 1 and X � 0 states. The model further assumes (v) the stochastic ATP hydrolysis in the CI (q ! 1) destabilizes the X � 1 state, and (vi) the stochastic ADP release from the CI and the subsequent ATP binding (q ! 0) destabilize the X � 0 state. phosphorylated, and D(k, t) = 0 when they are all unphosphorylated. Phosphorylation is promoted when KaiA binds on the CII [39,55], and dephosporylation proceeds when it unbinds [39]. We represent this tendency by writing d dt Dðk; tÞ ¼ k p H þ ðk; tÞ 1 À D k; t ð Þ ½ � À k dp H À ðk; tÞDðk; tÞ; ð1Þ where H + (k, t) = z/(1 + z) and H − (k, t) = 1/(1 + z) represent the effects of binding and unbinding of KaiA with z ¼ P C 6 A 2 ðk; tÞ=P 0 and a constant P 0 . For changing D between 0 and 1 in �12 h, k p and k dp should be of the order of 0.1 h −1 . ATPase reactions. Both CI and CII domains have ATPase activity, hydrolyzing about 10 ATP molecules in each CI domain and several ATP molecules in each CII domain in a day [23,24]. It is reasonable to consider that ATP is consumed in the CII for supplying a phosphate group in the P process [40], but the reason for the ATP consumption in the CI has been elusive. With the present treatment D(k, t) implicitly represents the ATPase reactions in the CII, and we more focus on the ATPase reactions in the CI. We consider the case ATP is abundant in the solution; therefore, the probability that the CI binds no nucleotide is small. Hence, we use the variable qðk; tÞ ¼ The ADP release and the subsequent ATP binding are the transition from q(i; k, t) = 1 to 0, and hydrolysis of the bound ATP is the transition from q(i; k, t) = 0 to 1. We simulate the stochastic ADP release and the ATP hydrolysis by treating q(i; k, t) as a stochastic variable changing with the lifetime of the ADP bound state, Δ ADP , and the frequency of hydrolysis, f hyd . The ATPase activity measured by the amount of the released ADP from KaiC is large in the P process [23], and the X � 0 (X � 1) structure binds ADP (ATP) [27]. These observations are consistent with the assumption that f hyd is a constant independent of X(k, t) and Δ ADP (k, t) is a decreasing function of X(k, t). In practice, we represent this tendency as where D 0 ADP is a constant determining the timescale and C X is a constant determining the seisitivity to the structure.
Feedback coupling through structural change. Allosteric transitions in protein oligomers typically have a timescale of 10 −3 � 10 −2 s [63], and we assume a similar timescale in the present problem. Because the other reactions in our system are much slower, we describe the KaiC structure as in quasi-equilibrium by treating the chemical states as the quasi-static constraints. Representing the constraints with a variable R(k, t), we have the expression, Xðk; tÞ ¼ Þ (See the Methods section). Here, the explicit form of R(k, t) constitutes major assumptions on the feedback couplings in the present model. Expanding R(k, t) up to the linear terms of P C 6 A 2 , P C 6 B i A 2j , D(k, t), and q(k, t), we have where d 0 is a constant to determine the average structure, and d 1 , d 2 , d 3 , and d 4 are constants defining the strength of the feedback coupling. We assume d 1 > 0, which stabilizes the X � 1 structure when KaiA binds on the CII, and d 2 > 0, which stabilizes the X � 0 structure when KaiB binds on the CI. Then, with the definitions of h A (k, t), f A (k, t), h B (k, t), and f B (k, t), we see that binding reactions constitute the positive feedback loops to stabilize the two states, the X � 1 and X � 0 states.
We use a constant d 3 > 0 in Eq 4, which destabilizes the X � 1 state when phosphorylated (D � 1) and destabilizes the X � 0 state when dephosphorylated (D � 0). This assumption is consistent with the X-ray crystallography analysis that Ser431 is phosphorylated in the X � 0 structure while dephosphorylated in the X � 1 structure [27]. Increase in D from D � 0 to � 1 is promoted by the KaiA binding to the CII of KaiC, which is promoted in the X � 1 state, which then destabilizes the X � 1 state through the d 3 term. Decrease in D from D � 1 to � 0 is promoted by the KaiA unbinding from the CII of KaiC, which is promoted in the X � 0 state, which then destabilizes the X � 0 state through the d 3 term. Therefore the d 3 > 0 in Eq 4 constitutes the negative feedback loops. Because k p and k dp in Eq 1 are small, this destabilization of X is a slow process gradually proceeding after the structural transition. Therefore, the d 3 term represents the time-delayed negative feedback loops to drive the oscillations between the X � 1 and X � 0 states. d 4 F(q(k, t)) in Eq 4 represents the coupling of structure with the ATPase reactions, which largely determines the balance between positive and negative feedback effects. This coupling was inferred from the observations that the ATP hydrolysis is necessary for binding KaiB to KaiC [51,54,64,65], and that the structure is modified upon ATP hydrolysis [24,27,60,66].
) and d 4 > 0, the ADP binding on the CI destabilizes the X � 1 state, while the ATP binding destabilizes the X � 0 state. Therefore, the stochastic ATP hydrolysis triggers the transition to the X � 0 state, and the stochastic release of ADP with the subsequent ATP binding triggers the transition to the X � 1 state (Fig 1). Thus, the model describes the positive feedback between the structure and the binding/ unbinding of KaiA and KaiB to/from KaiC (the d 1 and d 2 terms), the negative feedback between the structure and the P/dP process (the d 3 term), and the transition-triggering effects of the ATPase reactions (the d 4 term). These couplings generate the cooperative chemical and structural oscillations in individual KaiC molecules.

Communication among many KaiC molecules at the ensemble level
We simulated the ensemble of N = 1000 or 2000 KaiC hexamers. For N = 1000 and V = 3 × 10 −15 l, the concentration of KaiC is C T = 3.3 μM on a monomer basis, which is near to the value 3.5 μM often used in experiments. We assumed the ratio A T : B T : C T = 1 : 3 : 3 as used in many experiments [23,54,67], where A T and B T are total concentrations of KaiA and KaiB on a monomer basis. The system was described by variables, tÞ and D(k, t) were calculated by numerically integrating the kinetic equations, and P C 6 A 2 ðk; tÞ and X(k, t) were calculated with the quasi-equilibrium approximation at each time step. q(k, t) was calculated by simulating the stochastic transitions of q(i; k, t) between 0 and 1 with frequencies D À 1 ADP and f hyd . Concentrations of free unbound KaiA dimer and KaiB monomer, x A (t) and x B (t), were calculated at each time step from the following equations of conservation, ð5Þ Competition between the 2nd and 3rd terms of the l.h.s. of Eq 5 provides communication among KaiC molecules leading to the synchronization. See the Methods section for details.  (Fig 2A). The nucleotide-binding state in the CI ring q(k, t) also exhibits transitions between the state rapidly fluctuating around q(k, t) � 0.3 and the state around q(k, t) � 0.6. The phosphorylation level D(k, t) follows these switching transitions with the slower rates of the P/dP reactions; D increases in the X � 1 state and decreases in the X � 0 state, showing saw-tooth oscillations.

Single-molecule and ensemble-level oscillations
At the ensemble level, these fluctuating oscillations in individual molecules are averaged, resulting in the regular oscillations of � ). The ensemble-averaged rate of the ADP release from KaiC, � q r ðtÞ (Methods), is large during the P phase as observed experimentally [23], and the ADP binding probability � qðtÞ is large during the dP phase consistently with the experimental observations [27]. In this way, the model reproduces the stable circadian oscillations at the ensemble level averaging the synchronized individual KaiC oscillations.

Modifications of the feedback strength
The binding of KaiA or KaiB may induce the global change of each KaiC subunit, shifting the position and orientation of subunit, whose effects represented by d 1 and d 2 in Eq 4 should be insensitive to the single-residue substitution at the CI-CII interface. On the other hand, P/dP in the CII or the ATPase reaction in the CI is the local atomic change around the phosphate group, whose effects are transmitted through chains of electrostatic and volume-excluding interactions, producing the allosteric communication between the CI and CII [27,66]. This process should be sensitive to the atomic interactions at the CI-CII interface and hence sensitive to the single-residue substitution. Here, we assume that substituting a CI-CII interface residue to the smaller volume one is represented by the decrease in d 3 and d 4 in Eq 4 with a scaling factor s < 1 to sd 3 and sd 4 while d 1 and d 2 in Eq 4 being kept in their original values. With s < 1, the simulated phosphorylation level oscillations indeed show the large amplitude and long period as expected from the decrease in the negative feedback strength, while the amplitude is small and the period is short for s > 1 (Fig 3A). On the contrary, with positive feedback enhancement with a scaling factor s > 1, changing d 1 and d 2 to sd 1 and sd 2 enlarges the amplitude and period, while they are reduced when s < 1 (Fig 3B). These effects of the feedback strength modifications were systematically examined in Fig  3C and 3D. When d 3 is scaled to sd 3 with d 4 , d 1 , and d 2 kept constant, the amplitude and period are extended as s is decreased. When s ≲ 0.5, the oscillations disappear as the system is caught at the X � 1 state losing the negative feedback destabilization of the X � 1 state. Similar behaviors were found when the structure-ATPase coupling was changed to sd 4 with d 3 , d 1 , and d 2 kept constant. Thus, the ATPase reactions give similar effects to the negative feedback in the present model. With the combined change to sd 3 and sd 4 , the period change is more prominent, ranging from 9 h to 120 h (Fig 3C), which explains the observed ten-times period change induced by the single-residue substitution at the CI-CII interface [26]. Modifying the positive feedback strength to sd 1 and sd 2 shows that the oscillations disappear when the positive feedback is too weak or too strong (Fig 3D). The period and amplitude are enlarged as s increases in between these boundaries of positive feedback strength (Fig 3C).

Thermal loosening of the strcutural coupling explains temperature compensation
We propose that the increased structural fluctuations of KaiC in the higher temperature should weaken the interactions at the CI-CII interface to decrease d 3 and d 4 in Eq 4, producing similar effects to the single-residue substitution (Fig 3A). These effects should enlarge the amplitude and period, compensating for the accelerated reactions in high temperatures. Here, we examine our hypothesis by simulating the oscillations in various temperatures.
Assuming that the fluctuations are proportinal to exp(-ΔE f /(k B T)) with a constant ΔE f , the feedback coupling strength at temperature T, d 3 (T) and d 4 (T), should be modified from those We call this temperature dependence of d 3 (T) and d 4 (T) Rule 1. Another rule comes from the activation energy ΔE gs!fs > 0 and ΔE fs!gs > 0 for the KaiB fold transformations. We assume h B / exp[−ΔE gs!fs /(k B T)] and f B / exp½À DE 0 fs!gs =ðk B TÞ�, which affects the binding/ unbinding kinetics of KaiB. Here, DE 0 fs!gs ¼ DE fs!gs þ DE bound and ΔE bound is the KaiC-KaiB binding energy. We call this assumption on h B and f B Rule 2. The other assumption corresponds to the temperature insensitivity of the ATPase reactions as observed in experiments [23,68]. We assume the temperature-insensitive ATPase reactions by imposing D 0 is a constant to determine the lifetime of the ADP bound state (Eq 3), and f hyd (T) is the hydrolysis frequency of the bound ATP. We call this assumption Rule 3. These three rules are summarized in Table 1.
The model has six rate constants ( Table 2) other than the four rate constants, h B , f B , D 0 ADP À 1 , and f hyd , discussed in Table 1. Each of six rate constants should depend on temperature in each distinctive way. However, for examining the hypothesis transparently, we assume a straightforward case of the same activation energy ΔE 0 in six rate constants, which leads to the same temperature dependence as k dp (T) = s(ΔE 0 ; T, T 0 )k dp (T 0 ), etc., where s is defined in Eq 6. In the case we do not impose any of Rule 1, Rule 2, or Rule 3 by assuming D 0  Table 1 even in the case ΔE 0 = 10k B T 0 in other six rate constants. We also show that Rule 1 plays a dominant role, and temperature compensation is realized without imposing Rule 2 or 3; temperature compensation is realized if Rule 1 is satisfied even when ΔE 0 = 10k B T 0 is used for all ten rate constants.  In order to distinguish the roles of three rules, we examined four cases ( Table 1). All three Rules apply in Case A (Fig 4A), Rules 1 and 3 in Case B (Fig 4B), Rules 1 and 2 in Case C ( Fig  4C), and Rules 2 and 3 in Case D (Fig 4D). The period is temperature compensated in Case A (Q 10 = 1.01), Case B (Q 10 = 0.93), and Case C (Q 10 = 1.02) with a slight overcompensation in Case B, while the period shows a distinct temperature dependence in Case D (Q 10 = 1.51), showing that Rule 1 plays a dominant role in temperature compensation. The overcompensation in Case B shows that Rule 2 enhances the temperature dependence of the period, which was compensated for by Rule 1 in Cases A and C. In Cases A, B, and C, amplitude sharply decreases as temperature decreases below 20˚C as observed experimentally [69]. We should  note that Rule 3, temperature insensitivity of ATPase reactions, does not play a significant role in the present examination ( Fig 4C); we discuss more on this point in the "Effects of the ATPase reactions" subsection.
A determinant role of Rule 1 is evident also from calculations with the varied temperature dependence of the structural fluctuations. Q 10 decreases as ΔE f increases ( Fig 5A); Q 10 = 1.51 (ΔE f = 0), 1.28 (ΔE f = 3k B T 0 ), 1.08 (ΔE f = 6k B T 0 ), and 0.90 (ΔE f = 9k B T 0 ); showing the temperature compensation with 6k B T 0 , and the distinct overcompensation with 9k B T 0 . A drop of the amplitude in the low temperature regime takes place at the higher temperature as ΔE f increases (Fig 5B), consistently with the expected period-amplitude correlation. Fig 5 shows a sensitive dependence of Q 10 on ΔE f . In KaiC mutants showing the altered oscillation period, we can expect the intact Q 10 � 1 when the mutations perturb the P/dP-reaction rates in the CII or the ATPase-reaction rates in the CI, which alters the period, but does not alter the CI-CII interface, and hence does not affect ΔE f much. However, with a single-residue substitution at the CI-CII interface of KaiC, the substitution can change ΔE f . In the experimental report, some substitutions at the CI-CII interface showed the intact Q 10 , but the others showed the enlarged Q 10 [26]. A possible explanation for the former case is that the protein region around the substituted residue compensates for the structural fluctuations' temperature dependence, leading to the robust ΔE f against the mutation, while in the latter case, such compensation does not sufficiently work, leading to the enlarged Q 10 . It is crucial to examine whether such a difference in ΔE f indeed occurs in those mutants.

Phase shift caused by a stepping change in temperature
A significant feature of circadian rhythms is the response of their oscillation phase to the external stimuli. In particular, the phase is distinctively shifted by a temperature change in circadian oscillators while their period is temperature compensated [70,71]. Also, in the in vitro KaiABC system, Yoshida et al. [72] found the phase shift indued upon temperature change. A step-up of temperature from T = 30˚C to 45˚C in the dP process advanced the phase, while a step-up in the P process delayed the phase. On the other hand, a step-down of temperature from T = 45˚C to 30˚C in the P process advanced the phase, and a step-down in the dP process delayed the phase. Here, we compare the simulated phase shift in the present model with the experimental results in [72]. In order to avoid the confusion coming from the slight period change upon temperature stepping, we defined the circadian time (CT) in the same way as in Ref. [72]: timing of the peak in the simulated rhythm at each temperature was set to CT 16, and the period length was normalized to be 24 h in CT (Fig 6A).
In the present model, oscillations of individual KaiC molecules are synchronized by the KaiA binding to the KaiC-KaiB-KaiA complexes, which entrains oscillating molecules in the ensemble into the dP phase [44]. With Rule 2 in Table 1, the T step-down increases the ratio h B /f B , enhancing the entrainment and stabilizing the dP phase. Thus, the T step-down in the dP phase elongates the dP process and causes the phase delay. On the other hand, the T stepup in the same dP phase causes the opposite effect of the phase advance.
This expectation is confirmed in the results shown in Fig 6B-6G. In Fig 6B and 6C, the temperature was stepped up from T = 30˚C to 45˚C. Fig 6B shows Table 1 were used with ΔE 0 = 6k B T 0 and N = 1000. https://doi.org/10.1371/journal.pcbi.1010494.g006

PLOS COMPUTATIONAL BIOLOGY
Temperature compensation of the KaiABC circadian rhythm 6C shows the phase advance caused by a T step-up in the dP phase at CT 17.3 (IT 95). In Fig  6D and 6E, the temperature was stepped down from T = 45˚C to 30˚C. Fig 6D shows the phase advance caused by a simulated T step-down in the P phase at CT 7.9 (IT 80.4), and Fig 6E shows the phase delay caused by a T step-down in the dP phase at CT 20 (IT 92). These results can be quantified by plotting a phase response curve (PRC). Fig 6F and 6G show PRCs for T step-up and T step-down, respectively. These simulated PRCs qualitatively reproduce the experimentally observed PRCs in [72]: upon T step-up, the phase was advanced at CT 16-CT 22 and delayed at CT 8-CT 14, and upon T step-down, the phase was advanced at CT 6-CT 11 and delayed at CT 14-CT 24. Thus, our model reproduces the essential features of the experimentally observed phase shift upon the stepping change of temperature.

Effects of the ATPase reactions
In order to clarify the effects of the ATPase reactions on the oscillations, we scaled the rate constants of the ATPase reactions in the present model. Fig 7A shows the amplitude and period calculated by modifying the inverse lifetime of the ADP bound state, D 0 ADP À 1 , and the hydrolysis frequency of the bound ATP, f hyd . They were scaled by a factor s a in two different ways; in case I, they were scaled as s a ðD 0 ADP Þ À 1 and s a f hyd , and in case II, ðs a D 0 ADP Þ À 1 and s a f hyd . In case I, period and amplitude only slightly depend on s a for s a > 0.5. This insensitivity to the scaling corresponds to the temperature insensitivity shown in Fig 4C, where we did not use Rule 3 but scaled the constants as D 0 ADP ðTÞ À 1 ¼ sðDE 0 ; T; T 0 ÞD 0 ADP ðT 0 Þ À 1 and f hyd (T) = s(ΔE 0 ; T, T 0 )f hyd (T 0 ). The insensitivity in case I and Fig 4C indicates that the change in the ATPase reaction rates does not affect the period when they were changed with the constraint Δ ADP � f hyd = const. The ATPase reactions should affect the period in two ways: The probability that the CI binds ADP during the P phase is correlated to how the X � 1 state is destabilized, and the probability that the CI binds ATP during the dP phase is correlated to how the X � 0 state is destabilized. However, changing the rates under the constraint Δ ADP � f hyd = const does not affect these probabilities. Therefore, the period was insensitive to the changes in the rate constants under the constraint Δ ADP � f hyd = const. In case II, this constraint is not satisfied, and the change in the rate constants brings about a significant change in the period (Fig 7A). We examined the correlation between the oscillation frequency (i.e., 1/period) and the ATPase activity calculated in the nonoscillatory condition in the absence of KaiA and KaiB. Fig 7B shows a clear correlation between the oscillation frequency and the ATPase activity in case II as experimentally observed [23,24]. The period is insensitive to the change in the ATPase rates only when the specific constraint is satisfied. This result suggests that the strategy to keep the specific constraint has been avoided in evolutionary design. Instead, the more general strategy of the temperature-insensitive ATPase reactions might have been evolutionarily realized in KaiC [23,68] as in the TTO regulator, CKIε/δ, in mammals [22].
The present analysis showed that the oscillation period changes when Δ ADP � f hyd varies; therefore, the temperature dependence of the period may be enhanced in mutants in which the temperature dependence of Δ ADP and/or f hyd was modified, making Δ ADP � f hyd dependent on temperature. Indeed, when D 0 ADP ðTÞ T 0 )f hyd (T 0 ) with ΔE a 6 ¼ ΔE b , the temperature dependence of the period should be enhanced. However, the temperature dependence of the structure-reaction coupling with ΔE f > 0 at the CI-CII interface should compensate for this effect. For example, the present model calculation with ΔE a = 5k B T 0 , ΔE b = 0, and ΔE f = 7k B T 0 showed Q 10 = 1.04, showing only a mild increase of Q 10 from Q 10 = 1.01 in Fig 4A. However, with the larger temperature dependence of the ATPase activity as ΔE a = 10k B T 0 and ΔE b = 0 with ΔE f = 7k B T 0 , we have Q 10 = 1.06; therefore, the further analyses are necessary for elucidating the relationship between the temperature dependence of the period and the temperature dependence of the ATPase activity. It is important to analyze how the Q 10 of the oscillation period varies among mutants showing the larger Q 10 of the ATPase activity [68].

Phase shift caused by a pulse of the increased ADP concentration
Examining the response of the oscillation phase to the change in the ATP/ADP concentration ratio should further test the roles of the ATP hydrolyses in the KaiABC oscillations. Rust et al. [73] experimentally examined the phase response of the in vitro KaiABC system by adding ADP to the reaction buffer at various timing. After several hours of ADP addition, they reduced the ADP concentration to the original level by adding pyruvate kinase, which converts ADP to ATP. Rust et al. found that this pulse of ADP addition delayed the oscillation phase when the pulse was added at around the trough of the oscillating phosphorylation level, while the pulse advanced the phase when the pulse was added at around the peak of the phosphorylation level [73]. In our model, the increased concentration of ADP should elongate the ADPbound state's lifetime, enhancing the transition probability from the P phase to the dP phase and reducing the probability of the opposite transition. Therefore, we expect that adding the ADP pulse at the oscillation peak helps the transition to advance the phase, while adding the ADP pulse at the oscillation trough disturbs the transition to delay the phase. This expectation is consistent with the results observed in [73], and we confirmed it with simulations as in the following. In our model, the effects of an increase in the ADP concentration should correspond to the increase in the ADP-bound state's lifetime. Therefore, we simulated the ADP pulse by increasing D 0 ADP at a particular time and reducing D 0 ADP to the original value 6 hours after the increase. We defined CT in the same way as in Fig 6 (Fig 8A) and used the same parameters as in Fig  4A. Fig 8B and 8D show the example phase shift with the simulated addition of ADP pulse around the trough and the peak of oscillations, respectively. As expected, the addition of the pulse around the trough delayed the oscillation phase (Fig 8B), and the addition of the pulse around the peak advanced the phase (Fig 8D). It is intriguing that the simulated change from the phase-delay to phase-advancement behaviors was abrupt, inducing a jump in the PRC at around CT 5 (Fig 8E), reproducing the observed jump in the experimental PRC [73]. The model predicts that the addition of the ADP pulse at the jump point at CT 5 (IT 81.6) advanced the oscillation phase of almost half of KaiC molecules and delayed the phase of the rest half of molecules, which strongly desynchronized the oscillations of individual KaiC molecules, diminishing the ensemble-level oscillations ( Fig 8C); it is important to experimentally examine this desynchronization as a test of the oscillation mechanism proposed in the present model.

Discussion
Using the model describing the feedback coupling among reactions and structural transitions in individual KaiC molecules and synchronization of many KaiC molecules, we showed that weakening the negative feedback strength extends the amplitude and lengthens the period of oscillations in the KaiABC system. This weakening can explain the observed wide range of period modification induced by the single-residue substitution at the CI-CII interface of KaiC [26]. We hypothesized that thermal fluctuations induce similar effects to the substitution at the CI-CII interface, which explained the stable temperature compensation in the KaiABC system. The ATPase reactions also affect the period, but the period is insensitive when the ATPase rates are changed under the specific constraint.
A possible test of the thermal weakening of the CI-CII structural coupling is to measure the temperature dependence of the thermal fluctuations of the interface residue with NMR or other spectroscopic methods. The model predicts that the temperature dependence of fluctuations should be anti-correlated with the Q 10 of the period, and such an anti-correlation was suggested from the recent neutron scattering data [68].
Another possible test is to observe the single-molecular behavior. In the typical oscillatory condition, oscillations of individual KaiC molecules are synchronized, inducing coherent circadian oscillations at the ensemble level (Fig 9A). If the synchronization is realized through the sequestration of KaiA into the KaiC-KaiB-KaiA complexes as assumed in the present study, synchronization should be lost when the binding affinity between KaiA and KaiB is reduced. It would be possible to design a KaiB mutant with a low affinity to KaiA, and in the present model, such a mutant is represented by reducing the h AB value. With this reduction, synchronization is lost, and the ensemble oscillations disappear despite the oscillations with the large amplitude remaining in individual KaiC molecules (Fig 9B). In such a desynchronized condition, period of individual oscillations is temperature compensated with the temperature-dependent oscillation amplitude in the present model (Fig 9C). In contrast, when the negative feedback strength does not depend on temperature, individual oscillations are not temperature compensated without showing a significant temperature dependence of amplitude ( Fig 9D). The temperature compensation hypothesis in the present study contrasts with the hypothesis of the competitive KaiA binding reactions [25]. With the latter hypothesis, temperature compensation is induced only from the ensemble level mechanism, not from the individual molecular mechanism, which should be distinguishable with the single-molecule observation in the condition that the synchronization is lost.
With the larger KaiA concentration than the standard value of A T /C T = 1/3, Ito-Miwa et al. showed that Q 10 of the oscillation period becomes large (Fig S6 in [26]). In our model with the parameters used in Fig 6, with the increased KaiA concentration of A T /C T = 1/2, the oscillation period was temperature compensated with Q 10 = 0.98, which disagrees with the experimental report. We should note that the results depend on the parameterization in the model; for example, with ΔE gs!fs = 20k B T and DE 0 fs!gs ¼ 0, Q 10 = 1.33 when A T /C T = 1/2. Therefore, the further examinations should be needed to compare the model results with the experimental report by carefully examining the mechanisms of the binding/unbinding kinetics of KaiA and KaiB to/from KaiC. Our analyses showed the importance of the molecular features which determine the coupling between reactions and structural transitions in the oscillating molecules. Investigations of the KaiABC system should highlight the significance of the regulation through the atomic reaction-structure coupling and help provide innovative methods of design for regulating the system dynamics of an ensemble of molecules.

Multifold feedback coupling in KaiC
We describe individual KaiC hexamers with coarse-grained variables; the binding state of KaiA on the CII domain, y C 6 A 2 ðkÞ, the binding state of KaiB and KaiA on the CI domain, y C 6 B i A 2j ðkÞ, the phosphorylation level, 0 � DðkÞ � 1, the structural state, W(k), and the nucleotide-binding state, q(i; k, t). Here, We consider the system to consist of N KaiC hexamers. Then, the system state at time t is described by a set of vectors,ỹ C 6 A 2 ¼ fy C 6 A 2 ð1Þ; y C 6 A 2 ð2Þ; � � � ; y C 6 A 2 ðNÞg, etc. The stochastic evolution of the system state is described by the probability distribution, where x A (t) and x B (t) are concentrations of free KaiA dimer and KaiB monomer unbound from KaiC, respectively. The structural change takes place in milliseconds or so in usual protein oligomers, and we assume a similar timescale in the present problem. This timescale is much shorter than the timescales of other reactions; therefore, we treat the variableW as in quasi-equilibrium. Then, the quasi-equilibrium free energy G quasi should be expanded byW as Here, because W(k) 2 = 1, the second-order term of the expansion only consists of the sum with k 6 ¼ l. In the present model, the interaction between different KaiC hexamers is indirect through the KaiA sequestration; therefore, we can put G 2 (k, l) = 0. Then, in the expression G quasi = G 0 + ∑ k W(k)G 1 (k), W(k) behaves like an Ising spin under the external field R(k, t) = −G 1 (k). Therefore, the average of W(k) in quasi-equilibrium is � W ðk; tÞ ¼ tanhðÀ bRðk; tÞÞ. By intoducing the order parameter of structure, 0 � X(k, t) � 1, as Xðk; tÞ ¼ ð � Wðk; tÞ þ 1Þ=2, we have showing the same form as in Eq 4 in the main text. Here, β = 1/(k B T) is inverse temperature, and R(k, t) is a quasi-equilibrium constraint reflecting the chemical state of the KaiC hexamer at time t. The explicit form of R(k, t) represents how the feedback relations among reactions and structure work in the system and is defined after the other variables in Eq 12 are transformed to an expression suitable to be handled. We use the mean-field approximation; Eq 12 is approximated by factorizing P into each degree of freedom as where q(k, t) = {q(1; k, t), q(2; k, t), � � �, q(6; k, t)} is a six-dimensional vector representing the nucleotide-binding state of each of six CI domains in the kth hexamer. We represent the nonequilibrium consumption of ATP by treating q(i; k, t) as a stochastic variable taking the value either 1 or 0 depending on whether the CI domain binds ADP or ATP (Eq 11). X(k, t), x A (t), and x B (t) are calculated by solving the relations explained in the subsections "Reaction-structure feedback coupling" and "Coupling of multiple oscillators" in this Methods section. Thus, by dropping the variables, q(k, t), X(k, t), x A (t), and x B (t) from the expression, Eq 14 is Because the KaiA binding and unbinding reactions are faster than the other reactions in the present system, we can approximate Eq 16 with the quasi-equilibrium approximation as d dt P C 6 A 2 ðk; tÞ ¼ 0. Then, we have with g C:A (k, t) = h A (k, t)/f A (k, t). As discussed in the main text, h A should be an increasing function of X(k, t) and f A is a decreasing function of X(k, t). We use the form, where h A0 and f A0 are the rate cosnstants to determine the time scale and A X > 0 is a constant to determine the sensitvity to the structure. In a similar way, by writing P C 6 B i A 2j ðk; tÞ ¼ Pðy C 6 B i A 2j ðkÞ ¼ 1; tÞ, the binding and unbinding of KaiA to and from KaiB should be in quasi-equilibrium, leading to where P C 6 B 0 ðk; tÞ ¼ P C 6 A 2 ðk; tÞ þ P C 6 B 0 A 0 ðk; tÞ, and the rate constants are Here, h B0 and f B0 are the rate constants defining the time scale and B X > 0 is a constant determining the sensitivity to the structure. Because the bindable conformation of KaiB appears with the thermal activation with the energy ΔE fs!gs , we assume h B0 and f B0 at temperature T are as explained in Table 1 in the main text. with the rate constants k p and k dp . Here, H + (k, t) = z/(1 + z) and H − (k, t) = 1/(1 + z) are the effects of binding and unbinding of KaiA to and from the CII, respectively, and z ¼ P C 6 A 2 ðk; tÞ=P 0 with a constant P 0 . ATPase reactions. We describe the non-equilibrium ATPase reactions by using a stochastic variable q(i;k, t). When ATP is bound on the ith CI domain of the kth KaiC hexamer, we write q(i; k, t) = 0. The ATP hydrolysis is represented by the transition from q(i; k, t) = 0 to q(i; k, t) = 1, which takes place at a random timing with the frequency f hyd . We approximate that f hyd does not depend on X, as explained in the main text. The state q(i; k, t) = 1 represents the ADP-bound state. It is not known how the release of inorganic phosphate (P i ) impacts the KaiC structure. In the present study, for simplicity, we do not distinguish the ADP+P i bound state just after the hydrolysis and the ADP bound state after the P i release. We assume that the lifetime of the ADP bound stateD ADP ðk; tÞ is stochastically fluctuating aŝ D ADP ðk; tÞ ¼ D ADP ðk; tÞ þ xðk; tÞ; ð25Þ

PLOS COMPUTATIONAL BIOLOGY
where ξ(k, t) is a random number satisfying hξ(k, t)i = 0 and hxðk; tÞxðk 0 ; t 0 Þi ¼ d kk 0 dðt À t 0 ÞD ADP ðk; tÞ. After the ADP release, the ATP rebinds, which turns the nucleotide-binding state from q(i; k, t) = 1 to 0. As explained in the main text, we assume f hyd ¼ const: ð26Þ and D ADP ðk; tÞ ¼ D 0 ADP 1 À tanh where C X and D 0 ADP are constants. The ensemble-averaged rate of the ADP release, � q r ðtÞ, shown in Fig 2B,  where θ(x) = 1 for x > 0 and θ(x) = 0 for x � 0. δt = 10 −3 h is the width of the simulation time step, and Δt = 0.2 h is the time window for the data sampling. Reaction-structure feedback coupling. R(k, t) in Eq 13 represents the major assumptions on the feedback coupling in the present model. We use Rðk; tÞ ¼ d 0 þ d 1 P C 6 A 2 ðk; tÞ À Coupling of multiple oscillators through conservation of molecules. The constraints coming from the conservation of the total concentrations of KaiA (Eq 5 of the main text) and KaiB are simplified with the quasi-equilibrium treatment of KaiA (Eqs 17 and 19) as g C:A ðk; tÞ 1 þ x A ðtÞg C:A ðk; tÞ P C 6 B 0 ðk; tÞ Simulations We simulated the system containing N = 1000 or 2000 KaiC hexamers by numerically integrating the kinetics with a time step of δt = 10 −3 h. The variables describing the system are P C 6 A 2 ðk; tÞ, P C 6 B i A 2j ðk; tÞ (0 � j � i � 6), q(i; k, t) (1 � i � 6), D(k, t), X(k, t), for k = 1 � N, x A (t) and x B (t). From given values of these variables at time t, the values at t + δt were obtained by (i) stochastically updating q(i; k, t) using constants f hyd and D 0 ADP with Eqs 25-27, (ii) evaluating the binding and unbinding constants of Eqs 18 and 22, (iii) updating P C 6 A 2 ðk; tÞ and P C 6 B i A 2j ðk; tÞ with Eqs 17, 20 and 21, (iv) updating x A and x B by solving Eqs 31 and 32, (v) updating D(k, t) with Eq 24, and (vi) updating X(k, t) using Eq 30. The period and amplitude were calculated from the trajectories, each having the length 3276.8 h obtained after the initial warming-up trajectories of 100 h length.