A Compact Model for the Complex Plant Circadian Clock

The circadian clock is an endogenous timekeeper that allows organisms to anticipate and adapt to the daily variations of their environment. The plant clock is an intricate network of interlocked feedback loops, in which transcription factors regulate each other to generate oscillations with expression peaks at specific times of the day. Over the last decade, mathematical modeling approaches have been used to understand the inner workings of the clock in the model plant Arabidopsis thaliana. Those efforts have produced a number of models of ever increasing complexity. Here, we present an alternative model that combines a low number of equations and parameters, similar to the very earliest models, with the complex network structure found in more recent ones. This simple model describes the temporal evolution of the abundance of eight clock gene mRNA/protein and captures key features of the clock on a qualitative level, namely the entrained and free-running behaviors of the wild type clock, as well as the defects found in knockout mutants (such as altered free-running periods, lack of entrainment, or changes in the expression of other clock genes). Additionally, our model produces complex responses to various light cues, such as extreme photoperiods and non-24 h environmental cycles, and can describe the control of hypocotyl growth by the clock. Our model constitutes a useful tool to probe dynamical properties of the core clock as well as clock-dependent processes.


MODEL EQUATIONS
The time evolution of the mRNA and protein levels of the variables CL (CCA1/LHY), P97 (PRR9/PRR7), EL (ELF4/LUX), P51 (PRR5/TOC1), and of the activity of P (PIF3,PIL1) are governed by the following differential equations: The parameters L and D represent light and darkness, respectively. Their values are L = 1 and D = 0 during light phases, and L = 0 and D = 1 during dark phases.
Hypocotyl growth is controlled by the core clock via the protein PIF (PIF4/PIF5). The time evolution of the mRNA and protein levels of PIF, and the hypocotyl length are described by the following equations: The activations and inhibitions are modeled using Hill terms. The Hill coefficients are set to 2, based on multiple reports of clock genes forming both homodimers and heterodimers, with the evening complex being approximated as a heterodimer of LUX and ELF4 (Pokhilko et al. (2012);Fujiwara et al. (2008); Wang et al. (2010); Henriques and Mas (2013); Lu et al. (2009)). When multiple transcription factors act on the same target, several possible equations can be used to describe the resulting regulation.
Here, we used classical equations for mutually exclusive inhibitors, in which the effect of two inhibitors A and B is given by ) n , rather than the non-exclusive form used in previous models, in which A and B combine as Both forms were considered during the building of the model, and the results obtained were very similar, as can be seen in Supplementary Figure 1. We verified that the qualitative defects of the clock mutants and the light-and entrainment-related properties of the clock hold true regardless of the form of the equations. Although the distinction between competitive and non-competitive inhibition is important in the case of enzymatic kinetics, the use of the Hill equation in describing transcriptional regulation is empirical rather than based on kinetic laws, making the difference less relevant here.

SBML version
An SBML version of the model is available in the BioModels database (http://www.ebi.ac.uk/ biomodels) under the identifier number MODEL1601130000.

COST FUNCTION
In 8L:16D, the cost function required CCA1/LHY mRNA to peak around dawn in the wild type, assigning a score of 0 for a peak between ZT22 and ZT4. This condition was sufficient to ensure the proper phase of all the variables, therefore no additional constraints were added. No conditions were imposed on any of the mutants.
The wild type and all mutants were required to have a detectable rhythm in constant light, defined as all variables having a minimal value of 0.1, as well as a minimal difference of 10% between their minimum and maximum values. The wild type was also required to have a detectable rhythm in continuous darkness. Any solution that did not meet those criteria was considered to be arrhythmic. An arrhythmic wild type incurred a very large score penalty.
The free-running period was calculated by computing the average time between successive maxima of each variable. The wild type obtained a score of 0 for a period between 24 and 25 hours in continuous light, and between 25 and 28 hours in continuous darkness. The cca1/lhy, prr5/toc1, and weak elf4/lux mutants received a zero score if their period in continuous light was at least 1h shorter than the wild type, with no lower bound. The prr9/prr7 mutant received a zero score if its period was at least 1h longer than the wild type, with no upper bound.
Parameter optimization for the hypocotyl growth module was done independently of the main clock module. Only the wild type was used for this fitting. The system was simulated for 360 hours, then the hypocotyl length was reset to zero and another 120h simulation was perfomed. The PIF4/5 mRNA level was fitted to experimental data in 8L:16D, 12L:12D, and 16L:8D, obtained from the DIURNAL database. Relative hypocotyl lengths in the WT after five days in photoperiods ranging from 0 to 24h in increments of 3h were fitted to data from Niwa et al. (2009), shown in Figure 9.  It should be noted that the optimization process yielded several sets of parameters with a similar low score, none of which provide a perfect fit for every single aspect of the clock, which is not unexpected given the high degree of simplification of the model. Those parameter sets can be roughly divided into two groups depending on how well they reproduce the elf4/lux and prr5/toc1 mutants. One group correctly predicts the levels of other clock genes and the short period of prr5/toc1 but has low PRR9/PRR7 and high CCA1/LHY in elf4/lux mutants (whereas they should be high and low, respectively). The other group reproduces the clock defects of elf4/lux but predicts an almost unaffected prr5/toc1 mutant. Both are equally successful at modelling the cca1/lhy and prr9/prr7 mutants. The other clock properties presented here (entrainment, photoperiodic adaptation) are independent of the choice of parameter sets. The set presented here was chosen as a compromise between the two extremes.  Figure 1. Effect of the Hill term on the system. The experimental data (in black) is the same as presented in Figure 1 of the main text. The blue lines are simulated profiles using the "exclusive inhibitions" variant of the model, which is the one presented in the main text. The red curves are simulations of a "nonexclusive inhibitions" variant, using the same parameter values with no further optimization. (C) Sensitivity analysis. The bars represent the range of values for each parameter that produce a detectable free-running rhythm, expressed as a fold change from the optimal value provided in Supplementary Table  1.