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.

Introduction

How is the coherent oscillation in biological system generated and how is it regulated by structures and reactions of constituent molecules? The system comprised of three cyanobacterial proteins, KaiA, KaiB, and KaiC, provides a unique opportunity to examine these problems[1-3]. As illustrated in Fig. 1, KaiC consists of tandemly repeated N-terminal (CI) and C-terminal (CII) domains[4] and assembles into a hexamer[5,6] by forming the CI-CII double rings[7]. By incubating KaiA, KaiB, and KaiC in the presence of ATP in vitro, two specific sites in the CII of each subunit of hexamer are repeatedly phosphorylated and dephosphorylated to exhibit a stable circadian oscillation[1,8,9].

Figure 1

Scheme of interactions among Kai proteins. KaiC assembles into a hexamer C6 by forming CI-CII double rings. The CI ring (green) binds to KaiB (pink) to form C6B with 1 ≤ i ≤ 6 when the CII is phosphorylated (P), which further binds to KaiA dimers (blue) to form C6BA2 with 1 ≤ j ≤ i [13]. The CII ring (yellow) binds to a KaiA dimer when dephosphorylated (dP) to form C6A2. Two specific sites, Ser431 and Thr432, in the CII of each subunit are phosphorylated (filled circle) or dephosphorylated (white circle). The CI domain exhibits slow ATPase reactions.

KaiA forms a dimer[10] and binds to the CII ring of KaiC hexamer[11] to promote the autophosphorylation of KaiC[12]. We write this KaiCA complex as C6A2. KaiB binds to the CI domain of each KaiC subunit to form C6B6 [13]. Binding of KaiA to the CII is suppressed in C6B6, which promotes autodephosphorylation of KaiC[14,15]. Recent electron microscopic and mass spectrometric analyses[13] showed that KaiB further binds to KaiA to form KaiCBA complex, C6B6A2 with 1 ≤ j ≤ 6. When KaiA concentration is not too high, participation of KaiA in forming C6B6A2 should sequestrate KaiA to suppress the C6A2 formation in the other KaiC hexamers, which may provide ensemble-level communication among KaiC hexamers. ATP molecules can bind on CI-CI and CII-CII interfaces of KaiC subunits; with slow ATPase activity, about 10 molecules of ATP are hydrolyzed per each CI domain and several ATP molecules per each CII domain in 24 hours[16,17].

A striking feature is the strong correlation between the ATPase activity and the phosphorylation rhythm: When the KaiC sequence is mutated, the induced modulation of frequency of ATP hydrolysis observed in the absence of KaiA and KaiB, which is the non-oscillatory condition, is correlated to the modulation of frequency of the phosphorylation rhythm observed in the oscillatory condition[16,17]. The ATPase activity of the truncated CI ring also shows correlation to the oscillation frequency of the phosphorylation rhythm[17]. These observations suggested the intrinsic relation between ATP hydrolysis in the CI and phosphorylation rhythm in the CII[16,18]. In the present paper, we develop a theoretical model to examine the mechanism behind this relation.

Many theoretical models have been proposed to analyze the Kai system[19-27], but most of them have been devoted to understanding the macroscopic ensemble-level oscillation; by assuming cyclic phosphorylation/dephosphorylation (P/dP) of individual KaiC molecules, various nonlinear mechanisms were proposed to explain how individual oscillations are synchronized to give rise to an ensemble-level oscillation. Example mechanisms were regulation through KaiA sequestration at various phases[8,19-23] or through monomer shuffling[24-27]. However, the discovery of the correlation between the ATPase activity and the phosphorylation rhythm has indicated that not only the ensemble-level regulation but the regulation within individual KaiC molecules through the interplay between the ATP hydrolysis and P/dP reactions should be important for determining the oscillation features. A 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, C6B0A0 is identical to C6. Probability of the kth C6 to be bound in forms of C6A2 and C6BA2 with 0 ≤ i ≤ 6 and 0 ≤ j ≤ i at time t are denoted by and , 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, and . The similar factorization approximation was used to describe single-molecular fluctuation of binding reactions to DNA[29]. Then, equations for and become where A and B are concentrations of unbound free KaiA and KaiB, respectively. See Methods section for more details. In Eq. 5, the rate constants of binding, h and h , and the rates of unbinding, f and f , should depend on the structure of individual KaiC hexamers, but we assume that h or f does not depend on the KaiC structure because KaiA in C6BA2 binds to KaiB and does not directly binds to KaiC. In solving Eqs 4 and 5, we have the constraints,where V is the system volume, and A and B 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  ≈ 1 when the kth KaiC hexamer is in the gs-state and X  ≈ 0 when in the cs-state. Then, the rate constants in Eqs 4 and 5 should depend on the structural order parameter X aswhere h , f , h , and f 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 (t) ≈ 1 when 12 sites of the kth KaiC hexamer are fully phosphorylated, and D (t) ≈ 0 when they are fully dephosphorylated. D is increased when KaiA binds to KaiC to form C6A2 and is decreased otherwise. This transition is represented by the equation of soft-spin dynamics aswhere k and k 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 is represented as a quasi-equilibrium average under the mean-field generated by D and other quantities aswhere β = 1/k B T, and c 0, c 1, c 2, and c 3 are constants. and represent the extent of KaiA binding and that of KaiB binding, respectively, where n represents the level of saturation of the effects of KaiB binding. See Methods section. q 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 (t) and D (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 , 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 , small D ) state for and the (large X , large D ) state for . On the other hand, when c 1 > 0 in Eq. 9, the large D leads the KaiC structure to the cs-state with small X . From Eq. 7, this brings about the small affinity of KaiA to KaiC, reducing in Eq. 4, which in turn reduces D in Eq. 8. Therefore, the c 1 > 0 term gives a negative feedback effect, which should stabilize the state with (large X , small D ) or (small X , large D ). 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 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 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 with q(i; k, t) = q 0 > 0 for t 0 ≤ t ≤ t 0 + δ when Pi 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 δ . 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 + Pi is bound. Here, we use a simple assumption that q(i; k, t) = 0 when ATP or ADP + Pi is bound or no nucleotide is bound on the CI of the ith subunit. In this model, the Pi release is assumed to take place randomly with a frequency f 0 at individual subunits, and the lifetime of the ADP-bound state δ in a subunit of the kth hexamer is assumed to depend on the structure aswhere γ ≥ 1 is a constant. 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 , the ATPase activity is indirectly affected by the phosphorylation level and the protein binding state. The ADP release is more frequent with shorter δ in the gs state with larger X ; 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, and , and the time scales of ATPase reactions, 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 ≈, 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 , 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, .

Results

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 and are plotted as functions of time, together with the ensemble average of probabilities and . We found that and oscillate in opposite phases, showing that there is a competition between two types of complexes C6A2 and C6BA2; when C6A2 dominates, X decreases and D increases with the action of KaiA on KaiC hexamers, and when C6BA2 dominates, there is not sufficient KaiA to form C6A2, which decreases D and increases X. In the latter case, KaiA is sequestered in C6BA2, 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 oscillation, , is plotted. is large when the KaiA concentration is larger than a threshold of A /C 6 ≈ 1.4, but with too large KaiA concentration, the sequestration mechanism does not work and the synchronization is lost to make . 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.

Figure 2

Simulated oscillation of the Kai system. (a) The simulated trajectory is shown by plotting averages over the ensemble of KaiC molecules as functions of time; the level of phosphorylation (red), structure (black, real line), probability to form the C6A2 complex, (dashed line), and probability to form the C6BA2 complexes, (dotted line) are plotted with A /C 6 = 2 and B /C 6 = 20, where C 6 = N/V is the concentration of KaiC hexamer. (b) The amplitude of oscillation, , is plotted on the plane of A /C 6 and B /C 6. Using units of V = 1 and k B T = 1, parameters are N = 1000, a = 1.6 h−1, k  = 0.4 h−1, k  = 0.4 h−1, c 0 = 10, c 1 = 6, c 2 = 0, c 3 = 3, f 0 = 2.4 h−1, q 0 = 0.8, δ 0 = 2.5 h, γ = 1.5, n  = 1, h  = 1.6 × 10−5 h−1, h  = 7.2 × 10−6 h−1, h  = 1.0 × 10−6 h−1,  f  = 0.8 h−1,  f  = 0.32 h−1, and f  = 0.2 h−1.

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 perturbing the system from the steady state with large and large . However, when q 0 is too large, the ATPase reactions stabilize the (large , small ) state and prevents oscillation. This is summarized in Fig. 4a by plotting the oscillation amplitude, , which becomes large when q 0 exceeds a certain threshold; therefore, 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 ; τ is a decreasing function of q 0, showing that the small and quick oscillation around the (large , small ) state dominates when q 0 is large. Figure 4b shows that when the lifetime of ADP bound state, δ 0, exceeds a certain threshold. For the large enough δ 0 with , 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.

Figure 3

Example trajectories calculated with different values of q 0. The ensemble averages of KaiC phosphorylation level (red), KaiC hexamer structure (black), probability to form the C6A2 complex, (green), and probability to form the C6BA2 complexes, (pink) are plotted as functions of time. The value of q 0 is written in units of k B T in each panel. Parameters other than q 0 are the same as explained in the caption of Fig. 2.

Figure 4

Effects of ATPase reactions on the circadian oscillation of the Kai system. (a–c) The oscillation amplitude (red) and the oscillation period τ (blue) are plotted as functions of (a) the amplitude of impact of the ADP bound state on the structure, q 0, (b) the lifetime of the ADP bound state, δ 0, and (c) the frequency of the Pi release, f 0. (d) The correlation between ATPase activity and oscillation frequency is shown by simulating with various values of f 0, where the ATPase activity was calculated in the condition of A  = B  = 0 and the oscillation frequency, f p, was calculated with A /C 6 = 2 and B /C 6 = 20. Parameters used are the same as in Fig. 2 except for those varied in each graph. Each point is the average of the data calculated from 10 trajectories, each having 2 × 103 h length.

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  = B  = 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].

Figure 5

The oscillation period τ is calculated by varying f 0. Change of τ is much smaller when f 0 δ 0 is kept constant as f 0 δ 0 = 4 (blue) than when δ 0 is kept constant as δ 0 = 2.5 h (red). Other parameters are the same as explained in the caption of Fig. 2. Each point is the average of the data calculated from 10 trajectories, each having 2 × 103 h length.

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 , shows stochastic fluctuation even in the case that the ensemble average shows a coherent oscillation. In Fig. 6, the stochastic fluctuation of individual KaiC hexamers is shown by superposing D (t) of several example KaiC hexamers, which were arbitrarily chosen from the N-hexamer ensemble. When the ATPase frequency, f 0, is large enough, individual oscillatory trajectories D (t) are fluctuating around the coherent oscillation of the ensemble average (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 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.

Figure 6

Stochastic fluctuation of oscillation of individual KaiC hexamers. Frequency of the Pi release, f 0, is (a) f 0 = 2.4 h−1 and (b) f 0 = 1.2 h−1. In top panels, simulated oscillations, D (t), of five KaiC hexamers are superposed; these hexamers were arbitrarily chosen from the ensemble of N = 1000 hexamers. Oscillations of different hexamers are shown with different colors. In bottom panels, the simulated oscillation of the ensemble, , is shown. Parameter values other than f 0 are the same as explained in the caption of Fig. 2.

The present analyses showed that individual stochastic oscillations are synchronized through sequestration of KaiA into C6BA2, 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 (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-21,23,25,28], the direct explicit dependence of the binding affinity on the phosphorylation level D (t) of KaiC was not considered in the present model. Instead, it was assumed that the phosphorylation level D (t) only indirectly affects the binding affinity through the modulation of the phosphorylation-level dependent structure X (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 (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 (t) as shown in Eq. 10 by neglecting the direct dependence of the ATPase activity on D (t); the ATPase activity is only indirectly affected by D (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 , i , j ), where n  = i  = j  = 0 when KaiA or KaiB is not bound to the kth KaiC and n  = 1 and i  = j  = 0 when a KaiA dimer binds to the CII ring of the kth KaiC. i is the number of bound KaiB on the CI ring of the kth KaiC and j is the number of bound KaiA dimers on KaiB. Eqs 4 and 5 can be derived by using the Hartree-like approximation[29] as and writing , 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

Mean-field equation for the structural transition

The mean-field relation of Eq. 9 is derived by assuming that free energy G of the kth KaiC hexamer is represented as a function of the structural order parameter of the ith subunit, x (k), aswhere x (k) is an Ising-like variavle with x (k) = 1 when the ith subunit is in the gs state and x (k) = −1 when it is in the cs state. When the structure transition is fast enough compared to the other reactions, variables D , , , and q can be regarded as static fields acting on x (k). J is the strength of a coupling between neighboring subunits, representing the cooperativity of the allosteric structural transition. In the case that , we can regard x (k) does not depend on i as x (k) = x(k). Then, by regarding as an Ising-spin Hamiltonian, the quasi-equilibrium average of x(k) leads to

By defining X as X  = (〈x(k)〉 + 1)/2 and regarding D , , , and q 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.