Role of ATP Hydrolysis in Cyanobacterial Circadian Oscillator

A cyanobacterial protein KaiC shows a stable oscillation in its phosphorylation level with approximately one day period when three proteins, KaiA, KaiB, and KaiC, are incubated in the presence of ATP in vitro. During this oscillation, KaiC hydrolyzes more ATP molecules than required for phosphorylation. Here, in this report, a theoretical model of the KaiABC oscillator is developed to elucidate the role of this ATP consumption by assuming multifold feedback relations among reactions and structural transition in each KaiC molecule and the structure-dependent binding reactions among Kai proteins. Results of numerical simulation showed that ATP hydrolysis is a driving mechanism of the phosphorylation oscillation in the present model, and that the frequency of ATP hydrolysis in individual KaiC molecules is correlated to the frequency of oscillation in the ensemble of many Kai molecules, which indicates that the coherent oscillation is generated through the coupled microscopic intramolecular and ensemble-level many-molecular regulations.

detailed molecular-level model of this interplay has been recently developed by Paijmans et al. 28 , and the further development of a comprehensive model that bridges between the molecular and ensemble levels is necessary.

Multifold feedback model
We consider the following binding/unbinding processes among Kai proteins; Here, C 6 B 0 A 0 is identical to C 6 . Probability of the kth C 6 to be bound in forms of C 6 A 2 and C 6 B i A 2j with 0 ≤ i ≤ 6 and 0 ≤ j ≤ i at time t are denoted by P k t ( , ) C A 6 2 and P k t ( , ) C B A i j 6 2 , respectively with k = 1, …, N, where N is the total copy number of KaiC hexamers in the system. We use the approximation that factorizes an N-body probability distribution function into one-body probability distributions, P k t ( , ) C A 6 2 and P k t ( , ) C B A i j 6 2 . The similar factorization approximation was used to describe single-molecular fluctuation of binding reactions to DNA 29 . Then, equations for P k t ( , ) C A 6 2 and P k t ( , )   where V is the system volume, and A T and B T are total concentrations of KaiA and KaiB, respectively. The small-angle X-ray scattering 30 , the NMR spectroscopy 31,32 , and the biochemical analysis 33 showed that during the phase of phosphorylation (P), KaiC hexamer takes a structure different from that during the phase of dephosphorylation (dP). Oyama et al. referred to the structure in the P-phase as ground state (gs), and the one in the dP-phase as competent state (cs); cs-KaiC can dissociate into monomers and is less stable than gs-KaiC 33 . We assume six subunits in a KaiC hexamer cooperatively undergo the structural transition for the hexamer to exhibit two structural states as assumed in previous models 20,24,25 ; we represent this allosteric transition by a continuous variable X as X k ≈ 1 when the kth KaiC hexamer is in the gs-state and X k ≈ 0 when in the cs-state. Then, the rate constants in Eqs 4 and 5 should depend on the structural order parameter X k as , , and f B0 are constants. With this definition of rate constants, KaiA binds to KaiC with larger affinity in the gs-state, which promotes P-reactions, while KaiB binds to KaiC with larger affinity in the cs-state, which promotes dP-reactions. In order to facilitate the extensive calculation of the many-molecular system, we adopt a simplified representation of the phosphorylation level of 12 sites in a single KaiC hexamer by using a continuous variable D as D k (t) ≈ 1 when 12 sites of the kth KaiC hexamer are fully phosphorylated, and D k (t) ≈ 0 when they are fully dephosphorylated. D k is increased when KaiA binds to KaiC to form C 6 A 2 and is decreased otherwise. This transition is represented by the equation of soft-spin dynamics as where k p and k dp are rate constants of P and dP reactions, respectively, and g(D) = aD(D − 1/2) 2 (D − 1) with a > 0 represents a soft-spin constraint to confine D in a finite range. We assume that structure transition is much faster than the other reactions, so that the structure X k is represented as a quasi-equilibrium average under the mean-field generated by D k and other quantities as where β = 1/k B T, and c 0 , c 1 , c 2 , and c 3 are constants.
k t ( , ) represent the extent of KaiA binding and that of KaiB binding, respectively, where n B represents the level of saturation of the effects of KaiB binding. See Methods section. q k represents the effect of ATPase reactions as explained later.
Equations 4-9 constitute loops of feedback among P/dP reactions and structural transitions, which brings about nonlinear behaviors of X k (t) and D k (t). When c 3 > 0 in Eq. 9, for example, the larger probability of KaiB binding to KaiC brings the structure to the cs-state with a small X k , which further enhances the KaiB binding through Eqs 5 and 7. Therefore, the c 3 > 0 term gives a positive feedback effect. If we would neglect the effects of the c 1 term, this positive feedback should stabilize the (small X k , small D k ) state for . On the other hand, when c 1 > 0 in Eq. 9, the large D k leads the KaiC structure to the cs-state with small X k . From Eq. 7, this brings about the small affinity of KaiA to KaiC, reducing P k t ( , ) C A 6 2 in Eq. 4, which in turn reduces D k in Eq. 8. Therefore, the c 1 > 0 term gives a negative feedback effect, which should stabilize the state with (large X k , small D k ) or (small X k , large D k ). The system tends to stay at these various states, but because of the competition among multiple feedback interactions, the system is easily moved when small perturbations are added to generate oscillation among the states. The q k term representing the effect of ATPase reactions provides such a perturbation.
The ATPase reactions are dominated by the CI ring 17 , whereas the P/dP reactions take place in the CII ring; therefore, a coordinated structural change of the CI and CII domains, i.e., allosteric communication between CI and CII mediates the effects of ATPase reactions and P/dP reactions. In the present model, X k represents such communication. There are accumulating experimental data suggesting that the ATPase activity in the CI is necessary for the binding of KaiC to KaiB 31,33,34 , indicating that the ATP hydrolysis enhances the binding affinity of KaiC to KaiB. Therefore, a plausible assumption is that the ATP hydrolysis changes the KaiC structure from gs to cs. This effect is represented in Eq. 9 by setting = ∑ = q t q i k t ( ) ( ; , ) k i 1 6 with q(i; k, t) = q 0 > 0 for t 0 ≤ t ≤ t 0 + δ k when P i is released at time t = t 0 and ADP is kept bound on the CI of the ith subunit of kth KaiC hexamer for the time duration δ k . KaiC hexamer is prevented from disassembling into monomers when ATP is bound on the CI 5,6 , which suggests q(i; k, t) < q 0 when ATP or ADP + P i is bound. Here, we use a simple assumption that q(i; k, t) = 0 when ATP or ADP + P i is bound or no nucleotide is bound on the CI of the ith subunit. In this model, the P i release is assumed to take place randomly with a frequency f 0 at individual subunits, and the lifetime of the ADP-bound state δ k in a subunit of the kth hexamer is assumed to depend on the structure as With the definition of Eq. 10, the ATPase activity does not directly depend on the phosphorylation level or the protein binding state. However, through the structure modulation X k , the ATPase activity is indirectly affected by the phosphorylation level and the protein binding state. The ADP release is more frequent with shorter δ k in the gs state with larger X k ; therefore, ATP is hydrolyzed more frequently during the P-phase as observed in the experiment 16 .
Parameters were chosen so as to make the reactions slow enough; the time scales of P/dP reactions, − k p 1 and − k dp 1 , and the time scales of ATPase reactions, − f 0 1 and δ 0 , were chosen to be of the order of 0.5-1 hour. This estimated order of the slow rates is consistent with the large activation free-energy for those reactions suggested by the structural analyses 17 . The rate of KaiA binding, which is ≈h A A T 2 , was chosen to be fast with the order of minutes to show the strong affinity when X ≈ 1, but the rate of KaiB binding, ≈h B B T , was chosen to be slow with time scale of hours to achieve balance between KaiA binding on a single site of the CII ring and KaiB binding on multiple sites of the CI ring; the slowness of KaiB binding is consistent with the experimental observation 22,35 . Within this estimated order of parameter range, parameters were calibrated to show a stable coherent oscillation of the ensemble average of D,

Figure 2a
shows an example trajectory obtained by numerically solving Eqs 4-10 for an ensemble of N = 1000 KaiC hexamers. Here, the ensemble averages D t ( ) and are plotted as functions of time, together with the ensemble average of probabilities . We found that P C A 6 2 and P CBA oscillate in opposite phases, showing that there is a competition between two types of complexes C 6 A 2 and C 6 B i A 2j ; when C 6 A 2 dominates, X decreases and D increases with the action of KaiA on KaiC hexamers, and when C 6 B i A 2j dominates, there is not sufficient KaiA to form C 6 A 2 , which decreases D and increases X. In the latter case, KaiA is sequestered in C 6 B i A 2j , which is the mechanism of synchronization of many KaiC hexamers and is necessary for the ensemble-level oscillation. This effect is shown in Fig. 2b, where the amplitude of D oscillation, ∆D, is plotted. ∆D is large when the KaiA concentration is larger than a threshold of A T /C 6T ≈ 1.4, but with too large KaiA concentration, the sequestration mechanism does not work and the synchronization is lost to make ∆ ≈ D 0. On the other hand, the upper limit of the KaiB concentration for the proper oscillation was not found within the reasonable range because the sequestration mechanism is not working for KaiB. Existence of such concentration ranges of KaiA and KaiB agrees with the experimental observation 36 and is similar to the results of the previous Monte Carlo-type simulation 20 in which sequestration of KaiA at a different oscillatory phase from the present model was assumed.
Effects of ATPase reactions are analyzed by changing q 0 , δ 0 , and f 0 , which are the parameters to define the ATPase reactions in the model. Figure 3 shows how the simulated trajectory changes when the degree of influence of ATPase reaction to the structure, q 0 , is changed. When q 0 is too small, the ATPase reactions are insufficient for Scientific REPORTS | 7: 17469 | DOI:10.1038/s41598-017-17717-z perturbing the system from the steady state with large X and large D. However, when q 0 is too large, the ATPase reactions stabilize the (large X , small D) state and prevents oscillation. This is summarized in Fig. 4a by plotting the oscillation amplitude, ∆D, which becomes large when q 0 exceeds a certain threshold; therefore, ∆ ≠ D 0 only within a limited range of q 0 . Also plotted in Fig. 4a is the oscillation period τ, which is defined by τ = 1/f p with f p being the frequency of the peak of Fourier transformed spectrum of D t ( ); τ is a decreasing function of q 0 , showing that the small and quick oscillation around the (large X , small D ) state dominates when q 0 is large. Figure 4b shows that ∆ ≠ D 0 when the lifetime of ADP bound state, δ 0 , exceeds a certain threshold. For the large enough δ 0 with δ  f 1 0 0 , the CI of each subunit almost always binds ADP and the system behavior shows saturation. τ shows a complex behavior as δ 0 is varied because the form of oscillation as a function of time is varied by changing δ 0 . The similar threshold and saturation behaviors are found in Fig. 4c when the frequency of ATP hydrolysis, f 0 , is varied.
Thus, the steady state turns into oscillation when the quantities that define the ATPase reactions in the model, q 0 , δ 0 , and f 0 , are larger than certain thresholds. Therefore, in the present model, ATP hydrolysis, which is randomly processed in individual CI domains, is a driving mechanism of the ensemble-level oscillation.
Notably, the oscillation period is modulated by changes in q 0 , δ 0 , and f 0 . Among them, the frequency of ATP hydrolysis, or the ATPase activity of KaiC, is most sensitively dependent on f 0 . In Fig. 4d, the averaged number of released ADP from each CI in 24 h is calculated as a measure of the ATPase activity in the non-oscillatory condition of A T = B T = 0 for various values of f 0 , and the thus calculated ATPase activity is compared with the frequency f p of phosphorylation rhythm calculated in the oscillatory condition. We find a clear correlation between ATPase activity and f p , and the slope of the line fitted to the results in Fig. 4d is 0.036. Because it is plausible that f 0 should be modulated by mutations, the calculated slope is consistent with the observed slope of 0.03-0.04 17 and 0.055 16 obtained from various mutant data. Thus, the perturbation of structure by the ATPase reactions in individual molecules is a key determinant of the oscillation frequency in the ensemble of many molecules.
As temperature is increased, reactions should be accelerated, which should increase f 0 and decrease δ 0 . It is suggestive that in Fig. 4b,c, the oscillation period is shortened when f 0 is increased and prolonged when δ 0 is decreased, so that the period is not changed largely when f 0 δ 0 is kept constant as shown in Fig. 5. Note that τ only insensitively depends on f 0 and δ 0 when f 0 and δ 0 are large enough as shown in Fig. 4; therefore, the oscillation period becomes insensitive to temperature change when f 0 or δ 0 is regulated to be large enough even in the case f 0 δ 0 is not kept constant. Further careful analyses are needed to elucidate the temperature dependence of f 0 and δ 0 , but the results in the present model showed that the control of oscillation period with ATPase reactions may work as a mechanism for the observed insensitivity of the oscillation period to temperature, i.e., temperature compensation, together with the ensemble-level mechanism of regulation for temperature compensation 21 .

Discussion
It is interesting to see how the stochastic fluctuation of oscillation is regulated by the ATP hydrolysis reactions. In the present model, due to the stochastic timing of ATP hydrolysis reactions in each CI domain of individual KaiC hexamers, the simulated oscillation of individual KaiC hexamers, D k , shows stochastic fluctuation even in the case    (Fig. 6a). When f 0 is too small, on the other hand, synchronization among individual oscillations becomes weak and the amplitude of the ensemble average oscillation D becomes small (Fig. 6b). Therefore, the ATP hydrolysis reactions are necessary for synchronizing individual stochastic oscillations to give rise to a coherent ensemble-level oscillation.
The present analyses showed that individual stochastic oscillations are synchronized through sequestration of KaiA into C 6 B i A 2j , and the synchronization promoted by this sequestration is enhanced by the frequent ATPase reactions; the frequent structural modulation induced by the ATPase reactions should perturb individual molecules, and this perturbation is necessary for individual hexamers to adjust to the ensemble-level oscillation. It is necessary to further analyze this mechanism quantitatively by comparing the degree of synchronization and the free energy flow induced by the ATP consumption. Theoretical description with combined use of landscape and probability-flux representations 23,37 should provide a useful means for this analysis.
We developed a simplified model of the reconstituted circadian clock. In particular, it was assumed that the binding affinity of KaiA and KaiB to KaiC depends on the structure X k (t) of KaiC hexamer as defined in Eq. 7. In contrast to the assumption of phosphorylation-level dependent binding in many of the previous theoretical models 8,[19][20][21]23,25,28 , the direct explicit dependence of the binding affinity on the phosphorylation level D k (t) of KaiC was not considered in the present model. Instead, it was assumed that the phosphorylation level D k (t) only indirectly affects the binding affinity through the modulation of the phosphorylation-level dependent structure X k (t) as shown in Eqs 7 and 9. This assumption should be reasonable because the atomic positions of the phosphorylation sites, Ser431 and Thr432, are distant from the binding sites of KaiA 11 and KaiB 13 , so that the allosteric communication through the change in X k (t) should mediate the effects of P/dP reactions and binding of KaiA or KaiB. In a similar way, the ATPase activity was assumed to depend on X k (t) as shown in Eq. 10 by neglecting the direct dependence of the ATPase activity on D k (t); the ATPase activity is only indirectly affected by D k (t) through Eqs 7 and 10. Thus, in the present model, different types of reactions, P/dP reactions, ATP hydrolysis, and Kai protein binding reactions, induce the allosteric structural transition, and the allosteric structural transition regulates these reactions, which constitutes multifold feedback interactions among reactions and structural transition. In the present paper, we focused on the intra-molecular dynamics induced by these multifold feedback interactions and analyzed the correlation between intra-molecular dynamics and the inter-molecular dynamics which gives rise to the ensemble-level oscillation. The model has rooms for improvement in explaining the further detailed features of the system. For instance, two sites of phosphorylation should be distinguished by using variables for individual phosphorylation sites, and structures of CI and CII in each subunit should be described by using structural variables for individual domains. Temperature dependence of each parameter should be carefully examined to analyze the mechanism of temperature compensation. Extension of the present model to these directions is an important avenue for further research.
The Kai oscillator provides a clear example that the intra-molecular regulation through the feedback relations among reactions and structure and the inter-molecular system-level regulation through concentration change are coupled with each other to realize an integrated behavior of stable oscillation. Physical perspective obtained by the analyses of this coupling should give further insights on the multilevel aspects of this important system.

Methods
Kinetic equations for the Kai protein binding. For the system having N hexamers of KaiC, stochastic dynamics of the system should be described by the master equation, which is the kinetic differential equation for P(n 1 , i 1 , j 1 ; n 2 , i 2 , j 2 ; …; n N , i N , j N ), where n k = i k = j k = 0 when KaiA or KaiB is not bound to the kth KaiC and n k = 1 and i k = j k = 0 when a KaiA dimer binds to the CII ring of the kth KaiC. i k is the number of bound KaiB on the CI ring of the kth KaiC and j k is the number of bound KaiA dimers on KaiB. Eqs 4 and 5 can be derived by using the Hartree-like approximation 29 as … =∏ = P n i j n i j n i j P n i j ( , , ; , , ; ; , , ) ( , , ) N N N k N k k k 1 1 1 2 2 2 1 and writing = = P k t P n ( , ) ( 1, 0, 0) k C A 6 2 , and . We note that Eq. 5 is the expression for the case of 2 ≤ i ≤ 5 and 1 ≤ j ≤ i − 1. For the other values of i or j, the corresponding equations are   where x i (k) is an Ising-like variavle with x i (k) = 1 when the ith subunit is in the gs state and x i (k) = −1 when it is in the cs state. When the structure transition is fast enough compared to the other reactions, variables D k , p k A , p k B , and q k can be regarded as static fields acting on x i (k). J is the strength of a coupling between neighboring subunits, representing the cooperativity of the allosteric structural transition. In the case that  By defining X k as X k = (〈x(k)〉 + 1)/2 and regarding D k , p k A , p k B , and q k as slowly changing time-dependent variables, Eq. 18 becomes Eq. 9. It is an intriguing issue to examine whether the free energy representation of Eq. 17 is valid by using atomistic molecular dynamics simulation.