Single-molecular and ensemble-level oscillations of cyanobacterial circadian clock.

When three cyanobacterial proteins, KaiA, KaiB, and KaiC, are incubated with ATP in vitro, the phosphorylation level of KaiC hexamers shows stable oscillation with approximately 24 h period. In order to understand this KaiABC clockwork, we need to analyze both the macroscopic synchronization of a large number of KaiC hexamers and the microscopic reactions and structural changes in individual KaiC molecules. In the present paper, we explain two coarse-grained theoretical models, the many-molecule (MM) model and the single-molecule (SM) model, to bridge the gap between macroscopic and microscopic understandings. In the simulation results with these models, ATP hydrolysis in the CI domain of KaiC hexamers drives oscillation of individual KaiC hexamers and the ATP hydrolysis is necessary for synchronizing oscillations of a large number of KaiC hexamers. Sensitive temperature dependence of the lifetime of the ADP bound state in the CI domain makes the oscillation period temperature insensitive. ATPase activity is correlated to the frequency of phosphorylation oscillation in the single molecule of KaiC hexamer, which should be the origin of the observed ensemble-level correlation between the ATPase activity and the frequency of phosphorylation oscillation. Thus, the simulation results with the MM and SM models suggest that ATP hydrolysis stochastically occurring in each CI domain of individual KaiC hexamers is a key process for oscillatory behaviors of the ensemble of many KaiC hexamers.

Cyanobacterial proteins, KaiA, KaiB, and KaiC, can reconstitute a circadian clock when they are incubated with ATP in vitro. In order to understand this prototypical oscillator, we need to analyze both synchronization of a macroscopically large number of oscillating molecules and microscopic reactions in individual molecules. We introduced two theoretical models to unify macroscopic and microscopic viewpoints. Simulation results suggest that ATP hydrolysis is necessary for synchronization and temperature compensation and that ATPase activity is correlated to the oscillation frequency in individual molecules. Thus, ATP hydrolysis stochastically occurring in individual molecules should determine important features of the ensemble-level oscillation.

A cyanobacterial protein, KaiC, shows a stable oscillation in its phosphorylation level with approximately 24 h period when three proteins, KaiA, KaiB, and KaiC, are incubated with ATP in vitro [1]. Much attention has been focused on this prototypical oscillator and close investigations have been made on interactions among Kai proteins [2,3]. As shown in Figure 1, KaiC monomer consists of tandemly repeated N-terminal (CI) and C-terminal (CII) domains [4] and assembles into a hexamer [5,6] forming the CI–CII double rings [7]. We denote this KaiC hexamer as C6. KaiA forms a dimer [8] and binds to the CII ring of C6 to form a complex, C6A2 [9]. KaiB binds to the CI of each subunit of KaiC hexamer, and a KaiA dimer can further bind to each KaiB to form C6B6A2 with 0 ≤ j ≤ 6 [10]. KaiC shows auto-phosphorylation (P) and auto-dephosphorylation (dP) at two specific sites, Ser431 and Thr432, of each CII domain, and the P reactions are promoted in the C6A2 complex [11-13]. In C6B6A2, binding of KaiA to the CII is suppressed, which promotes the dP reactions of KaiC [14,15]. Thus, the temporal appearance of different kinds of complexes, accompanied by the switching between P-phase and dP-phase, is the important feature of the oscillation in the KaiABC system.

Figure 1

Coarse-grained models of the KaiABC system. a) Schematic of interactions among Kai proteins considered in the MM model. Two phosphorylation sites, Ser431 and Thr432, of each subunit in the CII domain (blue) are repeatedly phosphorylated (pS, pT) and dephosphorylated (S, T). KaiA dimer (red) binds to KaiC hexamer C6 to form C6A2. The CI (pink) of each subunit binds to KaiB (cyan) to form C6B with 1 ≤ i ≤ 6, which further binds to KaiA dimers to form C6BA2 with 1 ≤ j ≤ i. b) Reactions and states in a KaiC subunit in the SM model. ATP hydrolysis reactions and KaiB binding/dissociation reactions occur in the CI domain (pink). Phosphorylation (P)/dephosphorylation (dP) reactions and KaiA binding/dissociation reactions take place in the CII domain (blue). CI and CII domains are coupled through allosteric structural change.

In order to elucidate the mechanism of this KaiABC clockwork, we need to understand both the microscopic atomic-scale dynamics of reactions in individual molecules [16] and the macroscopic ensemble-level synchronization among many KaiC molecules; because the ensemble-level oscillation vanishes when individual molecules oscillate independently of each other, synchronization of a macroscopically large number of molecules is necessary for maintaining the coherent oscillation as observed in test tube. Therefore, it is important to understand mechanisms at both microscopic and macroscopic levels.

A key observation for solving the problem of the relationship between mechanisms at two levels is the slow ATPase reactions in KaiC. The CI domain of each KaiC subunit hydrolyzes approximately 10 ATP molecules and the CII domain of each subunit hydrolyzes several ATP molecules per one day [16,17]. In particular, by comparing KaiC molecules of various mutants, it was shown that the ATPase activity is correlated to the frequency of the ensemble-level P/dP rhythm [16,17], where the ATPase activity was measured by the frequency of ADP release in the condition in the absence of KaiA and KaiB, which is the non-oscillatory condition, and the frequency of the ensemble-level P/dP rhythm was measured in the oscillatory condition with appropriate concentrations of KaiA and KaiB. Furthermore, the ATPase activity of the truncated CI domain shows correlation to the frequency of the ensemble-level P/dP rhythm [16]. These observations suggested the intrinsic relationship between ATP hydrolysis occurring in the CI of individual molecules and the P/dP rhythm of the ensemble of many molecules [16-18].

Many theoretical models were developed by describing concentrations of various molecular species with continuous variables in kinetic differential equations to model the ensemble-level macroscopic dynamics of the oscillating system [11,19-24]. With these macroscopic models, various nonlinear mechanisms, including KaiA sequestration at various phases [11,19-23] and monomer shuffling [24], were proposed to explain how the oscillation is maintained. Other models described the system with the Monte Carlo methods by explicitly considering stochastic P/dP reactions in individual molecules and stochastic binding reactions among many molecules [23-27]. These Monte Carlo-type models should be useful for analyzing the relation between microscopic and macroscopic mechanisms. However, the ATP hydrolysis reactions and the activated intramolecular dynamics with ATP consumption were not considered in these models. The role of ATP consumption was analyzed by recent models [28,29], and the further comprehensive analyses are required for elucidating the relationship between ATP hydrolysis and the P/dP rhythm [29].

Here, in the present paper, we explain two “mesoscopic” models of oscillation for bridging the microscopic and macroscopic- level descriptions; the many-molecule (MM) model, and the single-molecule (SM) model. In both models, structure and reactions in individual KaiC molecules are represented by coarse-grained variables. For example, structure of the kth KaiC hexamer is represented by an order parameter, X (X in the SM model by dropping the suffix k); structural transition is described by the change between states with X ≈ 0 and X ≈ 1. In particular, the role of ATP hydrolysis is investigated with these models by considering the effects of ATPase reactions on the structural change; X is perturbed by the ATP hydrolysis in these models, and binding/dissociation reactions, P/dP reactions, and ATPase reactions are affected by X. With the MM model in the present paper, intermolecular interactions are analyzed to examine the necessary conditions for the stable ensemble-level oscillation. This model was introduced in our previous publication [29], and is used in the present paper to analyze the mechanism of synchronization and temperature compensation. In the SM model, more emphasis is placed on the intramolecular dynamics to analyze the origin of oscillation of individual molecules. In the present paper, we discuss the linkage between microscopic and macroscopic oscillations through the combined analyses with the MM and SM models.

Methods

The many-molecule (MM) model

In order to analyze the oscillation of ensemble of molecules, we use the many-molecule (MM) model introduced in our previous publication [29]. In the MM model, we describe structure and reactions of individual molecules in a coarse-grained manner and calculate the coupled dynamics of N = 1000 hexamers of KaiC by explicitly following the dynamics of individual hexamers. As represented in Figure 1a, we assume that a KaiA dimer binds to a KaiC hexamer to form C6A2. KaiB monomer binds to each subunit of KaiC to form C6B with 1 ≤ i ≤ 6, and KaiA can further bind to C6B to form C6BA2 with 1 ≤ j ≤ i [10]. We assume that the dimerization reaction 2A ⇌ A2 is fast enough and equilibrated in the solution, and the concentration of A2 is always proportional to the square of the concentration of A. Therefore, all the reactions involving A2 are represented in terms of A for the convenience of model description. Then, the binding/unbinding reactions in the model are represented as

We should note that C6B0A0 in Eq. 1 is identical to C6. We define P(k, t) and P(k, t) as the probabilities of the kth KaiC hexamer C6 to be bound in forms of C6A2 and C6BA2, respectively, with 0 ≤ i ≤ 6 and 0 ≤ j ≤ i at time t and k = 1, …, N. Thus, the kinetic equation for P is given as

and the kinetic equation for P with 2 ≤ i ≤ 5 and 1 ≤ j ≤ i – 1 is

See the Methods section of Ref. 29 for the equations in more general cases with 0 ≤ i ≤ 6 and 0 ≤ j ≤ i. Here, h, f, h, f, h and f are rate constants, and A and B are concentrations of unbound free KaiA and KaiB molecules, respectively. When we numerically solve Eqs. 2 and 3, we impose the following constraints;

where A and B are total concentrations of KaiA and KaiB, respectively, and V is the volume of the system. In this paper, we use units with V = 1.

Biochemical and structural observations [30-33] have suggested that the structure of KaiC hexamer at the phosphorylated state is different from the structure of KaiC hexamer at the unphosphorylated state. The structure at the unphosphorylated state is more stable than the structure at the phosphorylated state [33]. Assuming that a KaiC hexamer has a tight structure at the unphosphorylated state and a loose structure at the phosphorylated state, we can represent the allosteric transition between tight and loose structural states with a continuous variable X as X ≈ 1 when the structure of the kth KaiC hexamer is in the tight structural state, and X ≈ 0 when it is in the loose structural state. In Eqs. 2 and 3, we assume that the binding rate constants, h and h, and the dissociation rate constants, f and f, depend on the KaiC structure, X(t);

where , and are constants. With this assumption, the binding affinity of KaiA to KaiC increases when the kth KaiC is in the tight structural state. On the contrary, the loose structure of KaiC promotes KaiB binding to the KaiC hexamer. Because KaiA binds to KaiB in making C6BiA2, and does not bind directly to KaiC, we assume that the binding rate constant, h, and the dissociation rate constant, f, are independent of X(t).

In the MM model, the phosphorylation level of twelve sites (Ser431 and Thr432 in six subunits) of a single KaiC hexamer is represented by a continuous variable D as D(t) ≈ 1 when the 12 sites of the kth KaiC hexamer are fully phosphorylated, and D(t) ≈ 0 when the 12 sites are fully dephosphorylated. We describe the transition between D(t) ≈ 1 and D(t) ≈ 0 with the following kinetic equation,

Here, k and k are rate constants of P and dP reactions, respectively, and ɡ(D) = aD(D – 1)(D – 1/2)2 with a > 0. ɡ(D) is a double-well potential function, which has been used to describe two-state transition dynamics of spins or model neurons with kinetic differential equations [34], and here it was introduced to assume the bistability of the oscillating variable D. We also assume that the structural transition of KaiC between tight and loose states is much faster than other reactions in the system; therefore, the structure X can be represented as a quasi-equilibrium average. In this paper, this average is determined under the mean-field generated by D, the binding probability of KaiA and KaiB, and the effect of ATPase reaction as

where β = 1/kT with k being the Boltzmann constant, and c, c, c, and c3 are constants. represents the level of KaiA binding to the kth KaiC hexamer, and the term is the binding level of KaiB. We use n = 1 in this paper. The term q(t) represents the effect of ATPase reactions.

The term c1D(t) with c1 > 0 in Eq. 7 brings a negative feedback effect in the system. KaiC structure tends to stay in the loose structure with large D and small X, which reduces the binding affinity of KaiA in Eq. 5 and decreases P(k, t) in Eq. 2. As a result, D is reduced in Eq. 6. The term with c3 > 0 represents a positive feedback. Large value, along with small X, brings KaiC structure to the loose state, which in turn increases the binding affinity of KaiB to KaiC through Eq. 5 and increases P (k, t) in Eq. 3. Competition among multiple feedback interactions leads the system to have multiple stationary states as described above. However, the system is easily shifted among the states when small perturbations are added to produce oscillation in the system. The term q(t), representing the effect of ATPase reactions, works as a perturbation.

The precise mechanism of ATPase reactions in KaiC has not been elucidated. In the other ATPases such as myosins [35], kinesins [36], and chaperonins [37], reactions proceed as (i) binding of ATP → (ii) hydrolysis of ATP → (iii) release of inorganic phosphate (Pi) → (iv) release of ADP. In F1-ATPase, on the other hand, it has been claimed that the ADP release takes place prior to the Pi release [38], although this remains controversial [39]. Because KaiC has a P-loop motif, which is similar to the corresponding ATP-binding motifs of myosins, kinesins and chaperonins, we assume in the present paper that the Pi release is prior to the ADP release in KaiC. We assume that ATP is stochastically hydrolyzed with a frequency f0 in each CI domain of subunits of each KaiC hexamer: q(i; k, t) = 0 when no nucleotide is bound or ATP is bound at time t in the ith subunit of the kth KaiC hexamer, and when ATP is hydrolyzed to be ADP and Pi, q(i; k, t) turns to nonzero as q(i; k, t) = q0 ≠ 0. We consider that this switching from 0 to q0 stochastically takes place with the frequency f0 at each subunit. After q (i; k, t) turns to nonzero at time t = t0, the ADP is kept bound for a time duration δ, which is represented by setting q(i; k, t) = q0 for t0 ≤ t ≤ t0 + δ. Then, at t = t0 + δ, the ADP is released and q(i; k, t) is set back to q(i; k, t) = 0. Because the observed ATPase activity is larger in the P-phase [17], we assume that the overall ATPase reaction rate is increased with a shorter ADP binding duration time δ when X is larger as

The overall ATPase reaction rate decreases when δ is larger; with the parameterization in Ref. 29, δ was larger than in the present parameterization, which resulted in the smaller ATPase activity than in the present paper, but the qualitative behavior of the model was the same as the present work. Biochemical and NMR data suggested that the ATPase activity in the CI is necessary for the binding of KaiC to KaiB [31,33,40], which indicated that the ATP hydrolysis enhances the binding affinity of KaiC to KaiB. Therefore, we assume that the ATP hydrolysis changes the KaiC structure from tight to loose states. We represent this effect by imposing q0 > 0 and writing the effects of ATP hydrolysis on the structure as .

We should note that the following alternative assumption is also possible; q(i; k, t) = q0 > 0 when Pi is released from the KaiC subunit but ADP is kept bound and q(i; k, t) = 0 when ATP is bound, ADP and Pi are bound or no nucleotide is bound on the CI [29]. There is no mathematical difference between this assumption and the assumption made in the last paragraph with the present simplified description of the ATPase reactions in the MM model. More precise modeling is necessary for distinguishing these different reaction schemes, but we do not go into these details in the present paper. We also note that the model feature remains unchanged even when the alternative mechanism of ATPase reactions with the ADP release prior to the Pi release is assumed; the mathematical structure of the model is not altered if we assume that q(i; k, t) = q0 when Pi is kept bound.

Because ATP hydrolysis takes place predominantly in the CI domain and P/dP reactions occur in the CII domain, communication between CI and CII is needed to explain the observed effects of ATP hydrolysis on the P/dP reactions. In the present MM model, allosteric structural change represented by the change in X in Eqs. 5–8 enables such communication. The overall ATPase reaction rate is modulated by X through the change of the duration time of ADP bound state δ as in Eq. 8, but the timing of ATP hydrolysis is determined stochastically with the frequency f0 in each subunit of each KaiC hexamer.

Eqs. 2–7 are numerically solved to follow the temporal change of D, X, P(k, t), P(k, t), and by stochastically varying q(i; k, t) with the frequency f0 and the lifetime δ with i = 1, …, 6 and k = 1, …, N. Oscillations in the ensemble level are monitored by calculating

and

Parameters in the MM model

We assume that the P/dP reactions are slow enough with the rates k~0.5–1 h−1 and k~0.5–1 h−1 to realize the period ~24 hour of the P/dP rhythm as was suggested by Ref. [21]. The time scale of ATPase reactions, and δ, should be also as slow as f0~1–2 h−1 and δ~0.5–1 hour; these slow rates of ATPase reactions are consistent with the structural observation [16]. We assume that the binding/dissociation rates of KaiA and KaiB defined as hA2, hA2, hB, f, f, and f should be h−1–min−1; such slow binding/dissociation of KaiB is consistent with the observed data [21,41]. With this rough estimation of parameter values, the ensemble oscillation in the MM model appears to be stable for the wide range of parameter values. We should note that the simulated period of oscillation depends on the parameter choice as discussed in the previous paper [29]. Because different values of period from 22 h to 26 h [1,42,43] were observed in different experimental conditions, we do not calibrate the parameters to reproduce a particular experimental report but use a typical set of parameters to grasp the important features of the oscillation. An example parameter set is summarized in Table 1. Unless otherwise mentioned, we discuss the results calculated with this parameter set in Results and Discussion section.

Table 1

Parameters in the MM model

Rate constant of KaiA binding to KaiC	hA1	6.6 × 10−6 h−1
Rate constant of KaiB binding to KaiC	hB0	2 × 10−6 h−1
Rate constant of KaiA binding to KaiB	hBA	1.8 × 10−6 h−1
Rate constant of KaiA dissociation from KaiC	fA0	6 × 10−1 h−1
Rate constant of KaiB dissociation from KaiC	fB1	4 × 10−1 h−1
Rate constant of KaiA dissociation from KaiB	fBA	2 × 10−1 h−1
Constant for variable confinement	a	2.2 h−1
Rate constant of phosphorylation	kp	4.4 × 10−1 h−1
Rate constant of dephosphorylation	kdp	4.4 × 10−1 h−1
Frequency of ATP hydrolysis	f0	2.2 h−1
Max. lifetime of ADP binding after hydrolysis	δ0	2.6 h
Effect of ATP hydrolysis on structure	q0	kBT0*
Base-line temperature effect on structure	c0	10kBT0*
Effect of phosphorylation on structure	c1	8kBT0*
Effect of KaiA binding on structure	c2	0
Effect of KaiB binding on structure	c3	4kBT0*
Total copy number of KaiA monomers	AT	2N
Total copy number of KaiB monomers	BT	20N
Total copy number of KaiC hexamers	N	1000
System volume	V	1

T0 = 30°C

In Table 1, the total concentrations of KaiA and KaiB, A and B, are set to 2N/V and 20N/V, respectively. In monomer basis, this corresponds to the concentration ratio of KaiA:KaiB:KaiC = 1:10:3. In many in vitro experimental measurements the ratio was set to KaiA:KaiB:KaiC = 1:3:3 [13,17,33], and the effects of varying these concentrations were tested by Nakajima et al. [43]. They showed that the oscillation amplitude and period were not changed when the KaiB concentration was varied between KaiA:KaiB:KaiC = 1:3:3 and 1:15:3 [43]. A and B used in Table 1 fall in this range. Values in Table 1 are shown with units of V = 1. Alternatively, by using units of V = 3 × 10−15 l, concentration of KaiC can be converted to 6N/V ≈ 3.3 μM with N = 1000, which is close to 3.5 μM, the KaiC concentration often used in in vitro experiment. A typical value of KaiB binding rate constant, which is with units of V = 1, can be converted to with units of V = 3 × 10−15 l. This value should be directly compared with the observed data when the experimental data of KaiB monomer kinetics becomes available.

The single-molecule (SM) model

In the MM model, the structure and reactions of individual molecules were described in a simplified manner to facilitate the calculation of an ensemble of a large number of molecules. In order to further analyze the role of ATPase reactions, we introduce a single-molecule model of KaiC hexamer by describing the more details of intra-molecular reactions and structural changes.

In the SM model, each subunit of a hexamer follows two conformational states; x = 0 when the ith subunit is in the loose structural state and x = 1 when it is in the tight structural state with i = 1, …, 6. Though both CI and CII domains have nucleotide binding sites, we focus on the CI domain because ATPase activity is predominant in the CI domain. The four nucleotide binding states in the CI of the ith subunit are represented as a = 0 when no nucleotide is bound, a = 1 when ATP is bound, a = 2 when ADP and Pi are bound, and a = 3 when Pi is released and ADP is kept bound. As shown in Figurer 1b, we define variables, Ser and Thr, as Ser = 1 (Thr = 1) when Ser431 (Thr432) of the ith subunit is phosphorylated and Ser = 0 (Thr = 0) when Ser431 (Thr432) of the ith subunit is dephosphorylated. We write KaiA = 1 when a KaiA dimer binds to the CII ring and KaiA = 0 when a KaiA dimer dissociates from the CII ring. We write KaiB = 1 when a KaiB monomer binds to the ith subunit and KaiB = 0 when a KaiB monomer dissociates from the ith subunit. Thus, as shown in Figure 1b, the state of KaiC subunits in the SM model is described by the variable KaiA and others;

Because the atomic structure is rapidly changing with the time scale of milliseconds or shorter but other variables should change more slowly in minutes or longer time scales, we describe the free energy G of the KaiC hexamer as a function of {x} as G = G (x1, x2, …, x6);

where b0, , b2, b3, b4, and J are constants. In the last term of Eq. 10, x7 = x1. If we regard the fast varying quantities {x} as “spin variables”, G of Eq. 10 should correspond to the “effective Hamiltonian” for spins, whereas the slow variables, Thr, Ser, KaiA, KaiB and θ(a) correspond to the magnetic fields acting on the spins. In other words, phosphorylation states, binding states, and ATPase states give biases for the structural change of each KaiC subunit. J represents coupling between two neighbor subunits; J > 0 implies that subunits change their structures in a cooperative way to give rise to two-state allosteric structure transition of the KaiC hexamer. The term xb4θ(a) represents coupling between ATPase reactions and structure. Here, we assume that θ(a) = −1 when a = 2 or 3, θ(a) = 1 when a = 1, and θ(a) = 0 when a = 0. These assumptions imply that ATP binding increases the tendency of KaiC hexamer structure to be tight while the binding of ADP + Pi or ADP loosens the structure.

For the time scale of minutes, we can represent the structural state of each KaiC subunit with a quasi-equilibrium thermal average of x as 〈x〉= Tr x exp(−βG)/Z by approximately treating Thr, Ser, KaiA, KaiB, and θ(a) as stationary variables, where Tr represents a sum over all possible values of {x} and Z is a partition function, Z = Tr exp(−βG). Then, we represent the hexamer state by using the following variables;

where X and D correspond to X and D in the MM model.

By comparing Eqs. 7 and 10, we find that b0, and , b2, b3, and b4 in the SM model correspond to c0, c1, c2, c3, and q0 in the MM model, respectively. Thus, we find multiple feedback relations working through Eq.10 in the SM model. When b2 > 0, binding of KaiA brings KaiC structure to the tight state as it decreases the energy G, which increases the binding affinity of KaiA to KaiC and enhances phosphorylation. When b3 > 0, binding of KaiB brings the structure of KaiC into the loose state. This effect is consistent with the mass spectrometry observation [44] and enhances dephosphorylation, which further increases the binding affinity of KaiB. Therefore, b2 > 0 and b3 > 0 give positive feedback effects in the system. On the other hand, when , phosphorylation of Ser431 stabilizes the loose structure of the KaiC subunit by lowering the energy G and this reduces the binding affinity of KaiA and stimulates dephosphorylation of Ser431. Hence, the term provides a negative feedback effect in the system. We assume so as to make the phosphorylation of Thr432 faster than the phosphorylation of Ser431 as observed in experiment [12,13]. In the present paper, we assume . Therefore, the term gives a positive feedback effect. As a result of the coexistence of these multiple feedback relations, the system has multiple metastable states including the phosphorylated loose-structure state and the dephosphorylated tight-structure state. This is why the system tends to stay one of these multiple states and makes transitions among them. ATPase reactions add perturbation to the structure and stimulate the transition between them. The term xb4θ(a) in Eq. 10 represents such perturbation through the ATPase reactions.

The binding/dissociation reactions of KaiA occur in the CII ring of KaiC hexamer. The reactions are defined as

We consider that the binding affinity of KaiA is regulated by the allosteric structure change of KaiC hexamer; binding of KaiA is enhanced when the hexamer structure is tight whereas binding of KaiA is suppressed when the structure is loose. Therefore, we assume that the rates in Eq. 12 are regulated by X as

A constant should depend on the concentration of free unbound KaiA molecules. Here in the present SM model, we use a simplified assumption that the free KaiA concentration is constant, so that is kept constant. This assumption represents the situation that the KaiA concentration is much larger than the KaiC concentration, so that the free KaiA concentration is not affected by the oscillation of individual KaiC molecules. In the MM model, as will be shown in Results and Discussion section, sequestration of KaiA to C6BA2 and the resultant temporal depletion of free KaiA concentration is the mechanism of synchronization among multiple KaiC hexamers; synchronization is lost when the free KaiA concentration is kept constant. Because synchronization is not discussed with the SM model, we here adopt a simple assumption that is kept constant.

The binding/dissociation reactions of KaiB in the CI domain of KaiC hexamer are described as

We assume that the rates in Eq. 14 are regulated by the structure of subunits 〈x〉 in the following manner;

In this way, the binding affinity of KaiB to the ith subunit of KaiC is decreased when the structure of the subunit is in the tight structure, and the KaiB binding affinity is large when the subunit is in the loose structure. Again, we use a simplified assumption that is kept constant through the simulation time course.

The P/dP reactions are described by

The rates in Eq. 16 are regulated by the binding status of KaiA as

Here, , and are dynamical variables, while, , and are constant parameters. We assume the following six equations for the change of nucleotide binding state;

Here, both the forward and backward reactions are considered for binding/unbinding of ATP and ATP hydrolysis. Binding reaction of Pi or ADP is not considered by assuming that their concentrations are small. The rates in Eq. 18 are assumed to depend on the structure of each subunit 〈x〉 as in the following manner;

Here, , and are dynamical variables, while , and are constant parameters.

We perform numerical simulations of the stochastic dynamics with SM model: By using 〈x〉 and X calculated with Eq. 10, the rates of reactions are calculated with Eqs. 12–19. KaiA, KaiB, a, Ser, and Thr are updated with these rates by the Gillespie algorithm [45]. Then, by using thus updated KaiA, KaiB, a, Ser, and Thr, the structural variables 〈x〉 and X in the new time step is calculated with Eq. 10.

Parameters in the SM model

Parameters in the SM model were set to have the same order of magnitude as the corresponding parameters in the MM model. In the SM model, more detailed reaction schemes are considered, so that the larger number of parameters are defined than in the MM model. For example, the P/dP reactions are considered at two sites, Ser431 and Thr432, which are regulated by the KaiA binding status. Therefore, the rates of the P/dP reactions are defined by two sets of parameters, those in the KaiA bound state, , and those in the KaiA unbound state, . Here, superscripts 1 and 0 imply the rates in the KaiA bound and unbound states, respectively. We consider that the phosphorylation rates are larger in the KaiA bound state and the dephosphorylation rates are larger in the KaiA unbound state. Therefore, and . Because the P/dP reactions at Thr432 precedes the P/dP reactions at Ser431, we consider and . We assume ATP hydrolysis rates depend on the structure of each subunit; the reaction rates are larger when the structure is tight. Therefore, we assume , and are larger than , and . Here, superscripts 1 and 0 indicate the rates in the tight and loose structures, respectively. An example set of parameters of rate constants and parameters in the free energy to determine the structure is summarized in Tables 2 and 3.

Table 2

Rate constants in the SM model*

ATP hydrolysis reactions in the tight structure	k01a,1	1.1
	k10a,1	1.0
	k12a,1	1.1
	k21a,1	1.0
	k23a,1	1.1
	k30a,1	1.1
ATP hydrolysis reactions in the loose structure	k01a,0	0.55
	k10a,0	0.5
	k12a,0	0.55
	k21a,0	0.5
	k23a,0	0.55
	k30a,0	0.55
P/dP reactions with KaiA unbound from KaiC	k0→1Thr,0	0.3
	k1→0Thr,0	3.9
	k0→1Ser,0	0.01
	k1→0Ser,0	0.5
P/dP reactions with KaiA bound to KaiC	k0→1Thr,0	4.8
	k1→0Thr,1	0.2
	k0→1Ser,1	0.6
	k1→0Ser,1	0.1
Binding and dissociation rates of KaiA	ha1	1.5
	fa0	0.1
Binding and dissociation rates of KaiB	hb0	1.5
	fb1	0.2

Values are in units of h−1.

Table 3

Free energy parameters in the SM model*

Base-line temperature effect	b0	0.13
Effect of Thr432 phosphorylation	b1Thr	−1
Effect of Ser431 phosphorylation	b1Ser	6
Effect of KaiA binding	b2	1.5
Effect of KaiB binding	b3	1
Effect of ATP hydrolysis	b4	4
Coupling between subunits	J	4.5

Values are in units of kT0 with T0 = 30°C.

Results and Discussion

Oscillation and synchronization in the ensemble of KaiABC molecules

In Figure 2a, we show an example of the ensemble-level oscillation simulated with the MM model. Shown are the KaiC phosphorylation level, D̄ (t), structure of KaiC hexamers, X̄(t), probability of C6 of forming the C6A2 complex, P̄ , and probability of forming the C6BA2 complexes, P̄ (t), where overbars represent that averages were taken over the ensemble of N hexamers. Here, we find that P̄ (t) and P̄ (t) show counter-phased oscillation. When P̄ (t) dominates, binding of KaiA to KaiC hexamer enhances phosphorylation, which increases D̄ and decreases X̄. On the contrary, when P̄ (t) is large, KaiA is sequestered into C6BA2, which depletes the free unbound KaiA. This depletion of free KaiA concentration enhances dephosphorylation, which then decreases D̄ and increases X̄. This sequestration of KaiA is the mechanism through which individual KaiC hexamers communicate with each other.

Figure 2

Ensemble and single molecule oscillations in the KaiABC system simulated with the MM model. a) The ensemble oscillation of the phosphorylation level D̄ (t) (red line), structure of KaiC hexamer X̄(t) (black real line), probability to form complex C6A2, P̄ (dashed line), and probability to form complexes C6BA2, P̄ (dotted line), are plotted as functions of time. Averaged over the ensemble of N = 1000 hexamers. b) The single molecule oscillation of phosphorylation level D(t) (red line), structure X(t) (black real line), P(k, t) (dashed line), P(k, t) (dotted line) of a KaiC hexamer arbitrarily chosen from the ensemble are plotted as functions of time. Parameters in Table 1 were used.

Individual oscillation of single KaiC hexamer is shown in Figure 2b, where oscillation trajectories of D(t), X(t), P(k, t) and P(k, t) of a KaiC hexamer arbitrarily chosen from the ensemble are plotted. When X(t) is large, the lifetime of bound ADP, δ, is short and ADP is frequently released from the CI. This short residence of ADP at CI does not affect X(t) much. However, when X(t) begins to be small, δ becomes longer and ATP hydrolysis provides larger influence to X(t). This gives rise to spike-like changes in X(t), which induces noisy fluctuating decrease of D(t). Synchronization of a large number of KaiC hexamers smooths out this fluctuation to give rise to a coherent oscillation of D̄ .

To analyze this synchronization furthermore, we plot individual oscillations D(t) of five KaiC hexamers arbitrarily chosen from N hexamers along with the ensemble-averaged oscillation D̄ (t). In Figure 3a, individual and ensemble oscillations calculated with the same parameter values as used in Figure 2 are shown. Fluctuating individual oscillations are synchronized and entrained into the ensemble oscillation, which gives rise to a stable coherent oscillation of the ensemble. The essential role of KaiA sequestration in this synchronization is clearly shown in Figure 3b, in which free KaiA concentration is fixed to be A = 0.2A. When the free KaiA concentration is fixed to be constant as in Figure 3b, the sequestration mechanism does not work and individual KaiC hexamers lose synchronization among them. Then, the ensemble-level oscillation vanishes though individual molecules remain oscillating with the large amplitude. Thus, the temporal change of the binding affinity of KaiA to KaiC due to the dynamical change of free KaiA concentration is essential for maintaining synchronization. In contrast, as shown in Figure 3c, synchronization and the ensemble oscillation are maintained when the free KaiB concentration is fixed to be B = 0.9B. This maintenance of synchronization shows that the temporal change in the concentration of free KaiB is not necessary for synchronization.

Figure 3

Oscillations of individual molecules in the ensemble and the ensemble-averaged oscillation simulated with the MM model. The phosphorylation level D(t) (blue) of five individual KaiC hexamers arbitrarily chosen from the ensemble are superposed with the ensemble average D̄ (t) (red). a) Parameters in Table 1 were used with f0 = 2.2 h−1. b) Concentration of free KaiA molecules is fixed to be A = 0.2A, but concentration of free KaiB molecules remains dynamically changing. c) Concentration of free KaiA molecules remains dynamically changing, but concentration of free KaiB molecules is fixed to be B = 0.9B. Frequency of ATP hydrolysis was made small as d) f0 = 0.8 h−1 and e) f0 = 0.5 h−1.

Because ATP hydrolysis plays a significant role to generate individual oscillations as shown in Figure 2b, it is interesting to see how the ATP hydrolysis affects the synchronization. In Figure 3d, the frequency of ATP hydrolysis, f0, was made smaller from f0 = 2.2 h−1 in Figure 3a to f0 = 0.8 h−1 in Figure 3d. We find synchronization is lost with this small frequency of ATP hydrolysis. Therefore, ATP hydrolysis with sufficient frequency is necessary for synchronization. ATP hydrolysis occurring with random timing in individual KaiC subunits should perturb each KaiC hexamer to adjust to the dynamical oscillation of free KaiA concentration; without this adjustment, individual KaiC hexamers stay oscillating with their own individual phases. As shown in Figure 3e, by further reducing the ATP hydrolysis frequency to f0 = 0.5 h−1 many KaiC hexamers stop oscillating and show small amplitude fluctuation. Therefore, when f0 is decreased, KaiC oscillators first lose synchronization, and then lose individual oscillations.

Together, these results showed that in the MM model KaiA sequestration is a basis for synchronization of individual KaiC hexamers, and ATP hydrolysis helps individual KaiC hexamers to synchronize under the mechanism of KaiA sequestration.

Temperature compensation

An important feature of the KaiABC system is its insensitivity to temperature change, i.e., temperature compensation; the observed Q10 defined by the relative change in the oscillation period upon 10°C temperature change was about 1.06 [1] to 1.01±0.02 [46]. However, the physical mechanism of such temperature insensitivity has been elusive. Here, we analyze this problem by using the MM model, particularly by investigating the role of ATP hydrolysis in determining temperature dependence of the oscillation period and amplitude.

The MM model has nine parameters for rate constants, , h, , f , a, k, and k, and two parameters for the ATP hydrolysis reactions, f0 and δ0. We write their values at temperature T around T0 = 30°C as

We compare various cases for the values of eleven activation energies, ΔE, ΔE, ΔE, ΔE, ΔE, ΔE, ΔE, ΔE, ΔE, ΔE, and ΔE;

Here, case I is the simplest assumption to make all the relevant activation energies same. Case II is the assumption of the stronger temperature dependence of the binding affinity of KaiA to KaiB. Hatakeyama and Kaneko [20] discussed that enhancement of the KaiA sequestration at higher temperature reduces the phosphorylation rate, which should lead to the temperature compensation. In the present MM model, this enhancement of KaiA sequestration is represented by the larger ΔE than others. As shown in the previous paper [29], the important features of the MM model are the dependence of the oscillation period on the frequency of the ATP hydrolysis f0 and on the lifetime of the ADP binding state δ0. In case III, the stronger temperature dependence of δ0 is assumed with the larger ΔE0; as temperature increases, the rate of ADP release is enhanced and δ0 becomes small with the larger temperature dependence than other reactions. The frequency of ATP hydrolysis was shown to be temperature insensitive [17]. A possible explanation for this temperature insensitivity is the regulation of the reaction rate by the diffusive non-activation type binding process of ATP to the CI, but the precise mechanism of this temperature insensitivity has not been known. Here, to represent this temperature insensitivity, we assume ΔE 0 = 0 in case IV. Cases V and VI are calculated to test the effects of assuming ΔE 0 = 0.

In Figures 4a and 4b, the calculated temperature dependences of period and amplitude of the ensemble-averaged oscillation D̄ (t) are shown for cases I–VI. As shown in Figure 4a, there are only small differences between case I and case II and between case V and case VI; therefore, the enhancement of KaiA sequestration at higher temperature contributes only slightly for temperature compensation in the present model. On the other hand, period calculated in cases III and IV is much more insensitive to temperature than period calculated in cases I and II, showing that the enhancement of shortening of the lifetime of the ADP-bound state at higher temperature gives a significant effect on temperature compensation. Temperature insensitivity of ATP hydrolysis frequency gives some effect to make period temperature insensitive as shown in the differences between case III and case IV and between II and VI in Figure 4a. Q10 of the oscillation period calculated between 30°C and 40°C is Q10 = 1.17 (case I), 1.16 (case II), 1.07 (case III), 1.02 (case IV), 1.13 (case V), and 1.10 (case VI).

Figure 4

Temperature dependence of period and amplitude of the ensemble oscillation D̄ (t) simulated with the MM model. a) Period τ and b) amplitude ΔD̄ . Six cases defined in Eq. 21 are compared: case I (red), case II (green), case III (blue), case IV (magenta), case V (cyan), and case VI (purple). Parameters in Table 1 were used.

As shown in Figure 4b, amplitude of oscillation becomes only slightly small as temperature decreases in cases I , II, V, and VI, but the amplitude approaches rapidly to 0 as temperature is decreased in cases III and IV. In cases III and IV, the oscillation amplitude continuously decreases and vanishes at a certain threshold temperature; in other words, the oscillation disappears with the Hopf bifurcation mechanism in low temperature, which is the behavior consistent with the experimental observation [42].

Summarizing these results, we found that the accelerated release of ADP from the CI with increased temperature largely contributes to temperature compensation of oscillation period, and insensitivity of ATP hydrolysis frequency to temperature gives some help to temperature compensation. The prolonged release of ADP at low temperature is consistent with the observed disappearance [42] of oscillation through the Hopf bifurcation at low temperature. Meanwhile, enhancement of KaiA sequestration does not help much for temperature compensation.

Oscillation in individual KaiC hexamer

In the preceding subsections, roles of ATP hydrolysis in the KaiABC oscillation were analyzed with the MM model. Though ATPase reactions randomly occur in individual molecules, they have significant effects on the simulated ensemble-level oscillation. In the following subsections, we use the SM model for further investigating the roles of ATP hydrolysis.

In Figure 5, we plot the simulated single-molecular phosphorylation levels, D, Ser, and Thr, and structure X; these quantities are defined in Eq. 11 in Method section. As shown in Figure 5a, the phosphorylation level and structure oscillate with opposite phases in a similar way to the oscillation in the MM model. From the simulated oscillation of Ser and Thr in Figure 5b, we find that the phase of Thr precedes the phase of Ser. In this way, the simulated single-molecular oscillation is noisy stochastic oscillation, but it captures the important features observed in the ensemble-level oscillation. Therefore, the physical origin of features of the ensemble-level oscillation exists at the single-molecular level.

Figure 5

Oscillations of single KaiC hexamer simulated with the SM model. a) Simulated oscillation of the phosphorylation level D(t) (red) and structure X(t) (black). b) Simulated oscillation of phosphorylation at Ser432 (red) and Thr431 (blue). Parameters in Tables 2 and 3 were used.

The mechanism of this single-molecular oscillation can be more closely analyzed by changing the parameters, , b2, b3, and b4 in Eq. 10, which represent the strength of coupling between reactions and structural change. In Figure 6, we plot the oscillation amplitude, δD, by introducing scaling factors α1, α2, and α3 to scale the parameters as (Fig. 6a), , b2, , α2b2, α2b3 (Fig. 6b), and b4 → α3b4 (Fig. 6c). Here, δD = 2 × (standard deviation of the distribution of D), where the values of D were sampled from the trajectory D (t) calculated over 224 steps. As discussed in Method section, α1 modulates the strength of negative feedback interaction between P/dP reactions and structure, α2 modulates the strength of positive feedback interactions between P/dP reactions and structure or between binding and structure, and α3 modulates the effect of the ATPase reactions on structure. As shown in Figure 6a, δD drops to a small value when α1 becomes small. This appreciable decrease of δD can be confirmed by examining the calculated trajectories as shown in Figure 6d. In Figure 6d, oscillations of D and X are masked by intense fluctuations when the negative feedback is lost with α1 = 0, while the oscillations remain distinct when the negative feedback is strong enough with α1 = 3. The oscillation features with α1 = 3 are similar to those found in Figure 5a, which were calculated with α1 = 1. This weakening of oscillation for the small α1 explains the behavior of δD in Figure 6a. As shown in Figure 6b, δD mildly decreases when α2 becomes small. Trajectories in Figure 6e show that the oscillations of D and X become somewhat noisy when the positive feedback is lost with α2 = 0, while the oscillations turn to trains of pulses when the positive feedback is strong with α2 = 3. As shown in Figure 6c, δD shows a peak at α3 ≈ 0.5 and diminishes both at small and large α3. As shown in Figure 6f, the system tends to stay long at the small X and small D state when α3 = 0, which prevents oscillations, while the system tends to stay at large D with frequent noisy oscillation of X when large α3 = 3. Therefore, there is a suitable range of the coupling strength between ATPase reactions and structure to keep the oscillation amplitude large. To summarize, we find that the large amplitude oscillation in individual KaiC hexamers is supported by negative and positive feedback interactions in structure change as well as the suitable strength of coupling between ATPase reactions and structure change.

Figure 6

Effects of variation of the coupling strength between reactions and structure in the SM model. a) Oscillation amplitude δD is plotted as a function of the scaling factor α1 of coupling between P/dP reactions at Ser431 and structure. α2 = α3 = 1. b) δD as a function of the scaling factor α2 of coupling between P/dP reactions at Thr432 or Kai protein binding reactions and structure. α1 = α3 = 1. c) δD as a function of the scaling factor α3 of coupling between ATPase reactions and structure. α1 = α2 = 1. Example trajectories of D (red) and X (black) calculated with d) α1 = 0 (upper panel) and α1 = 3 (lower panel) with α2 = α3 = 1, e) α2 = 0 (upper panel) and α2 = 3 (lower panel) with α1 = α3 = 1, and f) α3 = 0 (upper panel) and α3 = 3 (lower panel) with α1 = α2 = 1.

Correlation between ATPase activity and oscillation frequency

The observed correlation between the frequency of ATP hydrolysis, or the ATPase activity, and the frequency of KaiC oscillation [16,17] is important for analyzing the relationship between the ATPase reactions occurring in individual KaiC hexamer and the ensemble-level P/dP rhythm. As discussed in previous subsections, many features of the ensemble-level oscillation are already found in the single-molecular level. Therefore, it is meaningful to analyze the relation between the ATPase activity and the frequency of the single-molecular oscillation.

In Figure 7a, the single-molecular ATPase activity was calculated with the SM model by introducing a scaling factor ɛ; the rate constants of ATPase reactions defined in Eq. 19 were uniformly changed from the values of Table 2 as , and with i = 1 and 0. With these scaled parameters with various values of ɛ, the number of released ADP molecules (i.e., the number of a = 3 → 0 transitions) was counted as a measure of the ATPase activity in the condition in the absence of KaiA and KaiB, which was represented in the model by imposing . Also shown in Figure 7a is the averaged peak frequency of the power spectrum of the Fourier transform of the D oscillation calculated with nonzero and defined in Table 2 for the corresponding values of ɛ. Because oscillations bear intense stochasticity, the averaged peak frequency and its error bar were derived from extensive calculations: 50 power spectra were calculated from the Fourier transform of 50 trajectories of D(t), with each trajectory having 218 data points. Then, the peak of the spectrum obtained by averaging 50 spectra was identified. This calculation was repeated 100 times by using different random numbers, and from the distribution of 100 peak frequencies, the averaged peak frequency was obtained and the error bar was calculated from the standard deviation. In spite of the intense stochasticity as seen from the large error bars in Figure 7a, we can see that the ATPase activity is correlated to the P/dP oscillation frequency in single KaiC hexamer.

Figure 7

Correlation between ATPase activity and frequency of phosphorylation oscillation. a) Single-molecular correlation calculated with the SM model. Each point represents the ATPase activity calculated in the non-oscillatory condition and the average frequency calculated in the oscillatory condition for various scaling values ɛ of the ATPase reaction rates: ɛ = 4 (blue), ɛ = 2 (red), ɛ = 1 (black), ɛ = 0.6 (magenta), and ɛ = 0.4 (cyan). Error bars are the standard deviation of the distribution of the peak frequency of the Fourier transform spectra. b) Ensemble-level correlation calculated with the MM model. Each point represents the ATPase activity calculated in the non-oscillatory condition and the average frequency calculated in the oscillatory condition for various values of f0. Each point is the average of 10 trajectories, each of which is 2 × 103 h long.

In Figure 7b, the results of the ensemble calculation with the MM model were plotted. Here, the number of released ADP in the ensemble were calculated in the condition of A = B = 0 with various values of ATP hydrolysis frequency f0, and the corresponding ensemble-level P/dP rhythm was calculated with the nonzero values of A and B as defined in Table 1. Figure 7b shows that there is a distinct correlation between ATPase activity and the P/dP rhythm in the ensemble.

Comparing the results of Figures 7a and 7b, we find that the important feature of the ensemble-level correlation between the ATPase activity and the P/dP rhythm already appears at the single-molecule level.

Conclusion

We constructed two models, the MM model, which explains oscillation dynamics of an ensemble of a large number of Kai molecules and the dynamics of constituent individual molecules, and the SM model, which explains oscillation dynamics of individual molecules in more details. In these two models, we have hypothesized that KaiC hexamers undergo allosteric transitions between tight and loose structural states. The binding affinity of KaiA and KaiB was assumed to depend on the structure of KaiC hexamer and the rates of P/dP reactions were defined to depend on the binding status of KaiA. The binding of KaiB to KaiC leads to the KaiABC complex formation, which sequesters KaiA to deplete the free KaiA concentration and affects the phosphorylation rate of the other KaiC molecules. ATPase activity of the CI domain of KaiC was assumed to be regulated by the structure of KaiC, and the ATP hydrolysis reactions were assumed to bring about the perturbation on the allosteric transition of the KaiC structure.

Simulations with the MM model showed that stochastic individual oscillations are synchronized through the communication among KaiC hexamers, which was realized by the sequestration of KaiA into the KaiABC complex, and the frequent ATPase reactions in the CI domain perturb individual KaiC hexamers to allow them to generate the ensemble-level oscillation. Thus, both KaiA sequestration and frequent ATPase reactions perform pivotal roles to maintain a coherent ensemble-level oscillation. Temperature dependence of the period and amplitude of the ensemble-level oscillation were examined with the MM model, and it was shown that the intense temperature dependence of the lifetime of the ADP bound state in the CI domain is important to realize the temperature compensation of the oscillation period, and such dependence is consistent with the observed disappearance of oscillation at low temperature [42]. It was also shown that temperature insensitivity of the frequency of ATP hydrolysis helps temperature compensation to some extent.

Simulations with the SM model showed that the many features of the ensemble-level oscillation shown with the MM model appear also in the single-molecular oscillation in the SM model though the single-molecular oscillation bears intense stochastic fluctuation. Results in the SM model showed that perturbation of the ATP hydrolysis reactions on the structure is necessary to maintain the large amplitude single-molecular oscillation and that ATPase activity and the P/dP oscillation frequency are correlated to each other at single-molecule level, which largely explains the correlation between ATPase activity and the P/dP oscillation frequency found in the ensemble level.

Thus, both MM and SM models suggest that ATP hydrolysis reactions stochastically occurring in individual KaiC hexamers play important roles in synchronization, temperature compensation, and frequency determination at the ensemble-level oscillation. In order to verify these suggestions, the present “mesoscopic” models should be further compared with the microscopic atomic-level observations of structure and reactions [16] and also with the macroscopic ensemble-level observations of phase shift, synchronization, and entrainment [46,47]. In this way, a comprehensive picture that unifies explanations of microscopic through macroscopic phenomena should open a new perspective of biomolecular systems working coherently in living cell systems.