PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of interfaceThe Royal Society PublishingInterfaceAboutBrowse by SubjectAlertsFree Trial
 
J R Soc Interface. 2012 September 7; 9(74): 2365–2382.
Published online 2012 April 12. doi:  10.1098/rsif.2012.0080
PMCID: PMC3405750

Digital clocks: simple Boolean models can quantitatively describe circadian systems

Abstract

The gene networks that comprise the circadian clock modulate biological function across a range of scales, from gene expression to performance and adaptive behaviour. The clock functions by generating endogenous rhythms that can be entrained to the external 24-h day–night cycle, enabling organisms to optimally time biochemical processes relative to dawn and dusk. In recent years, computational models based on differential equations have become useful tools for dissecting and quantifying the complex regulatory relationships underlying the clock's oscillatory dynamics. However, optimizing the large parameter sets characteristic of these models places intense demands on both computational and experimental resources, limiting the scope of in silico studies. Here, we develop an approach based on Boolean logic that dramatically reduces the parametrization, making the state and parameter spaces finite and tractable. We introduce efficient methods for fitting Boolean models to molecular data, successfully demonstrating their application to synthetic time courses generated by a number of established clock models, as well as experimental expression levels measured using luciferase imaging. Our results indicate that despite their relative simplicity, logic models can (i) simulate circadian oscillations with the correct, experimentally observed phase relationships among genes and (ii) flexibly entrain to light stimuli, reproducing the complex responses to variations in daylength generated by more detailed differential equation formulations. Our work also demonstrates that logic models have sufficient predictive power to identify optimal regulatory structures from experimental data. By presenting the first Boolean models of circadian circuits together with general techniques for their optimization, we hope to establish a new framework for the systematic modelling of more complex clocks, as well as other circuits with different qualitative dynamics. In particular, we anticipate that the ability of logic models to provide a computationally efficient representation of system behaviour could greatly facilitate the reverse-engineering of large-scale biochemical networks.

Keywords: systems biology, circadian gene networks, Boolean logic, photoperiodism, Arabidopsis thaliana

1. Introduction

Circadian rhythms are the fundamental daily oscillations in metabolism, physiology and behaviour that occur in almost all organisms, ranging from cyanobacteria to humans [1]. The gene regulatory networks (GRNs), or clocks, that generate these rhythms regulate the expression of associated genes in roughly 24-h cycles. Circadian networks have been studied in a variety of experimentally tractable model systems, revealing that different organisms share structurally similar circuits based around interlocking sets of positive and negative gene–protein feedback loops [24]. In mammals, circadian rhythms are being increasingly recognized as important to healthy phenotypes, playing a role in ageing [5], cancer [6], vascular disease [7] and psychiatric disorders [8], as well as modulating innate immunity [911].

Clocks synchronize to their environment by using light and temperature to regulate the levels of one or more components of the feedback loops. This ensures that key biological processes are optimized relative to dawn and dusk, benefiting growth and survival [12,13]. For the clock to provide such an adaptive advantage, the phase must change appropriately when the clock is subject to regular perturbations—particularly seasonal changes in daylength (the photoperiod). However, as well as exhibiting flexible responses to variations in the input light signal, the clock must also exhibit robustness to irregular perturbations, such as genetic mutations and the intrinsically stochastic environment of the cell.

Temperature also plays a critical role as an environmental time cue. Across different species, the clock is relatively insensitive to temperature in that the period of free-running oscillations typically has a Q10 value close to 1 [1416]. This latter phenomenon, known as temperature compensation, is generally considered to be one of the defining properties of the circadian clock and has been suggested to be a key requirement for stability of the clock's phase relationship under seasonal temperature variations [17,18].

The ability of ordinary and delay differential equations (DEs) to reproduce the underlying continuous dynamics of biochemical networks, and to parametrize individual reactions has led to the construction of DE models in a number of circadian organisms. These include the fungus Neurospora crassa [1925], the fly Drosophila melanogaster [2629], the mammal Mus musculus [3033] and the higher plant Arabidopsis thaliana [3438]. Such models have proved useful in uncovering the general design principles of circadian oscillators, as well as providing a quantitative framework within which to interpret experimental results [4,38]. In particular, novel insights have been gained into the mechanisms promoting robustness with respect to photoperiod changes [25], temperature fluctuations [18,23] and molecular noise [3941]. The DE models have also yielded experimentally testable predictions that have led to the discovery of novel circadian regulators [36].

However, a significant drawback of the DE approach is that the values of the kinetic parameters controlling each individual reaction have to be specified, and for clocks these are typically unknown. When constructing a DE model, it is therefore necessary to calculate the particular combination of parameter values giving an optimal fit to experimental data [18,25,3437]. For realistic systems involving large numbers of reactions, this optimization procedure is computationally very expensive making exhaustive parameter searches intractable. With increasing parameter numbers also comes a need for data with which to constrain the optimization, placing a greater demand on experiment in terms of finance, time and ethics. These concerns mean that there is a pressing need for modelling approaches that minimize the number of parameters required, while adequately capturing the essential dynamical behaviour of the system of interest.

Here, we develop just such an approach, based on Boolean logic. In Boolean models, the activity of each gene is described with a two-state variable taking the value ON (1) or OFF (0), meaning that its products are present or absent, respectively. Biochemical interactions are represented by simple, binary functions that calculate the state of a gene from the activation state of its upstream components [4250]. This approximation dramatically reduces the state–space of the system, mapping the infinite number of different continuous system states in a DE model to a finite number of discrete states in the Boolean equivalent.

An additional important advantage of using a logic approach is that the total number of parameters is substantially reduced. For a given gene, the full set of reactions determining its state through a particular interaction is parametrized by a single signalling delay, representing the net time taken for these reactions to cause a change in state. Fitting to experimental data introduces an associated discretization threshold, with expression levels above the threshold taken to correspond to the ON state of the gene and levels below it to the OFF state [4648,50]. In fitting a particular experimental dataset, each delay becomes a multiple of the sampling interval, while only a bounded subset of thresholds will yield distinct Boolean expression patterns. This means that the total number of parameter combinations is finite and can be enumerated. Thus, by building a logic version of a DE model, an infinite model can be converted into a finite one with fewer parameters to be optimized. This extends the scale and complexity of GRNs that can be studied by Boolean models far beyond the practical scope of DEs.

In this work, we introduce the first Boolean models of circadian networks. By constructing logic analogues of a number of established DE clock models, we demonstrate that in each case the Boolean models are capable of accurately reproducing the higher-order properties—particularly photoperiod responses—of their DE counterparts. This suggests that the complex, biological signal transduction simulated by the DE models can be captured in Boolean equivalents possessing significantly smaller parameter sets. We introduce a general method for optimizing Boolean models that avoids the qualitative and often subjective terms characteristic of the cost functions used to fit the parameters of large DE clock models. Furthermore, we show that our fitting algorithm is capable of determining the optimal Boolean model configuration associated with a given circuit topology and experimental dataset. In particular, our algorithm successfully predicts the recently discovered repressive action of the circadian gene TOC1 on LHY in the central feedback loop of the Arabidopsis clock [51].

Taken together, our results show that Boolean models can quantitatively distinguish between a range of putative regulatory structures on the basis of the system dynamics. This identifies Boolean logic as a viable technique for reverse-engineering circadian networks, complementing approaches based on DEs. Moreover, our work also suggests novel hybrid modelling approaches based on employing Boolean models as a first step towards the construction of more detailed DE formulations. More generally, we propose that our methodology provides an efficient way of systematically modelling complex signalling pathways, including other oscillatory circuits and systems characterized by steady-state dynamics.

2. Results

2.1. Logic models employ significantly fewer parameters

We selected four recent circadian oscillator models of increasing complexity with which to assess the suitability of a Boolean formulation. The simplest of these was a Neurospora model based on a single negative feedback loop with a single light input [19] (figure 1a). This is represented by three DEs parametrized by 13 kinetic constants. The second model was a modified version of the 1-loop Neurospora circuit in which there is a pair of negative feedback loops associated with different isoforms of the active protein [18] (figure 1b). The extra feedback loop results in five DEs parametrized by 18 kinetic constants. The third model was an Arabidopsis circuit based on a pair of interlocking feedback loops with three light inputs [35] (figure 1c). It is described by 13 DEs together with 64 kinetic constants. The final model considered was a 3-loop Arabidopsis circuit obtained by adding an extra feedback loop and light input to the 2-loop system [36] (figure 1d). This yields 16 DEs parametrized by 80 kinetic constants.

Figure 1.
Circuit diagrams for the clock models. Genes are boxed and arrows denote regulatory interactions. Diamonds represent light inputs. (a) The single-loop Neurospora model [19]. FRQ protein represses production of frq transcript. Light acts on the network ...

In a Boolean model, an interaction, j, between two components Xi and Xk is quantified by the corresponding signalling delay τj. This is the time taken for the biochemical processes represented by j to convert a change in the state of Xi into a change in the state of Xk (electronic supplementary material, figure S1). The signalling delays are thus parameters that determine the dynamics of the model, with different combinations of delays yielding different attractor states (e.g. steady-states or limit cycles) [43,44,52,53]. In order to relate the discrete dynamics to the continuous variations in expression level observed experimentally, it is necessary to introduce a discretization process. A hypothetical time course is shown in the electronic supplementary material, figure S2a; and the result of applying a particular discretization threshold is shown in the electronic supplementary material, figure S2b. In the general case, the choice of threshold is dependent on the effective range of activity and expression that is critical for signal propagation. As a result, each threshold also becomes a parameter of the logic model [46,47,50].

For logic descriptions of the circadian networks shown in figure 2, each edge is parametrized by a signalling delay τj, and each vertex by a threshold Ti. Thus, a network is characterized by the vectors τ = (τj) and T = (Ti). In figure 3, we compare the total number of parameters in the logic and DE versions of each network (see also the electronic supplementary material, table S1). We can see that dramatically fewer parameters are required in a logic description compared with the corresponding DE formulation. Indeed for the largest network considered, 3-loop Arabidopsis, the Boolean description reduces the number of parameters by a factor of 4.

Figure 2.
The abstract topologies for the logic representations of the clock circuits shown in figure 1. Genes are boxed and arrows denote regulatory interactions. (a) The single-loop Neurospora model. (b) The 2-loop Neurospora model. (c) The 2-loop Arabidopsis ...
Figure 3.
The number of parameters required for each clock configuration as DE models (white bars) and logic models (black bars).

2.2. Logic model configurations consistent with differential equation models are optimal

In identifying the logic models that best reproduce the corresponding DE dynamics, we considered variations to the structures of the networks shown in figure 1. Starting from the abstract topologies of figure 2, each edge is activating or inhibiting, and where two or more edges lead into a vertex, the corresponding inhibition/activation signals are combined to determine the state of the gene.

In the Boolean formalism, the manner in which regulatory signals affect a gene's expression level is represented by the corresponding logic gate. This is a binary function that specifies the current state of the gene, ON (1) or OFF (0), for each possible combination of input states [46,54]. For genes with a single input, there are two possible functions for which the output varies with the input: the identity gate, in which the output follows the input (0 → 0 and 1 → 1); and the NOT gate, in which the input is inverted (0 → 1 and 1 → 0). These represent activation and repression of the gene by its regulator, respectively. Genes with two inputs are commonly modelled using either an AND or an OR gate [46,49,54]. For simplicity, we considered each multi-input gate to be a composition of ANDs and ORs (see §4 for further details). The particular combination of logic gates used to model a network is referred to here as its logic configuration (LC): this encodes the regulatory structure of the network in a compact fashion.

It follows that each abstract topology of figure 2 gives rise to 2E+M possible LCs, where E is the number of edges, and M is the number of vertices with more than one input. The total number of possible LCs is 4 for 1-loop Neurospora, 32 for 2-loop Neurospora, 256 for 2-loop Arabidopsis and 2048 for 3-loop Arabidopsis. In each case, a subset of these LCs is consistent with the pattern of activation and inhibition in the corresponding DE model. That there can be more than one such LC in each case is due to the choice of AND or OR where a vertex has multiple edges leading into it. For Neurospora, the 1-loop circuit has a unique LC consistent with the DE model, whereas for the 2-loop circuit, there are two DE LCs. The 2- and 3-loop Arabidopsis networks yield four and eight DE LCs, respectively.

For a given logic model, we would expect the DE LCs to most closely reproduce the DE dynamics. To test this hypothesis, we optimized LCs to synthetic experimental data obtained from the DE system, and then compared their predictive performance. For each LC, the match to the continuous dynamics was quantified by finding the combinations of parameters (signalling delays and discretization thresholds) that minimized a quantitative cost function. In order to be able to objectively compare the ability of the Boolean models to reproduce experimental time courses against that of their DE counterparts, we employed a cost function that closely mirrored those commonly used to optimize continuous models [18,25,3436,38]. The cost function we used measured the goodness-of-fit of each logic model to synthetic data generated in both 24 h light–dark (LD) cycles and the appropriate free-running light regime (continuous dark, DD, for Neurospora; continuous light, LL, for Arabidopsis). At each vertex, the cost score was calculated as the correlation between the discretized time series for the downstream species and the time-delayed predicted output calculated from the discretized data for the upstream species. Scores were summed across all vertices and light regimes to give the final cost value (see §4 for details).

After ranking the optimized LCs by score, we then assessed whether the top-ranking LCs comprised viable clock circuits by checking that they were capable of (i) generating self-sustained oscillations with a circadian period in constant conditions and (ii) entraining to LD cycles over a realistic range of photoperiods [55].

For the Neurospora and 2-loop Arabidopsis logic models, the full set of LCs was fitted to synthetic data. The best-performing LCs for these networks are shown in figures 4 and and55a. It can be seen that the optimization method produces a clear separation of the LCs by score. Moreover, for each network, one of the DE LCs is uniquely identified as the optimal circuit yielding a viable clock. In the case of 3-loop Arabidopsis, we considered the subset of configurations obtained by setting the gates common with the 2-loop Arabidopsis circuit to their optimized values. This mirrored the construction of the 3-loop DE model which was derived from the 2-loop system by adding an additional feedback loop while fixing all other interactions, and then optimizing the parameters of the new loop [36]. Constraining the structure of the circuit in this fashion yielded eight possible LCs (the two edges in the LHY–PRR loop described as activation or inhibition and the AND/OR interaction at LHY; figure 2d). Figure 5b shows that in this system, as for the other models, a DE LC emerges as the optimal circuit.

Figure 4.
The results of exploring the logic configurations (LCs) belonging to the abstract topologies of the (a) 1-loop and (b) 2-loop Neurospora models. Cost scores are shown for the optimal fit of each LC to synthetic data. The LCs are indexed by their decimal ...
Figure 5.
The results of exploring the logic configurations (LCs) belonging to the abstract topologies of the (a) 2-loop and (b) 3-loop Arabidopsis models. Cost scores are shown for the optimal fit of each LC to synthetic data. As in figure 4, each LC is ...

2.3. Optimal Boolean models have biological time-series characteristics

The time series generated by the optimal configurations in LD cycles are shown in figure 6. The corresponding DE simulations are also plotted for comparison.

Figure 6.
Time series for the differential equation and Boolean versions of the clock models in 12:12 LD cycles. Two 24-h cycles are plotted for each model. (a,b) 1-loop Neurospora; (c,d) 2-loop Neurospora; (e,f) 2-loop Arabidopsis; (g,h) 3-loop Arabidopsis. Differential ...

In each case, it is clear that the Boolean models capture the same qualitative dynamics as their DE counterparts. Different species are switched on and off relative to one another with phases that match the patterns of rising and falling expression in the corresponding continuous time series. Moreover, the delays between the switching times are similar to the phase differences between the peaks and troughs of the DE solutions.

It should be noted that both Boolean Arabidopsis models reproduce the acute light response in the Y gene, as well as the circadian response in Y around dusk (cf. figure 6eh). This demonstrates the ability of the Boolean circuits to simulate biochemical processes that occur on different time scales within the same system.

The optimal LCs give equally good matches to the DE dynamics in simulated free-running conditions, as can be seen in the electronic supplementary material, figure S3.

2.4. Optimal Boolean models have biological photoperiodic behaviour

In order to assess the extent to which the Boolean models reproduce the DE dynamics in a more global and quantitative fashion, we compared phase changes over a range of biologically realistic photoperiods. In the Boolean framework, the times at which solutions switch between 0 and 1 emerge as natural phase measures. For continuous time courses—such as those generated by DE models—the point in the circadian cycle at which the expression level of a molecular species decreases below a set threshold can be employed as a phase marker [25,56,57]. This suggested using the time at which each species decreases below its discretization threshold as a phase measure for the DE simulations, and the time of the 1 → 0 transition as the equivalent marker in the Boolean models.

The phase–photoperiod relationships computed in this fashion are shown for the Neurospora models in figure 7a,b. For both networks, the photoperiodic behaviour of the Boolean and DE models is very close: indeed for the 1-loop network, they are exactly equivalent. Figure 7c,d plots the photoperiod simulations obtained with the Arabidopsis circuits. Here too, the phase–photoperiod profiles are very similar, with the addition of the LHY–PRR loop to the 2-loop model causing a transition from a predominately dusk-locked system to a dawn-locked one [57]. In particular, the Boolean 2-loop Arabidopsis circuit exactly reproduces the dual light response in the Y gene, in which the acute peak tracks dawn, and the circadian peak tracks dusk. This suggests that the logic circuits possess sufficient dynamic flexibility to perform the complex integration of environmental signals that is a central property of circadian systems.

Figure 7.
Comparing the photoperiodic behaviour of the Boolean and differential equation (DE) versions of each model. For Boolean models, the phase of each species is taken as the time within the LD cycle of the ON to OFF transition (downward triangles). Analogously, ...

2.5. Boolean models can determine circadian network structure from experimental data

The success of the logic models in recovering the correct DE configurations from synthetic data suggested that for a fixed abstract topology, our optimization procedure has the capacity to determine the logic network most consistent with a given dataset. We tested this finding further by optimizing the 3-loop Arabidopsis logic circuit to highly sampled experimental time series recorded using luciferase (LUC) imaging in constant light from a wild-type strain [57]. All possible LCs were considered, corresponding to a network inference carried out assuming no prior biological knowledge. The cost function optimized was the same as that used for fitting to synthetic LL data. As previously, viable clock circuits were taken to be those yielding autonomous limit cycles with circadian periods.

The results of fitting to experimental data are presented in figure 8a. It can be seen that the second highest-ranking LC giving a viable circuit is a DE configuration. This configuration, GDE, is in fact the same as that previously determined to be optimal from the synthetic Arabidopsis datasets. Moreover, figure 8b shows that GDE emerges as the top-ranking clock configuration if the regulatory structure is constrained to incorporate the central LHY–TOC1 negative feedback loop of the corresponding continuous model. Interestingly, the optimal configuration, GOPT, differs from GDE in the sign of the TOC1–LHY interaction, with TOC1 repressing LHY in GOPT as opposed to activating it (compare electronic supplementary material, figure S4a,b). This result is consistent with the recent experimental characterization of the LHY–TOC1 circuit as a double-negative feedback loop, rather than the single-negative loop assumed in the 2- and 3-loop DE models [51].

Figure 8.
Identifying the logic configurations (LCs) of the 3-loop Arabidopsis model giving the best fits to experimental data. Plotting conventions are as described in figures 4 and and5.5. (a) Top-ranking LCs. The black arrow indicates the optimal ...

Figure 9b,c shows the simulations of the free-running clock obtained from the optimal and DE configurations, respectively. For both LCs, the simulated oscillation period and timing of gene expression are close to that of the natural system (cf. figure 9a). This can be seen clearly by comparing the durations for which each gene is ON in the two Boolean models with the corresponding peaks in the continuous data.

Figure 9.
Simulations generated by the optimal fits of the 3-loop Arabidopsis model to experimental data. (a) Experimental expression profiles for the genes CCA1, TOC1, GI and PRR9 in free-running conditions (LL). Expression levels were determined using LUC reporter ...

3. Discussion

3.1. A new approach to quantitative circadian modelling based on Boolean logic

Circadian clocks have become popular systems for studying the relationship between gene–protein dynamics and phenotype. The molecular machinery underlying these networks has been relatively well characterized, and this has led to great interest in developing predictive computational models of the clock. Such models are usually formulated as sets of DEs. The high level of biochemical detail afforded by this approach has allowed DE models to successfully address a range of issues regarding the functional relationship between the architecture of the clock and circadian homeostasis. One homeostatic mechanism that has been comprehensively investigated from a theoretical perspective is temperature compensation. By assuming that certain subsets of a circuit's kinetic parameters are temperature-dependent, temperature has been incorporated into a number of circadian DE models, ranging from minimal circuits built on the Goodwin oscillator [15,23,58] to more detailed formulations that explicitly model the underlying biochemical reactions [18,24]. These models have provided new insights into the possible mechanisms that lead to circadian mutants with altered compensation properties, as well as suggesting generic motifs that may facilitate the tuning of the phase–temperature relationship.

However, clock models based on DEs are characterized by large numbers of kinetic parameters, the values of which are typically unknown. Optimizing these to experimental data is a major computational bottleneck, which imposes a hard bound on the maximum system size that can be studied. In addition, as the fitting problem is typically highly underdetermined and involves noisy, undersampled experimental data, robust optimization can require the construction of cost functions targeting specific qualitative features of the data, introducing a degree of arbitrariness to the fitting procedure [18,25,3436].

The need for minimal clock parametrizations to address these issues has been recognized elsewhere in the literature. For example, recent Neurospora work has shown that it is possible to use a simple two-parameter function to represent all the intermediate processes between the expression of a clock gene and its action on a downstream target, while still maintaining sufficient flexibility to accurately simulate biological temperature and photoperiod responses [18,25].

An alternative technique for modelling GRNs is provided by Boolean logic. By assuming discrete expression levels, Boolean models provide an even greater reduction in complexity, albeit at the cost of reduced biochemical precision. Previous studies have exploited this reduction to study GRN state–space structures [48,59,60], and to introduce probabilistic models that parametrize statistical transitions between states [45,61]. The different limit cycles enumerable by a single Boolean GRN have been interpreted in some models as examples of cell differentiation [62,63], making limit cycle properties of special interest. Elsewhere, Boolean studies have focused on the regulatory logic underlying transcription [6466], the immune response [44,53] and apoptosis [49].

Here, we have presented the first models of biological clocks based on the Boolean formalism. We constructed logic versions of four DE models simulating circadian rhythms in the organisms N. crassa and A. thaliana. To derive the best logic representations, both regulatory structure (LC) and parameters (signalling delays and discretization thresholds) were optimized to synthetic experimental data generated from the DEs. In each case, we found that the LC yielding the best fit to data was consistent with that of the corresponding DE system. This suggested that the Boolean models were capable of identifying the particular combination of activating and inhibiting elements best able to account for a set of experimental expression patterns (figures 4 and and5).5). We confirmed this hypothesis by successfully applying our optimization algorithm to real Arabidopsis LUC data (figure 8). These results demonstrate that the logic models possess sufficient predictive power to perform biological network inference, despite their relative simplicity. Moreover, the optimization to LUC data highlights the importance of network structure in determining dynamic behaviour [67,68]. This optimization gave rise to a pair of putative 3-loop network architectures with distinct parameter values that generated very similar expression time series (cf. figure 9b,c and electronic supplementary material, figure S4). One of these architectures matched the 3-loop Arabidopsis DE model. It therefore contains the DE model's core LHY–TOC1 negative feedback loop, which was based on an experimental study that demonstrated the repression of TOC1 by LHY while also inferring the activation of LHY by TOC1 [69]. Significantly, the other architecture—which gave the best fit to data—agreed with more recent biochemical work showing that TOC1 is in fact a repressor of LHY [51]. This optimal architecture is also consistent with a complementary computational study showing that an expanded version of the 3-loop DE model incorporating a negative TOC1–LHY interaction yields more accurate simulations of TOC1 knockout and overexpression mutants [70].

We also found that although the cost function used to fit synthetic data only assessed goodness-of-fit in simulated 12:12 LD cycles, the logic models gave very good fits to the DEs in long and short days also, closely reproducing the relationships between photoperiod and expression phase (figure 7). The coordination of biochemical activity with the timing of dawn and dusk is a key system-level property that can be used to assess the relationship between the structure of a circadian network and its evolutionary flexibility [57]. Our photoperiod simulations indicate, somewhat surprisingly, that this property can be accounted for by significantly reduced dynamic models possessing much smaller parameter sets. In particular, the combined dawn- and dusk-tracking observed in the 2-loop Arabidopsis DE model (figure 7c) is accurately replicated by its logic formulation which has less than one-quarter the number of parameters. A similar reduction was observed for 3-loop Arabidopsis, with the 80 parameters describing the DE models decreasing to 20 in its Boolean equivalent (electronic supplementary material, table S1). The ability of simple, discrete models to simulate biologically realistic photoentraiment may also have important implications for the nature of circadian signal processing, in addition to being interesting from a modelling perspective. Specifically, it is consistent with hypotheses based on experimental data suggesting that the mechanism by which daylength and light intensity information is transmitted to output pathways may be partly digital in nature [71].

It should be noted that despite the success of the logic models in reproducing light responses, the coarser representation of network dynamics they provide means that the reproduction of certain circadian properties poses problems for the Boolean approach. A particular restriction of note relates to the modelling of temperature compensation. Temperature can be readily introduced into any DE model by following the methodology originally introduced by Ruoff [72,73]. In this scheme, the temperature dependence of a kinetic parameter is assumed to be governed by the Arrhenius relation. This leads to a simple balance equation for the corresponding activation energies and control coefficients which, when satisfied, guarantees a near-zero period-temperature derivative at the balance point. In a Boolean model, the parameters that determine the free-running period τFR are the discrete delays τj, and so simple balance equations can be derived from the form of the function determining the dependence of τFR on the τjs. For example, in the 2-loop Neurospora model, it can be shown that τFR = τ1 + τ2 + τ3 + τ4. Thus, compensation can be achieved simply by choosing the temperature dependence of each delay τj, so that their sum is constant at each point over the range of interest. However, as the τjs are generic parameters which in each instance summarize several biochemical processes, any balance equation derived in this manner will not in general be grounded in physical chemistry, reducing its biological relevance.

3.2. Computational advantages of the Boolean formulation

In addition to providing a more compact parametrization, a significant strength of the logic modelling approach is that it greatly reduces the complexity of the optimization procedure itself. This enables optimal configurations and delay–threshold combinations to be distinguished with a greater objectivity. The critical simplification afforded by the Boolean formulation is that the cost score is computed from bitstrings, which take the value 0 or 1, rather than data that vary over a much broader range of values. This removes the need for ad hoc cost function terms for normalizing the data and robustly computing period and phase. Indeed, the cost function we employed in this work was extremely simple, computing the correlation between the predictions of the model and the corresponding discretized experimental data at each vertex.

A second important simplification relates to the structure of the parameter space. In DE models, each interaction delay is dependent on processes that occur over a range of time scales (e.g. transcription, translation, degradation, etc.). This means that as well as being uncountably infinite, the parameter space can cover several different orders of magnitude, making it difficult to establish a priori bounds. In contrast, each signalling delay in a logic model is constrained to be a multiple of the data sampling interval tS, and is bounded above by the maximum possible time tMAX over which the corresponding interaction can occur. tMAX is itself bounded by the duration of the experimental measurements and can be further restricted on the basis of biological knowledge (here, the free-running period was used to bound all delays—see §4 for details). The set of possible signalling delays is thus finite. The set of possible discretization thresholds is also finite as varying the threshold for a given species will generate a set number of different binary time series, meaning it is only necessary to consider thresholds for which distinct bitstrings are obtained.

For logic models, the parameter space as a whole is therefore finite, and can be objectively bounded. Furthermore, for a fixed abstract topology, the set of all possible LCs is also finite; it is simply equal to the set of possible Boolean functions. This means that it is, in principle, possible to comprehensively search across all possible regulatory structures and parameter combinations to determine the best fit to data. Such a search is impossible with DE systems. Furthermore, searching across different patterns of activation and inhibition for a DE model is often problematic as it requires some a priori assumptions to be made regarding the underlying biochemical mechanisms; e.g. specification of which reactions may be cooperative and the ranges of the corresponding Hill coefficients [18,25,3436]. In practice, however, the number of possible network and parameter combinations in a Boolean formulation can become too large for a complete search with the computational resources available, making it necessary to constrain the optimization.

3.3. Refining the optimization protocol

Here, in order to ensure computational tractability, we subsampled the delay–threshold space, while also fixing the parameters controlling the impulse light inputs. In addition, we restricted ourselves to a subset of logic circuits by assuming that (i) multiple light inputs to a gene are combined with an OR gate and (ii) the net light signal directly modulates the state of the gene through either an AND or an OR gate, depending on the free-running light regime (figure 2). Of course not all possible circuits will be biologically reasonable ones (e.g. any circuit for which gates with light inputs uniformly output 0 or 1 would be unviable). Nonetheless, it is reasonable to expect that some interactions may be more accurately modelled by gates that are not of the simple AND/OR type [64], meaning that better performing LCs may have been overlooked.

In view of these restrictions, the fact that our optimal Boolean circuits match the dynamics of their target datasets, both synthetic and experimental, is very encouraging. Indeed, we anticipate that it should be possible to find circuit configurations and parameter sets giving good fits to data over a broader range of genetic and environmental perturbations. For example, our Boolean 3-loop Arabidopsis model shows consistency with much of the photoperiodic behaviour of its DE counterpart (figure 7d). However, it does not reproduce the phase response observed in the DE model as photoperiod is decreased, for which some components switch between dawn- and dusk-dominance in a complex manner [57].

A probable contributing factor is that our current optimization method involves computing the score at every parameter combination over a fixed lattice. This can be inefficient, particularly where the stoichiometry of an interaction and/or its molecular dynamics require the density of interacting species to accrue beyond a certain value before the interaction is statistically likely to occur. In such cases, the optimal threshold choice for the discretization is likely to be found within a narrow band of possibilities. Furthermore, high threshold resolutions can be required to resolve topological degeneracies in the cost function which make optimal networks difficult to distinguish (see the electronic supplementary material, §S1.2). The accuracy of the optimization could thus be increased by employing logic variants of global search methods capable of providing a more computationally efficient exploration of the parameter space, such as evolutionary algorithms [74]. This would increase the predictive power of the models, better positioning us to address features such as the dawn and dusk switching in shorter photoperiods observed in the 3-loop model.

3.4. Future directions

Finally, we note that there are many promising avenues for further developing the approaches introduced in this study. From a theoretical perspective, there is scope for extending established methods for analysing Boolean models. Of particular interest are techniques for determining the simplest inequalities involving linear combinations of the time delays that are required to traverse a given path through the state–space [44,52,53]. In the context of circadian models, this would involve developing general methods for deriving linear inequalities that result in free-running cycles with the target period. As the computational demands of optimization grow exponentially with the number of parameters to be fitted, restricting the search to the corresponding solution set would be expected to dramatically reduce the computational load.

More generally, hybrid logic/DE algorithms would see Boolean methods used firstly to determine the optimal model configuration from data, and then to identify the regions of parameter space within which an equivalent DE model is likely to give a good fit. This would exploit the ability of logic models, demonstrated in this work, to generate a coarse, but quantitative representation of the system dynamics in a systematic and efficient manner. We anticipate that regulatory networks of much greater complexity than those considered here could be quantitively modelled using such an approach.

4. Material and methods

4.1. Datasets used for model fitting

In order to incorporate some of the variability in expression levels characteristic of real experimental data, synthetic datasets were generated using the variant of Gillespie's stochastic simulation algorithm (SSA) introduced by Gonze et al. [39]. Fits to deterministic data obtained from direct integration of the DEs gave very similar results (data not shown). However, the stochastic datasets yielded a clearer separation between LCs, most likely owing to the lifting of topological degeneracies. For all circuits, final analyses were therefore restricted to the results obtained with stochastic time courses.

In generating the synthetic time courses, the scaling (or extensivity) parameter Ω was set close to the minimum value yielding self-sustained, unforced oscillations in each case. The final values used for simulations were 1-loop Neurospora, Ω = 25; 2-loop Neurospora, Ω = 50; 2- and 3-loop Arabidopsis, Ω = 1000. Time series were generated for five cycles in both entrained (LD 12:12) and free-running conditions (DD for the Neurospora models; LL for Arabidopsis). This choice mirrors the combination of experimental datasets typically chosen for parameter optimization in computational circadian studies [18,25,3436,38]. Time series were then subsampled every 0.5 h to give the data used for model fitting. Plots of the synthetic LD datasets are shown in the electronic supplementary material, figure S5.

Experimental gene expression levels were measured using LUC imaging, carried out as previously described [16]. Images were recorded using Hamamatsu C4742-98 digital cameras operating at −75°C under the control of Wasabi software (Hamamatsu Photonics, Hamamatsu City, Japan) with a sampling interval of 1.5 h. Bioluminescence levels were quantified using Metamorph software (MDS, Toronto, Canada). The resulting expression profiles were detrended for amplitude and baseline damping using the mFourfit function of BRASSv3 [35]. The detrended time series can be seen in figure 9a.

4.2. Implementing circadian networks in Boolean logic

For n biochemical species, let Xi(t) [set membership] {0,1} denote the activity of species i at time t, where t is a multiple, ktS, of the sampling interval tS. The update rule is then expressed as

equation image
4.1

Here, the Lks represent external light inputs to the circuit and τ1, … , τN+m are the signalling delays, with An external file that holds a picture, illustration, etc.
Object name is rsif20120080-i1.jpg for all i,j. Assuming a 24 h day and writing tDAWN and tDUSK for the times of dawn and dusk, respectively, the Lks are described by the following function:

equation image
4.2

Setting pk = 0 creates a constant darkness signal (DD); setting pk = 24 corresponds to constant light (LL); setting pk = tDUSKtDAWN yields a continuous LD cycle. Parameter sets with pk < tDUSKtDAWN yield LD cycles with a pulse of length pk at dawn.

The functions si are Boolean functions si : {0,1}n+m → {0,1} representing the interactions between genes and light inputs that determine the state of species Xi. We introduce a logical dependency function (or logic gate) G(·,g) to describe these interactions. This function takes a single parameter g [set membership] {0,1}, the value of which determines the type of reaction modelled. Formally, we consider two types of operator. The first operator acts on a single Boolean input Y [set membership] {0,1}, implementing either the identity or NOT gate, modelling activation and repression, respectively:

equation image
4.3

The second type of operator acts on two Boolean inputs Y,Z [set membership] {0,1}, implementing either the AND or the (inclusive) OR dependency. If either species can fulfil the interaction, an OR dependency is used. If both species are required in the interaction, then an AND dependency is used. Thus, for species Y [set membership] {0,1} and Z [set membership] {0,1}

equation image
4.4

The functions si in (4.1) are formed as compositions of these dependencies. For example, the update rule for gene TOC1 in the optimal 2-loop Arabidopsis model has the form

equation image
4.5

Using (4.3) and (4.4), this can be written as a composition of logic functions

equation image
4.6

with g6 = 0 and g1 = g8 = 1. The resulting Boolean function models a reaction in which LHY must be downregulated in order for Y to activate TOC1 production.

A given set of interaction functions si has an associated adjacency matrix A = (Aij) defined by

equation image
4.7

where L = (L1, … , Lm). Thus, Aij = 1 if species Xj can change the state of Xi, and so matrix A describes the abstract topology of a model (these topologies can be seen in figure 2). It follows that a circadian clock model is parametrized by its adjacency matrix A, a set of delays τ = (τ1, … ,τN+m) and a set of gates G = (g1, … , gd). Writing X = (X1, … , Xn), this yields the following compact, vectorized form of the update rule:

equation image
4.8

The set of gates G is the model's logic configuration (LC). As each configuration G = (g1, … , gd) is a bitstring, it can be represented uniquely by the corresponding decimal expansion D(G) defined as

equation image
4.9

For the models considered, the LCs are enumerated in terms of their decimal expansions for simplicity in figures 45 and 8.

Of the LCs consistent with A, a subset matches the pattern of activation and inhibition in the corresponding DE model. These are referred to as the DE LCs. For example, for the 2-loop Neurospora model, frq activates both isoforms of FRQ, giving g1 = g2 = 0, and both isoforms repress transcription, giving g3 = g4 = 1 (figures 1b and and22b). There are thus a pair of DE LCs in this case, 00110 and 00111, depending on whether both isoforms are required for repression (g5 = 0, corresponding to OR) or a single isoform is sufficient (g5 = 1, corresponding to AND). These DE LCs are encoded by the integers 6 and 7, respectively.

Finally, in order to construct the simplest logic models consistent with the general form of circadian DE systems, each Boolean function si is assumed to have the form

equation image
4.10

such that (i) Hi : {0,1}2 → {0,1} implements either the AND or OR gate; (ii) siF : {0,1}n → {0,1} encodes the structure of the free-running system; and (iii) siL : {0,1}m → {0,1} determines how multiple light inputs are integrated with siL(0, … , 0) = 0 and siL(1, … , 1) = 1. For consistency, it is therefore necessary that setting each light input Lk to the relevant constant value L recovers the free-running system (L = 0 for DD; L = 1 for LL). This condition is equivalent to Hi(siF(X),L) = siF(X). For the Neurospora models, where the free-running condition is DD, the consistency condition is achieved using an OR gate because x OR 0 = x. In the case of the Arabidopsis circuits, where the free-running condition is LL, the AND gate is appropriate as x AND 1 = x.

For both Arabidopsis models, the multiple light inputs to Y are combined using an OR gate: this is because in the corresponding DE models, both inputs can independently upregulate transcription [35,36]. In addition, the 3-loop model incorporates only the pulsed light input to the PRR gene because removing the continuous light input from the DE system had a negligible effect on its photoperiodic behaviour.

Full details of the logic formulation of each clock network are given in the electronic supplementary material, §S2.

4.3. Optimization and constraints

In order to identify the optimal combination of signalling delays τ = (τj), thresholds T = (Ti) and logic gates G = (gl) for a given model and dataset, we must introduce a cost function to be minimized. We used a simple function based on a correlation between the predicted time courses generated by the model and the corresponding data. Further details can be found in the electronic supplementary material, §S1.1. Optimal LC–parameter combinations are given in tables 1 and and22.

Table 1.
Optimal parameter sets: synthetic data. The logic configurations, G, delays, τj, and discretization thresholds, Ti (0 < Ti < 1), yielding the best fit of each logic model to synthetic time series. The values used for photoperiod ...
Table 2.
Optimal parameter sets: experimental data. The logic configurations (LCs), G, delays, τj, and discretization thresholds, Ti (0 < Ti < 1), yielding the top two fits of the 3-loop Arabidopsis logic model to experimental LUC time ...

The most comprehensive strategy in seeking to minimize the cost function is to systematically calculate the cost for all possible configurations and parametrizations of the model; that is all delays 0 ≤ τjtMAX, where tMAX is the longest time scale considered in the system, and all thresholds 0 < Ti < 1. However, when the model becomes too complex to permit a global analysis, it becomes necessary to introduce further constraints to restrict the parametrization to be considered. Such constraints should be as objective as possible.

One viable, objective constraint is to limit all the signalling delays around a closed loop so that they sum to no more than the period τFR of free-running oscillations in the target dataset.

For 1-loop Neurospora, this gives the single delay bound

equation image
4.11

where τFR = 22 h.

For 2-loop Neurospora, the constraint results in two delay bounds

equation image
4.12

where again τFR = 22 h.

For 2-loop Arabidopsis, there are three delay bounds

equation image
4.13

In this case, τFR = 25 h.

For 3-loop Arabidopsis, there are four delay bounds: those for the 2-loop model above together with an extra bound introduced by the addition of the LHY–PRR loop:

equation image
4.14

For this circuit, τFR = 24 h.

Two further constraints were introduced for reasons of computational tractability. Firstly, each delay τj was restricted to integer multiples of a minimum delay resolution τR, itself a multiple kτtS of the data sampling interval tS. Secondly, each threshold Ti was bounded within a subinterval [TMIN,TMAX] of [0,1], and also restricted to integer multiples of a minimum threshold resolution TR.

4.3.1. Optimization to synthetic data

For all models, TMIN and TMAX were set to the values 0.2 and 0.8, respectively. For the Neurospora models, kτ was set to 1 and TR to 0.025; for the Arabidopsis models, kτ was initially set to 2 and TR to 0.05. Scores were then recalculated with kτ = 1 and TR = 0.025, within intervals [τi − 3τR,τi + 3τR] and [Tj − 5TR, Tj + 5TR] centred around the parameter combinations giving the best scores. For Arabidopsis optimizations, light parameters controlling the impulse inputs (L1 and L3 in the 2-loop circuit; L1, L3 and L4 in the 3-loop circuit) were fixed at the values shown in table 1. These were determined from discrete approximations to the corresponding continuous curves in the DE models.

Each parameter set was then checked for a functioning circadian clock. Firstly, the solution to the free-running model was generated under the appropriate continuous light conditions, using the discretized data as initial conditions. If a limit cycle was obtained with period within 20% of the DE model, then this was fed into the model under simulated 12:12 LD cycles. Periodic, entrained oscillations were taken to indicate a viable clock circuit.

For the Neurospora and 2-loop Arabidopsis models, optimizations included all possible LCs. The 3-loop Arabidopsis optimizations were restricted to the subset of LCs defined by G = (10011011xyz), x,y,z [set membership]{0,1}; these are the configurations that result from fixing gates g1, … , g7 to their optimized values in the 2-loop circuit.

The photoperiod simulations shown in figure 7 were obtained by locally re-optimizing the logic circuits yielding viable clocks to simulations of the DE models under 12:12 LD cycles, and then calculating the maximum symmetric photoperiod interval (12 − PMAX, 12 + PMAX) over which both model formulations generated stable, entrained solutions. The cost was minimized with a simulated annealing algorithm [18,34,75], using a cost function for which the data and predicted time courses were taken to be single cycles of the entrained solution in the DE and logic models, respectively.

4.3.2. Optimization to experimental luciferase data

In fitting the 3-loop Arabidopsis model to the LUC data shown in figure 9a, genes were matched to model variables in the following manner: (i) CCA1 was equated to LHY on the basis that the LHY variable in the equivalent DE model groups together the effect of the two genes [36]; (ii) TOC1 was matched to TOC1; (iii) GI was matched to Y owing to experimental results showing that this gene can account for some of the action of Y in the DE model [36]; and (iv) PRR9 was matched to PRR, as this variable combines the PRR7 and 9 genes in the DE formulation [36].

Because the exact biological correlate of variable X is currently unknown, it was not used for costing the fit. Consequently, g2 and τ2 were both fixed at 0 (note that fixing g2 reduces the number of possible LCs from 2048 to 1024). This preserves the dynamics of all components in the model excluding X, together with the delay bounds used for constraining the parameter space. It also means that TOC1 and the dummy variable X have identical time series. Thus, in practice, the discretized TOC1 expression time series were used as a proxy for X data when calculating the cost at the LHY vertex. In addition, as the LUC data were measured in LL, the parameters p1p4 controlling the light inputs were all set to 24 (cf. equation (4.2)). In LL conditions, the absence of dawn and dusk means that the choice of the delay parameters τ9τ12 is arbitrary. We therefore set the values of these delays to 0. TMIN and TMAX were fixed at 0.2 and 0.8, respectively. kτ was initially set to 2 and TR to 0.2. Scores were then recalculated with kτ = 1 and TR = 0.05, within intervals [τi − 2τR,τi + 2τR] and [Tj − 2TR,Tj + 2TR] centred around the best-scoring parameter sets. The optimal parameter set for each LC was assessed to determine whether it yielded a viable clock by first generating the solution to the model using the discretized LUC time series as initial conditions, and then checking that this gave a limit cycle with free-running period within 20% of 24 h.

4.4. Software

The numerical routines for parameter optimization and model simulation were initially developed in Matlab (Mathworks, Cambridge, UK) and C. The scoring algorithms used for global parameter sweeps were subsequently converted into Java and run on a task farm computer consisting of 118 Intel Harpertown quad-core processors. All software used is available on request from the corresponding author.

Acknowledgments

The authors thank Kieron Edwards for making available the luciferase data used for model fitting and Declan Bates and Jonathan Fieldsend for helpful comments on the manuscript. SynthSys Edinburgh is a Centre for Integrative Systems Biology funded by BBSRC and EPSRC award no. D019621. This work was supported by the Wellcome Trust programme grant no. WT066784 (to P.G.), and has made use of resources provided by the Edinburgh Compute and Data Facility (http://www.ecdf.ed.ac.uk/) through the eDIKT initiative (http://www.edikt.org.uk). A.P. was funded by a BBSRC Summer Studentship and a Moray Endowment small grant awarded by the University of Edinburgh. The authors also acknowledge the help of Konrad Paszkiewicz and Robin Batten in accessing HPC facilities provided by the University of Exeter under its Systems Biology initiative.

References

1. Dunlap J. C., Loros J. J., DeCoursey P. J. 2003. Chronobiology: biological timekeeping. Sunderland, MA: Sinauer Associates.
2. Young M. W., Kay S. A. 2001. Time zones: a comparative genetics of circadian clocks. Nat. Rev. Genet. 2, 702–715 (doi:10.1038/35088576) doi: 10.1038/35088576. [PubMed] [Cross Ref]
3. Bell-Pedersen D., Cassone V. M., Earnest D. J., Golden S. S., Hardin P. E., Thomas T. L., Zoran M. 2005. Circadian rhythms from multiple oscillators: lessons from diverse organisms. Nat. Rev. Genet. 6, 544–556 (doi:10.1038/nrg1633) doi: 10.1038/nrg1633. [PMC free article] [PubMed] [Cross Ref]
4. Zhang E. E., Kay S. A. 2010. Clocks not winding down: unravelling circadian networks. Nat. Rev. Mol. Cell Biol. 11, 764–776 (doi:10.1038/nrm2995) doi: 10.1038/nrm2995. [PubMed] [Cross Ref]
5. Khapre R. V., Samsa W. E., Kondatov R. V. 2010. Circadian regulation of cell cycle: molecular connections between aging and the circadian clock. Ann. Med. 42, 1695–1700 (doi:10.3109/07853890.2010.499134) doi: 10.3109/07853890.2010.499134. [PubMed] [Cross Ref]
6. Gery S., Koeffler H. P. 2010. Circadian rhythms and cancer. Cell Cycle 9, 1097–1103 (doi:10.4161/cc.9.6.11046) doi: 10.4161/cc.9.6.11046. [PubMed] [Cross Ref]
7. Takeda N., Maemura K. 2010. Circadian clock and vascular disease. Hypertens. Res. 33, 645–651 (doi:10.1038/hr.2010.68) doi: 10.1038/hr.2010.68. [PubMed] [Cross Ref]
8. Westrich L., Sprouse J. 2010. Circadian rhythm dysregulation in bipolar disorder. Curr. Opin. Invest. Drugs 11, 779–787. [PubMed]
9. Lange T., Dimitrov S., Fehm H. L., Westermann J., Born J. 2006. Shift of monocyte function toward cellular immunity during sleep. Arch. Intern. Med. 166, 1695–1700 (doi:10.1001/archinte.166.16.1695) doi: 10.1001/archinte.166.16.1695. [PubMed] [Cross Ref]
10. Liu J., Malkani G., Mankani G., Shi X., Meyer M., Cunningham-Runddles S., Sun Z. S. 2006. The circadian clock period 2 gene regulates gamma interferon production of NK cells in host response to lipopolysaccharide-induced endotoxic shock. Infect. Immun. 74, 4750–4756 (doi:10.1128/IAI.00287-06) doi: 10.1128/IAI.00287-06. [PMC free article] [PubMed] [Cross Ref]
11. Keller M., Mazuch J., Abraham U., Eom G. D., Herzog E. D., Volk H. D., Kramer A., Maier B. 2009. A circadian clock in macrophages controls inflammatory immune responses. Proc. Natl Acad. Sci. USA 106, 21 407–21 412 (doi:10.1073/pnas.0906361106) doi: 10.1073/pnas.0906361106. [PubMed] [Cross Ref]
12. Ouyang Y., Andersson C. R., Kondo T., Golden S. S., Johnson C. H. 1998. Resonating circadian clocks enhance fitness in cyanobacteria. Proc. Natl Acad. Sci. USA 95, 8660–8664 (doi:10.1073/pnas.95.15.8660) doi: 10.1073/pnas.95.15.8660. [PubMed] [Cross Ref]
13. Dodd A. N., Salathia N., Hall A., Kévei E., Tóth R., Nagy F., Hibberd J. M., Millar A. J., Webb A. A. 2005. Plant circadian clocks increase photosynthesis, growth, survival, and competitive advantage. Science 309, 630–633 (doi:10.1126/science.1115581) doi: 10.1126/science.1115581. [PubMed] [Cross Ref]
14. Gardner G. F., Feldman J. F. 1981. Temperature compensation of circadian period length in clock mutants of Neurospora crassa. Plant Physiol. 68, 1244–1248 (doi:10.1104/pp.68.6.1244) doi: 10.1104/pp.68.6.1244. [PubMed] [Cross Ref]
15. Ruoff P., Rensing L., Kommedal R., Mohsenzadeh S. 1997. Modeling temperature compensation in chemical and biological oscillators. Chronobiol. Int. 14, 499–510 (doi:10.3109/07420529709001471) doi: 10.3109/07420529709001471. [PubMed] [Cross Ref]
16. Gould P. D., et al. 2006. The molecular basis of temperature compensation in the Arabidopsis circadian clock. Plant Cell 18, 1177–1187 (doi:10.1105/tpc.105.039990) doi: 10.1105/tpc.105.039990. [PubMed] [Cross Ref]
17. Brunner M., Diernfellner A. 2006. How temperature affects the circadian clock of Neurospora crassa. Chronobiol. Int. 23, 81–90 (doi:10.1080/07420520500545805) doi: 10.1080/07420520500545805. [PubMed] [Cross Ref]
18. Akman O. E., Locke J. C. W., Tang S., Carré I., Millar A. J., Rand D. A. 2008. Isoform switching facilitates period control in the Neurospora crassa circadian clock. Mol. Syst. Biol. 4, 164. (doi:10.1038/msb.2008.5) doi: 10.1038/msb.2008.5. [PMC free article] [PubMed] [Cross Ref]
19. Leloup J.-C., Gonze D., Goldbeter A. 1999. Limit cycle models for circadian rhythms based on transcriptional regulation in Drosophila and Neurospora. J. Biol. Rhythms 14, 433–448 (doi:10.1177/074873099129000948) doi: 10.1177/074873099129000948. [PubMed] [Cross Ref]
20. Smolen P., Baxter D. A., Byrne J. H. 2001. Modeling circadian oscillations with interlocking positive and negative feedback loops. J. Neurosci. 21, 6644–6656. [PubMed]
21. Sriram K., Gopinathan M. S. 2004. A two variable delay model for the circadian rhythm of Neurospora crassa. J. Theor. Biol. 231, 23–38 (doi:10.1016/j.jtbi.2004.04.006) doi: 10.1016/j.jtbi.2004.04.006. [PubMed] [Cross Ref]
22. Francois P. 2005. A model for the Neurospora circadian clock. Biophys. J. 88, 2369–2383 (doi:10.1529/biophysj.104.053975) doi: 10.1529/biophysj.104.053975. [PubMed] [Cross Ref]
23. Ruoff P., Loros J. J., Dunlap J. C. 2005. The relationship between FRQ-protein stability and temperature compensation in the Neurospora circadian clock. Proc. Natl Acad. Sci. USA 102, 17 681–17 686 (doi:10.1073/pnas.0505137102) doi: 10.1073/pnas.0505137102. [PubMed] [Cross Ref]
24. Hong C. I., Jolma I. W., Loros J. J., Dunlap J. C., Ruoff P. 2008. Simulating dark expressions and interactions of frq and wc-1 in the Neurospora circadian clock. Biophys. J. 94, 1221–1232 (doi:10.1529/biophysj.107.115154) doi: 10.1529/biophysj.107.115154. [PubMed] [Cross Ref]
25. Akman O. E., Rand D. A., Brown P. E., Millar A. J. 2010. Robustness from flexibility in the fungal circadian clock. BMC Syst. Biol. 4, 88. (doi:10.1186/1752-0509-4-88) doi: 10.1186/1752-0509-4-88. [PMC free article] [PubMed] [Cross Ref]
26. Goldbeter A. 1995. A model for circadian oscillations in the Drosophila period protein PER. Proc. R. Soc. Lond. B 261, 319–324 (doi:10.1098/rspb.1995.0153) doi: 10.1098/rspb.1995.0153. [PubMed] [Cross Ref]
27. Leloup J. C., Goldbeter A. 1998. A model for circadian rhythms in Drosophila incorporating the formation of a complex between the PER and TIM proteins. J. Biol. Rhythms 13, 70–87 (doi:10.1177/074873098128999934) doi: 10.1177/074873098128999934. [PubMed] [Cross Ref]
28. Ueda H. R., Hagiwara M., Kitano H. 2001. Robust oscillations within the interlocked feedback model of Drosophila circadian rhythm. J. Theor. Biol. 210, 401–406 (doi:10.1006/jtbi.2000.2226) doi: 10.1006/jtbi.2000.2226. [PubMed] [Cross Ref]
29. Smolen P., Hardin P. E., Lo B. S., Baxter D. A., Byrne J. H. 2004. Simulation of Drosophila circadian oscillations, mutations, and light responses by a model with VRI, PDP-1, and CLK. Biophys. J. 86, 2786–2802 (doi:10.1016/S0006-3495(04)74332-5) doi: 10.1016/S0006-3495(04)74332-5. [PubMed] [Cross Ref]
30. Forger D. B., Peskin C. S. 2003. A detailed predictive model of the mammalian circadian clock. Proc. Natl Acad. Sci. USA 100, 14 806–14 811 (doi:10.1073/pnas.2036281100) doi: 10.1073/pnas.2036281100. [PubMed] [Cross Ref]
31. Leloup J. C., Goldbeter A. 2003. Toward a detailed computational model for the mammalian circadian clock. Proc. Natl Acad. Sci. USA 100, 7051–7056 (doi:10.1073/pnas.1132112100) doi: 10.1073/pnas.1132112100. [PubMed] [Cross Ref]
32. Becker-Weimann S., Wolf J., Herzel H., Kramer A. 2004. Modeling feedback loops of the mammalian circadian oscillator. Biophys. J. 87, 3023–3034 (doi:10.1529/biophysj.104.040824) doi: 10.1529/biophysj.104.040824. [PubMed] [Cross Ref]
33. Mirsky H. P., Liu A. C., Welsh D. K., Kay S. A., Doyle F. J., III 2009. A model of the cell-autonomous mammalian circadian clock. Proc. Natl Acad. Sci. USA 106, 11 107–11 112 (doi:10.1073/pnas.0904837106) doi: 10.1073/pnas.0904837106. [PubMed] [Cross Ref]
34. Locke J. C. W., Millar A. J., Turner M. S. 2005. Modelling genetic networks with noisy and varied experimental data: the circadian clock in Arabidopsis thaliana. J. Theor. Biol. 234, 383–393 (doi:10.1016/j.jtbi.2004.11.038) doi: 10.1016/j.jtbi.2004.11.038. [PubMed] [Cross Ref]
35. Locke J. C. W., Southern M. M., Kozma-Bognar L., Hibberd V., Brown P. E., Turner M. S., Millar A. J. 2005. Extension of a genetic network model by iterative experimentation and mathematical analysis. Mol. Syst. Biol. 1, 2005.0013. (doi:10.1038/msb4100018) doi: 10.1038/msb4100018. [PMC free article] [PubMed] [Cross Ref]
36. Locke J. C. W., Kozma-Bognar L., Gould P. D., Fehér B., Kevei E., Nagy F., Turner M. S., Hall A., Millar A. J. 2006. Experimental validation of a predicted feedback loop in the multi-oscillator clock of Arabidopsis thaliana. Mol. Syst. Biol. 2, 59. (doi:10.1038/msb4100102) doi: 10.1038/msb4100102. [PMC free article] [PubMed] [Cross Ref]
37. Zeilinger M. N., Farré E. M., Taylor S. R., Kay S. A., Doyle F. J., III 2006. A novel computational model of the circadian clock in Arabidopsis that incorporates prr7 and prr9. Mol. Syst. Biol. 2, 58. (doi:10.1038/msb4100101) doi: 10.1038/msb4100101. [PMC free article] [PubMed] [Cross Ref]
38. Pokhilko A., Hodge S. K., Stratford K., Knox K., Edwards K., Thomson A. W., Mizuno T., Millar A. J. 2010. Data assimilation constrains new connections and components in a complex, eukaryotic circadian clock model. Mol. Syst. Biol. 6, 416. (doi:10.1038/msb.2010.69) doi: 10.1038/msb.2010.69. [PMC free article] [PubMed] [Cross Ref]
39. Gonze D., Halloy J., Goldbeter A. 2002. Biochemical clocks and molecular noise: theoretical study of robustness factors. J. Chem. Phys. 116, 10 997–11 010 (doi:10.1063/1.1475765) doi: 10.1063/1.1475765. [Cross Ref]
40. Gonze D., Halloy J., Goldbeter A. 2002. Robustness of circadian rhythms with respect to molecular noise. Proc. Natl Acad. Sci. USA 99, 673–678 (doi:10.1073/pnas.022628299) doi: 10.1073/pnas.022628299. [PubMed] [Cross Ref]
41. Vilar J. M., Kueh H. Y., Barkai N., Leibler S. 2002. Mechanisms of noise-resistance in genetic oscillators. Proc. Natl Acad. Sci. USA 99, 5988–5992 (doi:10.1073/pnas.092133899) doi: 10.1073/pnas.092133899. [PubMed] [Cross Ref]
42. Kauffman S. A. 1969. Metabolic stability and epigenesis in randomly constructed genetic nets. J. Theor. Biol. 22, 437–467 (doi:10.1016/0022-5193(69)90015-0) doi: 10.1016/0022-5193(69)90015-0. [PubMed] [Cross Ref]
43. Thomas R. 1991. Regulatory networks seen as asynchronous automata: a logical description. J. Theor. Biol. 153, 1–23 (doi:10.1016/S0022-5193(05)80350-9) doi: 10.1016/S0022-5193(05)80350-9. [Cross Ref]
44. Kaufman M., Andris F., Leo O. 1999. A logical analysis of T cell activation and anergy. Proc. Natl Acad. Sci. USA 7, 3894–3899 (doi:10.1073/pnas.96.7.3894) doi: 10.1073/pnas.96.7.3894. [PubMed] [Cross Ref]
45. Shmulevich I., Dougherty E. R., Kim S., Zhang W. 2002. Probabilistic Boolean networks: a rule-based uncertainty model for gene regulatory networks. Bioinformatics 18, 261–274 (doi:10.1093/bioinformatics/18.2.261) doi: 10.1093/bioinformatics/18.2.261. [PubMed] [Cross Ref]
46. Watterson S., Marshall S., Ghazal P. 2008. Logic models of pathway biology. Drug. Discov. Today 13, 447–456 (doi:10.1016/j.drudis.2008.03.019) doi: 10.1016/j.drudis.2008.03.019. [PubMed] [Cross Ref]
47. Yu L., Watterson S., Marshall S., Ghazal P. 2008. Inferring Boolean networks with perturbation from sparse gene expression data: a general model applied to the interferon regulatory network. Mol. Biosyst. 4, 1024–1030 (doi:10.1039/b804649b) doi: 10.1039/b804649b. [PubMed] [Cross Ref]
48. Saez-Rodriguez J., Alexopoulos L. G., Epperlein J., Samaga R., Lauffenburger D. A., Klamt S., Sorger P. K. 2009. Discrete logic modelling as a means to link protein signalling networks with functional analysis of mammalian signal transduction. Mol. Syst. Biol. 5, 331. (doi:10.1038/msb.2009.87) doi: 10.1038/msb.2009.87. [PMC free article] [PubMed] [Cross Ref]
49. Schlatter R., Schmich K., Avalos Vizcarra I., Scheurich P., Sauter T., Borner C., Ederer M., Merfort I., Sawodny O. 2009. ON/OFF and beyond: a Boolean model of apoptosis. PLoS Comput. Biol. 5, e1000595. (doi:10.1371/journal.pcbi.1000595) doi: 10.1371/journal.pcbi.1000595. [PMC free article] [PubMed] [Cross Ref]
50. Watterson S., Ghazal P. 2010. Use of logic theory in understanding regulatory pathway signaling in response to infection. Future Microbiol. 5, 163–176 (doi:10.2217/fmb.10.8) doi: 10.2217/fmb.10.8. [PubMed] [Cross Ref]
51. Gendron J. M., Pruneda-Paz J. L., Doherty C. J., Gross A. M., Kang S. E., Kay S. A. 2012. Arabidopsis circadian clock protein, TOC1, is a DNA-binding transcription factor. Proc. Natl Acad. Sci. USA 109, 3167–3172 (doi:10.1073/pnas.1200355109) doi: 10.1073/pnas.1200355109. [PubMed] [Cross Ref]
52. Thomas R. 1983. Logical description, analysis and synthesis of biological and other networks comprising feedback loops. Adv. Chem. Phys. 55, 247–282 (doi:10.1002/9780470142790.ch20) doi: 10.1002/9780470142790.ch20. [Cross Ref]
53. Kaufman M., Urbain J., Thomas R. 1985. Towards a logical analysis of the immune response. J. Theor. Biol. 114, 527–561 (doi:10.1016/S0022-5193(85)80042-4) doi: 10.1016/S0022-5193(85)80042-4. [PubMed] [Cross Ref]
54. Shmulevich I., Dougherty E. R., Zhang W. 2002. From Boolean to probabilistic Boolean networks as models of genetic regulatory networks. Proc. IEEE 90, 1778–1792 (doi:10.1109/JPROC.2002.804686) doi: 10.1109/JPROC.2002.804686. [Cross Ref]
55. Troein C., Locke J. C., Turner M. S., Millar A. J. 2009. Weather and seasons together demand complex biological clocks. Curr. Biol. 19, 1961–1964 (doi:10.1016/j.cub.2009.09.024) doi: 10.1016/j.cub.2009.09.024. [PubMed] [Cross Ref]
56. Tan Y., Dragovic Z., Roenneberg T., Merrow M. 2004. Entrainment dissociates transcription and translation of a circadian clock gene in Neurospora. Curr. Biol. 14, 433–438 (doi:10.1016/j.cub.2004.02.035) doi: 10.1016/j.cub.2004.02.035. [PubMed] [Cross Ref]
57. Edwards K. D., et al. 2010. Quantitative analysis of regulatory flexibility under changing environmental conditions. Mol. Syst. Biol. 6, 424. (doi:10.1038/msb.2010.81) doi: 10.1038/msb.2010.81. [PMC free article] [PubMed] [Cross Ref]
58. Ruoff P., Rensing L. 1996. The temperature-compensated Goodwin model simulates many circadian clock properties. J. Theor. Biol. 179, 275–285 (doi:10.1006/jtbi.1996.0067) doi: 10.1006/jtbi.1996.0067. [Cross Ref]
59. Ay F., Xu F., Kahveci T. 2009. Scalable steady state analysis of Boolean biological regulatory networks. PLoS ONE 4, e7992. (doi:10.1371/journal.pone.0007992) doi: 10.1371/journal.pone.0007992. [PMC free article] [PubMed] [Cross Ref]
60. Hau L. D., Kwon Y. K. 2011. The effects of feedback loops on disease comorbidity in human signaling networks. Bioinformatics 27, 1113–1120 (doi:10.1093/bioinformatics/btr082) doi: 10.1093/bioinformatics/btr082. [PubMed] [Cross Ref]
61. Krawitz P., Shmulevich I. 2007. Basin entropy in Boolean network ensembles. Phys. Rev. Lett. 98, 158701. (doi:10.1103/PhysRevLett.98.158701) doi: 10.1103/PhysRevLett.98.158701. [PubMed] [Cross Ref]
62. Ribeiro A. S., Kauffman S. A. 2007. Noisy attractors and ergodic sets in models of gene regulatory networks. J. Theor. Biol. 247, 743–755 (doi:10.1016/j.jtbi.2007.04.020) doi: 10.1016/j.jtbi.2007.04.020. [PubMed] [Cross Ref]
63. Garg A., Cara A. D., Xenarios I., Mendoza L., Micheli G. D. 2008. Synchronous versus asynchronous modeling of gene regulatory networks. Bioinformatics 24, 1917–1925 (doi:10.1093/bioinformatics/btn336) doi: 10.1093/bioinformatics/btn336. [PubMed] [Cross Ref]
64. Buchler N., Gerland U., Hwa T. 2003. On schemes of combinatorial transcription logic. Proc. Natl Acad. Sci. USA 100, 5136–5141 (doi:10.1073/pnas.0930314100) doi: 10.1073/pnas.0930314100. [PubMed] [Cross Ref]
65. Gerstung M., Timmer J., Fleck C. 2009. Noisy signaling through promoter logic gates. Phys. Rev. E 79, 011923. (doi:10.1103/PhysRevE.79.011923) doi: 10.1103/PhysRevE.79.011923. [PubMed] [Cross Ref]
66. Hunziker A., Tuboly C., Horváth P., Krishna S., Semsey S. 2010. Genetic flexibility of regulatory networks. Proc. Natl Acad. Sci. USA 107, 12 998–13 003 (doi:10.1073/pnas.0915003107) doi: 10.1073/pnas.0915003107. [PubMed] [Cross Ref]
67. Milo R., Shen-Orr S., Itzkovitz S., Kashtan N., Chklovskii D., Alon U. 2002. Network motifs: simple building blocks of complex networks. Science 298, 824–827 (doi:10.1126/science.298.5594.824) doi: 10.1126/science.298.5594.824. [PubMed] [Cross Ref]
68. Prill R., Iglesias P., Levchenko A. 2005. Dynamic properties of network motifs contribute to biological network organization. PLoS Biol. 3, e343. (doi:10.1371/journal.pbio.0030343) doi: 10.1371/journal.pbio.0030343. [PMC free article] [PubMed] [Cross Ref]
69. Alabadí D., Oyama T., Yanovsky M. J., Harmon F. G., Más P., Kay S. A. 2001. Reciprocal regulation between TOC1 and LHY/CCA1 within the Arabidopsis circadian clock. Science 293, 880–883 (doi:10.1126/science.1061320) doi: 10.1126/science.1061320. [PubMed] [Cross Ref]
70. Pokhilko A., Fernandez A. P., Edwards K. D., Southern M. M., Halliday K. J., Millar A. J. 2012. The clock gene circuit in Arabidopsis includes a repressilator with additional feedback loops. Mol. Syst. Biol. 8, 574. (doi:10.1038/msb.2012.6) doi: 10.1038/msb.2012.6. [PMC free article] [PubMed] [Cross Ref]
71. Dodd A. N., Love J., Webb A. A. 2005. The plant clock shows its metal: circadian regulation of cytosolic free Ca2+. Trends Plant Sci. 10, 15–21 (doi:10.1016/j.tplants.2004.12.001) doi: 10.1016/j.tplants.2004.12.001. [PubMed] [Cross Ref]
72. Ruoff P. 1992. Introducing temperature-compensation in any reaction kinetic oscillator model. Interdiscip. Cycle Res. 23, 92–99 (doi:10.1080/09291019209360133) doi: 10.1080/09291019209360133. [Cross Ref]
73. Ruoff P. 1994. General homeostasis in period and temperature-compensated chemical clock mutants by random selection conditions. Naturwissenschaften 81, 456–459 (doi:10.1007/BF01136649) doi: 10.1007/BF01136649. [Cross Ref]
74. Mendes P., Kell D. 1998. Non-linear optimization of biochemical pathways: applications to metabolic engineering and parameter estimation. Bioinformatics 14, 869–883 (doi:10.1093/bioinformatics/14.10.869) doi: 10.1093/bioinformatics/14.10.869. [PubMed] [Cross Ref]
75. Kirkpatrick S., Gelatt C. D., Vecchi M. P. 1983. Optimization by simulated annealing. Science 220, 671–680 (doi:10.1126/science.220.4598.671) doi: 10.1126/science.220.4598.671. [PubMed] [Cross Ref]

Articles from Journal of the Royal Society Interface are provided here courtesy of The Royal Society