Skip to main content

Single molecules can operate as primitive biological sensors, switches and oscillators



Switch-like and oscillatory dynamical systems are widely observed in biology. We investigate the simplest biological switch that is composed of a single molecule that can be autocatalytically converted between two opposing activity forms. We test how this simple network can keep its switching behaviour under perturbations in the system.


We show that this molecule can work as a robust bistable system, even for alterations in the reactions that drive the switching between various conformations. We propose that this single molecule system could work as a primitive biological sensor and show by steady state analysis of a mathematical model of the system that it could switch between possible states for changes in environmental signals. Particularly, we show that a single molecule phosphorylation-dephosphorylation switch could work as a nucleotide or energy sensor. We also notice that a given set of reductions in the reaction network can lead to the emergence of oscillatory behaviour.


We propose that evolution could have converted this switch into a single molecule oscillator, which could have been used as a primitive timekeeper. We discuss how the structure of the simplest known circadian clock regulatory system, found in cyanobacteria, resembles the proposed single molecule oscillator. Besides, we speculate if such minimal systems could have existed in an RNA world.


Molecular networks exhibit a variety of dynamical behaviours based on the necessities of the cell [1,2,3]. Switch-like dynamics, such as toggle-switches, are usually found in regulatory processes of decision-making systems, as the system needs an all-or-none response. Toggle-switches transform a gradient input signal into an on/off output [4,5,6,7]. Changes in the input level cause that the responding molecule switches between low and high activities and/or concentrations (i.e. between off and on states). These changes cause shifts in the dynamical stability of the response, so the steady states of the system can be followed on a bifurcation diagram or signal-response curve [3, 8, 9]. Due to the nonlinear nature of the toggle-switch, if the signal is below a critical level (CL1), the response will be kept in one state (i.e. low concentration or off state); when the signal increases enough to go through this critical level, the response will abruptly change to the other state (i.e. high concentration or on state) (Fig. 1). For intermediary values of the signal, a bistable regime may arise. Over this region, the state of the response is not determined solely on the signal value, but it also depends on the earlier state of the response. Thus, an increase of the signal allows the system to pass the CL1 and reach the on state. At this stage, the system will not be able to go back to its previous state, off state, until the signal decreases enough to pass another, lower, critical level (CL2) [4,5,6, 10, 11] (Fig. 1).

Fig. 1

Approximate Majority (AM) system. a Wiring diagram of the AM network. PP form attracts molecules into this state, while OO form does the opposite, both via catalytic reactions. Background reactions happen between all conformations at a low rate (grey shades on arrows). As an illustration of the embedded positive feedback loops, striped grey arrows indicate the pure positive feedback loops. The striped dash-end indicate the double-negative positive feedback loop. b Bifurcation analysis of the AM system. The plot shows how the concentration of the catalytic molecules, OO and PP, is affected by the availability of phosphate donor (nt). The stable steady state (ss) is indicated by a solid line, while the unstable steady state (us) is defined by the dash-dotted line. CL1 and CL2 indicate the threshold points at which the jump between the off and the on states happen. Black striped arrows indicate the direction of the jumps. See networks topology in Methods for more information about the structure. Consult Table 2 in Methods for values of each parameter. c List of reactions. The left column shows the catalytic reactions, while the right column presents the spontaneous reactions. On top of each reaction arrows, the parameter names that affect the specific reaction are indicated

Biological systems present remarkably nonlinear dynamical features. A typical highly non-linear mechanism is the multi-step modification of proteins, like phosphorylation, methylation or ubiquitination [12, 13]. Consider a molecular population where most molecules are in the same off state (i.e. dephosphorylated, unmethylated, etc.) (below CL1, Fig. 1b). In this stage, the molecular population can repeatedly bind activators (i.e. phosphate donors in the case of phosphorylation) to be modified on a single site at each encounter (distributive multistep modification), without changing the overall off state response. The active or on state will be abruptly achieved when the molecules get modified at a critical ratio or critical site. This critical event leads to a major conformational change that will activate them (above CL1; Fig. 1b).

Responding molecules and their activators are embedded in positive feedback loops to coordinate a biological toggle-switch [2, 3, 7]. Positive feedback loops can contain only positive (activating) interactions between the involved molecules (pure positive feedback loop). They can also contain an even number of negative interactions (double-negative positive feedback loop) [10, 14]. Such interactions usually lead to bistability and the hysteresis effect. A bistable system offers two stable states: when the input level is below the critical level CL2, or above CL1, the system rests in a unique stable state, either the on or the off state. Between these critical levels, the system can rest in both stable states (Fig. 1b). The system will rest in either of the states depending on its previous history. Such effect is called hysteresis. Hysteresis can stabilise the system in either state even with a noisy input signal, making the system resilient to signal perturbations or noise [15].

The Approximate Majority system, a simple network first developed as an algorithm for distributed computing [16], presents these properties: nonlinearity, bistability and hysteresis [5, 6]. This network relies on three positive feedback loops (two pure positives and one double-negative) to generate both bistability and hysteresis (Fig. 1a; background grey dashed arrows indicate the implicit feedback loops of the system). The nonlinear response is obtained by successive multi-step modifications of molecules. Although these features are necessary properties for toggle-switch behaviour, the major advance of the Approximate Majority (AM) network is its efficient minimal topography. The AM system is considered as a minimal topology due to the number of molecular types that define it: a single kind of molecule that can be found in three different states. The two extreme forms hold catalytic activities for their own activation and inhibition, while an intermediate, undecided state separates them [5]. From the computing point of view, the AM systems is described as an efficient algorithm, since it is able produce a fast and robust toggle-switch by using a minimal number of components [16].

The dynamical features of the AM network have been associated with biological regulatory networks at various levels of complexity [6]; from the epigenetic cell memory system [17] to the G2/M transition of the cell cycle [5]. However, as far as our current knowledge goes, existing biological systems do not hold an exact organization of a single autocatalytic molecule that can work as an efficient toggle-switch. Earlier we have shown that AM can emulate the behaviour of much more complex biological switches and it stands as the smallest of such systems [6, 18]. In this work we would like to study how currently observed bistable biological regulatory networks could have emerged from a simple network, like AM.

The AM is based on a non-biological system; to get it closer to biological examples, we go into the molecular details of what AM can represent. AM presents a single intermediate state, which from a molecular perspective, is not precise. In the more realistic Two Intermediates (TI) system (Fig. 2), we split the single intermediate state of AM into two intermediate forms. These forms indicate the modification pattern of the two distinct sites. Thus, TI could be considered as an ancient biological system, where a single molecule can be modified in two separate sites. These modifications change its conformation in such way that it will hold different activities in the extreme forms. The TI network could work as a primitive molecular sensor due to its toggle-switch behaviour. This minimal system could have evolved by undergoing structural modifications, altering its topology and having an impact on its dynamical function.

Fig. 2

Two intermediates (TI) system. a Wiring diagram of the TI system. Single-modified forms (OP, PO) are physically different. b Bifurcation analysis of TI system. The graphic shows how the concentration of the catalytic molecules, OO and PP, is affected by the availability of phosphate donor (nt). We speculate that the TI system could have served as a primitive molecular sensor of phosphate donor level. Thus, TI could get into the PP form only above a critical threshold, where it remains until the phosphate level drops below a lower critical level. The stable steady state (ss) is indicated by a solid line, while the unstable steady state (us) is labelled by the dash-dotted line. See networks topology in Methods for more information about the structure. Table 2 in Methods for values of each parameter. c List of reactions. The left column shows the catalytic reactions, while the right column gives the spontaneous reactions. On top of each reaction arrow, the parameter names that affect the specific reaction are indicated

Under these premises, we derive reduced models from the TI system (with a basic set of parameters) and study their dynamical behaviours through steady-state analysis. We find that the toggle-switch behaviour, and therefore the function to work as a sensor, can be preserved for a large variety of topological changes. This highlights that regulatory networks derived from the AM class of systems are robust and efficient decision makers. Furthermore, we find that some topological alterations can transform the toggle-switch dynamics into an oscillatory one. This could allow the evolution of new biological functions originating from a minimal AM-like structure. We discuss two biological examples that could have emerged from a minimal system like AM. The circadian rhythm regulatory network of cyanobacteria shows high similarities to the oscillatory derivatives of the AM network. Dynamical features of riboswitches, that are also ribozymes, could also be related to the presented systems. These highlights that evolution might have found this system simple but robust, relying on these AM-like type of systems to build up some current regulatory networks.


The two intermediates (TI) system, a realistic AM-like system working as a molecular sensor

The Two Intermediates (TI) network is composed of a single kind of molecule with two antagonistic autocatalytic activities. This motif is a double-site modification system where the autocatalytic states act as enzymes. This network can generate a toggle-switch behaviour due to interactions between the different forms, embedded in several positive feedback loops (Fig. 2a). Molecules of the TI network can be found in three different forms, based on the number of modifications they hold: none (OO), single (OP/PO) or double (PP). They turn between conformations either spontaneously (slow) or by a catalyst (fast). OO and PP conformations hold catalytic activities, so they are not only substrates in the reactions but catalysts (see Methods). PP forms increase the number of modifications in the other forms (pure positive feedback loop), while the catalytic form OO removes these modifications (pure positive feedback loop). Thus, each catalytic conformation attracts molecules to that specific state, and away from the other (double-negative feedback loop).

Modifications usually induce structural changes in the conformations. The addition of a modification may affect the exposure of other sites, by burying or exposing them. Molecules of the TI system present two sensitive sites of modification, and depending on which site is modified first, a different single-modified form will be obtained (either OP or PO) (Fig. 2a). Thus, two separated patterns of modification are established in the TI network. For example, to reach a double-modified form (PP) through OP, the second site (_P) is always modified first, followed by the first site (P_). To come back to the unmodified form (OO) through the same single-modified form (OP), the first site (P_) needs to be unmodified first.

From a dynamical point of view, the TI system is a suitable candidate to work as a molecular sensor, since its toggle-switch dynamics create an all-or-none response (Fig. 2b). A molecular sensor informs downstream processes about changes in the environment. Particularly, such all-or-none switch can convert an analogue change in the input into a digital (none or all) output. Considering TI as a minimal topology for sensing, the simplest way for carrying out its function is responding to changes in the molecules that affect its modifications. When the sensed compound is in excess, the molecules of TI are fully modified; if the sensed compound level is low, molecules are fully unmodified.

The TI system, as many biological processes, needs an energy supply to accomplish its function. Phosphates, and more concretely highly energetic nucleotides like NTP, were early established as energy molecules [19, 20]. Thus, the simplest way for spontaneous modifications to take place is to use the cleaved phosphate from a phosphate donor, like NTPs, to modify itself. Using the same mechanism, low energetic nucleotides, such as NDPs or NMPs, can be utilised as phosphate acceptors to remove these modifications. Equivalently, the catalytic form PP uses NTPs to add phosphates into the other forms, attracting them to the double-modified state. OO form uses low energetic nucleotides to remove these phosphates, attracting the molecules to the unmodified state. In a way, the catalytic form PP works like kinases, and the catalytic form OO as a kind of phosphatase that reverses the kinase reactions by using the same catalytic core, as it happens at some bifunctional kinases [21,22,23].

The concentration of free NTP, or any phosphate donor (nt from now on), can be sensed by AM and TI networks, as all modification reaction rates depend on the nt level. Thus, we consider the TI network as a feasible primitive nucleotide or phosphate sensor. Molecules are fully modified (double-phosphorylated, PP form) when nt is above a critical level; they switch back to full unmodified (dephosphorylated or OO form) only at a much lower nt level.

The TI system preserves the toggle-switch dynamic for a large variety of topological alterations

The TI system shows robust switching, as the system presents a wide bistability region (Fig. 2b). It is well-known that double-site phosphorylation systems, which involve the action of kinases and phosphatases, can generate bistability [24,25,26]. A simple double-site modification motif on the MAPK activation process has been investigated in that respect [24]. The MAPK motif may look similar to the TI network at first sight, however, the mechanism and its implications greatly differ. The MAPK system relies on external kinases and phosphatases to produce a bistable response. This bistable regimen is the result of the sequestration of these external enzymes in substrate-enzyme complexes. On the contrary, the TI system achieves the bistability through autocatalytic conversion of the single molecule into the antagonistic enzymatic conformations. These conversions, and the produced bistable region, are controlled, but not driven by the external level of a phosphate donor.

The TI system has a flexible structure, as the two catalytic forms can be reached through two separate paths of phosphorylation (through OP or PO). Flexible structures, such as the topology of the TI network, may retain their dynamical behaviour despite detrimental alterations. That is, even if the conformations lose their catalytic and/or spontaneous activity, the network structure is robust enough to maintain the toggle-switch dynamics. Due to our interest in the properties of the network instead of the effects of individual parameters, a dummy set of parameters is chosen. All catalytic reaction rates are equal to the unit (1 AU); all spontaneous reaction rates are set at the same low value of 0.05 AU (see Methods for more details). When we systematically test TI for the loss of path reactions (a catalytic and its associated spontaneous reactions), we find that systems with the loss of any one of the 8 reaction paths preserve the toggle-switch dynamics (Additional file 1: Figure S1).

To show alternative dynamics, the TI network needs to suffer structural changes that affect at least two reaction paths (Fig. 3). There is a total of 16 different combinations when two reaction paths are lost (two catalytic and their spontaneous activities are removed). However, due to the symmetry of the TI network, all these systems can be represented by 8 topologies (Fig. 3a). Some of these networks show directionality for the conversions between OO and PP forms, but remarkably, the toggle-switch dynamic is preserved in 6 out of the 8 possible topologies (Fig. 3a). As we have previously mentioned, the symmetry of the topologies is extended to the parameter choice. Although the symmetry of the parameter set could potentially influence the bistable dynamics, random parameter perturbations in the AM system have only a minor effect on the switching dynamics [15].

Fig. 3

Analysis of removing reactions in pairs from the TI system. a Bifurcation diagrams showing the steady states of the systems where two catalytic reactions, and their corresponded spontaneous reactions, have been removed from the TI network. The arrows above each panel indicate which reactions have been dropped. The stable state (ss) is indicated by a solid line, the unstable steady state (us) is defined by the dash-dotted line, and the maxima and minima of oscillations (amp – for amplitude) are shown by dashed lines. b Wiring diagram of the All Intermediate (AI) network. The reactions that convert PO into OO and PP are both missing, so all molecules end up in this single-phosphorylated form. c Wiring diagram of the Broken Paths (BP) network. The reaction from PO to PP was dropped together with the reaction from OP to OO. d Time course diagram for the BP network in the oscillatory regime (phosphate donor nt = 2.2 AU). e Zoom in the bifurcation diagram of the BP system (panel “p2 d0 bp2 bd0”). See networks topology in Methods for more information about the structure. Table 2 in Methods contain the values of parameters

There are two topologies, All Intermediate (AI) and Broken Paths (BP), that do not exhibit this behaviour (Fig. 3c, b). The AI network does not show any dynamical behaviour for changes in nt (Fig. 3a “p2 d1 pb2 bd1” panel, and 3b). In this topology, the abilities to change one of the single-phosphorylated forms (PO) are lost. The PO form serves now as a sink, and eventually, all molecules end up in this single-modified form (Fig. 3b). The BP network loses the bistability, although it keeps a sigmoid response for changes in nt. Strikingly, at the inflaction of this curve, the system shows a small region of oscillations (Fig. 3a “p2 d0 bp2 bd0” panel, and 3c). In this topology, the ability to add phosphates to one single-phosphorylated form (PO) is lost, together with the ability to remove phosphates from the other single-phosphorylated form (OP) (Fig. 3c). Thus, BP is a topology that works with directionality: PP can be reached only via OP, and OO can be reached via PO. The OO and OP forms are in balance with each other at low nt levels. This produces a reduction of the steady-state concentration of OO, compared with other topologies (Fig. 3a). When the phosphate donor increases beyond the threshold, a balance between PP and PO forms take control. Thus, together with the loss of hysteresis, the catalytic forms find the equilibrium at lower concentrations. Despite the loss of these dynamical features, over a small range of phosphate donor values (nt = [1.65, 2.40]), the BP network introduces a new behaviour: the oscillatory dynamics.

Alterations may lead the toggle-switch behaviour into oscillatory dynamics

The emergence of oscillations in such a minimal system (BP system) could raise new, more complex responses. Oscillations enable biological systems to anticipate regular changes, either internal or external, and prepare the processes accordingly [27,28,29]. The simplest example is the circadian clock [22, 23, 30,31,32,33,34]. Biological networks produce two types of sustained oscillations, relaxation or delay oscillators [1, 26, 35,36,37,38]. Delay oscillators are composed of at least three components, contain at least a negative feedback loop (an odd number of negative interactions) and generate some sort of delay. Relaxation oscillators are composed of at least two components, containing negative and positive feedback loops. Based on these criteria the BP system works as a relaxation oscillator (Fig. 3d). In this type of oscillator, the concentration of molecules slowly changes over a long period of time, and abruptly varies in a short time [35, 36, 38].

Depending on the phosphate donor level, the BP network can show both switch-like and oscillatory behaviours. However, this system is not efficient in either of these. When it is working as an oscillator, this dynamic is limited to a small range of phosphate donor values (nt = [1.65, 2.40] AU) and it does not even achieve the maximal amplitude (Fig. 3d). This reduction occurs because the PP and OP forms are together at the same level, as the OO and PO forms are appearing at the same time (Fig. 3a, d). When we look at a wider range of phosphate donor (nt) values, we can observe that it works as a simple switch, without hysteresis. For the same reasons as above, the catalytic states cannot convert all molecules into their own form (Fig. 3a).

Thus, the BP network is not able to work as an efficient molecular sensor. If we assume that TI is an ancient molecular sensor, the modifications that lead to BP could have destroyed its original function, driving this system to disappear from ancient organisms. However, the appearance of the oscillatory behaviour might contribute to further explorations of new topologies in favour of keeping and improving this oscillatory dynamics, and exploit it for other functions (Fig. 4, Additional file 2: Figure S2). In this scenario, the BP network may have suffered further structural alterations, removing spontaneous and/or catalytic reactions to achieve a system with robust oscillatory dynamics (Fig. 4, Additional file 2: Figure S2). Strikingly, it is enough to drop the spontaneous dephosphorylation of PO into OO (bd1 parameter), to greatly enlarge the oscillatory region for changes in phosphate donor levels (nt = [1.60, 52.35]; Fig. 4b, Additional file 2: Figure S2).

Fig. 4

BD1 system and the behaviour of altered topologies diverged from it. a Wiring diagram of the BD1 network. b Bifurcation (up) and time course (down) analysis of the BD1 network. The bifurcation diagram shows the changes in the steady states of the catalytic forms. The stable steady state (ss) is indicated by a straight line, the unstable steady state (us) is defined by the dashed-dot line, and the maxima and minima of oscillations (amp for amplitude) are shown by dashed lines. The phosphate donor for the time course diagram is set to nt = 4 AU. c Bifurcation analysis of the effect of the loss of one extra reaction from the BD1 network. The legend above each panel indicates which parameters have been removed (by name of the parameter and the arrow that represents it). Labels same as on panel (b)

This system (BD1 network) shows prominent relaxation-type of oscillations (Fig. 4b). In the BD1 system, PO behaves as a substrate, which slowly increases through the slow conversion from the PP form (since OO is almost totally absent, no autocatalysis helps this). When PO reaches a critical level, enough OO form appears to turn the autocatalytic steps on, converting eventually everything into OO. This form is not stable, so it is quickly converted back to PP through OP, and the whole process of building up PO starts from the beginning. This behaviour is maintained until high phosphate donor levels (nt = 52.35) are achieved. From this level, the oscillations are lost, since OO is getting totally absent and all molecules end up in the PO form (Fig. 4b, Additional file 2: Figure S2).

As shown above, the TI network can serve as a robust biological switch. It is resistant to perturbations in the form of loss of reaction paths in the network (Fig. 3a). Similarly, the BD1 network can preserve the oscillatory dynamic despite of various removals of catalytic and spontaneous reactions (Fig. 4c). Remarkably, the removal of three specific reactions leads to loss of the oscillatory behaviour. Both, the loss of the catalytic dephosphorylation of PO into OO (d1) and the catalytic phosphorylation of OO into OP (p0) (Fig. 4c), cut the remaining paths that can covert PP to OO and vice-versa. The loss of the spontaneous dephosphorylation of PP into PO (bd2) leads to a tiny regimen of oscillations (nt = [1.6–1.8] (AU)) before reaching an equilibrium (Table 1), where most molecules are converted into the double phosphorylated form (PP), while some goes to the single-modified form PO.

Table 1 Dynamical regimens of phosphate donor (nt) for each model. Columns indicate the type of behaviour; rows the models

Minimal oscillatory topologies derived from the TI system

Interestingly, other reductions of the BD1 network can lead to minimal systems, preserving the essential reactions for oscillatory dynamics (Fig. 5). The smallest topologies maintain the key reactions, which could not be removed from BD1 (Fig. 4): [1] the catalytic phosphorylation of OO into OP (p0), [2] the catalytic dephosphorylation of PO into OO (d1), and [3] the spontaneous dephosphorylation of PP into PO (bd2). To close the loop and reach all conformations, the phosphorylation of OP into PP can be driven either by the catalytic (Catalytic Oscillator (CO) network) or by the spontaneous reactions (Spontaneous Oscillator (SO) network) (Fig. 5a, b).

Fig. 5

The minimal oscillatory systems: Catalytic Oscillator (CO) network and Spontaneous Oscillator (SO) networks. Wiring diagram of CO (a) and SO (b) networks. Both networks maintain the spontaneous dephosphorylation of PP into PO (grey arrow), OO has the control over the modification of PO into itself, and PP keeps its catalytic control over the transition from OO to OP. In the CO network PP also catalyses the phosphorylation of OP into PP, while in the SO network this reaction is spontaneous. c Bifurcation analysis. The graphic shows how the steady-state concentration of the catalytic molecules, OO and PP are affected by the availability of phosphate donor (nt). Stable steady states (ss) precedes the oscillatory solutions (unstable steady state (us) – dashed line, max and min of the amplitude of oscillation (amp) - striped lines). d Time course diagram. The phosphate donor is set to nt = 4 AU. e Maximal amplitude of the BD1, CO and SO networks. f Period of the BD1, CO and SO networks over a range of phosphate donor values. See networks topology in Methods for more information about the structure and Table 2 in Methods for values of parameters

CO and SO networks work as relaxation oscillators (Fig. 5d). Similar to the BD1 network, the PO form works as a substrate and antagonises the catalytic forms OO and PP. In contrast to the BD1 network, the amplitudes of the catalytic forms of the smallest topologies increase with the phosphate donor (Fig. 5e). The CO system reaches the upper bound of amplitude faster than the SO system, merely because in the CO system PP catalyses the phosphorylation of OP, instead of leaving it for spontaneous conversion. For most of the oscillatory period, the molecules are slowly converted back to PO from the PP form. Once PO is leaking through OO, OO form turns on its autocatalysis, and converts all PO into OO. Then, OO is quickly converted back to PP. In the SO system we see that OP appears temporally, while in the CO system, the autocatalysis by PP is too fast to allow the appearance of high levels of OP. Despite this difference, both systems run with similar periods of oscillation (Fig. 5f). These two networks can be considered as minimal single molecule oscillators. Below, we discuss their possible relevance as timekeeping regulatory networks. Note, that the CO and SO networks can be converted symmetrically to systems where OO and PP are reverted.


Various evolutionary theories try to explain the appearance of complex regulatory interactions [39,40,41]. But it is generally proposed that the emergence of self-replicators and autocatalytic proteins was a key step in this [42, 43]. Following this idea, we introduced the TI model as a biologically feasible version of the Approximate Majority algorithm [16] as a minimal autocatalytic biological switch.

We could not identify an exact biological example for the TI system. We would require a single molecule that can take two opposing catalytic activities, and modify itself by these catalytic reactions. There are several systems which show similar network topologies and dynamics [17], but in these cases the feedback loops are not direct, instead intermediate molecules are affecting the catalytic reactions. Still, it is tempting to consider TI as a topology that could have served as a molecular sensor in primitive chemical or biological systems. Most biological processes are controlled by the energy state of cells [44]. As the TI network uses high energy nucleotides to convert itself between various states, it could have been acting as a nucleotide sensor in ancient systems. Particularly, TI would be found in PP conformation when the energy in the system is high, and in OO conformation when that energy is low (Fig. 2a,b). Since TI can robustly keep this feature (Fig. 3), one of the derivatives of TI could have been enough to serve as such a molecular sensor. In fact, the AM network is a culmination of the reduction, as this topology just preserves a single path for phosphorylation and dephosphorylation (Fig. 1a).

The energy state of a biological system is closely linked to the environmental conditions that surround it. As an example of this, photosynthetic organisms have a higher energy level during daylight, while reducing it in dark conditions [31, 34]. If the TI network would be used by such a cell, then it would be in the PP conformation during the daytime, and in an OO conformation at night, following the changes in ATP availability. As a function of light conditions, these alternating states of the TI network could affect internal processes. Since further loss of reaction paths can convert TI into a minimal oscillators such as BD1 (Fig. 4), CO or SO (Fig. 5), this could have helped the evolution of primitive circadian clocks that can anticipate changes in light conditions [45].

The minimal biological oscillator, called the KaiABC system controls the circadian clock in cyanobacteria [30, 46]. Importantly, the KaiABC-system generates sustained autonomous oscillations in vitro, when the three purified proteins, KaiC, KaiB and KaiA, and adenosine triphosphate (ATP) are mixed together [30]. In vivo, specific proteins regulating molecular activities during the day or night are binding to distinct forms of KaiC, the central molecule of the system [32, 33]. KaiC contains two phosphorylation sites (T432 and S431; abbreviated as T and S). During a 24 h period, it goes through a specific path of phosphorylation status changes: ST - > S/pT - > pS/pT - > pS/T - > ST (where p labels if a site is phosphorylated). The phosphorylation status of KaiC is the main driver of the oscillatory behaviour, and also the indicator of which processes (day or night) will be activated [47]. During daylight, KaiC autophosphorylates, and during the night, KaiC dephosphorylates itself using the same active site. The autokinase activity of KaiC is enhanced by KaiA, and the autophosphatase activity enhanced by KaiB [48, 49].

Figure 6 shows that the topology of the KaiABC-system shows high similarity to the above proposed minimal systems derived from TI. All the minimal models contain a single type of molecule with two sensitive sites of modification, four different conformations, and a phosphate donor is acting as the driver of the modifications. However, the dynamical behaviour of TI and most of its single reaction loss derivatives show bistability, not oscillations (Fig. 3). Stable relaxation oscillatory dynamics appears only for further reductions of the system (Fig. 3).

Fig. 6

Similarities in the topologies of the CO oscillatory network and the KaiC system of cyanobacteria. a CO system (b) KaiABC-system of S. elongatus. KaiA (blue triangle) stimulates the autokinase activity of KaiC to covert the fully dephosphorylated form of KaiC (red) into its fully phosphorylated form (blue). KaiB (flat orange circle) binds to fully phosphorylated KaiC, this leads to the displacement of KaiA and enhances the autophosphatase activity of KaiC. c KaiBC-system of Prochlorococcus. KaiB (flat orange circle) binds KaiC in its fully phosphorylated form (blue) and enhances its autophosphatase activity, but KaiA is missing, so KaiC can phosphorylate itself without support from KaiA [61,62,63,64]

Although the topology of the KaiABC-system resembles the minimal oscillatory systems CO and SO, there is a major difference between them. The derivate systems from the TI network do not require auxiliary proteins to work as oscillators. In contrast, the autoregulation of KaiC requires KaiA and KaiB to run together the circadian clock of the cyanobacteria Synechococcus elongatus [48, 49]. The presence of external molecules has an impact on the type of oscillatory behaviour that the systems show. While the CO and SO models show relaxation-type of oscillations, experiments and proposed models of the KaiABC system show sinusoid-like oscillations [50,51,52]. Similarly, the small oscillatory regime of the more complex BP system above shows more sinusoid-like oscillation (Fig. 3d). It has also been proposed that the core of the cyanobacteria circadian clock works as a relaxation oscillator [53, 54], which better match to the results of the minimal CO and SO models. Sinusoid-like oscillators can arise from a relaxation oscillator by adding an extra negative feedback loop, causing the delay that leads to a more symmetric oscillation. This external negative feedback loop comes through the interaction of external molecules KaiB and KaiC in cyanobacteria, and if similar effects are added to our CO and SO models we can also convert them to show more symmetric oscillations (Additional file 3: Figure S3).

Homologs of KaiC and KaiB have been found in various bacteria and archaea species [55,56,57,58,59]. Where KaiA is missing, the KaiC homologs show alterations in their sequences, which might facilitate their autophosphorylation reactions. These systems are closer to the proposed CO system, although here KaiB is still required for proper dephosphorylation steps (Fig. 6). The KaiBC system can sustain the circadian oscillations in purple bacteria [60], but in the marine cyanobacteria Prochlorococcus species, this system shows an hourglass behaviour. In these systems, the phosphorylation autonomously runs in one direction, but it requires a daily reset induced by environmental triggers [61,62,63,64]. There are organisms where both KaiA and KaiB are missing, while KaiC is still present [65]; but it is still not clear what related dynamics, if any, KaiC shows in these species [66,67,68]. Following our models, KaiC alone could exhibit either switch-like dynamics and work as a molecular sensor, like TI or one of its divergent topologies, or it could also show oscillatory dynamics, like CO or SO. The presented transition from TI towards CO and SO shows how a biological switch can be converted into an oscillator. This might have been a path the KaiC system followed; it could have started as a sensor of nucleotide level, then turned into an oscillator, when it could predict in advance the time when environment changes are expected [45].

Single molecule sensors of energy levels might have been required even at earlier stages of evolution. In a possible RNA world, riboswitches could have sensed various metabolites [69]. Riboswitches that are also ribozymes could have adjusted their self-regulatory rates dependent on the level of surrounding nucleotides [70]. In that way, these RNA molecules could have worked like the proposed TI network (Fig. 2) or one of its derivatives (Fig. 3) to sense nucleotide levels. This could have helped the emergence of replicating RNA molecules [71], which could have adjusted the initiation of their replication to the presence or absence of free nucleotides. These ideas sound theoretical, and indeed evolution could have taken different paths, but TI-like, robust RNA based bio-switches could be readily built [72], and probably these could be further modified to work as CO-like oscillators. Such oscillators could have been usefully in predicting daily fluctuations in energy resources, leading to the emergence of primitive circadian clocks with similar network structure as the one currently running (by proteins) in cyanobacteria. Similar, ATP dependent regulation of the cell cycle oscillations have also been observed recently [73, 74] further highlighting that ATP sensing based control of periodic processes could be highly beneficial and evolution could have selected for the appearance of such systems.


We have shown that a single molecule with two autocatalytic conformations can work both as a biological switch and as an oscillator. Such an oscillator is driving the circadian clock in cyanobacteria. The proposed network could be the simplest biologically plausible system that could perform these functions. In a proposed realisation as a phosphorylation-dephosphorylation switch it could also serve as a nucleotide or energy sensor. Based on these findings we propose that such a molecule could have served as a cornerstone during the emergence of self-replicating complex systems.


Networks topology

The topologies (wiring diagrams) represent a chemical reaction network. The biological entities that interact to produce the different models are described as proteins. The AM system is composed by a single protein that can be found in three different conformations, called OO, OP and PP. The TI and its derived systems are also based on the same single molecule. However, the proteins in these systems can adopt an extra conformation, the PO conformation.

This single molecular species is a protein regulated through (de)phosphorylation events. Because these single molecules are modelled assuming only two phosphosites, the impact of the (de)phosphorylation events is interpreted as following: the conformation OO denotes that the protein has not been modified, or it has been double dephosphorylated; the OP and PO conformations indicate a single (de)phosphorylation event; and the PP state shows the effects of being double phosphorylated. Note that the letter P indicates which phosphosite is phosphorylated, while the letter O shows the phosphosite that is not modified.

In the graphical representations, the different conformations of the main molecule (OO, OP, PO and PP) are treated as distinct chemical species. Interactions between conformations are represented by directional arrows, which represent the direction of the conversion reactions. As the reactions can be spontaneous and/or catalysed, the colour of the arrows indicates the mechanism of conversion. Black arrows indicate that the reaction happens either through the spontaneous or the catalytic mechanism. Grey arrows present the reactions that only occur spontaneously, while coloured arrows (either red or blue) describe the reactions that only happen through catalysis.

Modelling framework

The biological scenario developed in the main text is modelled under the concept of distributive step-modifications [75]. The distributive step-modification denotes the process by which an entity, that presents several sensitive sites for modification, takes at most one modification during a single encounter. This mechanism is explicit in the mathematical description when it is formulated following the mass action law. For example, a protein in the double dephosphorylated conformation OO adopts the double phosphorylated conformation PP after two separated modification encounters. This two step-modifications are defined separately through the mass action law, which is translated into a Hill-type dynamics of coefficient n = 2 [55].

The phosphorylation events modify the structure of the protein, and they may also change the function [105, 106]. Under this biological premise, the different conformations show the roles that proteins take in the systems. Concretely, the double (de)phosphorylated conformations OO and PP hold opposite catalytic activities. The OP and PO conformations are intermediary conformations that serve as a bridge between the catalytic conformations OO and PP. The existence of these intermediary conformations is key for modelling distributive step-modification processes. These conformations are the ones that introduce nonlinearity on the models.

The conversion between conformations may occur through two different mechanisms. A conformation undergoes either a spontaneous transformation, or a catalytic driven reaction. The spontaneous transformation is referred to a (de)phosphorylation event that the own molecule promotes. In the catalytic reaction, this (de)phosphorylation event is driven by one of the catalytic conformations, either OO or PP. These different mechanisms are included in the mathematical models through different terms in the quantitative descriptions. The catalytic driven terms include the catalyst and the substrate, while the spontaneous terms only contain the molecular species that changes.

The AM and TI systems are composed by either three or four conformations, respectively. However, the mathematically definition of one of the conformation is omitted in the quantitative models. Particularly, the missed equation is the one that details the behaviour of the intermediary conformation OP, which can be algebraically calculated, following the assumption that the total concentration of all the forms (tot) is fixed. In the case of the AM system this is: OP = tot - PP - OO. The remaining differential equations are:

$$ \frac{\mathrm{d} OO}{\mathrm{d}t}=d0\ast nd\ast OO\ast \left( tot- PP- OO\right)-p0\ast nt\ast PP\ast OO- bp0\ast nt\ast OO+ bd0\ast nd\ast \left( tot- PP- OO\right) $$
$$ \frac{\mathrm{d} PP}{\mathrm{d}t}=p3\ast nt\ast PP\ast \left( tot- PP- OO\right)-d3\ast nd\ast OO\ast PP+ bp3\ast nt\ast \left( tot- PP- OO\right)- bd3\ast nd\ast P $$

This approach is taken in order to ensure that the numerical simulations will be within the limits of the desired molecular concentration. The mathematical system of equations that describes the TI system and its derivatives also replaces the equation of the OP form for the subtraction OP = tot - PP - OO – PO:

$$ \frac{\mathrm{d} OO}{\mathrm{d}t}=-p0\ast nt\ast PP\ast OO-p1\ast nt\ast PP\ast OO+d0\ast nd\ast OO\ast \left( tot- PP- OO- PO\right)+d1\ast nd\ast OO\ast PO- bp0\ast nt\ast OO- bp1\ast nt\ast OO+ bd0\ast nd\ast \left( tot- PP- OO- PO\right)+ bd1\ast nd\ast PO $$
$$ \frac{\mathrm{d} PP}{\mathrm{d}t}=p3\ast nt\ast PP\ast \left( tot- PP- OO- PO\right)+p2\ast nt\ast PP\ast PO-d2\ast nd\ast OO\ast PP-d3\ast nd\ast OO\ast PP+ bp3\ast nt\ast \left( tot- PP- OO- PO\right)+ bp2\ast nt\ast PO- bd2\ast nd\ast PP- bd3\ast nd\ast PP $$
$$ \frac{\mathrm{d} PO}{\mathrm{d}t}=p1\ast nt\ast PP\ast OO-p2\ast nt\ast PP\ast PO+d2\ast nd\ast OO\ast PP-d1\ast nd\ast OO\ast PO+ bp1\ast nt\ast OO- bp2\ast nt\ast PO+ bd2\ast nd\ast PP- bd1\ast nd\ast PO $$

In both systems of equations, all terms are preceded by two parameters, the reaction rate and the modifier rate that modulates the propensity of the reaction. The reaction rates identify the type of reaction, spontaneous or catalytic, and its direction, towards OO or PP. Thus, spontaneous reactions that move to OO forms are indicated by bd0–3, while if it moves to PP forms, the parameters bp0–3 are used. The catalytic reactions that go towards OO are described by d0–3, and the ones that go towards PP are defined by the parameters p0–3. As catalytic reactions are faster, the reaction rates of spontaneous reactions are lower than the rates of the reactions that are catalysed.

The modifier rate identifies if the modification event is either a dephosphorylation or phosphorylation case. Phosphorylation events add a phosphate group from a phosphate donor in a sensitive site, while dephosphorylation events remove these phosphate groups and add them to a phosphate acceptor. The level of phosphate donor and acceptor are represented by the nt and nd parameters, respectively. The nt parameter is included in all spontaneous and catalytic reactions that move towards PP, while the nd parameter is contained in all spontaneous and catalytic reactions that go towards OO. In a biological context, the availability of the phosphate donor is usually more limited than the acceptor. Thus, the nd parameter stays as a fixed threshold in all models, while the nt parameter varies over a determined range. The specific values for each parameter previously described can be found in Table 2.

Table 2 Values of the parameters to conduct bifurcation analysis. Columns are the networks, and rows are the parameters. Initial conditions do not matter for these calculations, for simulations of oscillations we set OP = PO = 0 and OO = 2.1 PP = 1.9 to start from an asymmetric state

Bifurcation and time-course analysis

The bifurcation diagrams presented here (Figs. 1b; 2b; 3a; 4c; 5c) show the steady-state behaviour of the enzymatic forms OO and PP as the phosphate donor (nt) varies. The time-course diagrams (Figs. 3d; 4b; 5d) present the concentrations of all molecules as a function of time for a certain value of phosphate donor (nt).

All time-course and bifurcation analysis are obtained through our methods developed in R (version 3.2.2) [76]. These methods exploit the functions provided in deSolve [77] and rootSolve packages [78, 79] to conduct the numerical solutions. The methods to conduct these analyses are available for use at The graphical representation of the analysis make use of the ggplot2 package [80]. Functions from the reshape2 package [81] allow to obtain an appropriate format for plotting the data with ggplot2.

The models presented here were previously tested on commercial software. These tools are Visual GEC ( and Oscill8 (



All Intermediate


Approximate Majority


Broken Paths


Critical level


Catalytic oscillator


Phosphate acceptor


Nucleotide diphosphate


Phosphate donor


Nucleotide triphosphate


Spontaneous oscillator


Two intermediates


  1. 1.

    Thomas R, Thieffry D, Kaufman M. Dynamical behaviour of biological regulatory networks—I. Biological role of feedback loops and practical use of the concept of the loop-characteristic state. Bull Math Biol. 1995;57(2):247–76.

    Article  PubMed  CAS  Google Scholar 

  2. 2.

    Thomas R. Laws for the dynamics of regulatory networks. Int J Dev Biol. 2002;42(3):479–85. Others

    Google Scholar 

  3. 3.

    Tyson JJ, Chen KC, Novak B. Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell. Curr Opin Cell Biol. 2003;15(2):221–31.

    Article  PubMed  CAS  Google Scholar 

  4. 4.

    Santos SDM, biology FJES. On the cell cycle and its switches. Nature. 2008;454(7202):288–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  5. 5.

    Cardelli L, Csikász-Nagy A. The cell cycle switch computes approximate majority. Sci Rep. 2012;2:656.

  6. 6.

    Cardelli L, Hernansaiz-Ballesteros RDRD, Dalchau N, Csikász-Nagy A. Efficient switches in biology and computer science. PLoS Comput Biol. 2017;13(1):e1005100. Available from:

  7. 7.

    Angeli D, Ferrell JE, Sontag ED. Detection of multistability, bifurcations, and hysteresis in a large class of biological positive-feedback systems. Proc Natl Acad Sci. 2004;101(7):1822–7.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  8. 8.

    Csikász-Nagy A, Battogtokh D, Chen KC, Novák B, Tyson JJ. Analysis of a generic model of eukaryotic cell-cycle regulation. Biophys J. 2006;90(12):4361–79. [cited 2016 Mar 15];Available from:

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  9. 9.

    Ingalls B. Mathematical modelling in systems biology: an introduction. Cambridge: MIT Press; 2013.

    Google Scholar 

  10. 10.

    Thomas R, d’Ari R. Biological feedback. Boca Raton: CRC press; 1990.

  11. 11.

    Walleczek J. Self-organized biological dynamics and nonlinear control: toward understanding complexity, chaos and emergent function in living systems. Cambridge: Cambridge University Press; 2006.

  12. 12.

    Gunawardena J. Multisite protein phosphorylation makes a good threshold but can be a poor switch. Proc Natl Acad Sci U S A. 2005;102(41):14617–22.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  13. 13.

    Kapuy O, Barik D, Sananes MRD, Tyson JJ, Novák B. Bistability by multiple phosphorylation of regulatory proteins. Prog Biophys Mol Biol. 2009;100(1):47–56.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  14. 14.

    Milo R, Shen-Orr S, Itzkovitz S, Kashtan N, Chklovskii D, Alon U. Network motifs: simple building blocks of complex networks. Science (80- ). 2002;298(5594):824–7.

    Article  CAS  Google Scholar 

  15. 15.

    Cardelli L, Csikász-Nagy A, Dalchau N, Tribastone M, Tschaikowski M. Noise reduction in complex biological switches. Sci Rep. 2016;6:20214.

  16. 16.

    Angluin D, Aspnes J, Eisenstat D. A simple population protocol for fast robust approximate majority. Distrib Comput. 2008;21(2):87–102.

    Article  Google Scholar 

  17. 17.

    Dodd IB, Micheelsen MA, Sneppen K, Thon G. Theoretical analysis of epigenetic cell memory by nucleosome modification. Cell. 2007;129(4):813–22.

    Article  PubMed  CAS  Google Scholar 

  18. 18.

    Cardelli L. Morphisms of reaction networks that couple structure to function. BMC Syst Biol. 2014;8(1):84.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  19. 19.

    Westheimer FH. Why nature chose phosphates. Science (80- ). 1987;235(4793):1173–8.

    Article  CAS  Google Scholar 

  20. 20.

    Hunter T. Why nature chose phosphate to modify proteins. Philos Trans R Soc B Biol Sci. 2012;367(1602):2513–6.

    Article  CAS  Google Scholar 

  21. 21.

    McLeish MJ, Kenyon GL. Relating structure to mechanism in creatine kinase. Crit Rev Biochem Mol Biol. 2005;40(1):1–20.

    Article  PubMed  CAS  Google Scholar 

  22. 22.

    Egli M, Mori T, Pattanayek R, Xu Y, Qin X, Johnson CH. Dephosphorylation of the core clock protein KaiC in the cyanobacterial KaiABC circadian oscillator proceeds via an ATP synthase mechanism. Biochemistry. 2012;51(8):1547–58.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  23. 23.

    Nishiwaki T, Kondo T. Circadian autodephosphorylation of cyanobacterial clock protein KaiC occurs via formation of ATP as intermediate. J Biol Chem. 2012;287(22):18030–5.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  24. 24.

    Markevich NI, Hoek JB, Kholodenko BN. Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades. J Cell Biol. 2004;164(3):353–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  25. 25.

    Ferrell JE. Feedback regulation of opposing enzymes generates robust, all-or-none bistable responses. Curr Biol. 2008;18(6):R244–5.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  26. 26.

    Ferrell JE, Ha SH. Ultrasensitivity part III: cascades, bistable switches, and oscillators. Trends Biochem Sci. 2014;39(12):612–8.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  27. 27.

    Elowitz MB, Leibler S. A synthetic oscillatory network of transcriptional regulators. Nature. 2000;403(6767):335–8.

    Article  PubMed  CAS  Google Scholar 

  28. 28.

    Sprinzak D, Elowitz MB. Reconstruction of genetic circuits. Nature. 2005;438(7067):443–8.

    Article  PubMed  CAS  Google Scholar 

  29. 29.

    Hasty J, McMillen D, Collins JJ. Engineered gene circuits. Nature. 2002;420(6912):224–30.

    Article  PubMed  CAS  Google Scholar 

  30. 30.

    Nakajima M, Imai K, Ito H, Nishiwaki T, Murayama Y, Iwasaki H, et al. Reconstitution of circadian oscillation of cyanobacterial KaiC phosphorylation in vitro. Science (80- ). 2005;308(5720):414–5.

    Article  CAS  Google Scholar 

  31. 31.

    Rust MJ, Golden SS, O’Shea EK. Light-driven changes in energy metabolism directly entrain the cyanobacterial circadian oscillator. Science (80- ). 2011;331(6014):220–3.

    Article  CAS  Google Scholar 

  32. 32.

    Cohen SE, Golden SS. Circadian rhythms in cyanobacteria. Microbiol Mol Biol Rev. 2015;79(4):373–85.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  33. 33.

    Tseng R, Goularte NF, Chavan A, Luu J, Cohen SE, Chang Y-G, et al. Structural basis of the day-night transition in a bacterial circadian clock. Science (80- ). 2017;355(6330):1174–80.

    Article  CAS  Google Scholar 

  34. 34.

    Puszynska AM, O’Shea EK. Switching of metabolic programs in response to light availability is an essential function of the cyanobacterial circadian output pathway. elife. 2017;6:e23210.

    Article  PubMed  PubMed Central  Google Scholar 

  35. 35.

    Strogatz S, Friedman M, Mallinckrodt AJ, McKay S. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering. Comput Phys. 1994;8(5):532.

    Article  Google Scholar 

  36. 36.

    Tsai TY-C, Choi YS, Ma W, Pomerening JR, Tang C, Ferrell JE. Robust, tunable biological oscillations from interlinked positive and negative feedback loops. Science (80- ). 2008;321(5885):126–9.

    Article  CAS  Google Scholar 

  37. 37.

    Ferrell JE, Tsai TY-C, Yang Q. Modeling the cell cycle: why do certain circuits oscillate? Cell. 2011;144(6):874–85.

    Article  PubMed  CAS  Google Scholar 

  38. 38.

    Csikász-Nagy A, Dalchau N. What makes a biological clock efficient? in Essays for the Luca Cardelli Fest, Eds. Abadi M, Gardner P, Gordon AD, Mardare R. 2014. p. 85-94.

  39. 39.

    Kauffman SA. The origins of order: self-organization and selection in evolution. In: Spin glasses and biology. Oxford: World Scientific; 1992. p. 61–100.

  40. 40.

    Eigen M. Selforganization of matter and the evolution of biological macromolecules. Naturwissenschaften. 1971;58(10):465–523.

    Article  PubMed  CAS  Google Scholar 

  41. 41.

    Gánti T. Organization of chemical reactions into dividing and metabolizing units: the chemotons. Biosystems. 1975;7(1):15–21.

    Article  PubMed  Google Scholar 

  42. 42.

    Szathmáry E, Demeter L. Group selection of early replicators and the origin of life. J Theor Biol. 1987;128(4):463–86.

    Article  PubMed  Google Scholar 

  43. 43.

    Kauffman SA. Autocatalytic sets of proteins. J Theor Biol. 1986;119(1):1–24.

    Article  PubMed  CAS  Google Scholar 

  44. 44.

    Atkinson DE. Energy charge of the adenylate pool as a regulatory parameter. Interaction with feedback modifiers. Biochemistry. 1968;7(11):4030–4.

    Article  PubMed  CAS  Google Scholar 

  45. 45.

    Simons MJP. The evolution of the cyanobacterial posttranslational clock from a primitive “phoscillator.” J Biol Rhythm 2009;24(3):175–182.

  46. 46.

    Tomita J, Nakajima M, Kondo T, Iwasaki H. No transcription-translation feedback in circadian rhythm of KaiC phosphorylation. Science (80- ). 2005;307(5707):251–4.

    Article  CAS  Google Scholar 

  47. 47.

    Nishiwaki T, Satomi Y, Kitayama Y, Terauchi K, Kiyohara R, Takao T, et al. A sequential program of dual phosphorylation of KaiC as a basis for circadian rhythm in cyanobacteria. EMBO J. 2007;26(17):4029–37.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  48. 48.

    Iwasaki H, Nishiwaki T, Kitayama Y, Nakajima M, Kondo T. KaiA-stimulated KaiC phosphorylation in circadian timing loops in cyanobacteria. Proc Natl Acad Sci. 2002;99(24):15788–93.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  49. 49.

    Kitayama Y, Iwasaki H, Nishiwaki T, Kondo T. KaiB functions as an attenuator of KaiC phosphorylation in the cyanobacterial circadian clock system. EMBO J. 2003;22(9):2127–34.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  50. 50.

    Rust MJ, Markson JS, Lane WS, Fisher DS, O’shea EK. Ordered phosphorylation governs oscillation of a three-protein circadian clock. Science (80- ). 2007;318(5851):809–12.

    Article  CAS  Google Scholar 

  51. 51.

    Zwicker D, Lubensky DK, ten Wolde PR. Robust circadian clocks from coupled protein-modification and transcription–translation cycles. Proc Natl Acad Sci. 2010;107(52):22540–5.

    Article  PubMed  PubMed Central  Google Scholar 

  52. 52.

    Chaves M, Preto M. Hierarchy of models: from qualitative to quantitative analysis of circadian rhythms in cyanobacteria. Chaos An Interdiscip J Nonlinear Sci. 2013;23(2):25113.

    Article  CAS  Google Scholar 

  53. 53.

    Kurosawa G, Aihara K, Iwasa Y. A model for the circadian rhythm of cyanobacteria that maintains oscillation without gene expression. Biophys J. 2006;91(6):2015–23.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  54. 54.

    Ma L, Ranganathan R. Kernel mechanism of the cyanobacterial circadian clock is a relaxation oscillator. In: Decision and control and European control conference (CDC-ECC), 2011 50th IEEE conference on. IEEE; 2011. p. 5850–5.

    Google Scholar 

  55. 55.

    Dvornyk V, Vinogradova O, Nevo E. Origin and evolution of circadian clock genes in prokaryotes. Proc Natl Acad Sci. 2003;100(5):2495–500.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  56. 56.

    Ming H, Miyazono K, Tanokura M. Cloning, expression, purification, crystallization and preliminary crystallographic analysis of selenomethionine-labelled KaiC-like protein PH0186 from Pyrococcus horikoshii OT3. Acta Crystallogr Sect F Struct Biol Cryst Commun. 2007;63(4):327–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  57. 57.

    Kang H, Kubota K, Ming H, Miyazono K, Tanokura M. Crystal structure of KaiC-like protein PH0186 from hyperthermophilic archaea Pyrococcus horikoshii OT3. Proteins Struct Funct Bioinforma. 2009;75(4):1035–9.

    Article  CAS  Google Scholar 

  58. 58.

    Maniscalco M, Nannen J, Sodi V, Silver G, Lowrey PL, Bidle KA. Light-dependent expression of four cryptic archaeal circadian gene homologs. Front Microbiol. 2014;5:79.

  59. 59.

    Schmelling NM, Lehmann R, Chaudhury P, Beck C, Albers SV, Axmann IM, Wiegard, A. Minimal tool set for a prokaryotic circadian clock. BMC Evol Biol. 2017;17(1):169.

  60. 60.

    Min H, Guo H, Xiong J. Rhythmic gene expression in a purple photosynthetic bacterium, Rhodobacter sphaeroides. FEBS Lett. 2005;579(3):808–12.

    Article  PubMed  CAS  Google Scholar 

  61. 61.

    Holtzendorff J, Partensky F, Mella D, Lennon J-F, Hess WR, Garczarek L. Genome streamlining results in loss of robustness of the circadian clock in the marine cyanobacterium Prochlorococcus marinus PCC 9511. J Biol Rhythm. 2008;23(3):187–99.

    Article  CAS  Google Scholar 

  62. 62.

    Axmann IM, Dühring U, Seeliger L, Arnold A, Vanselow JT, Kramer A, et al. Biochemical evidence for a timing mechanism in Prochlorococcus. J Bacteriol. 2009;191(17):5342–7.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  63. 63.

    Axmann IM, Hertel S, Wiegard A, Dörrich AK, Wilde A. Diversity of KaiC-based timing systems in marine cyanobacteria. Mar Genomics. 2014;14:3–16.

    Article  PubMed  Google Scholar 

  64. 64.

    Mullineaux CW, Stanewsky R. The rolex and the hourglass: a simplified circadian clock in prochlorococcus? J Bacteriol. 2009;191(17):5333–5.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  65. 65.

    Loza-Correa M, Gomez-Valero L, Buchrieser C. Circadian clock proteins in prokaryotes: hidden rhythms? Front Microbiol. 2010;1:130.

  66. 66.

    Dvornyk V, Knudsen B. Functional divergence of the circadian clock proteins in prokaryotes. Genetica. 2005;124(2–3):247–55.

    Article  PubMed  Google Scholar 

  67. 67.

    Loza-Correa M, Sahr T, Rolando M, Daniels C, Petit P, Skarina T, et al. The legionella pneumophila Kai operon is implicated in stress response and confers fitness in competitive environments. Environ Microbiol. 2014;16(2):359–81.

    Article  PubMed  CAS  Google Scholar 

  68. 68.

    Ma P, Mori T, Zhao C, Thiel T, Johnson CH. Evolution of KaiC-dependent timekeepers: a proto-circadian timing mechanism confers adaptive fitness in the purple bacterium Rhodopseudomonas palustris. PLoS Genet. 2016;12(3):e1005922.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  69. 69.

    Breaker RR. Riboswitches and the RNA world. Cold Spring Harb Perspect Biol. 2012;4(2):a003566.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  70. 70.

    Talini G, Gallori E, Maurel M-C. Natural and unnatural ribozymes: back to the primordial RNA world. Res Microbiol. 2009;160(7):457–65.

    Article  PubMed  CAS  Google Scholar 

  71. 71.

    Szathmáry E, Smith JM. From replicators to reproducers: the first major transitions leading to life. J Theor Biol. 1997;187(4):555–71.

    Article  PubMed  Google Scholar 

  72. 72.

    Etzel M, Mörl M. Synthetic riboswitches: from plug and pray toward plug and play. Biochemistry. 2017;56(9):1181–98.

    Article  PubMed  CAS  Google Scholar 

  73. 73.

    Wang T, Zhao J, Ouyang Q, Qian H, Fu Y V, Li F. Phosphorylation energy and nonlinear kinetics as key determinants for G2/M transition in fission yeast cell cycle. arXiv Prepr arXiv161009637. 2016;.

    Google Scholar 

  74. 74.

    Papagiannakis A, Niebel B, Wit EC, Heinemann M. Autonomous metabolic oscillations robustly gate the early and late cell cycle. Mol Cell. 2017;65(2):285–95.

    Article  PubMed  CAS  Google Scholar 

  75. 75.

    Salazar C, Höfer T. Multisite protein phosphorylation–from molecular mechanisms to kinetic models. FEBS J. 2009;276(12):3177–98.

    Article  PubMed  CAS  Google Scholar 

  76. 76.

    R Core Team. R: A Language and Environment for Statistical Computing. Vienna; 2015. Available from:

  77. 77.

    Setzer KS and TP and RW. Solving Differential Equations in R: Package deSolve. J Stat Sofware. 2010;33:1–25. Available from:

    Google Scholar 

  78. 78.

    Soetaert K. rootSolve: nonlinear root finding, Equilibrium and steady-state analysis of ordinary differential equations. 2009.

    Google Scholar 

  79. 79.

    Soetaert K, PMJ H. A Practical Guide to Ecological Modelling. Using R as a simulation platform. Dordrecht: Springer; 2009. p. 372.

  80. 80.

    Wickham H. ggplot2: elegant graphics for data analysis. New York: Springer; 2009.

  81. 81.

    Wickham H. Reshaping data with the reshape package. J Stat Softw. 2007;21(12):1–20. Available from:

    Article  Google Scholar 

Download references


We thank Chris Barnes, Neil Dalchau, Carl Johnson, James Locke, Peter Parker and Andrew Phillips for insightful discussions.


This work was supported by Microsoft Research through its PhD Scholarship Programme, the European Social Fund (EFOP-3.6.2–16–2017-00013) and the Royal Society (Royal Society Research Professorship RP120138).

Availability of data and materials

Codes of models and analysis methods are available at

Author information




RDHB performed research, RDHB, LC and ACN conceived the study, analysed results and wrote the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Attila Csikász-Nagy.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Figure S1. Removal of reactions from the TI network. a) Bifurcation diagrams showing the steady states of the systems where a reaction path is lost from the TI network. The arrows above each panel indicate which reactions have been dropped from the TI network. The name of the reaction rates that are dropped are also indicated above each panel. The stable state (ss) is indicated by a straight line, while the unstable state (us) is defined by the dashed-dot line. b) The wiring diagram of the Two Intermediates (TI) system. (PDF 285 kb)

Additional file 2:

Figure S2. Removal of reactions from the BP network. a) Bifurcation diagrams showing the steady states of the systems where either a catalytic or a spontaneous reaction have been removed from the BP network. The arrows above each panel indicate which reactions have been dropped from the BP network. The name of the reaction rates that are dropped are also indicated above each panel. The stable state (ss) is indicated by a straight line, the unstable state (us) is defined by the dashed-dot line, and the maximal and minimal amplitudes (ma) are shown by the striped lines. b) Broken Paths (BP) system. (PDF 313 kb)

Additional file 3:

Figure S3. Negative feedback driven oscillations of a TI derivative system. a) Wiring diagram of a TI system derivate that interacts with an external molecule that presents two conformations (A and B). b) List of reactions. The left column shows the catalytic reactions driven by the main TI system, while the right column presents the reactions established with the external conformations A and B. On top of each reaction arrow, the parameter names that affect the specific reactions are indicated. c) Time-course diagram showing the sinusoidal behaviour of the system. Initial concentrations of the molecules: OO = 2; PP = 1; B = 2 AU, all others are 0. Parameters: p 0  = 1, p 3  = 0.5, d 1  = 0.5, d 2  = 1, k 3  = 0.1, tot = 3. (PDF 390 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Hernansaiz-Ballesteros, R.D., Cardelli, L. & Csikász-Nagy, A. Single molecules can operate as primitive biological sensors, switches and oscillators. BMC Syst Biol 12, 70 (2018).

Download citation


  • Computational biology
  • Mathematical modelling
  • Multistability
  • Bistability
  • Oscillation
  • Networks
  • Circadian rhythm
  • RNA world
  • Evolution
  • Approximate majority