Abstract

Cyanobacteria are the simplest known cellular systems that regulate their biological activities in daily cycles. For the cyanobacterium Synechococcus elongatus, it has been shown by in vitro and in vivo experiments that the basic circadian timing process is based on rhythmic phosphorylation of KaiC hexamers. Despite the excellent experimental work, a full systems level understanding of the in vitro clock is still lacking. In this work, we provide a mathematical approach to scan different hypothetical mechanisms for the primary circadian oscillator, starting from experimentally established molecular properties of the clock proteins. Although optimised for highest performance, only one of the in silico‐generated reaction networks was able to reproduce the experimentally found high amplitude and robustness against perturbations. In this reaction network, a negative feedback synchronises the phosphorylation level of the individual hexamers and has indeed been realised in S. elongatus by KaiA sequestration as confirmed by experiments.

Visual Overview

Synopsis

In most eukaryotes, an internal clock drives numerous activities into daily cycles—a circadian clock, which is ticking with a periodicity of about 24 h. The period of this free‐running rhythm is highly robust against many changes in the natural environment, for example, in cyanobacteria the clock can compensate for variations in the ambient temperature. But for certain external stimuli (e.g. light, temperature, nutrients), the circadian rhythm is able to be entrained. The optimal temporal coordination of biological processes and the adaptation to daily fluctuations play a critical role in the survival of diverse organisms. Photoautotrophic organisms like plants and cyanobacteria are subjected to a daily light–dark rhythm during their photosynthetic activity and have been demonstrated to posses a free‐running circadian clock as well. In particular, for the cyanobacterium Synechococcus elongatus, a robust circadian rhythm has been observed under constant darkness conditions and even for complete suppression of the cellular transcription and translation activity (Tomita et al, 2005). Moreover, only three different cyanobacterial proteins (KaiA, KaiB, and KaiC) are sufficient to achieve a temperature‐compensated circadian rhythm of phosphorylation cycles in vitro (Nakajima et al, 2005). Despite numerous excellent experimental investigations, it remained unclear how a biochemical mechanism involving just three proteins can keep its timing so precisely over a long period, as biochemical events are known to be intrinsically stochastic.

We used the in vitro clock of S. elongatus as an instructive model system for the circadian mechanism in a unicellular organism. The phosphorylation–dephosphorylation cycle of KaiC functions as a circadian oscillator by mixing the three Kai proteins and ATP in a test tube. Thereby, KaiA and KaiB are modulating KaiC phosphorylation state by forming complexes (Kageyama et al, 2003; Garces et al, 2004) with so far unknown stoichiometry. To identify the protein complex formation during a 24 h cycle, we monitored the molecular composition and weight of complexes by 2D gel separation experiments. Consistent with the recent observations, we demonstrated that KaiC forms stable hexamers. After incubation, the hexamers start to autophosphorylate, a process that is enhanced by KaiA. At maximum KaiC phosphorylation, KaiB dimers begin to associate and from stable complexes with KaiC, initiating the dephosphorylation of KaiC. Interestingly, at low phosphorylation levels of KaiC, KaiA is found in a complex with KaiC and KaiB. However, the exact molecular mechanism generating temperature‐compensated 24‐h oscillations is yet unknown.

A novel mathematical approach was provided to generate different hypothetical mechanisms for a basic circadian oscillator. Starting from the biochemically defined reaction network known for the clock proteins and verified by our experiments, we introduced systematically different feedback interactions connecting different states of the KaiC phosphorylation cycle to achieve oscillatory behaviour (Tyson et al, 2003). We obtained 224 possible network topologies using a single feedback loop. A global optimisation of reaction constants for high amplitude was employed for each generated network topology based on experimental observations that in vivo phosphorylation is able to oscillate close to the maximum and minimum phosphorylation levels as found in vitro (Nishiwaki et al, 2004; Tomita et al, 2005). For the feedback mechanisms allowing oscillations, stability in phase and frequency was analysed as it has been shown that cyanobacterial cells possess a stable phase over several generations even under constant low‐light conditions (Mihalcescu et al, 2004).

Finally, only one of the in silico‐generated reaction networks (Figure 1) was able to reproduce the experimentally found high amplitude of oscillations and robustness against changes resulting from cell division, protein synthesis, and degradation. In this reaction network (Figure 1C), a negative feedback, which can be realised by KaiA sequestration to low‐phosphorylated KaiBC complexes (Figure 2), synchronises the phosphorylation level of the individual KaiC hexamers. Thus, our theoretical analysis suggests that the main oscillatory mechanism is a consequence of KaiA sequestration. Intriguingly, the experimentally observed behaviour can be simulated by our mathematical approach without the exact knowledge of the biochemical reaction constants. Moreover, our optimisation procedure revealed that for maximal oscillations of the found network only a small amount of KaiC hexamers is needed to form complexes with KaiA and KaiB and the remaining large fraction of phosphorylated KaiC never undergoes the complete cycle of our reaction network and is available as free hexamers at all times (Figure 2).

In our approach, the implementation of a systematic rewiring algorithm to identify oscillatory network topologies and a subsequent global optimisation procedure for high amplitudes to identify appropriate reaction constants roughly simulate a process comparable with biological evolutionary selection of a circadian clock. Also, there exists now strong experimental evidence that the theoretically found core mechanism of KaiA sequestration is indeed realised in S. elongatus (Kageyama et al, 2006). Thus, the cyanobacterial clock might have evolved by an optimisation process selecting for a minimal circadian pacemaker of high amplitude and robustness.

The cyanobacterial clock is designed to compensate for correlated variations of its components

Introduction

The coordination of biological activities into daily cycles provides an important advantage for the fitness of diverse organisms, from bacteria to humans (Ouyang et al, 1998; Johnson and Golden, 1999; Sharma, 2003; Johnson, 2004; Woelfle et al, 2004). Central to this coordination is an internal clock that drives gene expression in an approximate 24 h rhythm—a circadian clock found in most eukaryotes and, among prokaryotes, exclusively in cyanobacteria.

The cyanobacterium Synechococcus elongatus PCC 7942 is used as an instructive model system for circadian mechanisms in unicellular organisms (Iwasaki and Kondo, 2004). Three proteins KaiA, KaiB, KaiC, which share no sequence similarity to any known eukaryotic clock component, form the core oscillator of the internal clock in cyanobacteria. These three proteins together with ATP are sufficient to generate temperature‐compensated circadian oscillations of KaiC phosphorylation in a test tube (Nakajima et al, 2005). Thus, in contrast to eukaryotic clock models, the cyanobacterial core oscillator operates independently of transcription and translation processes (Tomita et al, 2005; Woelfle and Johnson, 2006). Detailed functional and structural studies on Kai proteins revealed the following main properties: (i) KaiC monomers bind ATP to form a stable hexamer (Nishiwaki et al, 2000; Mori et al, 2002; Hayashi et al, 2003, 2006); (ii) KaiC has both autophosphorylation and dephosphorylation activities; in the absence of KaiA, about 20% of KaiC is phosphorylated (Nishiwaki et al, 2000; Tomita et al, 2005); (iii) a KaiC monomer possesses two main phosphorylation sites (T432 and S341) (Xu et al, 2004), resulting in a total of 12 phosphorylation sites for each KaiC hexamer; (iv) KaiA forms a dimer and enhances autophosphorylation activity of KaiC (Kitayama et al, 2003; Xu et al, 2003; Uzumaki et al, 2004; Vakonakis and LiWang, 2004; Ye et al, 2004; Pattanayek et al, 2006); (v) KaiB attenuates KaiA‐enhanced phosphorylation of KaiC (Kitayama et al, 2003); and (vi) the three Kai proteins can form stable complexes with yet unknown stoichiometry during the subjective night (Kageyama et al, 2003; Garces et al, 2004). It has been further shown experimentally that there exists no intercellular communication that leads to synchronisation between individual cells (Mihalcescu et al, 2004) and temperature compensation of the circadian clock can be readily seen in the in vitro assay as a consequence of temperature‐independent phosphorylation kinetics of the KaiC hexamers (Tomita et al, 2005). It is still unclear if rhythmic phosphorylation of KaiC is the only oscillatory mechanism in cyanobacteria or if a second oscillator, based on transcriptional–translational feedbacks, plays an essential role in the in vivo system. There is some indication for this conjecture as it has been found that KaiA enhances transcription of the kaiBC operon whereas KaiC suppresses its own synthesis (Iwasaki et al, 2002; Nair et al, 2002; Kutsuna et al, 2005).

Although there exist detailed analyses of structure, binding affinities, and phosphorylation kinetics of Kai proteins, no consistent mathematical model has been proposed so far that assigns functional roles to the molecular mechanisms involved in the circadian phosphorylation process. Despite that just three proteins are involved in the in vitro assay, there is enough room for attributing hypothetical complexes of various stoichiometry effects on the phosphorylation kinetics that lead to sustained oscillatory behaviour of the KaiC phosphorylation level.

In a previous theoretical approach, monomer exchange was proposed as a possible mechanism to synchronise KaiA‐enhanced phosphorylation of the hexamers (Emberly and Wingreen, 2006). Fully phosphorylated hexamers should then undergo a structural change that allows for the creation of higher order clusters of KaiC hexamers. Within these clusters, enhanced phosphorylation activity by KaiA is inhibited. The clusters dephosphorylate consequently and are assumed to decay to free hexamers whenever the phosphorylation level of a cluster falls below a critical value. A second approach (Mehra et al, 2006) attributes phosphorylated KaiC an autocatalytic activity that induces a positive feedback on the KaiC phosphorylation process and results in synchronisation at the highest phosphorylation level of the hexamers. Phosphorylated KaiC binds KaiB and KaiA, and dephosphorylates under release of KaiA and KaiB. Finally, in a third approach (Kurosawa et al, 2006), it was conjectured that KaiB can exist in an active and inactive state. Phosphorylated KaiC enhances deactivation of KaiB and active KaiB enhances dephosphorylation of KaiC. This two‐fold inhibition results in an effective positive feedback, similar to the one described for the second approach. To arrive at strong oscillatory behaviour, the authors assume a switch‐like deactivation of KaiB by phosphorylated KaiC.

The first two approaches are in strong contrast to recent experiments (Kageyama et al, 2006) and the experimental results shown in this work. Furthermore, there exists so far no experimental support for the hypothesis that either KaiA or KaiB exists in more than one configurational state. However, it is possible to construct other hypothetical reaction schemes that are in agreement with experimental data but represent different biochemical mechanisms leading to oscillatory behaviour.

To track down the most likely clock mechanism, we propose a novel approach that allows to scan potential oscillatory clock mechanisms starting from the experimentally established molecular interactions among the Kai proteins. A necessary requirement for a reaction network to show oscillatory behaviour is the existence of at least one feedback loop (Tyson et al, 2003). We therefore generated different feedback loops that connect the different phosphorylation and binding states of KaiC in a systematic way. In our approach, we scanned potential oscillatory mechanisms that involve one feedback loop only. For each generated reaction system, a global optimisation of reaction constants was employed to find sustained oscillations of high amplitude. The reason to optimise for high amplitude arises from in vivo experiments. As shown in Tomita et al (2005), the relative amount of KaiC with at least one phosphorylated residue ranges between ∼20 and 80% under constant darkness condition. Thus KaiC phosphorylation oscillates close to the minimal and maximal achievable phosphorylation levels as observed in in vitro experiments (20–90%) (Nishiwaki et al, 2004). In the absence of light, synthesis of the Kai proteins is suppressed and their cellular abundance remains constant over time. Under these conditions, we expect the in vivo circadian clock to show similar behaviour as in the in vitro experiments. By scanning different feedback mechanisms, we arrive at a set of clock mechanisms that can be ranked with respect to their stability in phase and frequency. We find that the most robust oscillatory network structure uses a negative feedback that has been realised in cyanobacteria by KaiA sequestration to low‐phosphorylated KaiBC complexes. The underlying assumption within this optimisation approach is that the cyanobacterial core clock is the result of strong evolutionary selection for robust oscillations against environmental perturbations, such as temperature changes, gene expression noise, and dilution of protein levels by cell divisions. Latter has been evidenced by in vivo experiments, where cells showed phase stability over several generations under constant low‐light conditions (Mihalcescu et al, 2004).

Results

Formation of KaiABC complexes in vitro

To identify the assembly of the Kai protein complexes during a 24 h cycle, we monitored their molecular weight using Blue Native PAGE (BN‐PAGE) separation. The electrophoretic mobility of KaiC in the first dimension, at a position that corresponds to approximately 360 kDa, implies that it forms stable hexamers (theoretical weight of a KaiC hexamer: 348 kDa). KaiA and KaiB are found in a lower range of the BN‐PAGE (Figure 1A). The migration behaviour in the BN‐PAGE revealed that KaiB might be in its tetrameric form whereas KaiA forms a stable dimer, consistent with recent observations (Hayashi et al, 2004; Hitomi et al, 2005; Iwase et al, 2005). Six to eight hours after incubation of the three proteins, an additional complex at around 480 kDa appears that reaches its highest intensity after 10–12 h. The second dimension gel shows (Figure 1D) that this larger complex consists of KaiC and KaiB. The difference between the two complexes is about 115 kDa. Thus, at least 5–6 KaiB dimers should be associated with one KaiC hexamer. This finding is in contrast with the theoretical approach of Emberly and Wingreen (2006) where the authors assume the existence of stable oligomeric formations of hexamers. The absence of KaiBC complexes within the first 4 h suggests that only highly phosphorylated KaiC hexamers can associate with KaiB. After 12 h, this well‐resolved complex starts to dissociate and KaiC appears to be present also at a position in the gel where low molecular weight components are migrating. This result coincides with the predicted monomer exchange among KaiBC complexes (Kageyama et al, 2006). After 18 h, diffuse‐migrating KaiC‐containing complexes appear at around 550–650 kDa in the gel. The denaturating second dimension revealed that among KaiB and KaiC, KaiA is present, to some extent, in these complexes. It should be noted that the used silver staining of SDS gels is not equally sensitive for different proteins; thus, the intensity of the stained protein spots in the second dimension does not directly reflect the stoichiometry of these proteins in the complex. The lower amount of the KaiA in this largest protein complex population and its diffuse migration suggest that KaiA is only loosely associated with this complex. Note that KaiA is also expected to bind only weakly to the KaiC hexamers and thus BN‐PAGE separation will overestimate the amount of free KaiA. In order to check to what extent KaiC‐containing complexes become phosphorylated over the incubation time, radioactive [γ‐32P]ATP was used in the KaiABC reaction mixture. Figure 1B suggests that at time point 12 h, the highly phosphorylated KaiC hexamers are mainly included in the 480 kDa KaiBC complex. Nevertheless, a small part of the KaiC label present is shifted down to low molecular weight monomers, indicating that monomeric KaiC also exists in phosphorylated form. In addition, we used different stoichiometries between KaiC and KaiB proteins to demonstrate the impact of the KaiB concentration on the formation of KaiBC complexes. Figure 1C clearly shows that the amount of KaiB in the reaction mixture is critical for the formation of 480 kDa complex.

Formation of Kai protein complexes during the in vivo phosphorylation cycle. (A) Proteins (0.2 μg ml−1 KaiA, 0.2 μg ml−1 KaiB, and 0.8 μg ml−1KaiC; corresponding to a molar ratio of about 1:3:3) were incubated under standard conditions for different periods of time. The protein samples were immediately applied to a 4–16% BN‐PAGE gel. Molecular weight marker proteins are indicated on the left. (B) Several standard incubation mixtures were supplemented with [γ‐32P]ATP (10 μCi per reaction mixture), subjected to BN‐PAGE, and the dried gel was autoradiographed. (C) Formation of protein complexes after 12 h incubation of a standard KaiC and KaiA mixture with varying concentrations of KaiB (1 × corresponds to the standard reaction mixture). (D) Lanes of interest (2, 8, and 18 h) were excised from a 1D BN gel, placed on top of a resolving gel for Tricine–SDS–PAGE and subsequently separated by electrophoresis. Arrows indicate newly arising KaiA and KaiB protein spots associated with larger KaiC‐containing complexes.

The core reaction network

The experimentally found dynamics of phosphorylation and complex formation allows us to establish a core reaction network that is schematically drawn in Figure 2. Here, the KaiC hexamers in the ith reaction step, C6(i−1)⇌C6(i), can be phosphorylated and dephosphorylated with rates k+i−1(t) and k−i(t), respectively. These rates depend on the concentrations of free KaiA and KaiB dimers and are therefore time‐dependent. Under physiological conditions, the maximal achievable average phosphorylation level of KaiC is less than 45% for serine and about 80% for threonine phosphorylation sites, resulting in approximately 90% KaiC with at least one residue phosphorylated (Nishiwaki et al, 2004). In our reaction scheme (Figure 1), we therefore take N=6 phosphorylation steps between the lowest (20%) and highest physiological phosphorylation states as a working hypothesis. KaiC hexamers at the highest phosphorylation state, C6(N), are assumed to undergo a conformational change, C6(N)⇌C6*(N). The state C6*(N) is then stabilised by complex formation with m≈6 KaiB dimers, B2 (see discussion of Figure 1),

Core reaction scheme of KaiC phosphorylation. The KaiC hexamers, C6, undergo phosphorylation and dephosphorylation depending on the actual reaction rates. Higher (lower) level of KaiA increases (decreases) the phosphorylation rates. Hexamers in the highest phosphorylation state bind six KaiB dimers to form a stable complex. The [B12C6] complex is assumed to dephosphorylate and releases KaiB dimers in the lowest phosphorylation state. The superscript indicates the number of phosphate groups added, starting from the minimal physiological phosphorylation level.

The KaiBC complex is expected to be stable for high phosphorylation levels and gets increasingly unstable for low phosphorylation levels as confirmed by Figure 1 and by Kageyama et al (2006). In Figure 2, we assume that full dissociation of a complex is only possible for the lowest phosphorylation state, [B2mC6(0)] → C6(0)+mB2, but the complexes can loosely bind additional KaiA and KaiB depending on their phosphorylation level. Our experiments (Figure 1B) indicate that KaiA is somehow deactivated upon interaction with complexes, and thus elevated autophosphorylation activity is inhibited for both free hexamers and KaiBC complexes (Kageyama et al, 2006). We assume further that KaiA and KaiB exist in a single conformational state only.

The core reaction scheme (Figure 2), common to all reaction networks studied in this work, cannot show oscillations of KaiC phosphorylation but approaches a stable fixed point for all initial conditions. This dynamic behaviour is a consequence of the fact that each KaiC hexamer follows its own stochastic phosphorylation cycle and there exists so far no mechanism in our model that leads to sufficient synchronisation between the phosphorylation levels of the individual hexamers. However, synchronisation can be realised biochemically by introducing reactions that result in at least one positive or negative feedback loop.

Systematic scan of reaction networks involving one feedback loop

In our mathematical approach, we describe the dynamics of the reaction network on the level of Michaelis–Menten kinetics, where in the simplest form we can assign KaiA and KaiB the role of enzymes whereas free hexamers and KaiBC complexes play the role of substrates. The total concentrations of KaiA and KaiB appear on this level of description always in product form with their corresponding reaction constants (see Supplementary material). Because latter constants are optimised for highest performance, we can absorb the total concentrations of KaiA and KaiB in the reaction rates. The dynamics of the system then depends exclusively on the distribution of phosphorylation levels among hexamers and KaiBC complexes. The fact that KaiC hexamers and KaiBC complexes can bind KaiA and KaiB with different affinity, depending on their phosphorylation level, results in an effective dependence of kinetic constants on the distribution of phosphorylation states in the system. For example, attributing low‐phosphorylated KaiBC complexes a high affinity for KaiA results in a negative feedback on the phosphorylation kinetics for the free hexamers—the rates ki+(t) or all i∈{0,…,5} decrease with increasing concentration of complexes, [B2mC6(j)] for j∈{0,…,3}, as a consequence of the implicit dependence of ki+(t) on KaiA.

To scan the space of potential oscillatory reaction networks that include just one feedback loop in a systematic way, we developed a ‘rewiring’ scheme by inserting and removing appropriate links in the core reaction scheme (Figure 2). Examples of possible feedback loops are shown in Figure 3. Here, arrows indicate positive feedbacks and bar‐ends negative feedback loops. The feedback in Figure 3C shows the above‐mentioned repression arising from KaiA sequestration by low‐phosphorylated KaiBC complexes.

Three examples out of eight reaction networks that show oscillatory behaviour with high amplitude. In the reaction schemes (left panels), bar ends indicate repressions and arrows enzymatic activations. Feedbacks are allowed to work only between the groups that are indicated by dashed boxes in (A). Each link that targets a group of low‐ or high‐phosphorylated hexamers (C6) or KaiBC complexes ([B12C6]) is assigned one additional reaction constant to be optimised (see Supplementary information). (A) Left panel: mechanism of KaiABC clock as proposed by Emberly and Wingreen (2006). The negative feedback arises from the assumption that hexamers can exchange monomers. The middle graph shows the time course of the KaiC phosphorylation level for protein concentrations taken from in vitro experiments with N=6 (green line) and N=20 (blue line) phosphorylation steps. Changes in phase and frequency owing to a collective five‐fold increase or 10‐fold decrease in protein concentrations are shown by the solid black (only partially drawn) and dashed red line, respectively, for N=20. The right panel shows the fraction of KaiC in the C6 (blue line) and in the [B12C6] (green line) state (N=20). (B) Left panel: positive feedback by assigning the highest phosphorylated hexamer, C66, a kinase activity (Mehra et al, 2006). Middle and right panels: as in (A) but for N=6 phosphorylation steps. (C) Left panel: negative feedback induced by sequestration of KaiA proteins by the KaiBC complexes. Middle and right panels as in (B).

To reduce the space of possible topologies, we grouped the states of the reaction chain in high‐ and low‐phosphorylated states and assume equal reaction rates for the phosphorylation and binding processes within these groups. These restrictions do have influence on the maximal achievable amplitude for each topology but not on their principal ability to show significant oscillatory behaviour and will be relaxed after a set of oscillatory candidates has been found. The rewiring occurs now exclusively between these groups and the associated repressions or activations act equally on all phosphorylation reactions within the groups. In the example (Figure 3A), a feedback loop from C6(2) to C6(1) is therefore prohibited. Moreover, a feedback can target more than one group in a consecutive manner, indicated by two links in Figure 3B and C, resulting in additional reaction constants to be optimised for. Additionally, we allow only the states C6(0), C6(N), [B2mC6(0)], or [B2mC6(N)] to act as enzymes for the simple reason that declaring other states to enzymes does not lead to new oscillatory mechanisms but results in oscillations of lower amplitude. For simplicity, enzymatic reactions are treated strictly bimolecular, for example, the transitions C6(i) → C6(i+1) enhanced by a positive feedback loop scale bilinearly with the concentration of enzymes and hexamers in the ith phosphorylation state. The reaction constants are optimised to show maximal amplitude where search in parameter space ranges over four orders of magnitude around a physiological meaningful initial guess (see Supplementary information).

Robustness of potential clock candidates

The in vitro experiments (Kageyama et al, 2006) show that the time course of KaiC phosphorylation is unchanged under a collective five‐fold increase in protein concentrations and damped oscillations of significantly lower amplitude are observed if the concentrations are decreased by a factor of 10 (data reprinted in Figure 6A). These results indicate high robustness of the core oscillator against collective overexpression that can arise from stochastic fluctuations in gene expression (Elowitz et al, 2002; Kollmann et al, 2005) or from the genome‐wide circadian regulation of mRNA synthesis (Takai et al, 2006; Wijnen and Young, 2006). For the KaiABC clock, correlations in protein synthesis are expected to be a direct consequence of KaiA promoting kaiBC expression and KaiB and KaiC being part of the same operon. Because outstanding robustness of a network can generally not be achieved by fine‐tuning of reaction constants, we will rank in the following the generated oscillatory topologies with respect to their robustness against overexpression. Within the above‐stated rules of our rewiring approach, it was possible to construct 224 network topologies with one feedback loop that use one group as source and can target between one and three groups of high‐ or low‐phosphorylated states in a consecutive manner (see Supplementary information for details). The corresponding reactions (activation or repression), attributed with each link, are chosen such that they all increase or decrease the clockwise reaction flow upon changes in the amount of molecules in the source group. Among the generated networks, 103 candidates show oscillatory behaviour. However, they do not necessarily represent distinct topologies, as for small enough binding constants a two‐link feedback can reduce to a one‐link feedback. In total, we found eight distinct topologies that oscillate with an amplitude that changes KaiC phosphorylation levels more than 30%.

Figure 3 shows three examples of possible clock candidates. The reaction networks (Figure 3A and B) correspond to the first and second approach discussed in Introduction (Emberly and Wingreen, 2006; Mehra et al, 2006). The negative feedback from high‐ to low‐phosphorylated hexamers in Figure 3A describes effectively—although not mathematically equivalent—the monomer exchange among hexamers that has been proposed as possible synchronisation mechanism (Emberly and Wingreen, 2006). This topology, however, suffers the same shortcomings as the model proposed by Emberly and Wingreen (2006) namely that a sufficiently high amplitude can only be achieved by increasing artificially the number of residues, N=20 (this work) and N=96 (Emberly and Wingreen, 2006). Note that increasing the number of phosphorylation steps always implies a higher maximal amplitude. Also we observe a strong phase and frequency shift upon a concerted increase and decrease in the concentration of Kai proteins (Figure 3A, middle panel). The positive feedback in Figure 3B mimics the autocatalytic behaviour used by Mehra et al (2006) to generate oscillations of high amplitude. This topology has a deficit in robustness against common changes in protein levels as shown in Figure 3B, middle panel, where sustained oscillatory behaviour is completely lost in both cases. The topology depicted in Figure 3C shows sufficiently strong oscillations and also reproduces the experimentally found behaviour with respect to changes in the concentrations of the Kai proteins. Latter network is the most promising candidate for the cyanobacterial in vitro clock from all reaction networks generated by our rewiring scheme. The networks depicted in Figure 3B and C are among the eight found topologies that show strong oscillatory behaviour. The other network topologies with high amplitude (see Supplementary information) use either an autocatalytic feedback mechanism as discussed in Figure 3B or a negative feedback similar to Figure 3C. However, these network topologies show a higher shift in frequency or phase upon elevation of protein concentrations than experimentally observed. We emphasise that in the rewiring scheme, reaction constants have been exclusively optimised for highest amplitude and not for robustness against changes in abundance of its components. Thus, the robust behaviour of topology (Figure 3C) is solely a consequence of its network structure. That outstanding robustness of biological networks is a consequence of its structure and cannot be achieved by fine‐tuning of reaction constants has been shown by mathematical investigation for the Drosophila segmentation network (von Dassow et al, 2000) and the chemotaxis network in Escherichia coli (Barkai and Leibler, 1997).

Biochemical interpretation of the most robust oscillatory reaction network

We can give now a biochemical interpretation to the negative two‐link feedback loop shown in Figure 3C. A simple interpretation is that KaiA binds with high affinity to low‐phosphorylated KaiBC complexes. That this molecular mechanism is indeed mathematically equivalent to the functioning of the feedback loop in our rewiring scheme (Figure 3C) is shown in Supplementary information. We included the concentrations of KaiA and KaiB explicitly in the kinetic equations (see Supplementary information) and optimised all reaction constants separately with respect to highest amplitude. The optimisation is performed under the constraints to reproduce the experimentally found dynamical behaviour of the phosphorylation kinetics in the absence of KaiB (Supplementary Figure S2) and invariance of phase and frequency upon a concerted five‐fold elevation of all protein concentrations (Kitayama et al, 2003; Xu et al, 2003; Tomita et al, 2005; Kageyama et al, 2006). The BN‐PAGE analysis (Figure 1) also indicates the appearance of free monomers in the system after KaiBC complexes have been built. This suggests a possible reshuffling mechanism of KaiC monomers, similar to the one introduced by Emberly and Wingreen (2006) and experimentally verified by Kageyama et al (2006). We complete our mathematical model by accounting for monomer exchange among KaiBC complexes and obtain with the appearing shuffling rate another parameter to be optimised for. The resulting time courses are shown in Figures 4 and 5. To discriminate the contributions of the two synchronisation mechanisms, KaiA sequestration and monomer exchange, we investigated these mechanisms separately with respect to their contribution to amplitude and robustness. Figure 4 shows the percentage of phosphorylated KaiC over four periods for two different concentrations of Kai proteins, the one used in the in vitro assay (Figure 1) and a concerted five‐fold increase of these concentrations. The oscillatory behaviour shown in Figure 4A is the result of optimisation for high amplitude and stability of phase and frequency upon a five‐fold increase of all protein levels in the absence of monomer exchange. Figure 4B shows similar time courses as Figure 4A but includes also the monomer exchange rate in the optimisation process. Here, a slight improvement in amplitude is observed along with an increase in the maximal amount of KaiBC complexes during one cycle. The reason why monomer exchange is not rewarded better as a synchronisation mechanism can be seen from Figure 4C. Here, we included the monomer exchange rate in the optimisation process for high amplitude but imposed no constraints on robustness of phase or frequency. We find that the amplitude of free hexamer and KaiBC complex abundances has increased strongly (Figure 4C, lower panel), indicating a change in the oscillatory mechanism. Although the KaiC phosphorylation oscillates between the minimal (20%) and maximal (90%) achievable phosphorylation levels, we observe at the same time a shift in frequency upon elevation of protein concentration and thus the loss of robustness in the system. Given the biochemical constraints of Figure 2, our theoretical analysis reveals that under selective pressure for high amplitude and robustness against concerted changes in protein levels, the oscillatory mechanism with highest performance is given by the reaction network depicted schematically in Figure 3C, left panel.

Time course of phosphorylated KaiC (upper panels) and free hexamers and KaiBC complexes (lower panels) for the model using KaiA sequestration as a negative feedback. Red lines arise from a concerted five‐fold increase of the native protein concentrations (blue lines) after optimisation. Lower panels: solid and dashed lines indicate the amount of free hexamers and KaiBC complexes, respectively. (A) With reaction constants from optimisation for both highest amplitude and robustness against concerted five‐fold elevation of protein concentrations in the absence of monomer exchange. (B) Optimisation as in (A) but accounting also for monomer exchange. (C) Optimisation for highest amplitude and accounting for monomer exchange but without any constraints on robustness.

Upper panels: abundances of the complexes KaiAC, KaiBC, and KaiABC over time, relative to the total amount of hexamers in the system. Lower panels: relative concentrations of free KaiA and KaiB. (A) Optimisation as in Figure 4A. (B) Optimisation as in Figure 4B.

The striking result that the molecular mechanism in the cyanobacterial circadian clock indeed corresponds to this theoretical found oscillator mechanism is shown in Figures 5 and 6. Here, our conjecture that KaiA sequestration to low‐phosphorylated KaiBC complexes is the primary oscillator and monomer shuffling is just a secondary effect is in good agreement with recent experiments (Kageyama et al, 2006), whose results are reprinted in Figure 6. The figure shows the relative amount of KaiAC, KaiBC, and KaiABC complexes found during two cycles. The mean levels of the time courses are not quantitative, as our model does not account for the experimentally observed excess binding of KaiA and KaiB to free hexamers and KaiBC complexes and the two‐step pulldown assay used in experiments (Kageyama et al, 2006) might partially dissociate KaiA from KaiABC complexes. However, the phase relation of the different time courses is in agreement with experiments. The KaiA–KaiB interactions peak close to 26 h after incubation, followed by the formation of KaiBC complexes with a maximum observed after 30 h, and the sequestration of KaiA by KaiBC complexes peaking close to 36 h after incubation. The oscillatory mechanisms with (Figure 5B) and without (Figure 5A) monomer exchange are quite similar when optimised for robustness against upregulation of protein levels, although the amplitude for KaiBC complexes increases significantly.

Time course of KaiC phosphorylation and complex formation in comparison with experiments. Upper panels: results from the mathematical modelling using the optimisation approach from Figure 4A; lower panels: experimental results from Kageyama et al (2006). (A) Relative level of phosphorylated KaiC for different protein levels. (B) Complex formation of KaiAC, KaiBC, and KaiABC.

Discussion

We can now give a picture of the functioning of the in vitro circadian clock as shown in Figure 7. Our optimisation procedure revealed that the maximum oscillation for the reaction network (Figure 3C) is achieved when only about 10–20% of the hexamers form stable complexes with KaiB. Consequently, the largest fraction of phosphorylated KaiC (80–90%) never takes the circular route in the reaction network but remains part of free hexamer pool. This counter‐intuitive finding is supported by experiments, indicating that—after transients have declined—about 70% of the hexamers do not bind KaiB (Kageyama et al, 2006). This result reveals that in order to achieve highest amplitude, the phosphorylation dynamics oscillates only to small extent between free hexamers and KaiBC complexes. Temporal variations in the amount of KaiBC complexes drive the phosphorylation level of free hexameric KaiC by rhythmically changing its phosphorylation rates. This is achieved by recruiting free KaiA to low‐phosphorylated KaiBC complexes, as evidenced by Kageyama et al (2006) and our experiments (Figure 1). The shuffling mechanism is probably just a secondary mechanism that increases the amplitude of oscillations. If one forces monomer shuffling to be the primary cause of synchronisation, for example, by lowering the maximum binding affinity of the complexes for KaiA or increasing the shuffling rate, we observe that close to 100% of KaiC takes the circular route (Figure 4C, lower panel). This behaviour is obvious as only complexes reshuffle (Kageyama et al, 2006) and therefore the circular route via complex building is favoured. However, monomer exchange in the absence of KaiA sequestration does not lead to significant oscillations (Supplementary Figure S3). Although a common synchronisation via monomer exchange and KaiA sequestration leads to highest oscillations, the robustness against collective changes in protein concentrations is lost. Thus, the realisation of the cyanobacterial clock via a robust negative feedback of KaiA sequestration as the primary mechanism is likely the result of a strong selective pressure for stability of phase and frequency against perturbations.

Cartoon of the KaiC phosphorylation cycle. Phosphorylated KaiC monomers are marked by filled circles. The relative amount of molecules involved in the transitions between the indicated states is indicated by the line thickness of the arrows. After incubation, autophosphorylation of KaiC hexamers is enhanced by KaiA. After 6 h, KaiBC complexes start to build. The KaiBC complexes dephosphorylate, presumably by inhibiting access of KaiA to its active sites for KaiC phosphorylation. Low phosphorylated KaiBC can bind KaiA with high affinity and this leads to dephosphorylation also of the KaiC hexamers. The synchronisation of the KaiBC phosphorylation level by monomer exchange plays only a minor role.

The striking fact that the experimentally observed dynamics can be reproduced by our approach (Figures 5 and 6) without detailed knowledge about reaction constants gives evidence for our initial assumption that the cyanobacterial core clock is the result of an evolutionary optimisation process towards robust sustained oscillations with high amplitude under functional constraints (Figure 2). In turn, strong indications that a system is highly optimised for a specific function allows one to employ a systematic search among potential molecular mechanisms in order to isolate reaction networks that show optimal performance.

It is commonly believed that evolution of biological networks starts from primitive structures with links (e.g. high affinities) and nodes (e.g. proteins) added consecutively over time as long as they provide significant advantage in fitness to the organism. For most prokaryotes, it is also known that energy consumption for protein synthesis is kept at minimum. Thus, it is tempting to conjecture that under high selective pressure, the evolved networks result in the smallest sufficiently robust topology for a given function (Kollmann et al, 2005). The results shown in this work give strong support for this hypothesis.

Materials and methods

Experimental methods

Purification of recombinant Kai proteins.

The recombinant GST fusion proteins were produced in E. coli BL21 strains provided by T Kondo and purified using GSTrap™ 4B columns and PreScission Protease (GE Healthcare, Uppsala, Sweden) as described (Iwasaki et al, 2002; Nishiwaki et al, 2004). The reconstitution of KaiC phosphorylation cycle in vitro was performed as described (Nakajima et al, 2005). Protein concentrations were estimated according to Lowry et al (1951).

Mathematical modelling

The phosphorylation dynamics of the KaiABC clock system is described on the level of ordinary differential equations using Michaelis–Menten kinetics. We optimised the reaction constants simultaneously for maximal amplitude and robustness against a concerted five‐fold increase in concentrations of Kai proteins using an adaptive simulating annealing routine (Lester Ingber, http://www.ingber.com/). Monomer exchange among KaiBC complexes is modelled on the level of mass action kinetics with rate proportional to squared concentration of KaiBC complexes. All simulations were performed under MATLAB from the MathWorks group. A detailed description of the mathematical modelling is found in the Supplementary information.

Acknowledgements

We thank Takao Kondo for providing the E. coli BL21 strains and Jana Wolf and Florian Geier for carefully reading the manuscript. This work was financially supported by the DFG Emmy Noether Programme and the SFB 618.