Robustness of Circadian Clocks to Daylight Fluctuations: Hints from the Picoeucaryote Ostreococcus tauri

The development of systemic approaches in biology has put emphasis on identifying genetic modules whose behavior can be modeled accurately so as to gain insight into their structure and function. However, most gene circuits in a cell are under control of external signals and thus, quantitative agreement between experimental data and a mathematical model is difficult. Circadian biology has been one notable exception: quantitative models of the internal clock that orchestrates biological processes over the 24-hour diurnal cycle have been constructed for a few organisms, from cyanobacteria to plants and mammals. In most cases, a complex architecture with interlocked feedback loops has been evidenced. Here we present the first modeling results for the circadian clock of the green unicellular alga Ostreococcus tauri. Two plant-like clock genes have been shown to play a central role in the Ostreococcus clock. We find that their expression time profiles can be accurately reproduced by a minimal model of a two-gene transcriptional feedback loop. Remarkably, best adjustment of data recorded under light/dark alternation is obtained when assuming that the oscillator is not coupled to the diurnal cycle. This suggests that coupling to light is confined to specific time intervals and has no dynamical effect when the oscillator is entrained by the diurnal cycle. This intringuing property may reflect a strategy to minimize the impact of fluctuations in daylight intensity on the core circadian oscillator, a type of perturbation that has been rarely considered when assessing the robustness of circadian clocks.

of the internal clock that orchestrates biological processes over the 24-hour diurnal cycle have been constructed for a few organisms, from cyanobacteria to plants and mammals.
In most cases, a complex architecture with interlocked feedback loops has been evidenced.
Here we present first modeling results for the circadian clock of the green unicellular alga Ostreococcus tauri. Two plant-like clock genes have been shown to play a central role in Ostreococcus clock. We find that their expression time profiles can be accurately reproduced by a minimal model of a two-gene transcriptional feedback loop. Remarkably, best adjustment of data recorded under light/dark alternation is obtained when assuming that the oscillator is not coupled to the diurnal cycle. This suggests that coupling to light is confined to specific time intervals and has no dynamical effect when the oscillator is entrained by the diurnal cycle. This intringuing property may reflect a strategy to minimize the impact of fluctuations in daylight intensity on the core circadian oscillator, a type of perturbation that has been rarely considered when assessing the robustness of circadian clocks.

Author Summary
Circadian clocks keep time of day in many living organisms, allowing them to anticipate environmental changes induced by day/night alternation. They consist of networks of genes and proteins interacting so as to generate biochemical oscillations with a period close to 24 hours. Circadian clocks synchronize to the day/night cycle through the year principally by sensing ambient light. Depending on the weather, the perceived light intensity can display large fluctuations within the day and from day to day, potentially inducing unwanted resetting of the clock. Furthermore, marine organisms such as microalgae are subjected to dramatic changes in light intensities in the water column due to streams and

Introduction
Real-time monitoring of gene activity now allow us to unravel the complex dynamical behavior of regulatory networks underlying cell functions [1]. However, understanding the collective behavior of even a few molecular actors defies intuition, as it depends not only on the topology of the interaction network but also on strengths and response times of its links [2]. A mathematical description of a regulatory network is thus necessary to qualitatively and quantitatively understand its dynamical behavior, but obtaining it is challenging. State variables and parameters are subject to large fluctuations [3], which create artificial complexity and mask the actual network structure. Genetic modules are usually not isolated but coupled to a larger network, and a given gene can be involved in different modules and pathways [4]. It is thus important to identify gene circuits whose dynamical behavior can be modeled quantitatively, to serve as model circuits.
One strategy for obtaining such circuits has been to construct synthetic networks, which are isolated by design [5][6][7]. As recent experiments have shown, an excellent quantitative agreement can be obtained by incorporating when needed detailed descriptions of various biochemical processes (e.g., multimerization, transport, DNA looping, etc.) [7].
Another strategy is to study natural gene circuits whose function makes them rela-tively autonomous and stable. The circadian clocks that drive biological processes around the day/night cycle in many living organisms are natural candidates, as these genetic oscillators keep track of the most regular environmental constraint: the alternation of daylight and darkness caused by Earth rotation [8][9][10][11]. Informed by experiments, circadian clock models have progressively become more complex, evolving from single loops featuring a self-repressed gene [12,13] to networks of interlocked feedback loops [14][15][16][17].
Here we report surprisingly good agreement between the mathematical model of a single transcriptional feedback loop and expression profiles of two central clock genes of Ostreococcus tauri. This microscopic green alga is the smallest free-living eukaryote known to date and belongs to the Prasinophyceae, one of the most ancient groups of the green lineage. Ostreococcus displays a very simple cellular organization, with only one mitochondrion and one chloroplast [18,19]. Its small genome (12.6 Mbp) sequence revealed a high compaction (85% of coding DNA) and a very low gene redundancy [20] (e.g., most cyclins and CDK are present as a single copy gene [21]).The cell division cycle of Ostreococcus is under control of a circadian oscillator, with cell division occurring at the end of the day in light/dark cycles [21]. These daily rhythms in cell division meet the criteria characterizing a circadian clock, as they can be entrained to different photoperiods, persist under constant conditions and respond to light pulses by phase shifts that depend on internal time [21].
Very recently, some light has been shed on the molecular workings of Ostreococcus clock by Corellou et al. [22]. Since the clock of closely related Arabidopsis has been extensively studied, they searched Ostreococcus genome for orthologs of higher plant clock genes and found only two, similar to Arabidopsis central clock genes Toc1 and Cca1 [22]. These two genes display rhythmic expression both under light/dark alternation and in constant light conditions. A functional analysis by overexpression/antisense strategy showed that Toc1 and Cca1 are important clock genes in Ostreococcus. Overexpression of Toc1 led to increased levels of CCA1 while overexpression of Cca1 resulted in lower levels of TOC1.
Furthermore CCA1 was shown to bind to a conserved evening element sequence (EE) that is required for the circadian regulated activity of Toc1 promoter. Whether Toc1 and Cca1 work in a negative feedback loop could not be inferred from this study since Ostreococcus clock appeared to rely on more than a simple Toc1 /Cca1 negative feedback loop.
Interestingly, Arabidopsis genes Toc1 and Cca1 were the core actors of the first plant clock model, based on a transcriptional loop where TOC1 activates Cca1 and the similar gene Lhy, whose proteins dimerize to repress Toc1 [23,24]. However, this model did not reproduce well expression peaks of Toc1 and Cca1 in Arabidopsis [24] and was extended to adjust experimental data [25]. Current Arabidopsis clock models feature several interlocked feedback loops [15,16]. This led us to investigate whether the transcriptional feedback loop model where Toc1 activates Cca1 and is repressed by Cca1 would be relevant for Ostreococcus.
We not only found that this two-gene loop model reproduces perfectly transcript profiles of Ostreococcus Toc1 and Cca1 but that excellent adjustment of data recorded under light/dark alternation is obtained when no model parameter depends on light intensity. This counterintuitive finding suggests that the oscillator is not permanently coupled to light across the 24-hour cycle but only during specific time intervals, which is supported by numerical simulations. In this article, we propose that the invisibility of coupling in entrainment conditions reflects a strategy to shield the oscillator from natural fluctuations in daylight intensity.

Experimental data and model adjustment
To characterize the temporal pattern of Toc1 and Cca1 expression in Ostreococcus, we used microarray data acquired in triplicate under 12:12 light/dark cycle, as described in [21] (Fig. 1). One Toc1 and two Cca1 mRNA time courses had no aberrant point.
Here, we use as target profiles the complete Toc1 profile and the complete Cca1 profile whose samples are obtained from the same microarray data as the Toc1 profile. We checked that the results described in this work are robust to the biological variations observed. Corellou et al. have also carried out an extensive work of genetic transformation in Ostreococcus, leading to transcriptional and translational fusion lines allowing one to monitor transcriptional activity and protein dynamics in living cells [22]. However, luciferase kinetics in this organism is still not well known and we postpone the analysis of luminescence time series to a future work. Model adjustment has thus been carried out using microarray expression data, which reflect accurately the endogeneous levels of mRNA. Although seeking quantitative agreement with luminescence time series was premature at this stage, predicted protein concentration profiles were compared with data from translational fusion lines as an additional test. A minimal mathematical model of the two-gene feedback loop comprises four ordinary differential equations (Eq. (2), Methods) with 16 parameters. Since detailed models extending the basic 4-ODE model (2) could only have led to better adjustment, we purposely neglected here effects such as compartmentalisation or delays due to transcription or translation so as to minimize the risk of overfitting and reliably assess the validity of the two-gene loop hypothesis.
Experimental data are recorded under 12:12 Light/Dark (LD) alternation so that the coupling which synchronizes the clock to the diurnal cycle must be hypothesized.
Circadian models usually assume that some parameters depend on light intensity (e.g., a degradation rate increases in the dark), and thus take different values at day and night. Parameter space dimension then increases by the number of modulated parameters.
Various couplings to light were considered, with 1 to 16 parameters depending on light intensity. We also tested adjustment to model (2) with all parameters constant, which allowed us to quantify the relevance of coupling mechanisms by measuring the difference between best-fitting profiles in the coupled and uncoupled cases.
The free-running period (FRP) of the oscillator in constant day conditions was fixed at 24 hours, which was the mean value observed in experiments [22], but we checked that our main results remain valid for other values of the FRP. In fact, we found that when FRP was freely adjustable, it usually converged to values close to or slightly below 24 hours. Fixing the FRP at exactly 24 hours is interesting in that coupling mechanisms are selected by adjustment only if they improve goodness of fit and not merely to achieve frequency locking.
A free-running model adjusts experimental data The first result is that an excellent agreement between numerical and experimental profiles is obtained, with a root mean square (RMS) error of a few percent (Figs. 2(A)-(B)).
There is no point in extending model (2) to improve adjustment of microarray data, which are compatible with the hypothesis of a Toc1 -Cca1 feedback loop. Moreover, the corresponding protein profiles (not adjusted) correlate well with luminescence signals from CCA1:Luc and TOC1:Luc translational fusion lines (Figs. 2(C)-(F)).
But the more surprising is that a non-coupled model, where all parameters are kept constant, adjusts experimental data ( Fig. 2(B), RMS error 3.6%) essentially as well as a fully coupled model where all parameters are allowed to vary between day and night ( Fig. 2(A), RMS error 3.3%). The corresponding parameter values are given in Table 1. When only one or a few parameters were modulated, goodness of fit significantly degraded compared to the uncoupled and fully coupled cases. This indicates that besides being biologically unrealistic, the model with all parameters modulated fits data merely because of its large parameter space dimension, and cannot be considered seriously. Moreover we simulated the transition from LD alternation to constant light (LL) or constant darkness (DD) conditions for this model and found that it still adjusted experimental data well in LL while displaying strongly damped oscillations in DD (Fig. S1). This confirms that adjustment relies on time profiles being close to free-running oscillator profiles and that adjustment by a fully coupled model is in fact accidental.
On the other hand the uncoupled model is equally unrealistic because it cannot be entrained to the day/night cycle, whereas it is observed experimentally that upon a phase shift of the light/dark cycle, CCA1 and TOC1 expression peaks quickly recover their original timings in the cycle. To verify that adjustment by a free-running oscillator model does not depend on the target profile used, we generated a large number of synthetic profiles whose samples where randomly chosen inside the interval of variation observed in biological triplicates, and adjusted a free-running oscillator model to them. In each case, we found that although RMS error slightly degraded compared our target profile (where mCCA1 and mTOC1 samples for a given time always come from the same microarray), it remained on average near 10 %, with visually excellent adjustment (Fig. S2). Last, it should be noted that assuming a FRP of 24 hours allows frequency locking to occur without coupling, but cannot induce best adjustment in this limiting case by itself.
Thus the paradoxical result that data points fall almost perfectly on the temporal profiles of a free-running oscillator is counterintuitive but must nevertheless be viewed as a signature of the clock architecture. As we will see, this in fact does not imply that the oscillator is uncoupled but only that within the class of models considered so far, where parameters of the TOC1-CCA1 loop take day and night values, the uncoupled model is the one approaching experimental data best. Nothing precludes that there are more general coupling schemes that adjust data equally well.
Before unveiling such models, we discuss now whether the simple negative feedback loop described by model (2) is a plausible autonomous gene oscillator. With two transcriptional regulations, it is a simpler circuit than the Repressilator, where three genes repress themselves circularly [5]. It is known that in this topology, oscillations become more stable as the number of genes along the loop increases. The two-gene feedback loop described by (2) could therefore seem to be a less robust oscillator than the Repressilator, and thus a poor model for the core oscillator of a circadian clock.
To address this issue, we checked robustness of adjustment with respect to parameter variations. We found that the experimental profiles can be reproduced in a wide region of parameter space around the optimum, which is quite remarkable given the simplicity of the model (Fig. S3). Moreover, a distinctive feature of the best fitting parameter sets is a strongly saturated degradation, in particular for Cca1 mRNA, with an extremely low value of K M C equal to 0.6% of the maximal CCa1 mRNA concentration (see table 1). In this situation, the number of molecules degraded per unit time is essentially constant and does not depend on the concentration except at very small values. This is consistent with the characteristic sawtooth shape of our target profile drawn in linear scale ( Fig. 1(B)).
The role of post-translational interactions in gene oscillators and circadian clocks has been recently emphasized (see, e.g., [26,27]), and in particular saturated degradation has since long been known to favor oscillations [9,28,29]. Recently, it has been been shown to act as a delay [30,31] and to be essential for inducing robust oscillations in simple synthetic oscillators [7,32,33] (compare Fig. 1(B) with Fig. 5 of [33]). Thus, strongly saturated degradation is very likely also a key dynamical ingredient of the natural gene oscillator studied here.

Adjustment by a model with gated coupling
Circadian models are usually coupled to diurnal cycle by changing some parameter values between day and night [12][13][14][15][16][17]. This assumes that all molecular actors involved in light input pathways have been incorporated and that their properties (e.g., degradation rates) react directly to light. Such couplings act over the entire cycle except when light-sensitive actors are present only transiently. For example, models of Arabidopsis clock feature an intermediary protein PIF3 that is necessary for induction of CCA1 by light but is shortly degraded after dawn so that CCA1 transcription is only transiently activated [15,24,25].
Gating of light input has been observed in several circadian clocks and may be important for maintaining proper timing under different photoperiods [34].
In our case, light/dark alternation has no detectable signature in the dynamics of Toc1 and Cca1 mRNA when the clock is phase-locked to the diurnal cycle. This suggests that the actors of the two-gene loop do not sense light directly, and are driven via unknown mediators, which modify their properties inside specific temporal intervals. Since the input pathway can have complex structure and dynamics, possibly featuring separate feedback loops, the windows of active coupling may be located anywhere inside the diurnal cycle and reflect light level at other times of the cycle. Coupling activation should depend both on time of day and on the intrinsic dynamics of the light input pathway, notwithstanding a possible feedback from the circadian core oscillator [35][36][37].
For simplicity, we restrict ourselves to models in which some parameters of the TOC1-CCA1 feedback loop are modified between two times of the day, measured relatively to dawn (ZT0). The start and end times of coupling windows are then model parameters instead of being fixed at light/dark transitions. This assumes that the input pathway tracks diurnal cycle instantaneously, without loss of generality for understanding behavior in entrainment conditions. In this scheme, resetting of the two-gene oscillator can be studied by simply shifting the oscillator phase relatively to the coupling windows. The results so obtained will be sufficient to show that there exist coupling schemes which leave no signature on mRNA profiles, and to study their properties.
What makes our approach original is not the gated coupling to diurnal cycle, which can be found in other models, but the fact that we do not try to model the actors of the input pathway, which can be complex. This is because we focus here on the TOC1-CCA1 feedback loop, which mostly behaves as an autonomous oscillator. Thus we only need to know the action of the unknown mediators on TOC1 or CCA1, the details of their dynamics being irrelevant.
We systematically scanned the coupling window start and end times, adjusting model for each pair. This revealed that many coupling schemes are compatible with experimental data. For example, TOC1 degradation rate δ P T can be modified almost arbitrarily in a large temporal window between ZT22.5 and ZT6.5 without degrading adjustment. This is shown in Figs. 3(A)-(C), where δ P T = 3δ 0 P T inside this window (here and below, δ 0 X denotes the uncoupled degradation rate of variable X). Although the coupling is active for 8 hours, this coupling scheme generates mRNA and protein profiles which are indistinguishable from those of a free-running oscillator. Indeed, modifying TOC1 stability in a window where protein level is low, as is the case for any subinterval of the ZT22.5-ZT6.5 window, does not perturb the oscillator.
We also found a family of time windows of different lengths centered around ZT13.33, inside which the CCA1 degradation rate δ P C can be decreased without significantly mod-ifying goodness of fit. In Figs. 3(D)-(F), we show the effect of having δ P C = δ 0 P C /2 between ZT12.8 and ZT13.95. In this coupling scheme, mRNA profiles are not affected but coupling activation has a noticeable effect on CCA1 level, which rises faster than in the uncoupled case. After the window, however, CCA1 level relaxes in a few hours to the uncoupled profile, losing memory of the perturbation. Near this time of the day, the CCA1 protein level appears to be slaved by the other variables: the perturbation induced by modified degradation does not propagate to the other variables, and when coupling is switched off, the protein level relaxes to its value in the uncoupled solution. Thus, the effect of coupling is not only small but transient. An important consequence, which we will exploit later, is that the two coupling windows shown in Fig. 3 can be combined without modifying adjustment, provided the perturbation induced by one window has vanished when the other window begins.
In these examples, adjustment is sensitive to the timing of these coupling windows: when the start time is modified slightly, the end time must be changed simultaneously so as to recover good adjustment. On the other hand, we found that adjustment error depends little on the coupling strength (measured by the ratio between degradation rates outside and inside the window), especially for short coupling windows.  Fig. 3 as well as for two other windows inside which the CCA1 protein degradation is reduced, one shorter and the other longer than the window in Fig. 3(B). The window of accelated TOC1 degradation is totally insentitive to modifications of the TOC1 degradation rate, which is due to protein levels being very low in this window. Windows of CCA1 stabilization are all the more insensitive to variations in CCA1 degradation rate as they are shorter. To quantify the sensitivity of a given window we define r max as the largest value of the ratio r = δ 0 P C /δ P C such that adjustment RMS respectively.
To gain better insight into the effect of a coupling window, we must take into account the fact that the induced variation in the entrained oscillations can be decomposed as a displacement along the limit cycle (resulting in a phase shift) and a displacement transversely to the limit cycle (resulting in a deformation of the limit cycle). To this end, we apply a variable phase shift to the entrained time profile and optimize this phase shift so as to minimize the adjustment error. We define the waveform error as the minimal value of the latter, and the phase error the value of the phase shift for which it is obtained.
A small waveform error indicates that we are following the same limit cycle as in the free-running case, possibly with a different phase than is observed experimentally. Waveform and phase errors for the three windows of CCA1 protein stabilization considered in Thus it appears that for short enough windows, the effect of the light coupling mechanism can be entirely captured by studing the phase response induced by the mechanism and that a necessary property of a coupling window is that it induces a zero phase shift of the free-running limit cycle (or a phase shift corresponding to the mismatch between the natural and forcing periods in the general case that we will consider later).

Systematic characterization of gated coupling mechanisms
Besides the two specific examples shown in Fig. 3, other coupling schemes are compatible with experimental data. In this section, we undergo a systematic approach in order to determine those coupling schemes that do synchronize the free-running model to the day/night cycle, while leaving no signature on mRNA profiles when the phase-locking regime is achieved. To this aim, a preliminary step is to identify those coupling schemes that synchronize in the limit of weak forcing using the tools of infinitesimal phase response curve, which can be defined in the framework of perturbation theory in the vicinity of periodic orbits [38][39][40]. Computation of the parametric impulse phase response curve [41] (Z piP RC ) characterizing a light-coupling mechanism corresponding to parameter variation dp allows one to determine time intervals specified by duration τ and median position t m such that when the mechanism is applied in this time interval, it generates a zero phase shift and phase-locking is stable to small perturbations (see Supporting materials). Such intervals satisfy: As with the examples considered in the previous section, some coupling mechanisms have robust adjustment properties in that a good adjustment is obtained at the two different coupling strengths for the same timings, which coincide with the timings computed in the weak coupling limit. In these cases, adjustment is robust to variations in the coupling strength, which suggests that for these coupling mechanisms, the weak coupling approximation remains valid up to large coupling strengths. For instance, light coupling mechanisms that temporarily increase TOC protein degradation (δ P T ) or CCA1 activation threshold (P T 0 ) in windows located during the day the night appear to be robust couplings. Similarly, decreasing CCA1 protein degradation (δ P C ) or TOC repression threshold (P C0 ) in windows occuring during night are robust light-coupling mechanisms.
Some other mechanisms do not display the same robustness because either the window timings corresponding to good adjustment depend sensitively on coupling strength (e.g., for positive modulation of mTOC1 degradation rate) or because no good adjustment can be found except for very short windows (e.g., modulation of mCCA1 degradation rate).
Other robust coupling mechanisms can be identified in Fig. S4, in which the coupling mechanisms not considered in Fig. 5 are characterized. Our analysis shows that several coupling mechanisms are compatible with the experimental data and that discriminating them requires more experimental data. In particular, monitoring gene expression in transient conditions will probably be crucial since the coupling mechanism leaves apparently no signature in the experimental data in entrainement conditions. For simplicity, we restrict ourselves in the following to models in which half-lives of TOC1 or CCA1 proteins are modified during a specific time interval that is determined in Fig 5(D).

Resetting
One may wonder about the purpose of coupling schemes with almost no effect on the oscillator. The key point is that our data have been recorded when the clock was entrained by the diurnal cycle and phase-locked to it. A natural question then is: how do such couplings behave when clock is out of phase and resetting is needed? We found that while the two mechanisms shown in Fig. 3 have poor resetting properties when applied separately (Fig. S5), a combination of both can be very effective. In Fig. 7(A)-(B), we show how the two-gene oscillator recovers from a sudden phase-shift of 12 hours using a two-window coupling scheme. As described above, we assume for simplicity that the two coupling windows remain fixed with respect to the day/night cycle. To design this coupling, we utilized the fact that modifying coupling strengths inside windows hardly affects adjustment. We could therefore choose their values so as to minimize the maximal residual phase shift after three days for all possible initial lags ( Fig. 7(C)). Interestingly, we found that the best resetting behavior is obtained when the start time of the window of modified TOC degradation coincides with dawn. Phase locking in this example is globally stable. However, resetting becomes slow when the residual phase shift is under an hour and the residual phase shift is variable (RMS phase error after 5 days is 25 minutes and maximum phase error is 1 hour), and ( Fig. 7(C)).
This inefficiency results in fact from the limitations of a model where the two parameters are modulated by a rectangular profile with fixed timing. Indeed, we will see later that impressive adjustment and resetting behavior can be simultaneously obtained when parameters are modulated with smooth profiles. Our numerical results thus show that a coupling scheme can at the same time be almost invisible when the oscillator is in phase with its forcing cycle and effective enough to ensure resetting when the oscillator is out of phase. By invisible, we mean that the time profile remains in a close neighborhood of the uncoupled one, so that the only effect of coupling is to fix the phase of the oscillation with respect to the day/night cycle.

Robustness to daylight fluctuations
Why would it be beneficial for a circadian oscillator to be minimally affected by light/dark alternation in normal operation? A tempting hypothesis is that while daylight is essential for synchronizing the clock, its fluctuations can be detrimental to time keeping and that it is important to shield the oscillator from them. If the entrained temporal profile remains close to that of an uncoupled oscillator at different values of the coupling parameter, then it will be naturally insensitive to fluctuations in this parameter. To gain insight into this fundamental question, we subjected the fully coupled and occasionally coupled clock models to fluctuating daylight.
With the light input pathway unknown, we must allow for the fact that light fluctuations may be strongly attenuated upon reaching the Toc1 -Cca1 loop. For example, the light signal could be transmitted through an ultrasensitive signaling cascade with almost constant output above an input threshold close to daylight intensities at dawn. The core oscillator would then be subjected to a driving cycle much closer to a perfect square wave than the intensity profile. We thus considered varying modulation depths for the core oscillator parameters to reflect this possible attenuation.
Although the two types of model adjust experimental data equally well when subjected to a regular alternation, they have completely different responses to daylight fluctuations.
In Fig. 8, we assume that light intensity is constant throughout a given day but varies randomly from day to day. For almost zero modulation, the fully coupled model of Fig. 2(B) maintains relatively regular oscillations of varying amplitude (Fig. 8(B)). When parameter values are modulated by only a few percent, however, this model behaves erratically: oscillations stop for a few days, expression peaks occur a few hours in advance,... We also studied the effect of fluctuations at shorter time scales. When light intensity was varied randomly each hour, but with the same mean intensity each day, the permanently coupled model was still affected but much less than in Fig. 8 (Fig. S6).

Influence of free-running period
The results described above may seem to rely on the FRP being equal to 24 hours.
When the FRP is smaller or larger, coupling is required to achieve frequency locking and pull the oscillation period to 24 hours. To investigate this more general case, we scaled kinetic constants of the free-running model used in Fig. 2(B) to shift the FRP to 25 or 23.5 hours. In both cases (short FRP and long FRP), we could find models with gated coupling that adjust perfectly the experimental data with a period of 24 hours (Fig. 9).
These models are very similar to those shown in Fig. 3, the only notable difference being that coupling windows are shifted so that the induced resetting corrects for the period mismatch. Interestingly, the coupling windows for a FRP of 25 hours are located near the light/dark and dark/light transitions. We found that these coupling schemes were also very robust to daylight fluctuations (Fig. S7), indicating that the modulation ratio (equal to 3 for the two windows) is not critical. We also found that without taking adjustment into account, the free running oscillator is entrained by the coupling windows shown in Fig. 9) within a wide range of modulation ratios, from a lower threshold of 1.05 (resp.

Gating by smooth profiles
Gating of light input by rectangular profiles does not reflect the fact that the concentration of the mediators modulating the oscillator typically vary in a gradual way. The existence of nested coupling windows such that models with shorter windows can adjust data with larger parameter modulation (see Fig. 4) suggests investigating the action of smooth gating profiles, with maximal parameter modulation near the center of the window. To this end, we considered 24-hour periodic, Gaussian-shaped, modulation profiles defined , which are parameterized by the times of maximal modulation t C , t T , the coupling durations σ C , σ T and the modulation depths k C and k T . To assess whether good data adjustment and resetting behavior could be obtained simultaneously, these six parameters were chosen so as to minimize the RMS residual phase error 5 days after an initial random phase shift ranging from -12 to 12 hours (see Methods). Note that this naturally forces adjustment to experimental RNA profiles.
The behavior of the model using the optimized modulation profiles (Figs. 10(A)-(B)) confirms the findings obtained with rectangular profiles (Fig. 10). The entrained RNA and protein time profiles shadow that of the reference free-running oscillator, with little evidence of the coupling (Figs. 10(C)-(E)). Phase resetting in response to a phase shift is excellent (Fig. 10(F)): RMS (resp. maximum) residual phase shift after 5 days is 2.4 min (resp., 10 min). This is all the more remarkable as the Gaussian shape of the modulation profile is artificial, which shows that the dynamical mechanism exploited here is robust and relatively insensitive to the shape of the modulation profile. Moreover, the oscillator is extremely resistant to daylight fluctuations ( Fig. 10(F)). In spite of its simplicity, the two gene-oscillator studied here thus fulfills key requirements for a circadian oscillator when modulated with the right timing.

Discussion
Our findings illustrate how mathematical modeling can give insight into the architecture of a genetic module. Not only can expression profiles of two Ostreococcus clock genes be reproduced accurately by a simple two-gene transcriptional feedback loop model, but furthermore excellent adjustment of mRNA data is provided by a free-running model. This counterintuitive result can be explained if coupling to the diurnal cycle occurs during specific temporal windows, where unidentified mediators interact with the TOC1-CCA1 oscillator in such a way that it experiences negligible forcing when it is in phase with the day/night cycle, and strong resetting when it is out of phase. We could exhibit many coupling schemes compatible with experimental mRNA temporal profiles, differing by the coupling mechanism or by the window timing. This indicates that identification of the actual light input pathway will require additional experimental data. Our analysis strongly supports the conjecture that Ostreococcus genes Cca1 and Toc1 are the molecular components of an oscillator at the core of Ostreococcus clock but does not exclude that other coupled oscillators or feedback loops exist.
Why would a circadian oscillator decouple from the day/night cycle when in phase with it so as to generate quasi-autonomous oscillations? A natural hypothesis is that this protects the clock against daylight fluctuations, which can be important in natural conditions [42]. In a vast majority of numerical simulations and experiments on circadian clocks reported in the literature, the day/night cycle is taken into account through a perfect alternation of constant light intensity and darkness. However, this is somehow idealized, as the primary channel through which clocks get information about Earth rotation, namely daylight, is variable.
In nature, the daylight intensity sensed by an organism depends not only on time of day but also on various factors such as sky cover or, for marine organisms such as Ostreococcus, the distance to sea surface and water turbidity, which can affect perceived intensity much more than atmosphere. Therefore, the light intensity reaching a circadian clock can vary several-fold not only from one day to the next but also between different times of the day. A clock permanently coupled to light is also permanently subjected to These results lead to enquire whether similar designs exist in other circadian clocks.
Although the importance of this problem was noted some time ago [42], the robustness of circadian clocks to daylight fluctuations and how this constraint shapes their molecular architecture have been little studied until very recently [43,44]. The discussion on how genetic oscillators can keep daytime has essentially focused on the most important sources of noise under constant conditions : temperature variations [40,45,46] or fluctuations in concentration due to small numbers of molecules [47,48]. However, an operating clock is naturally subjected to an external forcing cycle, which is yet another source of fluctuations.
We thus conjecture that a circadian clock must be built so as to be insensitive to daylight intensity fluctuations when entrained by the day/night cycle, just as it is insentitive to molecular or temperature fluctuations, and that this can be achieved by keeping the oscillator as close to the free-running limit cycle as possible, scheduling coupling at a time when the oscillator is not responsive. An important consequence of this principle is that it allows us to discriminate between different possible coupling mechanisms for a given model, as our analysis revealed dramatic differences in the ability of different parametric modulations to buffer fluctuations. It also allows us to determine the preferred timing for a given coupling mechanism, which may prove very helpful when trying to identify the molecular actors which mediate the light information to the clock.
When the FRP is close to 24 hours, as in much of our analysis it is easy to understand why robustness to daylight fluctuations requires that the forced oscillation shadows the free-running solution. Robustness manifests itself in the time profile remaining constant when subjected to random sequences of daylight intensity. This includes strongly fluctuating sequences as well as sequences of constant daylight intensity at different levels.
Thus, the oscillator response should be the same at high and low daylight intensities, which implies that the solution must remain close to the free-running one as forcing is increased from zero. Note that this only holds in entrainment conditions, where coupling is not needed. When the clock is out of phase, strong responses to forcing are expected, with resetting being faster as forcing is stronger.
When the natural and external periods are significantly different, the problem may seem more complex as coupling is required to correct the period mismatch. There is a minimal coupling strength under which the oscillator is not frequency-locked and entrainment cannot occur. Nevertheless, we showed that properly timing the coupling windows is as effective for oscillators with FRP of 23.5 and 25 hours as for the 24-hour example we had considered. Again, the forced solution remains close to the free-running limit cycle even if proceeding at a different speed to correct the period mismatch. This also shows that FRP is not a critical parameter for adjustment of the experimental data used here.
A consequence of the small deviation of the limit cycle from the free-running one when coupling strength is varied is that oscillations should vary little upon a transition from LD to LL or DD conditions (see, e.g., Figs. 9(G)-(H)). We searched the litterature for examples of such behavior. Ref. [13] provides a interesting comparison of models for the Drosophila and Neurospora circadian clocks which is illustrative for our discussion. In this study, the variation in amplitude is much less pronounced for the Drosophila model than for the Neurospora one (see Fig. 2 of [13]). Concurrently, the sensitivity of the phase of the entrained oscillations to variations in the light-controlled parameter is much smaller for the Drosophila model (see Fig. 3 of [13]), which is a necessary condition for robustness to daylight fluctuations. Another interesting comparison involves the one-loop and two-loop models of Arabidopsis clock [24,25]. The one-loop model clearly modifies its behavior upon entering DD conditions from LD (see Fig. 5 of [24]) while the twoloop model preserves its average waveform when transiting from LD to LL, except for the disappearance of the acute response to light at dawn (see Fig. 6 of [25]). Thus, the two-loop model not only reproduces experimental data better but also seems more robust.
The Drosophila and Neurospora clock models analyzed in [13] also differ in their response to forcing when their FRP is close to 24 hours [49]. A number of circadian models cannot be entrained when their FRP is too close to 24 hours because complex oscillations, period-doubled or chaotic ones, are observed easily for moderate to strong forcing. Indeed, it is expected that near resonance between the forcing and natural periods, the strong response exalts nonlinearities and favors complex behavior. Again, the Drosophila clock model appears to be more robust in this respect [49]. We stress that making the coupling invisible in entrainment conditions naturally addresses this issue. Dynamically uncoupling the oscillator from the diurnal cycle in entrainment conditions makes it immune both to fluctuations in daylight intensity and to destabilization in the face of strong forcing.
An important problem is how a clock with occasional coupling can adjust to different photoperiods so as to anticipate daily events all along the year. We can only touch briefly this question here as it requires understanding how the temporal profile of the coupling windows changes with photoperiod and thus a detailed description of the unknown light input pathways and additional feedback loops that control the timing of these windows.
The key point is that the phase of the entrained oscillations is controlled by the position of the coupling windows. Thus the role of light input pathways and additional feedback loops, whose internal dynamics will typically be affected by input from photoreceptors and feedback from the TOC1-CCA1 oscillator, is to time the coupling windows as needed for each photoperiod so that the correct oscillation timing is generated [35][36][37]. This question will be addressed in a future work, together with the analysis of the luminescence time series recorded for differents photoperiods.
Our results also bring some insight into the recent observation that a circadian clock may require multiple feedback loops to maintain proper timing of expression peaks in response to noisy light input across the year [43]. We have shown here that a single twogene loop can display impressive robustness to daylight fluctuations when its parameters are modulated with the right timing. As noted when discussing the response to different photoperiods, this requires the presence of additional feedback loops to generate the biochemical signal needed to drive the core oscillator appropriately, and which we have not yet identified and modeled in Ostreococcus. Robustness to fluctuations thus implies a minimal level of complexity.
Finally, robustness to intensity fluctuations may explain why it is important to have a self-sustained oscillator at the core of the clock, as a forced damped oscillator permanently needs forcing to maintain its amplitude, and is thereby vulnerable to amplitude fluctuations. Confining the dynamics near the free-running limit cycle allows to have a pure phase dynamics for the core oscillator, uncoupled from intensity fluctuations. Understanding how to construct it will require taking into account the sensitivity of the free-running oscillator to perturbations across its cycle [50].
A simple organism as Ostreococcus can apparently combine mathematical simplicity with the complexity of any cell. The low genomic redundancy of Ostreococcus is certainly crucial for allowing accurate mathematical modeling, leading to better insight into the clock workings. Ostreococcus therefore stands as a very promising model for circadian biology, but also more generally for systems biology.

Materials and Methods
A minimal mathematical model of the transcriptional loop where Toc1 activates Cca1 which represses Toc1, consists of the following four differential equations: Eqs (2)  and cooperativity n C . Similarly, Cca1 transcription rate is µ C (resp., µ C + λ C ) at zero (resp., infinite) TOC1 concentration, with threshold P T 0 and cooperativity n T . Translation of TOC1 and CCA1 occurs at rates β T and β C , respectively. For each species Y , the Michaelis-Menten degradation term is written so that δ Y is the low-concentration degradation rate and K Y is the saturation threshold.
Model (2)  performed with the SEULEX algorithm [52]. Adjustment was carried out with 14 (resp. 2) Quad-Core Intel Xeon processors at 2.83 GHz during 72 hours for the 28-dimensional (resp. 12-dimensional) parameter space. Convergence was checked by verifying that the vicinity of the optimum was well sampled. In the uncoupled case, the ODE system is invariant under time translation so that its solutions are defined up to an arbitrary phase.
An additional routine was then used to select the best-fitting phase.
To study the effect of daylight fluctuations, parameters were modulated as follows.
L(t) ∈ [0, 1] is the randomly varying light intensity, with L ref = 0.5 the reference level.
We define the reference modulation depth of the Y parameter taking value Y L at standard light level and Y D in dark as k ref   , where the shaded area illustrates the sawtooth shape of the Cca1 mRNA profile, which will be used later as evidence of a strongly saturated enzymatic degradation. This area has been obtained by fitting a straight line through Cca1 data points at ZT12, ZT15 and ZT18 on one hand and at ZT21, ZT0 and ZT3 on the other hand.   Fig 2(B)) and with coupling (solid lines). Gray areas indicate coupling activation. In the left (resp. right) column, TOC1 (resp. CCA1) degradation rate is multiplied by 3 (resp. divided by 2) from ZT22.5 to ZT6.5 (resp. from ZT12. 8        . Adjustment by models with gated coupling when FRP is different from 24 hours. Gated coupling can also synchronize free-running clock models with a FRP of 23.5h or 25h without leaving any signature in mRNA profiles. Top left, (A)-(C): numerical solutions of model (2) for a FRP of 23.5h, subjected to coupling windows shown as shaded areas. TOC1 (resp. CCA1) protein degradation rate is multiplied (resp. divided) by three from ZT3 to ZT5.5 (resp. ZT15 to ZT17.5). Top right, (D)-(F): numerical solutions of model (2) for a FRP of 25h, subjected to coupling windows shown as shaded areas. TOC1 (resp. CCA1) protein degradation rate is multiplied (resp. divided) by three from ZT22.75 to ZT24 (resp. ZT11.75 to ZT12).    Fig. 2(B)).

Supporting Information
Gated-coupling design in the weak modulation limit The free oscillator model was shown to adjust remarkably well the RNAmicroarray data from LD12:12 experiments. A tempting hypothesis is that the synchronization of the free running oscillator to the day-night cycle involves a light-dependent gated-coupling mechanism that has restricted effect on the RNA traces when phase locked. We develop here a systematic method to repertoriate the coupling schemes that synchronize the free oscillator to the diurnal cycle while preserving the adjustment score obtained in the absence of coupling. For enough weak coupling strength, any coupling schemes that achieve the correct locking phase preserve the adjustement score. Those coupling schemes can be found in the framework of perturbation theory in the vicinity of a periodic orbit [1,2,3], assuming that the driving force period is enough close to the internal clock period. We consider the state vector of a nonlinear oscillator, which represents the concentration of the molecular clock components. In constant dark conditions, the concentration vector X evolves according to: Eq. 3 has a periodic solution X γ (t) corresponding to a stable limit cycle of period T close to 24 hours. We assume that the coupling between the light and the circadian oscillator is mediated by a set of N components (k is the index), which modulate the parameter vector in the direction of dp k : where the 24h-periodic scalar function L k (t, τ k , (t m ) k ) represents the temporal profile of activation (rectangular-or gaussian-shaped profiles in the present paper) of the lightdependent component k with τ k and (t m ) k characterizing the effective coupling window duration and center (t = 0 correspond to the night-day transition or CT0).
A small enough parametric impulse perturbation applied at phase u induces an infinitesimal change of the circadian oscillator phase defining a T -periodic scalar function Z piP RC (u, dp) called infinitesimal phase response curve [2] or, to be more precise, parametric impulse phase response curve [3]: Z piP RC (u, dp k ) = (dp k ) T .Z p (u) where Z p (u) = [ ∂F(X γ (u)) ∂p ] T ∂φ(X γ (u)) ∂X Then, the phase change induced by the light sensed during the daytime can be derived from the convolution of the temporal profile of the light-sensing components with the piPRC: ∆φ = k=1,N T 0 L k (u, τ k , (t m ) k )Z piP RC (u + φ, dp k )du (7) where φ is the phase of the oscillator at CT0. A stable entrainment state requires that the scalar functions L and Z ipP RC satisfies: T 0 L k (u, τ k , (t m ) k ) Z piP RC (u + φ * , dp k )du = δφ * k=1,N T 0 L k (u, τ k , (t m ) k ) Z ′ piP RC (u + φ * , dp k )du < 0 where φ * is the locked phase (relative to CT0) and δφ * is the phase change induced by the period mismatch between the free oscillator and the day-night period, which is assumed to be small with respect to T .
Phase RMS error after 5 day/night cycles is 25 min while the maximum error is 1 hour.