article: The principal pacemaker of the circadian clock of the cyanobacterium S. elongatus is a protein phosphorylation cycle consisting of three proteins, KaiA, KaiB and KaiC. KaiC forms a homohexamer, with each monomer consisting of two domains, CI and CII. Both domains can bind and hydrolyze ATP, but only the CII domain can be phosphorylated, at two residues, in a well-defined sequence. While this system has been studied extensively, how the clock is driven thermodynamically has remained elusive. Inspired by recent experimental observations and building on ideas from previous mathematical models, we present a new, thermodynamically consistent, statistical-mechanical model of the clock. At its heart are two main ideas: i) ATP hydrolysis in the CI domain provides the thermodynamic driving force for the clock, switching KaiC between an active conformational state in which its phosphorylation level tends to rise and an inactive one in which it tends to fall; ii) phosphorylation of the CII domain provides the timer for the hydrolysis in the CI domain. The model also naturally explains how KaiA, by acting as a nucleotide exchange factor, can stimulate phosphorylation of KaiC, and how the differential affinity of KaiA for the different KaiC phosphoforms generates the characteristic temporal order of KaiC phosphorylation. As the phosphorylation level in the CII domain rises, the release of ADP from CI slows down, making the inactive conformational state of KaiC more stable. In the inactive state, KaiC binds KaiB, which not only stabilizes this state further, but also leads to the sequestration of KaiA, and hence to KaiC dephosphorylation. Using a dedicated kinetic Monte Carlo algorithm, which makes it possible to efficiently simulate this system consisting of more than a billion reactions, we show that the model can describe a wealth of experimental data.
  Circadian clocks, which allow organisms to anticipate changes in day and night, are a fascinating example of biological rhythms. One of the most studied and best characterized models of circadian oscillations is the cyanobacterium Synechococcus elongatus. It is now known that it combines a protein synthesis cycle [1–3] with a protein phosphorylation cycle [4], and in 2005 the latter was reconstituted in the test tube [5]. This stimulated a detailed characterization of its design principles, in a fruitful collaboration between experiments and modeling [6–17]. Yet, it has remained an open question how the clock is driven thermodynamically. The central components of the protein phosphorylation cycle are three proteins, KaiA, KaiB, and KaiC. KaiC forms a homohexamer with two phosphorylation sites per monomer, serine 431 and threonine 432 [3], which are phosphorylated and dephosphorylated in a well-defined temporal order [8,18]. KaiA stimulates phosphorylation [19–21], while KaiB negates the effect of KaiA [20–23]. Each individual hexamer can be thought of as a tiny oscillator, going through a phosphorylation cycle in roughly 24 hours. To generate macroscopic oscillations, the individual oscillators have to tick in phase with each other. Due to the inherent noise in biochemical reactions, however, the phosphorylation cycles of the respective hexamers will inevitably run out of phase when there is no mechanism to synchronize them. Modeling in combination with experiments indicate that the oscillations of the individual hexamers are synchronized via the mechanism of differential affinity: When the majority of the KaiC hexamers is in the dephosphorylation phase of the cycle, there are front runners that have already completed their cycle, and are now ready to be phosphorylated again. Phosphorylation requires KaiA, however, and the limited supply of KaiA is sequestered by the laggards that are still in the dephosphorylation phase of the cycle, where they bind KaiA very strongly together with KaiB. By taking KaiA away from the front runners, the laggards thus force the front runners to slow down, narrowing the distribution of phosphorylation states across hexamers [8,9]. This is the essence of the mechanism of differential affinity and it appears to be active not only during the dephosphorylation phase of the cycle [8,9], but also during the phoshorylation phase [9,24]. Monomer exchange between hexamers [6], observed in experiments [11], is an alternative synchronization mechanism. However, theoretical studies by us and others suggest that monomer exchange is not critical for stable oscillations [8–10]. While it is clear that the clock is driven by the turnover of ATP [25,26], how fuel turnover drives the phosphorylation oscillations is still unclear. In previous models [9,24,27], phosphorylation is driven by ATP hydrolysis, while dephosphoryation proceeds via the spontaneous release of the phosphate groups from the threonine and serine residues. Intriguingly, however, recent experiments have revealed that during the dephosphorylation phase of the clock, ATP is regenerated [28,29]: the phosphate groups on the serine and threonine residues are transferred back to ADP. If phosphorylation and dephosphorylation do not cause any net turnover of ATP, what then drives the clock? Clocks are necessarily dissipative, entailing a net turnover of fuel molecules per cycle. KaiC consists of two highly homologous domains, called the CI and the CII domain [30]. Both domains can bind and hydrolyze ATP [31,32], but only the CII domain can be phosphorylated [33]. The ATP regeneration experiments indicate that the phosphorylation and dephosphorylation of CII proceeds via the transfer of phosphate groups between the threonine/serine residues and the nucleotide bound to CII [28,29], leaving open the possibility that there is no net turnover of ATP on the CII domain. Here, we argue that the hydrolysis of ATP in the CI domain is the principal energetic driver of the clock. We present a new mathematical model of the post-translational Kai circadian clock in S. elongatus, which is based on the idea that ATP hydrolysis in CI drives a conformational switch of KaiC. Previously, it has been predicted that ATP hydrolysis plays an important role in driving conformational transitions [27,34] and that these transitions are vital to generating the oscillations [9,24], predictions that have found experimental support [27,34–38]. Our model, however, goes farther, predicting that ATP hydrolysis in the CI domain is the place where detailed balance must be broken in order to generate sustained oscillations. Our model is inspired by that of Van Zon et al. [9]. KaiC switches between an active conformation in which the phosphorylation level tends to rise, and an inactive one in which it tends to fall [9]. The model describes how KaiA binds to the CII domain of KaiC in the active conformation, and how KaiA can then drive phosphorylation by acting as a nucleotide exchange factor [39], stimulating the exchange of ADP for ATP. The model predicts that as the phosphorylation level in the CII domain rises, the release of ADP from CI slows down. The ADP-bound state makes the inactive conformation of KaiC more stable, causing the hexamer to flip to the inactive state and triggering dephosphorylation. In our model, ATP hydrolysis in the CI domain thus provides the thermodynamic driving force for the oscillations, while the phosphorylation in the CII domain provides the timer for the hydrolysis in the CI domain. While the coupling between ATP hydrolysis in the CI domain and phosphorylation in the CII domain is the main feature of the new model, leading to novel testable predictions, the model can also describe a wealth of additional experimental data. The differential affinity of KaiA for the different KaiC phosphoforms naturally explains the characteristic sequence in which the threonine and serine sites are phosphorylated [8,18]. In addition, while retention of ADP in the CI domain triggers a switch between the active and inactive KaiC conformations, the binding of KaiB to CI stabilizes the inactive state further, and leads, as in previous models [8,9, 24], to the sequestration of KaiA, necessary for synchronizing the oscillations. The model predicts that the slow binding of KaiB, as observed experimentally [40], introduces a delay between the moment that a given KaiC hexamer reaches its point of maximum phosphorylation, and hence no longer needs KaiA to progress along the phosphorylation cycle, and the moment that the same KaiC actually sequesters KaiA. In our model, this delay is essential because it allows the laggards to reach the top of the cycle before the front runners take away KaiA. Although our model predicts that ATP hydrolysis in CII is sufficient for generating cycles in individual hexamers, hydrolysis in CI is necessary for creating macroscopic oscillations: ATP hydrolysis in CI drives the conformational switch to the inactive state, which leads to the synchronization of the hexamers through the sequestration of KaiA. Lastly, the model can explain the experimental observation that the oscillation period is robust to variations in steady-state ATP/ADP levels, while the system can be entrained by transient changes in this ratio, which is one of the important mechanisms for coupling the clock to light [41].
  Our model of the in-vitro Kai circadian clock [5] builds on the hexamer model developed by Van Zon and coworkers [9]. But in contrast to that model, and following the models developed by Rust et al. [8,24,27], it explicitly keeps track of the two phosphorylation sites on each of the monomers, as well as their nucleotide binding states. The purpose of this section is to give an overview of the new model and its state variables, and to provide background information on ideas from previous models and their experimental justification. The new ingredients of the model, as well as their experimental movitation, are discussed only briefly; they are discussed in much more detail in the sections below. A cartoon of the model, illustrating the states KaiC can be in and how the system progresses through the cycle, is shown in Fig 1. Here we give a detailed description of how the phosphotranfer reactions, the ratio of ATP to ADP in the binding pockets and differential affinity of KaiA together give rise to the ordered phosphorylation of the serine and threonine sites in the CII domain. In steps, we present the foundations of our model together with the experimental results that underlie it and give a detailed mathematical description of the resulting free energies and reaction rates. In the previus section we discussed the phosphorylation dynamics of KaiC and the interaction between its CII domain and KaiA, and how these effects combine to generate the ordered phosphorylation of the threonine and serine sites in KaiC. The CI domain does not seem to play a crucial role here, in particular since the phosphorylation dynamics in the presence of KaiA only, is unaffected in a KaiC mutant where hydrolysis in the CI domain is deactivated [27]. This raises the question of what role the CI domain fulfills in the Kai oscillator. Here we describe how hydrolysis in CI together with the binding of KaiB to the CI domain drives the conformational switch of the hexamer and how the slow binding of KaiB, together with the subsequent sequestration of KaiA, helps to synchronize the ensemble of KaiC hexamers.
  To test the validity of our model of the CII domain and to find the correct parameter values shown in Table 2, we compare with the rich body of quantitative experimental results on the in-vitro Kai system. As is done in these experiments, we will study the behavior of different combinations of the main actors: KaiA, KaiB and KaiC and the ATP fraction, αATP. First, we will study the dephosphorylation dynamics of phosphorylated KaiC in the absence of KaiA and KaiB. We will investigate the dependence on the ATP fraction and compare to experimental results. Furthermore, we test our hypothesis that phosphotransfer from the threonine and serine sites to ADP is the only possible pathway for dephosphorylation, by comparing to a model where release of the phosphate into the bulk is possible. Next, we study the effects of KaiA on the phosphorylation dynamics and the influence of the bulk ATP fraction on the speed and steady state level of phosphorylation. Lastly, we study if the ordered phosphorylation dynamics of the serine (S) and threonine (T) residues persists when the system of KaiA and KaiC has reached steady state. In the subsequent section, we address the role of KaiB and the oscillatory dynamics. All simulations were performed at the experimental standard concentrations of 0. 6μM for KaiA and KaiC, which corresponds to simulating 720 KaiC hexamers and 720 KaiA dimers in a volume of 2 cubic micron. To test our model of the CI domain and to find the correct parameter values, shown in Table 3, we compare with the experiments on the evolution of the ATP fraction in the binding pockets of KaiC and the ATPase rate of KaiC [25,26,28,39]. As is done in these experiments, we will study the behavior of different combinations of the main actors: KaiA, KaiB and KaiC and the ATP fraction αATP. First we will study a system containing only KaiC, which allows us to constrain the parameters setting the rate of the hydrolysis cycle in the CI domain. Next, we look at the effect of dephosphorylation on the transient ATP fractions in the CI and CII domains, which provides an informative and testable prediction for how the CII domain regulates the ATP fraction in CI, which is indeed one of the key characteristics of our model. Then we study the effect of KaiA on the steady state ATPase rate and the dynamics of the ATP fraction in phosphorylating KaiC. For the full oscillating system, we will study the phase difference between the phosphorylation level and ATPase rate. Then we study whether the oscillations are robust to changes in the ATP fraction in solution, and whether our model can reproduce the experimental observation that the period of the oscillation is insensitive to these changes. Again, all simulations were performed with 720 KaiC hexamers and 720 KaiA dimers corresponding to the experimental standard condition of a 0. 6μM concentration in a volume of 2 cubic micron.
  Here we give an overview of the most important conclusions we can draw from our model. Although our model is able to reproduce most of the available experimental data, there are a few observations that it cannot replicate in its current form. First, in our model, KaiB barely interacts with unphosphorylated KaiC, in contrast with the observation that KaiB lowers the ATPase activity of a solution with only KaiC [25], and that KaiB can bind unphosphorylated WT KaiC [56] at micromolar concentrations. In our model, it is essential that unphosphorylated KaiC is predominantly in the active conformation: Monomers in the U state increase the dissociation rate from ADP in the CI domain, thereby stabilizing the active state which leads to the dissociation of the sequestered KaiA and KaiB at the end of the cycle. The consequence is that unphosphorylated KaiC is predominantly in the active conformation, which has a very low affinity for KaiB. Making the inactive state of unphosphorylated KaiC more stable would remedy this shortcoming of the model, because in the inactive state KaiC can bind KaiB. Consistent with this idea, very recent experiments suggest that the inactive conformational state is indeed more stable: About half of the unphosphorylated KaiC hexamers in a system without KaiA or KaiB, are in the inactive conformational state [59]. However, increasing the affinity of KaiB for unphosphorylated KaiC also increases, in the current model, the capacity of unphosphorylated KaiC to sequester KaiA, which impedes the release of KaiA at the end of the cycle. For future research, it will be interesting to see whether by amending the model these experimental observations can be reproduced. Secondly, the oscillations in our model come to a standstill when the bulk ATP fraction drops below 40%, while in experiments they continue to exist until the ATP fraction reaches 25% [27,53]. There are two reasons why the oscillations come to a standstill when the ATP fraction is reduced sufficiently. The first is that a lower ATP fraction lowers the rate of phosphorylation, while the rate of dephosphorylation remains constant. This reduces the number of KaiC hexamers that make it to the top of the cycle during the delay between the moment KaiC reaches the top of the phosphorylation cycle (and no longer needs KaiA to progress along the cycle) and the moment KaiC sequesters KaiA (and the laggards that are still in the phosphorylation phase can no longer progress to the top). Because there are fewer hexamers that make it to the top of the cycle, fewer hexamers will later participate in sequestering KaiA, which shortens the time window in which KaiA is sequestered, weakening the synchronization mechanism. In addition, the interplay of a lower phosphorylation rate and a constant dephosphorylation rate causes the distribution of KaiC phosphoforms during the dephosphorylation phase to widen. This increases the likelihood that the front runners reach the bottom of the cycle before all KaiA is sequestered; the front runners will then be phosphorylated again, causing the synchronization mechanism to break down. These effects are exacerbated by the fact that at lower ATP fraction the individual hexamers traverse a smaller phosphorylation cycle, as discussed in a forthcoming publication [60]. While these are the principal mechanisms by which oscillations cease to exist at sufficiently low ATP concentrations, a question is whether there are mechanisms that can mitigate these effects. One possibility is to make the time window of sequestration more deterministic. In the model presented here the binding of KaiB and the subsequent sequestration of KaiA is not directly coupled to the phosphorylation state of KaiC, in contrast to the dynamics in the Van Zon model [9]. As a result, both the distribution of the delay till, and the window of, KaiA sequestration is rather broad. This decreases the likelihood that during dephosphorylation there will be a time window in which there are enough hexamers to sequester KaiA. It thus conceivable that making the window of KaiA sequestration more deterministic, e. g. by making dephosphorylation of the respective monomers within a hexamer more concerted or by more tightly coupling KaiB-KaiC binding to the KaiC phosphorylation state, extends the range of ATP concentration over which the model exhibits oscillations. We have not investigated the effects of the concerted phosphorylation of hexamers, because it deviates too much from our current model where monomers phosphorylate independently. Another way to enhance the sequestration capacity might be to increase the affinity of unphosphorylated KaiC for KaiB (and thus, indirectly, KaiA), but as mentioned above, this impairs the release of KaiA at the end of the cycle. Furthermore, it is conceivable that adding an explicit cooperativity term in Eq 10, describing the binding of ADP and KaiB to the CI domain, might make the oscillations more robust. Finally, it is also possible that including monomer exchange might improve synchronization of the KaiC hexamers and thus allow oscillations to persist to lower ATP fractions, but such an effect likewise cannot readily be included in the current model. Lastly, the ATP consumption in our model is slightly higher than observed. We hypothesized that this can be attributed to the ATP hydrolysis in the CII domain. To provide support for this idea, we looked at an alternative model where the binding of KaiA suppresses ATP hydrolysis in CII. This model did show ATPase rates very similar to experimentally observed values. Here we explore the possibilities for experimentally verifying the predictions from our model. Our model is based on two important ingredients: 1) The relative stability of the two conformations is determined by the ATP fraction in the binding pockets of the CI domain and 2) this fraction is set by the number of phosphorylated serine and threonine sites in the CII domain of the hexamer. The dependence of the conformation on the nucleotide binding state of CI can be tested by measuring the ATP fraction in the binding pockets of KaiC [25,39] while at the same time probing the conformational state as in the study of [59], where the authors track the fractions of active and inactive KaiC hexamers over time. Our analysis predicts a positive correlation between the fraction of ADP bound to CI and the fraction of inactive KaiC, including KaiC mutants with a different hydrolysis rate constant in CI [25] and in an assay where KaiC is in the presence of KaiA and KaiB, and oscillates over time. The dependence of the ATP fraction in the binding pockets of the CI domain on the phosphorylation state of the CII domain can be tested by measuring the concentrations of monomers in the U, T, D and S phosphorylated state [8,54] and again the ATP fraction in the nucleotide binding pockets [25,39]. We predict that serine-phosphorylated KaiC slows down the dissociation of ADP from CI, decreasing its ATP fraction in the binding pocket, and threonine-phosphorylated KaiC should antagonize this effect. Starting with different ratios and levels of monomers phosphorylated at their threonine and serine sites, either using KaiC phosphomimics or aliquots from an oscillating system, the ADP fraction in the CI binding pockets should show a positive correlation with the difference between the number of phosphorylated serine sites and threonine sites. Adding KaiB should enhance the effect, because it cooperatively stabilizes the inactive state with ADP via the MWC mechanism. Our model can also explain the observation described in [61], where they found that adding ATP to unphosphorylated KaiC leads to a transient dip in the ATPase rate: The transient phosphorylation of KaiC temporarily lowers the CI-ADP dissociation rate. Clearly, it would be of interest to repeat these experiments starting with KaiC in different phosphorylation states, and in the presence and absence of KaiB. Our model predicts that, starting from fully phosphorylated KaiC, during dephosphorylation the ATP fraction in the binding pockets will exhibit a dip (see Fig 9). Related to this, and more specifically, our model predicts that the ADP dissociation rate in the CI domain is set by the relative number of phosphorylated serine and threonine sites in the CII domain, and not by their absolute levels. This implies that not all the residues have to be phosphorylated before a hexamer can switch to the inactive conformation and complete its cycle. Indeed as we show in a forthcoming publication, at lower ATP fractions of the buffer, hexamers go through a smaller phosphorylation cycle. This could in principle be tested experimentally if it were possible to track individual hexamers as they go through their phosphorylation cycle. While these experiments test our predictions on the connection between the CI and CII domain, our analysis also predicts an interesting consequence of the idea that phosphotransfer is the major pathway for dephosphorylation [28,29]. This could be tested by revisiting the experiments on dephosphorylation of radioactively labeled KaiC [28], but now in solution with non-hydrolyzable ATP. Since the ATP can not be hydrolyzed, there will be no ADP in the CII binding pockets, and the phosphate groups on the S and T sites in KaiC can not be transferred to ADP in the CII domain. This should significantly slow down the dephosphorylation speed if indeed phosphotransfer is the dominant pathway. Our model of the interaction of KaiA with the CI and CII domains of KaiC also allows us to make predictions for how long KaiA is bound to one of the domains during an oscillation. In our model, the dissociation rate of KaiA from the CII domain is much higher than the frequency of the oscillation. This is essential for differential affinity, where KaiA continually binds different KaiC hexamers during the phosphorylation phase to promote the phosphorylation of hexamers that are lagging behind. Our model therefore predicts a sharp peak in the distribution of times for which the CII domain of a certain KaiC hexamer is bound to KaiA during an oscillation period, ΔtCII·KaiA. Indeed, Fig 12A shows a clear single peak in this distribution. When we lower the dissociation rate of KaiA from the CII domain, this distribution broadens, Fig 12B, indicating differential affinity is hampered. One might think that when the dissociation rate is decreased even further, the distribution in ΔtCII·KaiA becomes bimodal, because once KaiA is bound to the CII domain of a certain hexamer, it will continue to stay bound to this hexamer during the whole phosphorylation phase of an oscillation cycle; because KaiA is limiting, this means that other KaiC hexamers will not, via their CII domain, bind KaiA during that cycle. In this case one fraction of hexamers does not or only very briefly binds KaiA via the CII domain while the other fraction is bound to KaiA for most of the time during the phosphorylation phase. However, in our model oscillations stop when we decrease the dissociation rate of KaiA to such a low level to allow for a bimodal distribution. Nevertheless, bi-modality could arise in experiments, which would indicate that the role of KaiA during the phosphorylation phase is very different from what we predict in our model. For KaiA bound to the CI domain, our simulations show exponentially distributed bound times, Fig 12B. This distribution is bimodal, as only 70% of the hexamers sequester KaiA, while the other 30% do not make it to the inactive state and bind six KaiB during an oscillation. Techniques to follow protein complex formation at the single molecule level have been developed [62], suggesting that future experiments might be able to reveal how long KaiA is bound to KaiC during an oscillation cycle.
  We model the post-translational Kai oscillator using kinetic Monte Carlo [63,64]. All code, written in C++, is available at https: //kmc-kaic. sourceforge. io. Each monomer has two nucleotide binding states for both the CI and CII domain and 4 phosphorylation states. Additionally, a hexamer can be in an active or inactive state. This results in (2 ⋅ 2 ⋅ 4) 6 ⋅ 2 = 225 different states a single hexamer can be in. Furthermore, there are the binding reactions with the KaiC hexamer: 2 for KaiA with the CII domain, 6 for KaiB monomers to the CI domain and another 6 reactions for KaiA binding to the CI domain of KaiC with 6 KaiB monomers. The total number of reactions therefore exceeds a billion. Clearly, the reaction combinatorics makes a straightforward ODE or Gillespie simulation unfeasible. Also the Quasi-Steady-State Approximation (QSSA) cannot be used, not only because of this combinatorial explosion of reactions, but also because most reactions do not exhibit a clear separation of timescales [65]. We thus set out to design a dedicated kinetic Monte Carlo algorithm to simulate the Kai system, where we do not have to write down all possible reactions explicitly. To this end, we designed an algorithm that keeps track of Nhex KaiC hexamers, which in turn consist of six explicitly simulated monomers, and N tot KaiA KaiA dimers. Since we have defined individual hexamers and monomers, we can calculate the propensity that a reaction occurs in a particular hexamer, and, in turn, the propensity that a reaction occurs in a monomer that is part of this hexamer. This allows us to create a layered version of the original algorithm by Gillespie [63], where we first determine the next reaction time and in which hexamer this reaction will take place. Then, in the next step, we determine which reaction or monomer of this hexamer fires. If a monomer fires, we choose which reaction of this monomer happens. This layered approach allows us to separate state changes in a hexamer that only modify the reaction propensities of that specific hexamer from state changes that influence the KaiA concentration in solution, which affects all hexamers. This greatly reduces the reaction combinatorics and the computational cost of the algorithm. Specifically, the state of the whole system, stot, consists of the state of the hexamers with index h, s h hex, and the number of KaiA dimers in solution, n sol KaiA: s tot = { { s h hex } h = 1 N hex, n sol KaiA }. The hexamer state contains the state vectors of its six monomers with index m, s h, m mon, the conformational state, C, the number of KaiA bound to CI, nCI·KaiA, the number of KaiB bound to CI, nCII·KaiA, and the number of KaiA bound to the CII domain, nCII·KaiA: s h hex = { { s h, m mon } m = 1 6, C, n CI · KaiA, n CI · KaiB, n CII · KaiA } h. The state vector of a monomer consists of the threonine and serine phosphorylation site and the nucleotide binding pockets of the CI and CII domain: s h, m mon = { S, T, n nucl CI, n nucl CII } h, m. Given the system state, we can calculate the firing propensity of reaction μ that changes the state vector of monomer m, which is part of hexamer h, q h, m μ (s h hex). Note that this reaction propensity can depend on the state of the whole hexamer, and therefore is a function of the hexamer state vector, and not only of the state vector of monomer m. The propensity for firing a reaction with index ν, which changes the state variables of hexamer h, denoted q h ν (s h hex, n sol KaiA), depends on the state vector of hexamer h only, and the number of KaiA in solution. Given these reaction propensities, we can calculate the accumulated propensities, denoted by q ˜, of firing a single monomer m in hexamer h, a single hexamer h and the total propensity as q ˜ h, m mon (s h hex) = ∑ μ q h, m μ, (19) q ˜ h hex (s h hex, n sol KaiA) = ∑ m = 1 6 q ˜ h, m mon + ∑ ν q h ν, (20) q ˜ tot ({ s hex } h = 1 N hex, n sol KaiA) = ∑ h = 1 N hex q ˜ h hex, (21) respectively. Since the firing of a reaction is a Markov process, the probability of hexamer h firing in the infinitesimal time interval [t + τ, t + τ + dτ] is P (τ, h | s tot, t) d τ = q ˜ h hex exp - q ˜ tot τ d τ. (22) Now, given two random numbers, ρ1, ρ2, drawn from a uniform distribution with domain [0,1), we calculate the next event time, τ, and the hexamer to fire, h, as τ = 1 q ˜ tot ln 1 ρ 1, (23) h = the smallest integer satisfying ∑ h ′ = 1 h q ˜ h hex > ρ 2 q ˜ tot, (24) respectively. Having defined the important propensities, our dedicated kinetic Monte Carlo algorithm becomes
summary: Circadian clocks are biological timekeeping devices with a rhythm of 24 hours in living cells pertaining to all kingdoms of life. They help organisms to coordinate their behavior with the day-night cycle. The circadian clock of the cyanobacterium Synechococcus elongatus is one of the simplest and best characterized clocks in biology. The central clock component is the protein KaiC, which is phosphorylated and dephosphorylated in a cyclical manner with a 24 hr period. While we know from elementary thermodynamics that oscillations require a net turnover of fuel molecules, in this case ATP, how ATP hydrolysis drives the clock has remained elusive. Based on recent experimental observations and building on ideas from existing models, we construct the most detailed mathematical model of this system to date. KaiC consists of two domains, CI and CII, which each can bind ATP, yet only CII can be phosphorylated. Moreover, KaiC can exist in two conformational states, an active one in which the phosphorylation level tends to rise, and an inactive one in which it tends to fall. Our model predicts that ATP hydrolysis in the CI domain is the principal energetic driver of the clock, driving the switching between the two conformational states, while phosphorylation in the CII domain provides the timer for the conformational switch. The coupling between ATP hydrolysis in the CI domain and phosphorylation in the CII domain leads to novel testable predictions.
keywords: phosphorylation
  chemical compounds
  phosphates
  enzymes
  enzymology
  threonine
  organic compounds
  phosphatases
  serine
  amino acids
  chemical dissociation
  proteins
  chemistry
  adenosine triphosphatase
  biochemistry
  hydrolysis
  organic chemistry
  post-translational modification
  hydroxyl amino acids
  biology and life sciences
  chemical reactions
  physical sciences
  atp hydrolysis
title: A thermodynamically consistent model of the post-translational Kai circadian clock