In plants, as in animals, the core mechanism to retain rhythmic gene expression relies on the interaction of multiple feedback loops. In recent years, molecular genetic techniques have revealed a complex network of clock components in Arabidopsis. To gain insight into the dynamics of these interactions, new components need to be integrated into the mathematical model of the plant clock. Our approach accelerates the iterative process of model identification, to incorporate new components, and to systematically test different proposed structural hypotheses. Recent studies indicate that the pseudo‐response regulators PRR7 and PRR9 play a key role in the core clock of Arabidopsis. We incorporate PRR7 and PRR9 into an existing model involving the transcription factors TIMING OF CAB (TOC1), LATE ELONGATED HYPOCOTYL (LHY) and CIRCADIAN CLOCK ASSOCIATED (CCA1). We propose candidate models based on experimental hypotheses and identify the computational models with the application of an optimization routine. Validation is accomplished through systematic analysis of various mutant phenotypes. We introduce and apply sensitivity analysis as a novel tool for analyzing and distinguishing the characteristics of proposed architectures, which also allows for further validation of the hypothesized structures.
In recent years, molecular genetic techniques have revealed a complex network of components in the Arabidopsis circadian clock. Mathematical models allow for a detailed study of the dynamics and architecture of such complex gene networks leading to a better understanding of the genetic interactions. It is important to maintain a constant iteration with experimentation, to include novel components as they are discovered and use the updated model to design new experiments. This study develops a framework to introduce new components into the mathematical model of the Arabidopsis circadian clock accelerating the iterative model development process and gaining insight into the system's properties.
We used the interlocked feedback loop model published in Locke et al (2005) as the base model. In Arabidopsis, the first suggested regulatory loop involves the morning expressed transcription factors CIRCADIAN CLOCK‐ASSOCIATED 1 (CCA1) and LATE ELONGATED HYPOCOTYL (LHY), and the evening expressed pseudo‐response regulator TIMING OF CAB EXPRESSION (TOC1). The hypothetical component X had been introduced to realize a longer delay between gene expression of CCA1/LHY and TOC1. The introduction of Y was motivated by the need for a mechanism to reproduce the dampening short period rhythms of the cca1/lhy double mutant and to include an additional light input at the end of the day.
In this study, the new components pseudo‐response regulators PRR7 and PRR9 were added in negative feedback loops based on the biological hypothesis that they are activated by LHY and in turn repress LHY transcription (Farré et al, 2005; Figure 1). We present three iterations steps of model development (Figure 1A–C).
A wide range of tools was used to establish and analyze new model structures. One of the challenges facing mathematical modeling of biological processes is parameter identification; they are notoriously difficult to determine experimentally. We established an optimization procedure based on an evolutionary strategy with a cost function mainly derived from wild‐type characteristics. This ensured that the model was not restricted by a specific set of parameters and enabled us to use a large set of biological mutant information to assess the predictive capability of the model structure. Models were evaluated by means of an extended phenotype catalogue, allowing for an easy and fair comparison of the structures. We also carried out detailed simulation analysis of component interactions to identify weak points in the structure and suggest further modifications. Finally, we applied sensitivity analysis in a novel manner, using it to direct the model development. Sensitivity analysis provides quantitative measures of robustness; the two measures in this study were the traces of component concentrations over time (classical state sensitivities) and phase behavior (measured by the phase response curve). Three major results emerged from the model development process.
First, the iteration process helped us to learn about general characteristics of the system. We observed that the timing of Y expression is critical for evening light entrainment, which enables the system to respond to changes in day length. This is important for our understanding of the mechanism of light input to the clock and will add in the identification of biological candidates for this function. In addition, our results suggest that a detailed description of the mechanisms of genetic interactions is important for the systems behavior. We observed that the introduction of an experimentally based precise light regulation mechanism on PRR9 expression had a significant effect on the systems behavior.
Second, the final model structure (Figure 1C) was capable of predicting a wide range of mutant phenotypes, such as a reduction of TOC1 expression by RNAi (toc1RNAi), mutations in PRR7 and PRR9 and the novel mutant combinations prr9toc1RNAi and prr7prr9toc1RNAi. However, it was unable to predict the mutations in CCA1 and LHY.
Finally, sensitivity analysis identified the weak points of the system. The developed model structure was heavily based on the TOC1/Y feedback loop. This could explain the model's failure to predict the cca1lhy double mutant phenotype. More detailed information on the regulation of CCA1 and LHY expression will be important to achieve the right balance between the different regulatory loops in the mathematical model. This is in accordance with genetic studies that have identified several genes involved in the regulation of LHY and CCA1 expression. The identification of their mechanism of action will be necessary for the next model development.
We developed a mathematical model of the Arabidopsis circadian clock, including PRR7 and PRR9, which is able to predict several single, double and triple mutant phenotypes.
Sensitivity Analysis was used to identify the properties and time sensing mechanisms of model structures.
PRR7 and CCA1/LHY were identified as weak points of the mathematical model indicating where more experimental data is needed for further model development.
Detailed dynamical studies showed that the timing of an evening light sensing element is essential for day length responsiveness
Circadian clocks are present in many organisms (Harmer et al, 2001; Ditty et al, 2003; Stanewsky, 2003) allowing them to predict cyclic changes in the environment. These clocks generate oscillating rhythms with a period length of approximately 24 h, the time needed for one rotation of the earth around its axis. These rhythms can be entrained by external stimuli such as light and temperature changes, but persist even in their absence. Significant progress has been made recently in identifying key genes involved in the plant circadian clock using the model plant Arabidopsis thaliana.
Although the elucidation of the molecular mechanisms of the circadian oscillator in plants is still in its infancy, recent studies indicate that it is likely to be formed of multiple interacting feedback loops (Salome and McClung, 2004). In Arabidopsis, the first suggested regulatory loop involves the transcription factors CIRCADIAN CLOCK‐ASSOCIATED 1 (CCA1) and LATE ELONGATED HYPOCOTYL (LHY), and the pseudo‐response regulator TIMING OF CAB EXPRESSION (TOC1). CCA1 and LHY are partially redundant genes expressed in the morning that have been shown to inhibit TOC1 expression by binding to a specific site in its promoter (Harmer et al, 2000; Alabadi et al, 2001). In turn TOC1, which peaks at dusk, is necessary for the correct pattern of CCA1 and LHY expression (Alabadi et al, 2002; Mas et al, 2003a). The important effects of the pseudo‐response regulators PRR7 and PRR9 on the clock have been described only recently (Farre et al, 2005; Salome and McClung, 2005). While single mutations in either PRR7 or PRR9 have only small effects, the prr7prr9 double mutant results in an extremely long phenotype; the free running period is longer than 30 h. Furthermore, PRR7 and PRR9 are potentially involved in light input to the clock, either in a direct or an indirect way through LHY/CCA1 (Farre et al, 2005). Although these two genes appear to play partially redundant roles in the generation of stable oscillations, PRR7 and PRR9 have distinct roles in light perception by the circadian system.
Together with experimentation, modeling generates insight into the genetic interactions. Mathematical models allow for a detailed study of the dynamics and architecture of complex gene networks (Savageau, 1971; Glass and Kauffman, 1973; Novak and Tyson, 1997; Smolen et al, 2000). Computational models of circadian oscillators capable of reproducing basic characteristics of the experimentally observed rhythms have been described for various organisms in the past (Neurospora (Leloup et al, 1999; Ruoff et al, 2001); mammals (Forger and Peskin, 2003; Leloup and Goldbeter, 2003); Drosophila (Tyson et al, 1999; Smolen et al, 2004)). Models have been constructed initially with structures that account for desired clock characteristics with a minimal number of components. It is however important to include novel components that are later discovered into these first generation oscillators in the next modeling step. To accomplish this process and keep models consistent with experimental understanding, constant iteration with experimentation is needed and plays a significant role in the validity of mathematical models. Another challenge facing modelers is parameter estimation. Single parameters do not always have a direct correlation to a biochemical process and the measurement of reaction rates is either noisy or not possible. Therefore, broad ranges of values had to be derived from the literature with little supporting data available to the biological system of interest. This makes modeling highly dependent on the choice of an initial working parameter set. Recently, the first mathematical models for Arabidopsis have been proposed (Locke et al, 2005a, 2005b). A key contribution of that work was the application of optimization methods to estimate parameters that best account for a selection of clock characteristics.
There are several features that can be used to assess model fidelity. Plants show oscillations with a free‐running period of approximately 24 h. The specific molecular oscillations occur with distinct phase relationships to each other and to the light–dark cycle, and also maintain stable phase relationships in free‐running conditions. Their clocks are shown to be entrained to 24 h oscillations by input intervals of 24 h and to be phase responsive. The phase is altered through single external stimuli, depending on the current phase at the time of interference. In plants, the only modeled input for the clock so far is light, although the inclusion of temperature influence will be essential for modeling the plant circadian system in the future. Another important attribute is the behavior exhibited by mutations that affect the level or activity of genes involved in the generation of rhythms.
Our approach presents an iterative process of model identification, to incorporate new components and to systematically test different structural hypotheses. A model structure is built using biological hypotheses and parameters are identified via an optimization routine. Parameters are chosen to minimize a cost function, which consists of the minimal number of terms necessary to achieve primary clock characteristics. The optimization procedure searches for a broad minimum of the cost function by iteratively combining the best solutions found in parameter space and subsequently decreasing the searching range. Because this search covers global parameter space, we ensure that model failure results from its structure rather than from improperly chosen parameters. A working model will account not only for optimized characteristics but should also be predictive for the remaining model metrics, in particular for several wild type and mutant phenotypes (Tomlin and Axelrod, 2005). A model is validated systematically by simulating the desired features for which biological data are available. This procedure reveals the strengths and weaknesses of the analyzed structure and assists in the development of new hypotheses. Finally, we introduce and apply sensitivity analysis as a novel tool for analyzing and distinguishing the characteristics of the proposed model architectures. To illustrate the iteration process, we incorporate PRR7 and PRR9 into the existing model of the circadian oscillator of Arabidopsis.
This study highlights the importance of details in genetic interactions and light input for the main characteristics of the plant circadian rhythm. Insight into system dynamics identifies important functionalities of the hypothetical component Y and results in a model that predicts a significant number of clock features. We motivate further examination of PRR7 and Y mechanisms along with the exact pathways of light influence on the system.
The Prr7‐Prr9‐Y network
We used the interlocked feedback loop model published in Locke et al (2005a) as the basis for the model extension. Since there is still insufficient experimental data on CCA1 and LHY to differentiate between their modes of regulation and their activity, CCA1 and LHY were considered as one component, called LHY. PRR7 and PRR9 were added in negative feedback loops based on the biological hypothesis that they are activated by LHY and in turn repress LHY transcription (Farre et al, 2005), giving rise to the extended PRR7‐PRR9‐Y model in Figure 1A. The new model structure thus comprises four loops: the first developed or core loop consisting of TOC1, X and LHY, the second interlocked loop formed by TOC1, LHY protein and the component Y and two additional negative feedback loops with PRR7 and PRR9 interacting with LHY.
Model equations were set up as mass balances using nonlinear ordinary differential equations in the form of Michaelis–Menten and Hill kinetics (Supplementary information). An optimal parameter set was identified using the described optimization procedure. The applied evolutionary strategy (ES) resulted in broad minimal optima of the cost function and was not trapped by narrow ‘fox‐holes’, which assured high performance of the algorithm in the neighborhood of the optimal solution. The cost function was designed to fit the wild‐type characteristics of period in constant light and constant darkness, and the phases of circadian RNA expression of the system components. The strong period lengthening effect of the prr7prr9 double mutant was the only non‐wild‐type characteristic included in the optimization procedure. This allows the evaluation of numerous mutant phenotypes for validation and comparison of the models. Correct phenotype behavior, therefore, resulted from the model structure rather than a particular parameter choice. The developed model structure was able to reproduce the optimized characteristics and accounts for the overall wild‐type behavior (Table I). LHY and TOC1 RNA levels peak at their desired target ZT times and the period length in the prr7prr9 double mutant exhibits the experimentally observed long period phenotype. However, in contrast to experimental data, the period was longer in constant light conditions than in constant darkness (Table I).
To further investigate the properties of the model structure, the system was perturbed and expanded on multiple structural locations. Biological data on mutations of all the model components, with the exception of the hypothetical components X and Y, are available (Table I) and can be directly compared to in silico disturbances. This procedure is analogous to the analysis of technical systems in engineering sciences where parameter perturbations are used to classify the performance of control structures (Doyle et al, 1992; Levine, 1996). Through the mutant analysis, model properties are systematically evaluated and different structures can be easily compared.
Since the removal of TOC1 would open two of the four loops in the studied model structures, we simulated the toc1 RNAi line #65 (Mas et al, 2003a), which shows a daily average of 85% reduction in RNA expression levels in comparison to the wild‐type (data not shown). In addition to the double mutations cca1lhy (Alabadi et al, 2002; Mizoguchi et al, 2002) and prr7prr9 (Farre et al, 2005; Salome and McClung, 2005) two novel mutant combinations were introduced. The prr9toc1RNAi mutant perturbs three different loops at a time, the PRR9/LHY loop is removed, while the drastic reduction of the TOC1 RNA level simultaneously affects X, LHY and Y. The triple mutant prr7prr9toc1RNAi causes the disturbance of four distinct feedback loops: PRR7/LHY, PRR9/LHY, TOC1/X/LHY and TOC1/Y/LHY. Single and double mutants were simulated as total knockouts by removing the corresponding RNA and starting the protein levels at zero initial condition. The cca1/lhy single mutants were modeled by halving the translation rate of LHY since one of the two components is still present. The toc1RNAi line was simulated by reducing the RNA level to approximately 86% of its value under light–dark cycles. This was achieved by gradually reducing the parameter for the TOC1 transcription rate. ztl mutants were simulated by halving the two parameters for light‐dependent degradation of TOC1 protein.
The PRR7‐PRR9‐Y model accounted for the novel triple mutation prr7prr9toc1RNAi (Table I) and it correctly showed elevated levels of LHY in the prr7prr9 double mutant (Supplementary Figure 1; Farre et al, 2005). However, this model did not correctly predict any of the single mutant phenotypes (Table I). In addition, PRR7 peaked slightly earlier and PRR9 slightly later than observed experimentally leading to almost overlapping peaks although experimentally a 3–4 h difference between the circadian peaks is observed under light/dark cycles. Since LHY was the only activating factor of both PRR7 and PRR9, different transcription rates alone cannot result in the observed long delay. We hypothesized that one way of achieving an earlier rise of PRR9 RNA levels would be to incorporate the observation that PRR9 expression is strongly regulated by light (Farre et al, 2005; Ito et al, 2005).
Light influence on PRR9—the PRR7‐PRR9Light‐Y network
We introduced light‐dependent regulation of PRR9 expression to generate a new model structure (Figure 1B). Experimental data indicate a complex regulation of PRR9 transcription. First, the expression of PRR9 under light/dark cycles and the light induction of PRR9 expression seem to be almost exclusively dependent on LHY (Farre et al, 2005). Second, there is a very low expression of PRR9 in ethyolated seedlings (Makino et al, 2001). Third, the decrease in amplitude of PRR9 RNA oscillations between light/dark cycles and constant light conditions is much larger than for PRR7 (Figure 2A), indicating that an acute light induction mechanism may regulate PRR9 transcription. The PRR7‐PRR9Light‐Y model therefore uses three terms to describe PRR9 transcription:
The first term reflects the acute light responsiveness of PRR9, whereas constant light activation during the subjective day is realized through the second term. As described by Locke et al (2005a), acute light response is mediated by a component P that is both active and unstable in light. The last term allows for light‐independent activation through LHY to enable RNA oscillations in constant darkness.
The parameters identified in the optimization were capable of reproducing most of the optimized characteristics (Table I). For PRR9 expression, the model predicted oscillations with strongly reduced amplitude in constant light in comparison to light/dark cycles (Figure 2A). It also resulted in the correct phase and period for PRR9 RNA oscillations in constant light (Figure 2B). When simulating the shift from light/dark to free running conditions, the model missed one cycle; this reflected a longer dynamic transient than desired due to large amplitude changes, but did not indicate the model was incapable of showing the correct phase response to the change in light input. When the system had reached its new limit cycle (cyclic trajectory) in constant light it sustained oscillations again, which was confirmed by continuing the simulation over a longer time period (Figure 2B). Although the oscillation amplitude was not specified in the optimization, the model predicted the low amplitude oscillations in constant darkness shown by experimental data (Nakamichi et al, 2004).
Despite its failure to exhibit an earlier peak of PRR9 expression in light/dark cycles (Table I), the new model had profound implications for understanding the system's behavior. Whereas all previous models resulted in an extremely long period phenotype for toc1RNAi mutants (Table I), this model correctly predicted a short period in constant light conditions when TOC1 expression is significantly reduced (Table I). In addition, the model correctly predicted the phenotype of the novel triple mutant prr7prr9toc1RNAi (Table I, Figure 3). These results strongly enforce our assumption that light is highly capable of modulating the system's response. We conclude that this iteration was an important step in the model evolution and, consequently, in advancing our understanding of the Arabidopsis clock network.
Interestingly, whereas the interlocked model responds with a clear peak shift to changes in day length, neither the PRR7‐PRR9‐Y nor the PRR7‐PRR9Light‐Y model showed phase delays under long days. One function of a circadian clock is the phasing of physiological processes to particular times of the day or night. It also enables the adjustment of the phase of specific processes in response to the length of light or dark periods. This is described as the phase angle of entrainment (Dunlap et al, 2004). The detection of light signals in the evening likely mediates the phase responses caused by the extension or shortening of the photoperiod. In the interlocked model, evening light information is received by TOC1 via the hypothetical component Y. Using the knowledge gained through mathematical modeling we analyzed the behavior of Y and TOC1 in the extended model in detail.
Detailed analysis of Gene Y
Both TOC1 and Y displayed differences in their expression pattern in the PRR7‐PRR9Light‐Y model in comparison to the interlocked model. Since no restrictions or characteristics of Y were included in the optimization procedure, the resulting expression pattern of Y RNA remained flexible and was likely to change between the different model structures. In the PRR7‐PRR9Light‐Y model Y was expressed in a single morning peak (Figure 4B), indicating the fusion between the acute light induced peak and the circadian peaks observed in the interlocked model (Figure 4A). Interestingly, the peak of TOC1 expression occurred 2 h earlier than in the interlocked model (Table I), peaking just before dusk in a 12 h light:12 h dark photoperiod. To gain a deeper understanding of the influence of Y on TOC1, the rate of TOC1 transcriptional activation was plotted over time for both models. Figure 5 shows that both Y expression patterns resulted in a very similar activation curve, indicating that in both model structures the morning peak is the key determinate for the system response of TOC1 by counteracting the morning repression through LHY.
One hypothesis for the lack of the photoperiodic response is based on the observation that in the PRR7‐PRR9Light‐Y model, the TOC1 RNA peak preceded dusk, whereas in the interlocked model, TOC1 RNA levels kept rising until light is off. TOC1 activation already decreased before darkness set in and therefore might not be responsive to the lengthening of the photoperiod. To test this hypothesis, we optimized for a parameter set that provided a circadian peak of TOC1 at ZT 13.8. This modification did not achieve any changes in the model's photoperiodic response. However, detailed analysis of Y transcriptional activation terms showed an evolving separation of the acute and the circadian peak of Y RNA with a later peak of TOC1 (Figure 6). Thus, the timing of the TOC1 RNA peak influenced the circadian Y RNA peak but not the phase angle under longer or shorter photoperiods.
Another possible explanation for the lack of change in the phase of entrainment in the PRR7 and PRR9 containing models is the phase of Y itself. The expression of Y was strongly influenced by the dynamics of TOC1 and LHY, which were different between the two models (Figure 6). Phase responsiveness was observed when Y peaked in the later day. This indicated that the acute light induced peak and the circadian Y expression peaks have distinct functions. The acute peak controlled most of the TOC1 expression and the circadian peak modulated the fine‐tuning response to day length. The hypothetical Y component realizes a mechanism to enable entrainment by the light to dark transition via its light induction term. However, the expansion of the model to include PRR7 and PRR9 modified the structure in such a way that the expression of Y became dominated by its acute light response term. This led to a single morning peak of Y RNA determined by the dark to light transition; we hypothesized that this caused the lack of phase responsiveness to photoperiod. To test this hypothesis, we developed a new model structure by modifying the regulation of Y expression.
The PRR7‐PRR9Light‐Y′ network
In this new model structure, we modified the influence of light on Y (Figure 1C). In the previous models the expression pattern and function of Y were mostly determined by its acute light activation term. In order to give more weight to the later circadian light induced term of Y, we removed its acute light activation and renamed it Y′.
Indeed, this model structure led to a later peak of Y′ RNA expression at approximately ZT 6, a circadian peak matching the function of Y′ as an activator of the evening gene TOC1. It retained all wild‐type basic characteristics, and also correctly reproduced the effects of the toc1RNAi, the prr9toc1RNAi and prr7prr9toc1RNAi mutations (Table I). Whereas the previous models were unable to predict prr7 and prr9 single mutant phenotypes, the PRR7‐PRR9Light‐Y′ structure led to the correct long period phenotypes. This model was still unable to predict the phenotypes of cca1/lhy single and double mutants and loss of ZTL (in the model described as dark‐dependent degradation of TOC1 (Mas et al, 2003b), motivating further modifications.
Interestingly, in contrast to the extended model iterations analyzed so far, this model was capable of reproducing the correct periods in constant light conditions. It predicts a significantly longer period length (1 h) under free‐running conditions in constant darkness than in constant light. This is a characteristic of the Arabidopsis circadian clock that is likely caused by action of light on the system, indicating that the new model better reflects the in vivo mechanisms.
In addition, as hypothesized, the changes in the model structure led to the recovery of the photoperiodic phase response. This occurred although the TOC1 RNA peak preceded dusk (Figure 7A and Supplementary Figure 2). The basis of this response was the modulation of Y′ expression by day length. Y′ peaked in the middle of the day and shifted its phase angle depending on the hours of light (Figure 7B). There is no experimental evidence that TOC1 is directly activated by light but it still responds to the photoperiod, indicating that it is likely to be regulated by a factor modulated by light. Interestingly, Y′ was necessary for the maintenance of rhythms in this model; a Y′ mutation lost all oscillations after one cycle in constant light.
The simulations of knockout and RNAi experiments were numerical experiments introducing large perturbations to the system. Sensitivity analyses supplement these experiments by studying the effects of small parametric perturbations. At a given time, the state of the system is defined by the concentration of each of the components. Mathematically, the set of concentrations is referred to as the vector of state variables, or simply, states. Classical state sensitivity analysis measures the effects of infinitesimal perturbations upon the state values (Savageau, 1971; Varma et al, 1999). In mathematical modeling, it has been used to assist the parameter identification process to fit models to experimental data.
The complexity of biological systems allows them to maintain characteristic behaviors in spite of environmental disturbances, such as temperature fluctuations (Csete and Doyle, 2002). For example, the Arabidopsis plant maintains the proper phase relationship with the light/dark cycle despite the lengthening and shortening of daytime. Mathematical models faithfully representing the biological system, display a similar robustness. We used sensitivity analysis as a local, quantitative, measure of robustness; the system's behavior was robust if it showed relative insensitivity to most parametric perturbations (Kitano, 2002).
The two measures in this study were the traces of component concentrations over time and phase behavior (measured by the phase response curve). Because parameters with low sensitivity supported the argument that the system was robust, careful attention was paid only to parameters with high sensitivity. High sensitivity may indicate that the parameter was meant to convey input, that it was (or should have been) regulated by other system components, or simply that our model structure was incorrect (Morohashi et al, 2002). For example, we expected phase behavior to have moderately high sensitivity with respect to the parameter representing light because light was a control input. Another possibility is that high sensitivity reflects a strong phenotype in the biological system. For example, we expect genes with arrhythmic loss of function phenotypes to be sensitive components in the mathematical model. In addition to making statements about the robustness of a system to one particular set of parameters and a single network structure, it has been shown that sensitivity rankings have been preserved across varying parameter sets and similar network structures (Stelling et al, 2004). We therefore performed a systematic sensitivity analysis of our three model iterations to study the influence of the structural changes.
Sensitivity analysis of model structure
Classical, or state, sensitivities were calculated to examine the general robustness properties, capturing the combined response of shape, phase, period and amplitude of the oscillations. In addition to the classical sensitivities, we computed state variable responses with more global information by computing sensitivities to large parametric perturbations, and to a large collection of parameters sets (for details see the Materials and methods section). The three methods produced like results (data not shown), strengthening statements about overall parameter sensitivities. We rank‐ordered the parameters from those with the greatest effect upon the system to those with the least. Grouping parameters according to their biochemical function provided a suitable means for examining the clock (Stelling et al, 2004). In addition to classifying parameters by function, we grouped parameters by component (state). To isolate the analysis of model structure from its interaction with its environment (via light input), we disregarded light as a parameter.
The functional analysis for all three models identified RNA degradation as the most sensitive process of the system (Supplementary Figure 3). One explanation used the relation of these kinetic parameters to housekeeping cell regulatory processes (Stelling et al, 2004). Whereas parameters specific to the circadian system were designed to tolerate disturbances (such as changes in input signals), parameters shared with other processes were not expected to undergo significant variations. However, recent reports counter that explanation by suggesting the existence of specific mechanisms for the degradation of circadian regulated RNAs (Gutierrez et al, 2002; Lidder et al, 2005). In this case, the high sensitivity of these parameters could be attributed to lack of regulatory detail in the model.
Parameters influencing X RNA but not X protein consistently showed the highest sensitivity in all three models (Figure 8). This was also observed for the interlocked model (data not shown). X had been introduced in the interlocked model as an activator of LHY expression (Locke et al, 2005a) to produce the delay between the peaks of expression of TOC1 and LHY observed in the experimental data. One interpretation for this result is the lack of detail in the description of X RNA regulation. This is consistent with the observation that X is not essential for oscillations in any of the three studied models. An essential function would lead to a high sensitivity for X protein, as observed for Y. This, however, is inconsistent with experimental data showing that genes involved in the regulation of CCA1 and LHY expression, such as ELF3, ELF4 and LUX are essential for rhythmicity (Hicks et al, 2001; Doyle et al, 2002; Hazen et al, 2005). Since X was not behaving like the known regulators of LHY, it indicates that the model could be improved through more detailed treatment of the regulation of LHY/CCA1 expression.
We observed that in all three models, PRR7 was more sensitive than PRR9. Since both single mutants showed only small effects on the system behavior, it is likely that this high sensitivity, particularly that of PRR7 RNA, is due to a lack of information about PRR7 RNA regulation. It has been shown experimentally that the expression of PRR7 is not as strongly influenced by CCA1 and LHY as PRR9 transcription (Farre et al, 2005), indicating that PRR7 is regulated by an additional mechanism yet to be discovered and introduced into the model.
Comparison of the parameter sensitivity rankings revealed distinct differences between the analyzed model structures. The PRR7‐PRR9‐Y and PRR7‐PRR9Light‐Y models displayed low sensitivity to the parameters associated with Y RNA levels. This is consistent with the observation that in these structures Y RNA did not respond to external light stimuli, such as changes in day length. The higher sensitivity of Y protein in the model, however, indicates that Y played an essential role. In the PRR7‐PRR9Light‐Y′ model, Y′ RNA was sensitive, reflecting its ability to be modified by light input. This high sensitivity together with its essential role for rhythmicity in constant light for all three models, suggests that Y/Y′ reflects a more complex regulatory mechanism. No activators of TOC1 have been identified so far using mutant screens, indicating that several redundant or partially redundant factors may be involved in this process.
Parameters associated with TOC1 displayed an increased sensitivity in the PRR7‐PRR9Light‐Y′ model compared to those in the previous models. This is likely to be the result of the combination of two effects. First, the high sensitivity of TOC1 RNA coincided with its ability to change phase depending on the light conditions. Second, loss of TOC1 led to a strong circadian phenotype (Table I) consistent with the experimental results (Strayer et al, 2000; Mas et al, 2003a).
Parameters determining the interactions between TOC1 and Y′ were among the five most sensitive parameters. This reflects a system highly dependent on the TOC1/Y′ feedback loop that is able to sense changes in the environment. The observation that the simulation of a cca1lhy double mutant does not lead to changes in the oscillations is consistent with a system structure based on TOC1/Y′. Since experimental data show that cca1lhy mutants have severely impaired rhythms under constant conditions (Alabadi et al, 2002; Mizoguchi et al, 2002), we expect the model to reflect that behavior with a relatively high sensitivity for LHY. This indicates that further details about LHY regulation and activity are necessary for a more accurate representation of the system.
Sensitivity analysis of model response to light
Light was introduced to the system via the parameter ld. This section studies the influence of ld on the behavior of the PRR7‐PRR9Light‐Y and PRR7‐PRR9Light ‐Y′ models, elucidating the fundamental differences between their respective phase behaviors. The former is strongly controlled by the light transition at dawn, while the latter is controlled by both dawn and dusk.
We analyzed the classical state sensitivities to perturbations of ld. The PRR7‐PRR9light‐Y model showed an extremely sensitive response to perturbations of ld, making ld the most sensitive parameter for all states (data not shown). This was not true for PRR7‐PRR9Light‐Y′, whose components could be categorized as those sensitive to light (LHY, PRR7, PRR9 and X protein) and those relatively insensitive to light (TOC1, Y′ and X RNA) (data not shown). Since, in contrast to TOC1 and Y, LHY is highly affected by the external light stimulus, the effects of light are mainly channeled through LHY.
To compare the different timing effects in detail, we turned from classical state sensitivity measures to those specific to oscillatory systems (Kramer et al, 1984; Ingalls, 2004; Rand et al, 2004). An impulse phase response curve (IPRC), derived from phase sensitivity (Kramer et al, 1984), represents the size and direction of the phase shift that occurs when the system undergoes an impulse perturbation. A curve describing phase shifts due to parametric perturbations is a parametric IPRC, while one describing phase shifts due to perturbations in state values is a state IPRC. Both are plotted as relative measures to ensure fair comparison between different parameters/states (e.g. a state IPRC reflects the phase shift due to a change in state value measured as a percentage of its amplitude).
We computed the parametric IPRC's for both models in constant darkness. As expected from the results of the state sensitivity analysis, changes in ld caused significant phase shifts in the PRR7‐PRR9Light‐Y model. These shifts were much larger than those caused by pulses in other parameters, and were outside the region of experimentally observed behavior (data not shown). In fact, the system restarted at ZT 0 whenever a light impulse was introduced (Supplementary Figure 4). In contrast, the PRR7‐PRR9Light‐Y′ model's range of phase resetting was in accordance with the observed behavior (Figure 9A). The ld IPRC for the PRR7‐PRR9Light‐Y′ model shared the basic characteristic of experimentally obtained PRCs (Covington et al, 2001), showing phase delays in the subjective day and phase advances in the subjective night. This effectively produces a shift towards the closest light‐on transition. It is interesting to notice that this computed PRC more closely resembled a Type I (weak) resetting pattern than a Type 0 (strong) resetting pattern, similar to the one reported for the interlocked model (Locke et al, 2005a).
The parameters with the largest influence on phase in the PRR7‐PRR9Light‐Y′ model were parameters characterizing the interaction of TOC1 and Y′ (Figure 8C). We compute the state IPRC in constant darkness in order to further investigate the indicated influence of TOC1/Y′ (Figure 9B). The state IPRC highlights the states with the largest effect on the circadian phase. It reveals a striking sensitivity to TOC1 nuclear protein, and a relatively high sensitivity to TOC1 and Y′, in general. In contrast to the ld parametric IPRCs (Figure 9A), state IPRCs related to TOC1 and Y′ (Figure 9B) describe phase advances in the subjective day and phase delays in the subjective night indicating that these states are sensing the light‐to‐dark transition. The modulation of Y′ and therefore TOC1 by light leads to a flexible evening sensor regulated by the circadian clock, whereas morning sensing is mediated by external light input on LHY. The computed parametric IPRC mirrors these two effects and their relative importance, while the state IPRC confirms the circadian regulated mechanism depending on the circadian peaks of TOC1 and Y′ rather than parameters.
The presented model results in a stronger evening entrainment than morning entrainment, although the Arabidopsis plant as a light‐sensing organism is observed to show a stronger resetting to the morning (Millar, 2003). This may be due to the uneven balance between the LHY/TOC1 and the TOC1/Y′ loop described above.
We used mathematical modeling tools to improve the iterative model development process to include recently identified genes into the model of the plant clock. In particular, the model was extended to include the two pseudo‐response regulators PRR7 and PRR9. Three iteration steps were performed to identify a description of the network structure leading to the best reproduction of experimental results. A phenotype catalogue was used to systematically validate the model development. We show that detailed description of the interaction between network components and between components and light stimuli had profound effects on the system's performance. In addition, sensitivity analysis was used to highlight the impact of disturbances on the robustness properties of the system.
All three developed models resulted in oscillations and circadian peaks for LHY, TOC1, PRR7 and PRR9 genes that correlate well with the experimental data. The parameter optimization procedure additionally succeeded in reproducing the long period phenotype of the prr7prr9 double mutation. The final PRR7‐PRR9Light‐Y′ model predicted a wide range of phenotypes including a novel prr7prr9toc1RNAi triple mutation, light entrainment and correct phase responsiveness to photoperiodic changes. Although a few phenotypes were still not correctly modeled, detailed analysis of the system response and parameter sensitivity analysis enabled the identification of the weak points of the model structure. These correlate well with gaps in the description of the biological system indicating where new experimental data are needed most.
Two specific results indicate that a better description of the regulation of PRR7 expression is necessary. First, the desired delay between PRR7 and PRR9 RNA expression observed under light/dark cycles could not be achieved solely by introducing light induction of PRR9, indicating that it might be caused by an effect on PRR7 expression. Second, although in the studied models the loss of function of PRR7 leads only to a weak phenotype, parameters associated with PRR7 transcription are highly sensitive. This suggests that a more detailed description of the regulation of PRR7 expression is necessary. As observed for PRR9, the incorporation of these details in the model is likely to affect the overall system response.
The simulation of a ztl mutation by halving the dark‐dependent TOC1 protein degradation rates did not lead to any phenotype in constant light (Table I), although the period in constant darkness was lengthened (period in constant darkness +0.8 h). Interestingly, although the prediction was not quantitatively precise, it reflects the experimental observation of ztl mutants having a stronger phenotype in constant darkness than in constant light (Somers et al, 2004). There are two compatible explanations for these effects. First, the model fails to weight the dark‐dependent term sufficiently. In agreement with this hypothesis, although the model correctly simulated TOC1 protein oscillations in both constant light and darkness, it did not reflect the lower average level in constant dark and the higher average level in constant light observed experimentally (Mas et al, 2003b). As noted by Mas et al (2003b), under light/dark cycles, TOC1 degradation begins at the end of the dark period suggesting a more complex regulation mechanism than a simple dark induced protein degradation.
Second, the simulation of a ztl mutation by reducing only the dark‐dependent degradation term of TOC1 protein represents an oversimplification. Experimental data indicate that in the ztl mutant, background TOC1 protein fails to cycle in both constant light and constant darkness (Mas et al, 2003b) suggesting that ZTL is not only involved in the dark‐dependent degradation of TOC1 but also in the regulation of TOC1 degradation under constant light. Since the simulation of halving the dark‐independent degradation rate of TOC1 protein and reducing the dark‐dependent degradation by one‐fourth, however, led to the correct tendency of a longer period in constant light conditions (period in constant light +0.4 h, data not shown), the model endorses this dual role for ZTL.
In spite of correctly predicting several phenotypes, none of the developed models was able to predict the cca1/lhy mutant phenotype (Table I). This contrasts with the capacity of the interlocked model to reproduce this phenotype (Table I; Locke et al, 2005a). This discrepancy reflects the differences between the modeling strategies presented in this work and the one used by Locke et al (2005a). In order to increase our prediction capacity, we chose a minimal number of optimization characteristics and did not include the double mutant phenotype in the optimization procedure as had been done for the development of the interlocked model. Since the proposed model structures failed to reproduce this phenotype, we carried out a detailed analysis in order to determine the reasons for this failure and suggest further modifications.
The basic structure common in all the studied models results in a system that is primarily regulated by the TOC1/Y(Y′) feedback loop, and, therefore, does not account for the cca1/lhy single or double mutants (Table I). The long period phenotype in the prr7prr9 double mutant in the model as in the experimental data was caused by the elevated levels of LHY/CCA1 expression (Farre et al, 2005; Supplementary Figure 1). This shows that the system was able to react to changes in LHY levels and therefore expected to achieve at least a qualitative prediction of the cca1lhy double mutant. Sensitivity analyses (Figure 8) highlighted the fact that the model relies heavily on the TOC1/Y′ feedback loop, which is able to cycle in the absence of LHY. One possible explanation might be that CCA1 and LHY are not as completely redundant as they have been modeled, thus a double mutation in planta is likely to eliminate a more complex structure than what has been modeled so far. It has also been recently reported that the CCA1 and LHY expression is differentially regulated (Gould et al, 2006). This additional level of complexity still needs to be better defined experimentally to be added to the model.
The large difference in amplitude observed for LHY and the PRR's between light–dark cycles and constant light conditions is a further indicator for a lack of information on LHY activation. The induction of LHY under entrained conditions relies heavily on the acute light response at the light to dark transition. This leads to a more severe attenuation in amplitude when transferring the system from light–dark cycles (acute effect) to constant light (no acute effect) than experimentally observed. PRR7 and PRR9, relying on LHY as their only activator, are affected by this amplitude attenuation and result in low amplitude oscillations themselves. The balance between the acute induction of LHY and other factors will change once a better description of LHY regulation is available. Together with further progress in the light signaling pathway, this will significantly affect the amplitude effects between entrained and free‐running conditions.
The component X as the activator of LHY was poorly described by the model leading to the high sensitivity value for parameters associated with X expression. There are several candidate genes that peak in the evening/night that are good candidates for X such as ELF3, ELF4 and LUX. Presently, it is difficult to hypothesize a mechanism for their activity because they all share the basic characteristics of being necessary for LHY/CCA1 expression and causing arrhythmicity under constant light conditions when mutated (Hicks et al, 1996; Covington et al, 2001; Hicks et al, 2001; Doyle et al, 2002; Hazen et al, 2005; Onai and Ishiura, 2005). Since all candidates are essential for rhythmicity, this might indicate that the current model lacks a significant part of the necessary structure. The essential role of X, however, will be tightly linked to the achievement of the correct prediction for a cca1lhy double mutation.
In summary, like most circadian models studied so far (Leloup et al, 1999; Tyson et al, 1999; Ruoff et al, 2001; Forger and Peskin, 2003; Leloup and Goldbeter, 2003; Smolen et al, 2004), our model does not predict every single characteristic of the system, but it created a framework where new components can be easily added and their effects analyzed. It is interesting to note that the correct phenotypic characteristics are kept through several model iterations indicating that the model development is proceeding in the correct direction.
We have also shown how system dynamical analysis provides insight into particular functionalities and guides model development as exemplified by the analysis of Y/Y′. Modeling allows for the visualization of particular rates and timing, for example, of transcriptional activation or inhibition, exploiting the actual advantage offered through mathematical modeling. This type of analysis led to the revelation of distinct functions for Y. Y/Y′ is essential for oscillations as an activator of TOC1 and transmitter of light signals. Mutations in Y/Y′ quickly lead to arrhythmicity in constant light, but not in constant darkness. Y/Y′ expression in the later day was the key for realizing a photoperiodic response and reproducing the correct periods in free running conditions, two important characteristics of the plant circadian clock achieved by the final model.
Our results suggest that two light‐sensing mechanisms are necessary to achieve the phase entrainment and the free‐running periods by regulating the effect of light on phase advances and delays. Evening sensing was not realized by a mechanism that detects the light to dark transition in a way analogous to the acute light responsive protein, but by components that track continuous light at the end of the day, in our model TOC1 and Y′. This is consistent with the idea that plants share parametric and nonparametric entrainment mechanisms (Millar, 2003). It is known that a complex interaction of photoreceptor signaling pathways regulate the light input to the clock (Millar et al, 1995; Devlin and Kay, 2000; Millar, 2003). In addition, the circadian clock in plants is able to gate the light input, it receives and it has been shown that ELF3 plays a role in this mechanism (McWatters et al, 2000; Covington et al, 2001).The Arabidopsis circadian clock models, built so far, do not include a gating factor. The PRR7‐PRR9light‐Y′ model was still able to achieve a PRC similar to the ones measured experimentally (Covington et al, 2001; Locke et al, 2005a). This might be due to the fact that light input on Y is limited by TOC1 repression at the end of the day and the unrepressed light input on PRR9 has little influence on the overall system since it modulates the weaker LHY feedback loop. The resulting model, thus, does not need additional gating yet. More and detailed experimental data on the hypothesized light mechanisms, however, will be indispensable in order to assess this part of the model in the future.
Sensitivity analysis in the Arabidopsis system has confirmed the relation between structure and robustness properties developed in the models for Drosophila circadian clocks (Stelling et al, 2004). The resulting sensitivity rankings reflect applied changes in the model structure in the shape of TOC1 and influence of Y, but also conserve structural similarities, the sensitivity of X and PRR7 RNA for instance. By validating models and, more important, by identifying model weaknesses, sensitivity analysis represents a new tool to guide iterative model development. PRC analysis identified phase‐sensing structures and will be essential in approaching a better representation of the entrainment mechanisms in future modeling.
The model of the circadian clock in Arabidopsis thaliana is still far from incorporating the entire data set provided by genetics. The goal is to be able to readily incorporate new experimental data in order to achieve a good predictive model of the studied system. The method and tools proposed herein improve this procedure by efficiently using the possibilities of mathematical modeling in a systematic approach. We have determined the strengths and weaknesses of the current model giving indication to further experiments necessary for its improvement. Finally, we show how mathematical modeling helps visualize genetic mechanisms that can explain general characteristics of circadian clocks such as entrainment.
Materials and methods
Models were implemented in MATLAB (Mathworks, Cambridge, UK) using Michaelis–Menten and Hill kinetics for transcription and degradation rates (see Supplementary information). As there is little data available to determine reaction rates explicitly, parameter fitting has become an important issue in mathematical modeling. Our model parameters were estimated applying an ES that minimizes a cost function, mainly based on wild‐type characteristics. As the modeling goal of this study was to incorporate PRR7 and PRR9, the long period effect of the prr7prr9 double mutant was included in the cost function. The 24‐h LD period is the most weighted term, then phase relationships of identified genes to light–dark cycles with more effort on the circadian peaks of TOC1 and LHY. Periods in constant light, constant darkness and the double mutant are equally weighted. The broadness of a peak is used as an indicator for degradation rates, the amplitude size, to assure detectable oscillations. The least weighted terms form standard deviations for oscillation sizes and peak phases. By mainly basing the optimization on wild‐type behavior and accounting for the correct period dimension in constant light conditions, the systematic phenotype analysis in constant light can be used to effectively evaluate the proposed model structure.
The optimization algorithm was initialized with oscillatory solutions resulting from a random search of 10 000 parameter sets in parameter space. Rhythms were demanded for free‐running conditions and light/dark cycles. The ES is an established optimization technique that is based on the principles of evolution and adaptation (Weicker, 2002). It is distinguished from the evolutionary algorithm, by using strategy parameters to adapt the step size or mutation strength during the routine, the so‐called self‐adaptation. By using a population‐based algorithm, attraction to local minima is reduced. The optimal parameter set found in the ES is further refined using a Hill procedure, a local optimizer that tracks the minimum in the environment of the starting point.
The random search was carried out by running the MATLAB code on a cluster, consisting of 47 × 2.8 GHz CPUs, using the Distributed Computing Toolbox. The MATLAB code for the optimization was compiled into C and executed on a 2.39 GHz Intel Pentium 4 CPU.
Parameter sensitivities measure the quantitative impact of perturbations in system parameters on characteristic system properties. A single parameter often captures a system of underlying reactions and in this case the sensitivity reflects the properties towards the underlying control mechanisms. The circadian clock model is a dynamical system that is represented by ordinary differential equations of the form
Classical state sensitivities along a specific state trajectory are defined by
Three different kinds of sensitivities are calculated:
Classical sensitivities, where S(t) contains the sensitivity coefficients
One‐dimensional sensitivities: scalar perturbationsFor each trial, one parameter is perturbed up by 40%. After the system has reached the new limit cycle, classical sensitivity coefficients are computed. The process is repeated with the parameter perturbed down by 40%.
Monte Carlo multi‐dimensional sensitivities: vector perturbationsA new parameter set is chosen by allowing each of the nominal values to vary 5% up, 5% down or to remain unchanged. After the limit cycle determined by this set parameter set is reached, classical sensitivity coefficients are computed. This is repeated 10 000 times, excluding all non‐oscillatory cases.
Classical sensitivities measure local properties, whereas one‐dimensional or multidimensional sensitivities allow for more globally valid conclusions. The three measures were well‐correlated, indicating that the robustness properties of the chosen parameter set are not locally unique, but rather relate to the model structure. All sensitivity rankings are computed using the BioSens toolkit (www.chemengr.ucsb.edu/~ceweb/faculty/doyle/biosens/BioSens.htm).
IPRCs predict the shift in phase due to a perturbation in either a parameter or a state (Taylor, Doyle III, Petzold, unpublished material). Phase is labeled ϕ and measures time. In the Y′ model LHY RNA peaks 1 h after dawn, so we associate ϕ=1 with the time at which LHY RNA peaks. If a system undergoes a phase advance of size Δϕ, then, once the system has settled, the peak of LHY RNA will occur 1−Δϕ hours after dawn.
The state IPRC at time t predicts a phase shift due to an instantaneous change in the value of state xi at time t:
and the relative IPRC (shown in the figures) is
The parametric IPRC measures a phase shift due to an infinitesimally small perturbation in parameter p that lasts an infinitesimally short amount of time:
and the relative IPRC (shown in the figures) is
unless p=0 (in which case we use the IPRC).
We computed the IPRC's using Matlab code.
The Arabidopsis thaliana prr7‐3 and prr9‐1 T‐DNA insertion lines (http://signal.salk.edu; prr7‐3 is SALK_030430, prr9‐1 is SALK_07551, here called prr7 and prr9, respectively), and the homozygous line TOC1RNAi #65 with the CCR2::luc reporter (Mas et al, 2003a) were used to generate the double mutant prr9TOC1RNAi and the triple mutant prr7prr9TOC1RNAi. The triple mutant was generated by crossing the double‐mutant prr7‐3prr9‐1 (Farre et al, 2005) to the TOC1RNAi line. All lines are in the Colombia‐0 (Col) background. Homozygous F4 lines were used for the circadian rhythms analysis. Seedlings were grown on Murashige and Skoog medium (Murashige and Skoog, 1962) with 0.8% agar and 3% sucrose for 6 days in light/dark cycles (12 h light, 12 h dark; 70 μmol m−2 s−1). Bioluminescence rhythms of single seedlings under constant light conditions were analyzed as previously described (Millar et al, 1995). Period length was estimated using FFT‐NLLS program (Millar et al, 1995; Plautz et al, 1997).
Expression analysis by reverse transcriptase‐mediated PCR
RNA extraction, cDNA synthesis and real‐time PCR detection were performed as described in Farre et al (2005). The gene IPP2 (isopentenyl:pyrophosphate:dimethyllallyl pyrophosphate isomerase, AT3G02780) was used as normalization control. The primers for PRR9 detection were described in Farre et al (2005). The following primers and probe were used for the detection of PRR7 expression:
Note added in proof
In a simultaneous publication in Molecular Systems Biology, Locke et al also describe an extension of the interlocked model to a three‐loop network by adding PRR7/PRR9 and further investigate their suggestion of GIGANTEA as a candidate for Y (Locke et al, 2006).
We thank Linda Petzold, Neda Bagheri, Ghislain Breton, Rudiyanto Gunawan, Sam Hazen and the Kay lab members for helpful comments and useful discussions. This work was supported by the Institute for Collaborative Biotechnologies through Grant DAAD19‐03‐D‐0004 from the US Army Research Office (FJD), IGERT NSF grant DGE02‐21715 (FJD), by NIH Grant #GM56006 (SAK), by NIH Grant #GM067837 (SAK), and by a postdoctoral fellowship from Human Frontier Science Program to EMF. This is manuscript # 18214‐BC of The Scripps Research Institute.
Supplementary Text [msb4100101-sup-0001.pdf]
Supplementary Figures [msb4100101-sup-0002.pdf]
Supplementary Table [msb4100101-sup-0003.xls]
Supplementary model 1 [msb4100101-sup-0004.xml.txt]
Supplementary model 2 [msb4100101-sup-0005.xml.txt]
Supplementary model 3 [msb4100101-sup-0006.xml.txt]
- Copyright © 2006 EMBO and Nature Publishing Group