Aurore Woller

Hélène Duez

Bart Staels

Marc Lefranc
email: marc.lefranc@univ-lille1.fr

A mathematical model of the liver circadian clock linking feeding/fasting cycles to clock function

Keywords: mammalian circadian clock, metabolism, mathematical model, SIRT1, AMPK, NAMPT, NAD+, REV-ERB, metabolic disorders

To maintain energy homeostasis despite variable energy supply and consumption along the diurnal cycle, the liver relies on a circadian clock synchronized to food timing. Perturbed feeding/fasting cycles have been associated with clock disruption and metabolic diseases, however the mechanisms are unclear. To address this question, we have constructed a mathematical model of the mammalian circadian clock incorporating the metabolic sensors SIRT1 and AMPK. The clock response to various temporal patterns of AMPK activation was simulated numerically, mimicking the effects of a normal diet, fasting and a high-fat diet. The model reproduces the dampened clock gene expression and NAD+ rhythms reported for mice on a high-fat diet, and predicts that this effect may be pharmacologically rescued by timed REV-ERB agonist administration. Our model thus identifies altered AMPK signaling as a mechanism leading to clock disruption and its associated metabolic effects, and a pharmacological approach to reset the clock in obesity.

Introduction

Life on Earth is subjected to periodic changes in its environment induced by Earth rotation. Most organisms anticipate these daily variations and orchestrate biological processes accordingly by relying on a circadian clock, a network of molecular interactions generating biochemical oscillations with a period close to 24 hours [START_REF] Dibner | The mammalian circadian timing system: organization and coordination of central and peripheral clocks[END_REF]. To synchronize with the diurnal cycle, circadian clocks incorporate sensors informing them of its progression. As the alternation of day and night is the primary environmental signal, light is generally the dominant circadian cue at the organismal level. In multicellular organisms, most cells contain a self-sustained circadian oscillator [START_REF] Dibner | The mammalian circadian timing system: organization and coordination of central and peripheral clocks[END_REF], which however does not receive the light signal directly. The mammalian circadian system thus relies on a central synchronizer, the suprachiasmatic nucleus (SCN), a group of neurons which receives photic inputs and drives other circadian oscillators in the brain and in other organs through various channels [START_REF] Dibner | The mammalian circadian timing system: organization and coordination of central and peripheral clocks[END_REF]].

Yet, the rhythms of circadian clocks throughout our bodies are not governed only by the light/dark cycle. The function of some organs exposes them to oscillating inputs which superimpose on the time of the day. In particular, the liver plays a central role in maintaining energy homeostasis while coping with large temporal variations of energy income, storage and utilization over the diurnal cycle. Accordingly, the feeding/fasting cycle is a more potent zeitgeber for the liver clock than systemic cues controlled by the SCN [START_REF] Damiola | Restricted feeding uncouples circadian oscillators in peripheral tissues from the central pacemaker in the suprachiasmatic nucleus[END_REF], Saini et al., 2013]. Understanding how the liver clock is entrained by metabolic cycles is thus crucial, especially as metabolism and the circadian clock are tightly integrated [Bass, 2012, Asher andSchibler, 2011].

On the one hand, clock disruption [Turek et al., 2005, Marcheva et al., 2010] or perturbations of the daily pattern of food intake [START_REF] Arble | Circadian timing of food intake contributes to weight gain[END_REF] have been associated with metabolic diseases. On the other hand, nutritional challenges or metabolic stress affect the clock period length and the expression profiles of core clock genes [Bass, 2012, Kohsaka et al., 2007, Hatori et al., 2012, Eckel-Mahan et al., 2013].

The implication of most liver-specific cycling transcripts in diverse aspects of metabolism further reinforces the link between liver circadian rhythms and metabolism [Panda et al., 2002, Zhang et al., 2014].

Here, we use mathematical modeling to address several important questions: What are the key factors that inform the liver clock about the cellular energetic state, so that it can be entrained to daily variations of food intake? Can we explain how perturbations of these cycles or nutritional stress disrupt the clock? Can we design t he r ap e ut i c strategies to counter such deleterious effects? Such a quantitative approach is needed to integrate the large number of experimental observations from the literature. Existing mathematical models of the mammalian circadian clock [e.g., [START_REF] Leloup | Toward a detailed computational model for the mammalian circadian clock[END_REF], Becker-Weimann et al., 2004, Forger and Peskin, 2003, Mirsky et al., 2009, Relógio et al., 2011, Korencic et al., 2012, St. John et al., 2014, Jolley et al., 2014] have focused on describing clocks that are freerunning or driven by light/dark cycles, and thus are not adequate to answer these questions. To understand how the clock is entrained by metabolism , we have t h u s constructed a model of the mammalian clockwork that tracks v a r i a t i o n s i n NAD+ and AMP, two metabolites which are central to biochemical reactions involved in energy production, storage and utilization [START_REF] Berg | Biochemistry[END_REF].

Indeed, it has been recently shown that molecular sensors of NAD+ and AMP regulate core clock genes, thus serving as metabolic inputs to the clock [START_REF] Huang | Circadian rhythms, sleep, and metabolism[END_REF].

The first sensor is the NAD+-dependent histone and protein deacetylase SIRT1, which regulates many metabolic pathways [START_REF] Chang | SIRT1 and other sirtuins in metabolism[END_REF]. The second sensor is the AMPK kinase, which is activated in situations of energy shortage [Hardie et al., 2012]. Systemic hormones such as insulin and glucagon have been shown to contribute also to clock resetting [START_REF] Mukherji | Shifting the feeding of mice to the rest phase creates metabolic alterations, which, on their own, shift the peripheral circadian clocks by 12 hours[END_REF], however we only consider intracellular factors in this first approach.

Our model reproduces temporal profiles of clock gene expression and NAD+ levels accurately, assuming AMPK activation patterns which are consistent with experimental profiles of AMP level [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF] and AMPK expression [Barnea et al., 2012], and predicts correctly the effect of SIRT1 loss-of-function and impaired AMPK activity through LKB1 deletion. This allows us to study how the clock is affected by various temporal patterns of AMPK activity, used as surrogates of cell metabolic state variations upon normal diet feeding, fasting and ad libitum high-fat diet feeding. Finally, we illustrate how the disruptive effect of an obesogene diet may be corrected by a pharmacological intervention.

The architecture of the mammalian clock network is well known [START_REF] Dibner | The mammalian circadian timing system: organization and coordination of central and peripheral clocks[END_REF]. The activators CLOCK and BMAL1 dimerize to induce transcription of several target genes, including the Period (Per 1, 2, 3) and Cryptochrome (Cry 1, 2) genes. Accumulating PERs and CRYs associate in a complex which enters the nucleus to inhibit CLOCK:BMAL1 transcriptional activity, thus repressing their own expression and forming the main feedback loop [START_REF] Dibner | The mammalian circadian timing system: organization and coordination of central and peripheral clocks[END_REF]. Additional loops are associated with other targets of CLOCK:BMAL1, the nuclear receptors REV-ERBα,β [START_REF] Preitner | The orphan nuclear receptor REV-ERBalpha controls circadian transcription within the positive limb of the mammalian circadian oscillator[END_REF], Guillaumond et al., 2005, Ueda et al., 2005, Liu et al., 2008, Cho et al., 2012, Bugge et al., 2012], and RORα,β, γ [START_REF] Akashi | The orphan nuclear receptor RORalpha regulates circadian transcription of the mammalian coreclock Bmal1[END_REF], Guillaumond et al., 2005, Ueda et al., 2005, Liu et al., 2008], which respectively repress and activate Bmal1 transcription. The REV-ERBs, which also repress Cry1 transcription [START_REF] Liu | Redundant function of REV-ERBalpha and beta and non-essential role for Bmal1 cycling in transcriptional regulation of intracellular circadian rhythms[END_REF], are essential for robust oscillations [START_REF] Relógio | Tuning the mammalian circadian clock: robust synergy of two loops[END_REF], Cho et al., 2012, Bugge et al., 2012]. Together with the RORs, they control important metabolic and physiological functions [START_REF] Duez | The nuclear receptors Reverbs and RORs integrate circadian rhythms and metabolism[END_REF], Duez and Staels, 2010, Cho et al., 2012, Bugge et al., 2012, Solt et al., 2012, Woldt et al., 2013, Everett and Lazar, 2014]. Conversely, several mechanisms by which metabolism drives the clock were recently described. Indeed, the NAD+-dependent deacetylase SIRT1 and the AMP sensor AMPK directly target core clock genes [START_REF] Huang | Circadian rhythms, sleep, and metabolism[END_REF]. SIRT1 complexes with CLOCK:BMAL1, deacetylates BMAL1 and PER2 as well as histones, modulating the transcriptional activity of CLOCK:BMAL1 and destabilizing PER2 [START_REF] Nakahata | The NAD+dependent deacetylase SIRT1 modulates CLOCK-mediated chromatin remodeling and circadian control[END_REF], Asher et al., 2008]. Additionally, SIRT1 activates PGC1α, which enhances Bmal1 expression by co-activating ROR [START_REF] Liu | Transcriptional coactivator PGC-1alpha integrates the mammalian clock and energy metabolism[END_REF]. In fact, NAD+ and SIRT1 are not only inputs but also core clock players, as the NAD+ salvage pathway, the main mammalian NAD+ biosynthesis pathway, is under circadian control [START_REF] Ramsey | Circadian clock feedback cycle through NAMPT-mediated NAD+ biosynthesis[END_REF], Nakahata et al., 2009]. Many cellular reactions, such as deacetylation by sirtuins, simultaneously hydrolyze NAD+ into NAM [START_REF] Cantó | NAD+ as a signaling molecule modulating metabolism[END_REF]. The salvage pathway recycles NAM into NMN and NMN into NAD+. NAMPT, the ratelimiting enzyme in the first step of the NAD+ salvage pathway, is directly regulated by CLOCK:BMAL1, resulting in circadian oscillations of the cellular NAD+ level [START_REF] Ramsey | Circadian clock feedback cycle through NAMPT-mediated NAD+ biosynthesis[END_REF], Nakahata et al., 2009]. Thus, CLOCK:BMAL1, NAMPT, NAD+ and SIRT1 form a negative feedback loop which contributes to the clock dynamics, as demonstrated by inhibiting NAMPT or SIRT1 [START_REF] Ramsey | Circadian clock feedback cycle through NAMPT-mediated NAD+ biosynthesis[END_REF], Nakahata et al., 2009, Bellet et al., 2013]. On the other hand, activated AMPK destabilizes CRY [Lamia et al., 2009], as well as indirectly PER by activating CKIε/δ [START_REF] Jordan | AMPK at the crossroads of circadian clocks and metabolism[END_REF]. Moreover, phosphorylation of PGC1α by AMPK is required before it can be deacetylated by SIRT1 [START_REF] Cantó | Interdependence of AMPK and SIRT1 for metabolic adaptation to fasting and exercise in skeletal muscle[END_REF].

The SIRT1 and AMPK signaling pathways are tightly coordinated [START_REF] Ruderman | AMPK and SIRT1: a long-standing partnership?[END_REF]. SIRT1 controls AMPK activation via deacetylation of liver kinase B1 (LKB1) [Lan et al., 2008, Hou et al., 2008], which phosphorylates and activates AMPK.

Conversely, AMPK activation leads to increased NAD+ levels and subsequent SIRT1 activation [Ruderman et al., 2010, Froy andMiskin, 2010]. It is unclear whether this effect occurs through up-regulation of NAMPT activity [Fulco et al., 2008, Brandauer et al., 2013] or via a modification of the NAD+/NADH ratio by metabolic processes induced by AMPK [START_REF] Cantó | AMPK regulates energy expenditure by modulating NAD+ metabolism and SIRT1 activity[END_REF]. However, the induction of NAMPT protein expression by nutrient restriction [START_REF] Yang | Nutrient-sensitive mitochondrial NAD+ levels dictate cell survival[END_REF] is consistent with the former hypothesis. Moreover, it has been reported recently that AMPK may regulate NAMPT activity post-translationally in skeletal muscle [START_REF] Brandauer | AMP-activated protein kinase regulates nicotinamide phosphoribosyl transferase expression in skeletal muscle[END_REF]. Importantly, such a regulation could account for the two peaks in NAMPT protein and NAD+ profiles observed in liver, one during the day and one during the night, in contrast with the single nighttime peak of Nampt mRNA (see Fig. 1 in [START_REF] Ramsey | Circadian clock feedback cycle through NAMPT-mediated NAD+ biosynthesis[END_REF]).

In conclusion, metabolism drives the circadian clock through many intertwined molecular interactions and feedback loops, whose coordination and relative importances can only be assessed mathematically

Results

Construction of the mathematical model

Our mathematical model describes the dynamics of the regulatory network shown in Fig. 1. To identify the core actors and interactions and to avoid overfitting, we kept the mathematical model as simple as possible. We grouped the three Period homologs as a single Per gene and the two cryptochromes as a single Cry gene.

Similarly, the Rev-Erb and Ror genes represent the two isoforms Rev-Erbα, β and the three isoforms Rorα, β, γ, respectively. The CLOCK protein was assumed to be constitutively expressed. Dbp served as an example of a clock-controlled gene, although other theoretical studies assumed that it belongs to the clock network [Korencic et al., 2012]. We did not take into account post-translational protein modifications except those induced by SIRT1 and AMPK, nor compartmentalization, considering that transport between cytoplasm and nucleus is fast on a circadian time scale.

Both the feeding/fasting and light/dark cycles influence the liver clock [Saini et al., 2013]. However, when these two cues are discordant, the entrainment phase is determined by food timing [START_REF] Damiola | Restricted feeding uncouples circadian oscillators in peripheral tissues from the central pacemaker in the suprachiasmatic nucleus[END_REF], Saini et al., 2013]. Thus, we neglected the influence of systemic signals from the SCN, and considered only how variations in AMPK and SIRT1 activities, induced by fluctuating AMP and NAD+ levels, drive the clock. This is legitimate as long as only the entrained regime is studied, not the resetting dynamics. Experimental data (Table S1 of [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF]) indicate that AMP levels cycle throughout the 24 hours, with one peak during the rest phase (around Zeitgeber time (ZT) 5) and a smaller one during the active period (around ZT17). To simplify the analysis, AMPK activity time course was described by a parametric profile with two pulses whose widths, amplitudes and timings are adjustable parameters, constrained to be consistent with experimental AMP level and AMPK abundance profiles (Figs. S1-S2, Supplemental Experimental Procedures and Supplemental Results). In this scheme, the action of SIRT1 on AMPK via LKB1 does not appear explicitly, but can be taken into account through the shape of the AMPK profile. SIRT1 activity was given by a function of NAD+ level, and PGC1α activity by a function of SIRT1 and AMPK activities (Supplemental Table S2). To take into account variations of nuclear PGC1α abundance without specifying the mechanisms regulating it, we assumed that it alternated between a low level from ZT3 to ZT11 and a high level for the remaining time, consistent with experimental observations [START_REF] Liu | Transcriptional coactivator PGC-1alpha integrates the mammalian clock and energy metabolism[END_REF].

NAD+ levels are influenced by (1) the variations of NAMPT activity due to the circadian regulation of Nampt [START_REF] Ramsey | Circadian clock feedback cycle through NAMPT-mediated NAD+ biosynthesis[END_REF], Nakahata et al., 2009] or to a hypothetical regulation of NAMPT by AMPK; and (2) metabolic processes such as glycolysis or beta-oxidation which modify the NAD+/NADH ratio [START_REF] Cantó | AMPK regulates energy expenditure by modulating NAD+ metabolism and SIRT1 activity[END_REF].

In a first approximation, the latter mechanism was neglected. This non-intuitive assumption still allowed us to reproduce experimental NAD+ profiles, and is consistent with observations that changes in NAD+ biosynthesis suffice to hinder SIRT3 deacetylase function [Peek et al., 2013], and that Bmal1 deficiency leads to drastically reduced Nampt expression and NAD+ levels [START_REF] Ramsey | Circadian clock feedback cycle through NAMPT-mediated NAD+ biosynthesis[END_REF]. A feature of our model is the up-regulation of NAMPT protein activity by AMPK through inhibition of NAMPT degradation. This is indeed a plausible way to account for the two-peak structure shown by the NAMPT protein profile together with the NAD+ profile (Fig. 1 in [START_REF] Ramsey | Circadian clock feedback cycle through NAMPT-mediated NAD+ biosynthesis[END_REF]), especially as the NAMPT protein and NAD+ peaks at ZT5 [START_REF] Ramsey | Circadian clock feedback cycle through NAMPT-mediated NAD+ biosynthesis[END_REF], Hatori et al., 2012] coincide with an AMP peak [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF] (Supplemental Information). The hypothesis that AMPK regulates NAMPT translation rate led to very similar results.

The resulting mathematical model consists of 16 ordinary differential equations describing the time evolutions of the mRNA and protein concentrations for the clock genes Bmal1, Per, Cry, Rev-Erb, Ror, the metabolic gene Nampt, of the mRNA concentration for the clock output gene Dbp and of the NAD+ level (Tables S1 andS2, Supplemental Experimental Procedures). It depends on 96 kinetic constants, most of which are unknown and must be estimated from experimental data (Table S3, Supplemental Experimental Procedures).

Expression time profiles of the main clock actors in mouse livers are accurately reproduced by the model

To validate the assumptions and ensure that the model adequately addresses the biological questions raised, we adjusted it to published gene expression and NAD+ temporal profiles (Table S3). We reused mRNA level profiles obtained by Hughes et al. [Hughes et al., 2009] from livers of mice entrained to a 12:12 light/dark (LD) cycle and then put in constant darkness and fed ad libitum. Taking advantage of the high temporal resolution (1 point/hour), expression profiles were adjusted by Fourier series with four harmonics (Fig. S3), and the latter were matched to in silico profiles. The numerical Per, Cry, Ror and Rev-Erb profiles were compared to the experimental Per2, Cry1, Rorγ and Rev-Erbα profiles, respectively. The NAD+ and AMP data were those obtained by Hatori et al. (Table S1 of [Hatori et al., 2012]) from livers of mice entrained to a 12:12 LD cycle and fed a normal chow diet ad libitum.

An excellent agreement between the predicted and experimental time profiles was obtained, which validates the simplifying assumptions used in the model (Fig. 2). In particular, the NAD+ profile, which is here slaved to NAMPT activity, displays the twopeak structure observed experimentally [START_REF] Ramsey | Circadian clock feedback cycle through NAMPT-mediated NAD+ biosynthesis[END_REF], Hatori et al., 2012].

The peak during the activity period (night) is linked to the circadian peak of NAMPT, while the one during the rest period (day) results from the up-regulation of NAMPT protein by AMPK (Fig. 2C). Thus, our mathematical model describes simultaneously the daily oscillations in the expression of clock and Nampt genes as well as in NAD+ levels, entrained by daily variations in AMP levels induced by the feeding/fasting cycle. Although NAD+ and AMP may seem to play different roles in this simplifying approach, both metabolites are essential in our model, as is clearly demonstrated by the disruptions observed when the metabolic sensors SIRT1 and AMPK are inactivated, as described in the next section.

SIRT1 and AMPK loss-of-function phenotypes are qualitatively reproduced

To validate our model further, we aimed to reproduce SIRT1 and LKB1 KO phenotypes in silico. Ablation of SIRT1 results in a sharp amplitude increase of Nampt, Dbp and Per mRNA [Bellet et al., 2013], consistent with the down-regulation of CLOCK:BMAL1 activity by SIRT1. Conversely, LKB1 knock-out, which impairs AMPK activation, generally decreases the amplitude of clock gene expression [Lamia et al., 2009]. To simulate the effect of LKB1 KO, AMPK activity was reduced to 3.75 % of its baseline value. With this approach, we were able to further constrain parameter values to reproduce the marked increase in the amplitude of gene expression observed in SIRT1 deficient mice [Bellet et al., 2013] as well as the qualitative changes in expression reported for LKB1 deficiency [Lamia et al., 2009] (Fig. 3). Interestingly, comparison of Nampt gene expression profiles between WT and SIRT1 KO mice allowed to estimate the strength of the modulation of CLOCK:BMAL1 activity by SIRT1, which was poorly constrained by the data from WT mice taken alone, and thus to assess the relative strength of the NAMPT-NAD+-SIRT1-CLOCK:BMAL1 feedback loop within the clock. Together, the adjustments in Figs. 2 and3 support the idea that the molecular network of Fig. 1 captures the main interactions needed to understand the entrainment of the liver clock by metabolism.

Effect of perturbations of AMP rhythms on the clock

Modern life style is characterized by a diversity of feeding patterns and correlates with an increased prevalence of metabolic disorders. Perturbations of the feedingfasting cycle have been associated with modifications in physiology and in clock gene expression [START_REF] Kohsaka | High-fat diet disrupts behavioral and molecular circadian rhythms in mice[END_REF], Hatori et al., 2012, Eckel-Mahan et al., 2013].

Interestingly, [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF] reported that mice under ad libitum high-fat diet (AL HFD) developed obesity, while those fed under time-restricted HFD (TR HFD) showed increased energy expenditure and nutrient utilization. AL HFD feeding led to dampening in clock gene expression, but clock gene expression in TR HFD fed mouse livers was similar in amplitude to those on a normal diet.

Additionally, Eckel-Mahan et al. [Eckel-Mahan et al., 2013] observed a strong reduction in the amplitude of Nampt expression and NAD+ level oscillations in AL HFD [Eckel-Mahan et al., 2013], potentially affecting the many metabolic roles of sirtuins [START_REF] Chang | SIRT1 and other sirtuins in metabolism[END_REF], including in oxidative metabolism [Peek et al., 2013]. However, the exact molecular mechanisms leading to these dysregulations remain poorly understood. Since it was noted that AMP levels are depressed in AL HFD fed mouse livers [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF], thus reducing AMPK activity, we sought to explore the putative effect of AMP rhythm modifications on this phenotype using our mathematical model. Specifically, two conditions were simulated in addition to the normal feeding/fasting rhythm corresponding to the AMPK activity pattern used in Fig. 2: a fasted-like state where AMPK activity remains high at all times, and a fed-like state with low AMPK activity, mimicking high-fat feeding, the amplitude of oscillations being significantly reduced in both cases. Remarkably, a decreased AMPK activity leads to a significant reduction in Nampt mRNA and protein expression, as well as in NAD+ levels (Figs. 4 and S4), as shown upon AL HFD feeding (Figs 3B and 4D-E in [Eckel-Mahan et al., 2013]). Conversely, constitutively active AMPK enhances Nampt mRNA, NAMPT protein and NAD+ levels (Figs. 4 and S4), and is also associated with a significant phase shift in NAD+ level and Nampt expression profiles, which is consistent with experimental observations in fasted mice [Peek et al., 2013]. The good agreement between numerical predictions and experimental observations strongly suggests that dampened AMPK activity rhythms contribute significantly to the modifications of clock gene and NAD+ rhythms observed under HFD conditions [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF], Eckel-Mahan et al., 2013], even though other additional nutrient sensors, such as mTOR and CREB, or other pathways are also likely involved. This is consistent with the fact that reduced amplitudes in Nampt and NAD+ profiles have been observed within 3 days of HFD [Eckel-Mahan et al., 2013] and thus cannot be caused by the concurrent development of obesity. While the loss of the circadian NAD+ peak has been shown to impair mitochondrial oxidative metabolism [Peek et al., 2013], it remains unclear whether obesity results from the loss of amplitude of gene expression, as hypothesized in [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF].

The model predicts the rescue of amplitude of gene expression by Cry1 deficiency in HFD-induced obese mice

To gain insight into this question, we searched the literature for other mechanisms besides time-restricted feeding which could protect animals against HFD-induced weight gain. A recent study reported that Cry1 -/-mice are protected against HFDinduced obesity [START_REF] Griebel | Mice deficient in Cryptochrome 1 (Cry1 (-/-)) exhibit resistance to obesity induced by a high-fat diet[END_REF]. We thus asked whether our model would predict a rescue of clock gene expression amplitude upon Cry1 deficiency, similar to what is observed during TR HFD [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF].

To mimic a selective knockout of Cry1 in our model, the transcription rate of Cry, which accounts for both Cry1 and Cry2, was halved. This mutant phenotype was then subjected to reduced AMPK rhythms. Remarkably, our mathematical model predicts that the amplitude of oscillations in NAD+ levels and expression of most clock genes is restored to physiological levels (Fig. 5). This result, together with the findings of [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF], provides additional support to the hypothesis that maintaining the amplitude of clock gene expression protects against metabolic alterations induced by HFD and prompts for further investigations on the pathophysiological effects of dampened clock gene rhythms.

Designing a pharmacological rescue of physiological expression profiles upon disruption of AMPK rhythms

If deleterious effects result from dampened clock gene expression, an important question is whether they can be reversed pharmacologically, and whether our mathematical model could help us to design such a strategy. One possible approach is to modulate transiently the activity of a clock protein through drug administration. A good target is REV-ERB because, on the one hand, REV-ERB reduces Cry transcription [START_REF] Liu | Redundant function of REV-ERBalpha and beta and non-essential role for Bmal1 cycling in transcriptional regulation of intracellular circadian rhythms[END_REF] and therefore could mimic a Cry1 deficiency, and, on the other hand, synthetic REV-ERB agonists are available [Meng et al., 2008, Solt et al., 2012, Woldt et al., 2013].

The administration of a hypothetical REV-ERB agonist pulse at distinct times of the day was simulated in conditions of dampened AMPK rhythms (Figs. 6 et S5),

showing that the time of administration is critical. When the timing of the agonist pulse was optimally chosen, the physiological amplitude and phase could be restored for most genes as well as for NAD+ (Figs. 6 A , S 5 B , S 6 ), thus recovering the normal activity profile of sirtuins [Peek et al., 2013]. At other times of the day, however, the oscillation amplitude was not restored and additional phase shifts were induced, with potential adverse effects (Figs. 6B, S5A, S5C, S6). Timed REV-ERB agonist administration may thus correct the potentially harmful effect of misaligned clock profiles [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF] or abolished NAD+ rhythms [Peek et al., 2013] in situations where AMPK activity rhythms are dampened, due to a direct action of the diet [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF] or to systemic effects associated with metabolic diseases [START_REF] Coughlan | AMPK activation: a therapeutic target for type 2 diabetes?[END_REF].

Discussion

The circadian clock keeps internal physiological processes synchronized with external periodic cues such as the light-dark cycle or the alternation between feeding and fasting phases [Bass, 2012]. To investigate how the feeding/fasting cycle entrains the hepatic circadian clock, we built a mathematical model incorporating the metabolic sensors SIRT1 and AMPK as regulators of the canonical clock actors. SIRT1 and AMPK provide readouts of NAD+ and AMP levels, which reflect the energetic status of the cell.

The response of NAD+ and AMP levels to feeding or fasting combines many simultaneous biochemical processes [START_REF] Berg | Biochemistry[END_REF] and is difficult to describe.

In particular, an open problem is to understand the two-peak structure of NAD+ profiles observed in ad libitum normal chow diet fed mouse livers [START_REF] Ramsey | Circadian clock feedback cycle through NAMPT-mediated NAD+ biosynthesis[END_REF], Hatori et al., 2012]. The larger peak near ZT13 is due to the circadian regulation of Nampt, but the origin of the smaller one near ZT5 has remained elusive. Here, we kept the description of the dynamics of NAD+ and AMP as simple as possible. We assumed that NAD+ levels are mostly governed by the salvage pathway and slaved to the activity of the rate-limiting enzyme NAMPT. This neglects the influence of processes such as glycolysis, beta-oxidation [START_REF] Berg | Biochemistry[END_REF], or de novo biogenesis.

A key hypothesis was that NAMPT activity varies not only due to the circadian regulation of Nampt, but also because of the regulation of NAMPT stability by AMPK. It was motivated by the observations that the two NAD+ peaks are synchronized with peaks in NAMPT protein expression (Fig. 1 in [START_REF] Ramsey | Circadian clock feedback cycle through NAMPT-mediated NAD+ biosynthesis[END_REF]), and that the ZT5 peak also coincides with a peak in AMP level (Table S1 in [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF]).

Moreover, it is generally agreed on that AMPK activation increases NAD+ levels even though different mechanisms for this effect have been proposed [START_REF] Cantó | AMPK regulates energy expenditure by modulating NAD+ metabolism and SIRT1 activity[END_REF], Fulco et al., 2008, Brandauer et al., 2013]. Our model is the simplest one that is consistent with all these facts. The entrainment of the clock by the feeding/fasting cycle is thus described by how it responds to a cycling AMPK activity profile which is consistent with existing data for AMP levels, which display two peaks near ZT5 and ZT17 [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF]. However, both AMPK and SIRT1 are important to drive the clock, with NAD+ oscillations amplifying and relaying AMP oscillations.

Remarkably, our mathematical model reproduces simultaneously experimental clock gene expression data [Hughes et al., 2009] and NAD+ level data [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF], accounting for the two-peak pattern in NAD+ [START_REF] Ramsey | Circadian clock feedback cycle through NAMPT-mediated NAD+ biosynthesis[END_REF], Hatori et al., 2012] and NAMPT protein [START_REF] Ramsey | Circadian clock feedback cycle through NAMPT-mediated NAD+ biosynthesis[END_REF] profiles, and reproduces several mutant phenotypes. Despite its limitations, our model thus seems to capture interactions that are essential for coupling the circadian clock to metabolism, showing that behaviors observed in various experimental conditions are consistent with the simple molecular network of Fig. 1.

Our model also predicts that the smaller, non-circadian, NAD+ peak vanishes with both constitutively high and low AMPK activity profiles, which mimic fasting and high-fat diet conditions, respectively [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF], Peek et al., 2013]. In the experimental data we used, the time intervals between the two AMP peaks and the two NAD+ peaks are 12 and 8-9 hours, respectively [Hughes et al., 2009, Hatori et al., 2012]. Rhythms with period lengths of 12 and 8 hours have been reported for hundreds of transcripts in liver, but also in other tissues [Hughes et al., 2009, Chiang et al., 2014, Vollmers et al., 2009]. However, the mechanisms generating these rhythms remain unclear, although some mechanisms have been proposed [START_REF] Westermark | Mechanism for 12 hr rhythm generation by the circadian clock[END_REF]. Harmonics are generally lost in vitro, suggesting that they are due to systemic cues not reproduced in vitro, but also in vivo when food access is restricted to the light period [Hughes et al., 2009]. Our study reinforces the idea that some ultradian rhythms, here those in NAD+ and NAMPT protein profiles, can be generated by the interplay between a circadian rhythm and a feeding rhythm.

The effects of AMPK and SIRT1 loss of function [Lamia et al., 2009, Bellet et al., 2013] suggest that perturbations of the fasting/feeding cycle disrupt the clock. We used our model to predict how changes in the AMP rhythms would affect the clock, motivated by the fact that a general decrease in AMP levels has been associated with dampened clock gene expression oscillations in mice subjected to AL HFD [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF].

AL HFD was also linked to drastically reduced oscillations in Nampt expression and NAD+ levels [Eckel-Mahan and Sassone-Corsi, 2013] and to obesity [Hatori et al., 2012, Eckel-Mahan andSassone-Corsi, 2013]. The loss of NAD+ oscillations has been shown to impair mitochondrial function [Peek et al., 2013], a deleterious effect of the HFD. Decreased amplitudes of clock gene expression have also been described in aged mammals [Froy, 2011], raising the question whether similar mechanisms are at work.

We found that our model indeed predicts that depressed AMP rhythms significantly reduce the amplitude of Nampt, NAD+, and clock gene expression oscillations, supporting the conclusion that these effects are associated [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF], and explaining why they can be observed within a few days of HFD feeding [Eckel-Mahan et al., 2013]. This is consistent with the reports that AMPK is required for the circadian oscillations of NAD+ and Nampt gene expression in the heart, skeletal muscle and fat [START_REF] Um | AMPK regulates circadian rhythms in a tissue-and isoform-specific manner[END_REF], and that mice on a time-restricted HFD, which restores the cycling of AMP, display normal clock rhythms [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF]. Such mice are also resistant to HFD-induced obesity, raising the question of a link between dampened oscillations and obesity. Interestingly, our model predicts that Cry1-/-mutants, which have recently been reported as resistant to HFD-induced obesity [START_REF] Griebel | Mice deficient in Cryptochrome 1 (Cry1 (-/-)) exhibit resistance to obesity induced by a high-fat diet[END_REF], maintain high amplitude oscillations when AMP rhythms are low. Thus, an important question is whether the adverse effects of dampened AMPK rhythms can be alleviated by a pharmacological approach.

The possibility of modulating REV-ERB activity has been recently reported [Meng et al., 2008, Solt et al., 2012, Woldt et al., 2013]. We thus simulated the administration of a hypothetical REV-ERB agonist. When it was delivered in the early night, the amplitude of rhythms in clock gene expression was restored with the correct phase, but administration at other time points led to strong phase shifts without changing the amplitude. This suggests that oscillations in clock gene expression can be rescued pharmacologically, but also highlights the danger of mistimed delivery of drugs with circadian targets [Meng et al., 2008], which is the case of many commonly prescribed drugs [START_REF] Zhang | A circadian gene expression atlas in mammals: Implications for biology and medicine[END_REF]. The optimal administration time obtained here must be refined, as it depends on a predicted REV-ERB profile which was not experimentally constrained. We believe that further studies will provide additional insights into the general principles of such pharmacological approaches and confirm that it is a robust strategy to enhance clock gene oscillations. Whether or not it would be effective in treating HFD-induced obesity remains to be shown, however rescuing the circadian NAD+ peak (Fig. 6 ) is likely very important since this peak is required for normal mitochondrial oxidative metabolism [Peek et al., 2013]. This result is all the more promising as a systematically decreased AMPK activity is observed in metabolic diseases such as type 2 diabetes [START_REF] Coughlan | AMPK activation: a therapeutic target for type 2 diabetes?[END_REF].

These results call for further investigations. Here, we studied only the entrained liver clock. The next step is to study how it resets when AMP rhythms are phase-shifted from day to day, mimicking the effect of fluctuating food timings. F o r e x a m p l e , i s the hepatic clock robust to fluctuations in the cycle that entrains it? This question was recently studied in an algal clock subjected to daylight fluctuations [Thommen et al., 2010, Morant et al., 2010, Thommen et al., 2012, Pfeuty et al., 2012]. To address this question, systemic cues from the SCN, which influence the resetting dynamics [Saini et al., 2013], must also be taken into account. A better modeling of PGC1α regulation is also required. Finally, the model should incorporate other actors mediating feeding/fasting cycles to clock, most notably insulin and glucagon [START_REF] Mukherji | Shifting the feeding of mice to the rest phase creates metabolic alterations, which, on their own, shift the peripheral circadian clocks by 12 hours[END_REF], allowing us to study how they cooperate with AMPK and SIRT1.

In conclusion, our results show that a simple, but quantitative, description of the entrainment of the circadian clock by metabolism is possible. They suggest that AMPK activity plays a key role in the metabolic entrainment of peripheral clocks and that an important clock input may be a post-transcriptional regulation of NAMPT protein abundance by AMPK. This supports the idea that strong AMP rhythms are essential for maintaining normal clock operation, emphasizing the importance of a sufficient long fasting period inside the diurnal cycle. Remarkably, we show how clock perturbations resulting from a challenged feeding/fasting cycle can be rescued by a well-timed pharmacological intervention. More generally, this modeling study raises the question of the importance of food and exercise timing to avoid metabolic disorders.

Experimental procedures

Experimental data for clock genes and Nampt expression in normal mouse livers [Hughes et al., 2009] were obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo/, accession number GSE11923). Experimental data for the NAD+ and AMP levels in livers of mouse fed a normal chow ad libitum were obtained from Table S1 of [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF]. Parameter estimation and numerical simulation were carried out using the COPASI pathway simulator [Hoops et al., 2006] and MATLAB [START_REF] Mat | Matlab release 2011b[END_REF]. although the smaller NAD+ peak disappears, and decreased in (B). To assess the accuracy of the prediction, we computed a reference target profile for each mutant phenotype. This profile was constructed by computing the fold change in expression between the WT and mutant phenotypes in [Bellet et al., 2013] and [Lamia et al., 2009],

Figure legends

and applying the same ratio to the WT profile used shown here and in Fig. 2. Table S1 List of variables of the mathematical model. Relates to Fig. 1.

Table S2 List of the differential equations and mathematical relations defining the mathematical model. Relates to Fig. 1.

Table S3 List of the kinetic constants obtained by parameter identification. For the convenience of the reader, Table S3 is split into Tables S3A-S3L. Experimental data for NAD+ level (green dots), AMP level (red dots) [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF], and AMPK protein abundance (magenta dots) [Barnea et al., 2012] are shown together with the putative AMPK profile used in our model (blue) and the Fourier series matching Nampt mRNA profile data (magenta line) from [Hughes et al., 2009]. The NAD+, AMP and AMPK abundance data are drawn with a interpolating periodic cubic spline to guide the eye. The key observations are the coincidences of (1) the AMP and NAD+ peaks near ZT5, and (2) of the Nampt mRNA and NAD+ peaks near ZT13. To help the reader appreciate these coincidences, the red and green vertical lines indicate the timings of estimated maxima of the AMP and NAD+ profile, respectively. The model AMPK profile is consistent with AMP level, AMPK abundance and NAD+ level profiles.

Supplemental Experimental Procedures

Supplemental Tables

+     Vmax_per • 1 + fold_per • [CB] Ka_per_cb • (1 + Act_SIRT) hill_per_cb 1 + [CB] Ka_per_cb • (1 + Act_SIRT) hill_per_cb • 1 + [PC] Ki_per_pc hill_per_pc     d[cry] dt = -dm_cry • [cry] (S2) +     Vmax_cry • 1 + fold_cry • [CB] Ka_cry_cb • (1 + Act_SIRT) hill_cry_cb 1 + [CB] Ka_cry_cb • (1 + Act_SIRT) hill_cry_cb • 1 + [PC] Ki_cry_pc hill_cry_pc • 1 1 + [Prot_rev] Ki_cry_rev hill_cry_rev     d[rev] dt = -dm_rev • [rev] (S3) +     Vmax_rev • 1 + fold_rev • [CB] Ka_rev_cb • (1 + Act_SIRT) hill_rev_cb 1 + [CB] Ka_rev_cb • (1 + Act_SIRT) hill_rev_cb • 1 + [PC] Ki_rev_pc hill_rev_pc     d[ror] dt = -dm_ror • [ror] (S4) +     Vmax_ror • 1 + fold_ror • [CB] Ka_ror_cb • (1 + Act_SIRT) hill_ror_cb 1 + [CB] Ka_ror_cb • (1 + Act_SIRT) hill_ror_cb • 1 + [PC] Ki_ror_pc hill_ror_pc     d[bmal] dt = -dm_bmal • [bmal] (S5) +     Vmax_bmal • 1 + fold_bmal • (1 + Act_PGC1a) • [Prot_ror] Ka_bmal_ror hill_bmal_ror 1 + [Prot_rev] Ki_bmal_rev hill_bmal_rev + [Prot_ror] Ka_bmal_ror hill_bmal_ror     d[Prot_per] dt = -dp_per • (1 + m_per_sirt • Act_SIRT + m_per_ampk • Act_AMPK) • [Prot_per] (S6) + kp_per • [per] -(kass_pc • [Prot_cry] • [Prot_per] -kdiss_pc • [PC]) d[Prot_cry] dt = -dp_cry • (1 + m_cry_ampk • Act_AMPK) • [Prot_cry] (S7) + kp_cry • [cry] -(kass_pc • [Prot_cry] • [Prot_per] -kdiss_pc • [PC]) d[Prot_rev] dt = -dp_rev • [Prot_rev] (S8) + kp_rev • [rev] d[Prot_ror] dt = -dp_ror • [Prot_ror] (S9) + kp_ror • [ror] d[Prot_bmal] dt = -dp_bmal • [Prot_bmal] (S10) + (kp_bmal • [bmal]) -(kass_cb • [Prot_bmal] -kdiss_cb • [CB]) d[PC] dt = + (kass_pc • [Prot_cry] • [Prot_per] -kdiss_pc • [PC]) (S11) -d_pc • [PC] d[CB] dt = + (kass_cb • [Prot_bmal] -kdiss_cb • [CB]) (S12) -d_cb • [CB] d[nampt] dt = -dm_nampt • [nampt] (S13)
Table S3A: mRNA and protein degradatation rates (15 parameters). ( 1) Most values correspond to lifetimes between one and a few hours, which is reasonable.

(2) The high value of NAMPT protein degradation rate may be linked to the fact that the peak of the Nampt mRNA profile reconstructed from experimental data of Hughes et al. lags the NAD+ peak in Hatori data [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF] by 0.5-0.75 hours, in contradiction with the strict dependence of NAD+ level on NAMPT protein abundance assumed in the model. Indeed, the NAMPT protein peak, itself lagging the Nampt mRNA peak, should be synchronized with the NAD+ peak. This temporal discrepancy can be explained by experimentalthe uncertainty due to the limited time resolution and the reconstruction procedure, as well as by the data coming from two independent experiments. By increasing NAMPT degradation rate as much as possible, the optimization procedure basically seeks to reduce the lag between Nampt mRNA and NAMPT protein peaks, in order to keep the lag between NAD+ and NAMPT protein peaks to a minimum. (3) The value of the ROR protein degradation rate is very low, indicating that its abundance in our model is almost constant. Thus, ROR plays a role in clockwork dynamics mainly via its co-activator PGC1a. (4)

The relatively high value of the Rev-Erb mRNA degradation rate compared to other mRNAs is consistent with the sharp non-sinusoidal shape of the experimental profile which indicates fast variations not compatible with a long lifetime.

Parameter Value (h S3D: Activation ratios (7 parameters). The activation ratio is the ratio between the transcription rates at full and zero activator concentrations, in the absence of any inhibitor (e.g., PER-CRY) or repressor (e.g., REV-ERB). The ratio between the maximal and minimum transcription rates can thus be much higher in the case of strong inhibition or repression. Low values of the activation ratio typically indicate that the transcription rate is more regulated by the inhibitor or the repressor than by the activator. This is for example the case of Cry, whose expression is known to increase in Bmal1 KO because the repression by REV-ERB is down-regulated more than the activation by CLOCK-BMAL. (1) The regulation threshold of gene X by protein Y is the concentration of Y at which transcription of X is half of its maximum value (which is reached in the absence of any inhibitor or repressor for an activator, and with the activator in excess for a repressor). ( 2) The regulation threshold indicates the range of values in which variations of the regulator concentration influence transcription the most, and thus should be compared to the maximum value of the regulator. For reference, the maximum values of the four regulator concentrations are :

[CB] max ∼ 0.36, [PC] max ∼ 10.2, [Prot_rev] max ∼ 0.51, [Prot_ror] max ∼ 1.83.
The larger the Hill coefficient, the narrower the range of concentrations in which variations of Y have a significant effect. (3) The activation ratio, the regulation threshold, the Hill coefficient ratio, etc., are closely related parameters which together determine the variation of regulation over one cycle. Thus, it may be difficult to interpret the value of an isolated parameter without taking into account the others. S3F: Hill coefficients (15 parameters). It is often considered that the Hill coefficient should be equal to the number of boxes in the promoter of the target gene (e.g., Rev-Erbα is activated through 3 E-boxes and thus hill_rev_cb should be 3) [Korencic et al., 2012]. We tested adjustment with hill_rev_cb equal to 3, and found that we could harldy reproduce the very sharp profile of Rev-Erba. This seems to indicate that the effective cooperativity is significantly higher than the number of E-boxes, probably due to the influence of additional regulation mechanisms beyond those taken into account in the model. Therefore, we did not constrain Hill coefficients, so that the value reported may not be the smallest one leading to a good adjustment. Supplemental Experimental Procedures

Mathematical model

The list of variables of the model is given in Table S1. The differential equations and mathematical relations governing their time evolution are liste in Table S2.

Eqs (S1)-( S4) and (S13) assume that the transcription rates of Per, Cry, Rev-Erb, Ror and Nampt are enhanced upon CLOCK-BMAL1 binding to their promoters, that this activation is inhibited by PER-CRY binding to CLOCK-BMAL1 on DNA, and that binding of CLOCK-BMAL1 to DNA is inhibited by SIRT1. Eq.( S2) also takes into account the repression of Cry by REV-ERB.

The formula for the transcription rate then follows by assuming thermodynamic equilibrium and by considering the fractions of time spent by the gene in its different states (free, bound by CLOCK-BMAL1, bound by CLOCK-BMAL1 and PER-CRY).

Eq. (S6) assumes that the PER protein is destabilized directly by SIRT1 and indirectly by AMPK (via CKI). Eq. (S7) assumes that CRY is destabilized by AMPK. Eq. (S14) assumes that NAMPT is stabilized by AMPK.

NAD kinetics, descrbibed by Eq. (S15) assumes that the total pool of NAD+, NAM, and NMN is constant (and is given by NAD total , with NAM and NMN grouped into a single non-NAD species since only the step from NAM to NMN is rate-limiting. We describe the conversion from NAD+ to NAM as an enzymatic process with constant rate, taking into account the approximately constant SIRT1 abundance, but neglecting the variations in the quantity of proteins being deacetylated. Moreover, we assume that this conversion does not operate below some threshold value NAD basal , consistent with the fact that many experimental NAD+ profiles diplay long plateaus of almost constant finite concentration (in particular in the fed state). On the contrary, the recycling of NAM/NMN is assumed to be catalyzed by NAMPT.

The values of the kinetic constants are given in Tables S3A-S3L. Many of these values are not uniquely determined, because there exist sets of parameters which may be changed in a coordinated way without degrading goodness of fit. Moreover, some parameter values depend on the mRNA and protein absolute levels which are not known and have been arbitrarily fixed here. As a matter of fact, the only parameters which have an absolute meaning are the degradation rates because they are given in inverse time units. Since mRNA and protein concentrations are only determined up to a scale factor, no effort was made to give them realistic values.

Mathematical description of pulses

The AMPK activity pulses, the PGC1a protein pulse and the agonist pulse were described mathematically by adapting the Input Signal Step Function (ISSF) introduced in [START_REF] Adams | The Input Signal Step Function (ISSF), a standard method to encode input signals in SBML models with software support, applied to circadian clock models[END_REF].

The ISSF is a multi-parametric periodic function which can reproduce a wide range of waveforms and is easily encoded in the Systems Biology Markup Language (SBML), which has been designed to describe models in systems biology. The use of a parametric profile allowed us to subject the AMPK profile to parameter optimization without a priori knowledge of the profile shape.

More precisely, we described the time evolution of a single pulse by function P (t, t c , T d , T w , T c ), whose arguments are the time t, the peak timing t c , the duration T d of the pulse (time during which the pulse is above half of the maximum), the time scale T w of the rising and decreasing fronts, and the period T c of the external cycles, and which is given by

P (t, t c , T d , T w , T c ) =S(t, t c -T d /2, 0, T w , T c ) -S(t, t c -T d /2, T d , T w , T c ) (S24) + S(t, t c -T d /2, T c , T w , T c ) -S(t, t c -T d /2, T c + T d , T w , T c )
where S is the step function:

S(t, t f , t s , T w , T c ) = 1 2 1 + tanh τ (t -t f , T c ) -t s T w (S25)
with t the time, t f the time location of the center of the step, t s a time shift, and T w is the time scale of the step. In (S25), τ (t, T c ) is the time elapsed since the last beginning of the day (ZT0), given the time t and cycle period T c :

τ (t, T c ) =t -T c • t T c (S26)
where x is the largest integer smaller than or equal to x. Expression (S25) guarantees that S gradually rises from 0 to 1 in a time interval given by 2T w .

A typical waveform described by the P function is shown in Fig. S1, which also clarifies the meaning of the t c , T d , and T w parameters. The structure of expression (S24) can be understood as follows. The first line is the difference of two step functions like those shown in red in Fig. S1, shifted by T d . This difference becomes high between the two rising fronts but vanishes elsewhere.

The second line in (S24) is the same expression as in the first line, but shifted by 24 hours so that the sum of the two lines is 24h-periodic given the dependence on τ (t, T c ). The function P in (S24) oscillates between 0 and 1 (asymptotically), but an expression oscillating between a and b, if needed, can be obtained as a + (b -a) • P .

Note that expression (S24) differs from that given in [START_REF] Adams | The Input Signal Step Function (ISSF), a standard method to encode input signals in SBML models with software support, applied to circadian clock models[END_REF] by the last term in the second line. This term, which was missing in [START_REF] Adams | The Input Signal Step Function (ISSF), a standard method to encode input signals in SBML models with software support, applied to circadian clock models[END_REF], is needed to make the function truly 24h-periodic.

Fourier fitting of the mRNA expression data

The utilization of the gene expression data from Hughes et al. [Hughes et al., 2009] was motivated by their high time resolution and their relative good reproducibility from the first day to the second day (gene expression was sampled from mouse livers each hour during 48 hours). Yet, the expression profiles were not perfectly periodic, unlike the solution curves of the mathematical model. Therefore, we fitted the expression time profiles to Fourier series in order to construct target profiles that are compatible with the numerical profiles.

Since the waveform of some genes (e.g., Rev-Erba, Bmal1) was clearly non-sinusoidal, we used Fourier series with 4 frequencies of the form:

S = a 0 + 4 k=1 s k sin 2kπ T + c k cos 2kπ T (S27)
All expression profiles were adjusted simultaneously, assuming a common period T . We found that the best global fit was obtained when T = 24.01 hours, which was rounded to 24 hours.

The resulting Fourier series are compared to the original raw time series in Fig. S3. They were used as target profiles in the subsequent parameter identification procedure.

Parameter identification

Parameter identification was carried out with the "Parameter estimation" module of the Copasi pathway simulator [Hoops et al., 2006], which minimizes an error function quantifying the discrepancies between the target profiles and the numerical solution for a given parameter set. We tested several of the algorithms proposed and found that the Hooke and Jeeves method yielded the best results. It is a derivative-free pattern search algorithm, similar to the simplex method.

Although relatively slow, this algorithm was much less prone to be blocked in a local minimum. Still, obtaining our final result required a series of educated guesses, correcting manually the most obvious discrepancy before restarting the optimization. We found that parameter values are not uniquely constrained, as different parameter sets which provide almost the same goodness of fit. A more detailed analysis of the relative parameter uncertainties will be given elsewhere.

Supplemental Results

Analysis of the NAD+, AMP and Nampt mRNA profiles and relation to the model AMPK profile

Two key hypotheses in our model are that (1) NAD+ levels are essentially determined by the salvage pathway, hence by NAMPT protein abundance; (2) NAMPT protein stability is regulated by AMPK (or more generally that NAMPT protein abundance is enhanced by AMPK). In this section, we provide support to these hypotheses by examining the NAD+ and AMP data from Hatori et al. [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF] together with the Fourier series profile for Nampt, reconstructed from the gene expression data of Hughes et al. [Hughes et al., 2009] (Fig. S2).

Regarding hypothesis (1), Fig. S2 shows that the big NAD+ peak near ZT13 is well synchronized with the reconstructed Nampt expression profile, which is all the more remarkable that they come from independent experiments. Curiously, the Nampt peak is lagging slightly behind the NAD+ peak, whereas the latter should occur after or close to the NAMPT protein peak, which should itself be preceded by the mRNA peak. This may be due to either the uncertainty in the reconstruction of the profile or to a small temporal misalignment between the Hughes and Hatori experiments.

Hypothesis (2) is supported by the almost perfect coincidence of the AMP and NAD+ peaks near ZT5, by the fact that AMPK activation upregulates NAMPT [Fulco et al., 2008, Yang et al., 2007, Brandauer et al., 2013] and, most importantly, by the fact that the NAMPT protein profile displays two peaks, including one near ZT5, as shown in Fig. 1 of [START_REF] Ramsey | Circadian clock feedback cycle through NAMPT-mediated NAD+ biosynthesis[END_REF]. The latter observation also supports hypothesis (1) since the two NAD+ peaks coincide with the two NAMPT protein peaks [START_REF] Ramsey | Circadian clock feedback cycle through NAMPT-mediated NAD+ biosynthesis[END_REF]. It should be mentioned that we also tested the hypothesis of an up-regulation of Nampt translation in our model, and that it was also compatible with the data. Therefore, our work points clearly at a post-translational regulation of Nampt, but more work is needed to identify the exact mechanism.

To summarize, the structure of the NAD+ temporal profile appears to be tightly correlated to the variations of Nampt expression and of AMP, providing support to our hypotheses.

The action of AMPK on the clock plays a critical role in our analysis. In this first approach, this action is described through an activity profile fixed externally. It is therefore important to check that the profile used in the analysis is consistent with experimental data available in the litterature.

The AMPK profile used in our analysis (Fig. S2), was obtained by model adjustment (described in Sec ). However this procedure generally requires a good initial guess because of the existence of many local minima. Since AMPK is primarily activated by AMP [Hardie et al., 2012], we hypothesized that each peak in the AMP level of [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF] would trigger a peak of AMPK activity, and therefore decided to describe AMPK activity by a two-pulse mathematical function (Sec. ), each pulse being synchronized with a peak in AMP (Fig. S2). We found that if the amplitudes of the two peaks were chosen in the same ratio as the two AMP peaks, then the action of AMPK would be much too strong near the second peak and distort the profiles, or it would be too weak near the first peak and the small NAD+ peak would not be reproduced. Thus we manually adjusted the ratio to be in the order of 4 to 6, and obtained the profile shown in Fig. S2 after adjustment.

However, AMPK activity depends not only on AMP level but also on AMPK expression and on SIRT1 activity since deacetylation of LKB1 increases AMPK phosphorylation [Lan et al., 2008, Hou et al., 2008]. Therefore, we compared our AMPK activity profile to an AMPK expression profile in liver from Barnea et al. [Barnea et al., 2012] and to the NAD+ profile from Hatori et al. [START_REF] Hatori | Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet[END_REF], which are shown in Fig. S2. The 5-fold decrease in AMPK expression between ZT3 and ZT17 naturally explains most of the difference between the amplitudes of the two AMPK activity peaks. Moreover, the slight advance of the AMPK activity peak relative to the AMP peaks can also explained by the fact that the peaks in AMPK expression and in NAD+ level (hence SIRT1 activity) do not coincide with AMP level peaks. More precisely, the first putative AMPK activity peak is about halfway between an AMPK expression peak and an AMP level peak, which is expected since AMPK activity should roughly be proportional to both. Similarly, the second putative AMPK peak is about halfway between peaks in NAD+ level and AMPK expression on the one hand, and the AMP peak on the other hand.

Globally, the agreement between the adjusted AMPK profile and the experimental data is thus extremely satisfying, given that we have considered together data from different experiments, and that we have not taken into account the possible role of AMPK nuclear localization.

Note that the model is not sensitive to absolute AMPK activity levels, only to the variations of some kinetic constants between the low and high activity levels. This makes our analysis largely independent of an accurate description of the relation between AMP levels and AMPK activity levels.

More precisely, assume that the modulation of a degradation rate d by AMPK is described by d = d 0 (1 + cA), where A is the AMPK activity level. Then the same degradation rate can be rewritten as d = d 0 (1 + c A ) where d 0 = d 0 (1 -cA 0 ), A = A + A 0 , and c = c/(1 -cA 0 ), essentially shifting the zero level of AMPK activity. Actually, the only quantity which is biologically observable is the variation of the degradation rate between low and high AMPK activity. Thus, the model should get this variation right, the rest is arbitrary. In this perspective, AL HFD is characterized by the absence of an alternation between high activity and low activity and thus between different degradation rates. It is this oscillation that counts.

Figure 1 :

 1 Figure 1: Network of molecular interactions taken into account in our mathematical model to describe the driving of the core mammalian liver clock (in the grey box) by metabolism via the two metabolites AMP and NAD+ (golden yellow ovals, with NAM, the inactive form of NAD+, in light yellow). The two metabolic sensors AMPK and SIRT1 appear as green boxes and their actions on the clock are indicated with green arrows. mRNAs are represented by salmon slanted boxes, proteins by red square boxes and protein complexes by orange ovals. Acetylation and phosphorylation are indicated with small red and green circles.

Figure 2 :

 2 Figure 2: Adjustment of the mathematical model to experimental data from in vivo experiments using WT mice. (A) the predicted profiles for clock gene expression and NAD+ level (solid blue lines) are compared to experimental data (red dots). (B) the time profiles of clock genes are plotted together to show the relative phases of maximum expression. (C) AMPK activity (blue), NAMPT protein (red) and NAD+ (green) profiles are shown, highlighting their synchronization. See also Figs. S1, S2 and S3.

Figure 3 :

 3 Figure 3: Effect of SIRT1 KO (A) and impaired AMPK activation via LKB1 KO (B) on clock gene expression and NAD+ profiles. WT and mutant phenotypes are shown in blue and red, respectively. The amplitude of oscillations is typically increased in (A),

Figure 4 :

 4 Figure 4: Influence of AMPK activity rhythms on the clock dynamics for a normal state (blue), a fed-like state with constantly low AMPK activity (red) and a fasted-like state

Figure 5 :

 5 Figure 5: Oscillations in the expression of the Nampt, Per and Dbp genes in Cry1 KO mutants exposed to dampened AMPK rhythms (orange), compared to WT cells facing the same challenge (red) or under normal AMPK rhythms (blue).

Figure 6 :

 6 Figure 6: Rescue of oscillations in the Nampt (top) and Rev-Erb (bottom) expression profiles of cells subjected to dampened AMPK rhythms, using a REV-ERB agonist. The profiles corresponding to treated and non treated cells are shown in orange and red,respectively, and compared to those corresponding to WT cells under normal AMPK rhythms (dashed blue). In the A panel, the agonist pulse (top inset) is optimally timed around ZT 13.7, leading to restored profiles, unlike when the agonist pulse is shifted by

  Figure 4

Figure S2

 S2 Figure S2 Experimental data from the literature for NAD+, AMP and AMPK protein abundance time profiles, target profile for Nampt mRNA and model AMPK activity profile. Relates to Fig. 2.

Figure S3

 S3 Figure S3 Fitting of experimental clock gene expression time series by Fourier series. Relates to Fig. 2.

Figure S4

 S4 Figure S4 Gene expression and NAD+ level time profiles predicted by the model for the fed, fasted and normal diet states. Relates to Fig. 4.

Figure S5

 S5 Figure S5 Clock gene expression and NAD+ level time profiles predicted by the model when an Rev-Erb agonist is delivered in the fed state at three different times of the day. Relates to Fig. 6.

Figure S6

 S6 Figure S6 Residual error between some time profiles generated by the agonist-treated clock model and the target experimental profiles. Relates to Fig. 6.

Figure S1 :

 S1 Figure S1: Structure of the AMPK activity pulse in the model and meaning of the pulse parameters t c , T d , and T w . A pulse profile described by (S24) is shown in blue, which is essentially obtained as the difference between the two step functions shown in red, modified to be 24hperiodic.

Figure S4 :Figure S5 :Figure S6 :

 S4S5S6 Figure S4: Gene expression and NAD+ level profiles in the fed (red), fasted (green) and normal diet (blue) states. The corresponding AMPK activity profiles are shown in the top left plot.

Core clock Cry Ror PGC1α P PGC1α Ac PGC1α Ac P

	Figure 1 0 0 1 2 Per2 mRNA A 0 0 1 2 Ror γ mRNA 0 0 1 2 3 Dbp mRNA Figure 2 A Figure 3 A 0 0 1 2 3 Nampt mRNA B 1 2 AMPK activity 4 B 1 Nampt 0 0 1 2 3 2 mRNA Rev-Erbα mRNA 0 0 3 B	12 12 12 WT Time (h) 24 24 Time (h) 24 Time (h) Normal diet 36 36 36 SIRT1 KO 48 48 48 6 12 18 24 30 36 42 48 0 12 0 1 Cry1 mRNA 0 12 0 1 2 Ref. Fed Time (h) 12 Time (h) 24 24 Time (h) Fasted 36 36 0 6 12 18 24 30 36 42 48 48 48 0 12 24 Time (h) 0 1 C Nampt mRNA 2.2 Model Experimental data 36 48 Time (h) 0 3 6 9 Dbp mRNA 0 6 12 18 24 30 36 42 48 Time (h) 1 2 3 NAD level WT LKB1 KO Ref. 2 3 NAD+ level Bmal1 mRNA 0 12 24 36 48 Time (h) 0 1 2 3 Rev-Erbα mRNA 2.5 Per2 Ror γ Dbp Cry1 Bmal1 Rev-Erbα Nampt 0 12 24 36 Time (h) C AMPK NAMPT NAD+ 48 1.0 1.2 1.4 1.6 1.8 2.0 NAD+ 6 12 18 24 30 36 42 48 Time (h) 0 6 12 18 24 30 36 42 48 Time (h) 0 0 6 12 18 24 30 36 42 48 Time (h) 0 1 2 3 Dbp mRNA 1 2 6 12 18 24 30 36 42 48 Time (h) 0 6 12 18 24 30 36 42 1 48 Time (h) Per2 mRNA
		3 0 0	6	12	18	24	30	36	42	48	2.0 0	6	12	18	24	30	36	42	48
	mRNA	2				Time (h)					1.5				Time (h)
											1.0
		1									0.5
		0 0	6	12	18	24	30	36	42	48	0 0.0	12	24	36	48
						Time (h)							Time (h)

18 24 30 36 42 48 54 60 Bmal1 18 24 30 36 42 48 54 60 Cry1 18 24 30 36 42 48 54 60 Per2 18 24 30 36 42 48 54 60 Rev-Erbα 18 24 30 36 42 48 54 60 Rorγ 18 24 30 36 42 48 54 60 Nampt Figure

  S3: For each of the genes taken into account in the model, the gene expression profiles obtained by Hughes et al.[Hughes et al., 2009] (in red) are compared to the Fourier series best fitting them (in black).

		Normal diet	Fed		Fasted
	AMPK activity	0 6 12 18 24 30 36 42 48 Time (h) 0 1 2 Nampt mRNA	0 6 12 18 24 30 36 42 48 Time (h) 0 1 2 3	Time (h) 0 6 12 18 24 30 36 42 48 1 2 3 NAD+ level
	0 6 12 18 24 30 36 42 48 Time (h) 0 1 2 Bmal1 mRNA	0 6 12 18 24 30 36 42 48 Time (h) 0 2 4 6 8 Rev-Erb mRNA	Time (h) 0 6 12 18 24 30 36 42 48 0 1 2 Ror 3 4 5 6 mRNA
	0 6 12 18 24 30 36 42 48 Time (h) 0 Per Cry 1 2 mRNA mRNA	0 6 12 18 24 30 36 42 48 Time (h) 0 Dbp 1 mRNA	Time (h) 0 6 12 18 24 30 36 42 48 0 2 4 6

Table S1 :

 S1 List of variables of the mathematical model. The time evolution of the variable is specified either by an ordinary differential equation (ODE), by a function of other variables, or by an externally imposed profile.

	Variable name Type	Description
	[per]	ODE	Concentration of Per mRNA
	[cry]	ODE	Concentration of Cry mRNA
	[rev]	ODE	Concentration of Rev-Erb mRNA
	[ror]	ODE	Concentration of Ror mRNA
	[bmal]	ODE	Concentration of Bmal1 mRNA
	[Prot_per]	ODE	Concentration of PER protein
	[Prot_cry]	ODE	Concentration of CRY protein
	[Prot_rev]	ODE	Concentration of REV-ERB protein
	[Prot_ror]	ODE	Concentration of ROR protein
	[Prot_bmal]	ODE	Concentration of BMAL1 protein
	[PC]	ODE	Concentration of PER-CRY protein complex
	[CB]	ODE	Concentration of CLOCK-BMAL1 protein complex
	[nampt]	ODE	Concentration of Nampt mRNA
	[Prot_nampt] ODE	Concentration of NAMPT protein
	[dbp]	ODE	Concentration of Dbp mRNA
	[NAD]	ODE	NAD level
	Act_SIRT	function Activity of SIRT1
	Act_PGC1a	function Activity of PGC1a
	Act_AMPK	external Activity of AMPK
	Prot_PGC1a external Nuclear abundance of PGC1a
	agonist_rev	external Concentration of REV-ERB agonist

Table S2 :

 S2 Differential equations and mathematical relations defining the mathematical model.

	d[per] dt	=	-dm_per • [per]	(S1)

Table S3B :

 S3B Complexation kinetic rates (4 parameters)

	Parameter Value	Unit	Description
	kass_cb	0.0162923358655	nmol -1 • l • h -1 CLOCK-BMAL association rate
	kass_pc	12.302005485	nmol -1 • l • h -1 PER-CRY association rate
	kdiss_cb	0.00613502224231 h -1	CLOCK-BMAL dissociation rate
	kdiss_pc	0.0365867175408	h -1	PER-CRY dissociation rate
	Table S3C: Maximal transcription rates (7 parameters)
	Parameter	Value (nmol • h -1 ) Description
	Vmax_bmal	0.0916862770193	Bmal maximal transcription rate
	Vmax_cry	0.702216664807		Cry maximal transcription rate
	Vmax_dbp	0.0802009609453	Dbp maximal transcription rate
	Vmax_nampt 3.49035201578		Nampt maximal transcription rate
	Vmax_per	0.730201742662		Per maximal transcription rate
	Vmax_rev	1.12297601784		Rev-Erb maximal transcription rate
	Vmax_ror	-1 ) 6.9843472736	Description Ror maximal transcription rate
	dm_bmal	0.827333442085	Bmal mRNA degradation rate
	dm_cry	0.319848706181	Cry mRNA degradation rate
	dm_dbp Table	0.384486532062	Dbp mRNA degradation rate
	dm_nampt 0.814311309051	Nampt mRNA degradation rate
	dm_per	0.323114598647	Per mRNA degradation rate
	dm_rev	4.0479072846	Rev-Erb mRNA degradation rate
	dm_ror	0.246575760727	Ror mRNA degradation rate
	dp_bmal	0.186413466797	BMAL protein degradation rate
	dp_cry	0.599026119971	CRY protein degradation rate
	dp_nampt 49.8841023982	NAMPT protein degradation rate
	dp_per	10.9446515392	PER protein degradation rate
	dp_rev	0.281564063057	REV-ERB protein degradation rate
	dp_ror	0.0340112196281 ROR degradation rate
	d_cb	0.197714012552	CLOCK-BMAL complex degradation rate
	d_pc	0.609290138329	PER-CRY complex degradation rate

Table S3E :

 S3E Regulation thresholds (15 parameters).

	Parameter	Value (dimensionless) Description
	fold_bmal	15.9258093373	activation ratio of Bmal by ROR
	fold_cry	1.1604489571	activation ratio of Cry by CLOCK-BMAL
	fold_dbp	400.0	activation ratio of Dbp by CLOCK-BMAL
	fold_nampt 1.57880681573	activation ratio of Nampt by CLOCK-BMAL
	fold_per	12.977351123	activation ratio of Per by CLOCK-BMAL
	fold_rev	73.2342431701	activation ratio of Rev-Erb by CLOCK-BMAL
	fold_ror	335.923333883	activation ratio of Ror by CLOCK-BMAL

Table S3G :

 S3G Translation rates (6 parameters). Again, the high value of kp_nampt may be forced by the possible temporal misalignment between the reconstructed Nampt expression profile and the NAD+ profile.

	Parameter	Value (dimensionless) Description
	hill_bmal_rev 4.32985205032	Hill coeff., regulation of Bmal by REV-ERB
	hill_bmal_ror 1.83992599578	Hill coeff., regulation of Bmal by ROR
	hill_cry_cb	9.1109447538	Hill coeff., regulation of Cry by CLOCK-BMAL
	hill_cry_pc	2.43715119318	Hill coeff., regulation of Cry by PER-CRY
	hill_cry_rev	4.20952050286	Hill coeff., regulation of Cry by REV-ERB
	hill_dbp_cb	7.32066818222	Hill coeff., regulation of Dbp by CLOCK-BMAL
	hill_dbp_pc	10.4312927466	Hill coeff., regulation of Dbp by PER-CRY
	hill_nampt_cb 1.91470474775	Hill coeff., regulation of Nampt by CLOCK-BMAL
	hill_nampt_pc 1.34080593157	Hill coeff., regulation of Nampt by PER-CRY
	hill_per_cb	8.52414053707	Hill coeff., regulation of Per by CLOCK-BMAL
	hill_per_pc	8.53897990872	Hill coeff., regulation of Per by PER-CRY
	hill_rev_cb	9.83701536835	Hill coeff., regulation of Rev-Erb by CLOCK-
			BMAL
	hill_rev_pc	3.31257899336	Hill coeff., regulation of Rev-Erb by PER-CRY
	hill_ror_cb	9.36456505724	Hill coeff., regulation of Ror by CLOCK-BMAL
	hill_ror_pc	1.84102439743	Hill coeff., regulation of Ror by PER-CRY
	Parameter Value (molecules per hour per mRNA) Description
	kp_bmal	0.628507384997	Bmal translation rate
	kp_cry	3.76912711677	Cry translation rate
	kp_nampt 58.9647983852	Nampt translation rate
	kp_per	13.2909782781	Per translation rate
	kp_rev	0.0513221194724	Rev-Erb translation rate
	kp_ror	0.0412765888526	Ror translation rate

Table S3H :

 S3H Protein stability modulation constants (4 parameters). (1) In our model, the activities of SIRT1 and AMPK have been defined to be dimensionless, hence the modulation constants are also dimensionless. (2) We included a possible modulation of PER by AMPK to take into account that CKI phosphorylated by AMPK may destabilize PER. (3) Given that AMPK activity peaks at about 2.5 in our model, the value of m_nampt_ampk indicates that NAMPT degradation rate is approximately divided by (1.63 × 2.5) ∼ 2.5 at peak AMPK activity.

	Parameter	Value	Description
	m_cry_ampk	0.07940386211	modulation of CRY stability by AMPK
	m_nampt_ampk 0.6352243791	modulation of NAMPT stability by AMPK
	m_per_ampk	0.005243953731 modulation of PER stability by AMPK
	m_per_sirt	0.005452322816 modulation of PER stability by SIRT

Table S3I :

 S3I NAD kinetics, Sirt1 and PGC1a activity (11 parameters)

	Parameter	Value	Unit		Description
	Vsirt	0.915854846427 N/A		Maximum SIRT1 activity
	Ksirt	0.75	nmol • l -1	Value of [NAD] at which SIRT1 activity is
					half of maximum
	d_nad	378.009130749	h -1		Rate of transformation of NAD into NAM
	Knad	0.321746039086 nmol • l -1	Value of [NAD] at which transformation
					into NAM is at half of maximum rate
	NAD_basal 0.9116166306	nmol • l -1	Value of [NAD] below which transforma-
					tion of NAD into NAM is inactivated
	Vnad	296.3933869	molecule		Maximum regeneration rate of NAD
			per	hour
			per NAMPT
			protein
	NAD_tot	4.166901679	nmol • l -1	Total concentration of NAD and NAM
	Knam	2.76496	nmol • l -1	Value of NAM at which NAD salvage rate
					is half of maximum
	Vpg	24.06372088	nmol -1 • h	Maximum activity of PGC1a
	Kpg1	0.046630145542 N/A		Michaelis-Menten constant for phosphory-
					lation of PGC1a by AMPK
	Kpg2	12.3526351747	N/A		Michaelis-Menten constant for deacetyla-
					tion of PGC1a by SIRT1

Table S3J :

 S3J Pulse parameters (12 parameters). The times tc1, tc2, and tc3 are given relative to ZT0.

	Parameter Value	Unit	Description
	tc1	4.38900149 h	Timing of the first AMPK pulse
	tc2	15.75	h	Timing of the second AMPK pulse
	tc3	18.875	h	Time of maximal nuclear PGC1a abundance
	Td1	2.25	h	Duration of the first AMPK pulse
	Td2	1.5	h	Duration of the first AMPK pulse
	Td3		h	Duration of the nuclear PGC1a presence
	Tr1	2.6	h	Rise time of the first AMPK pulse
	Tr2	1.8	h	Rise time of the second AMPK pulse
	Tr3	0.5	h	Rise time of nuclear PGC1a
	amp1	6.0	N/A	Amplitude of the first AMPK pulse
	amp2	0.9778008	N/A	Amplitude of the second AMPK pulse
	amp3	0.803062	nmol -1 • l Amplitude of the nuclear PGC1a abundance pulse

Table S3K :

 S3K Chronotherapy timings (4 parameters). The time tc4 is given relative to ZT0.

	Parameter Value	Description
	tc4	13.664	Timing of the agonist pulse
	Td4	2.83718	Duration of the agonist pulse
	Tr4	1.86794	Rise time of the agonist pulse
	amp4	0.465852 Amplitude of the agonist pulse

Table S3L :

 S3L Miscellaneous constants used to describe perturbations

	Name	WT SIRT1 KO LKB1 KO HFD fasting
	Csirt	1	0	1	1	1
	Campk 1	1	0.0375	0.05 0.05
	offs	0.02 0.02	0.02	0.02 2.6

AcknowledgmentsWe thank B. Pfeuty and Q. Thommen for discussions and K. Beuke for reading the manuscript. The authors acknowledge funding supports from INSERM, CNRS, Contrat de Plan Etat-Région (CPER) via the project "Photonics4Society", the Région Nord Pas-de-Calais, the European Regional Development Fund (ERDF), the "European Genomic Institute for Diabetes" (E.G.I.D., ANR-10-LABX-46), the "Centre Européen pour les Mathématiques, la Physique, et leurs Interactions" (CEMPI, ANR-11-LABX-0007) and European Commission. Research grants from the European Foundation for the Study of Diabetes (EFSD), The Fondation Francophone pour la Recherche sur le Diabète (FFRD), the European Commission (FP7) via consortium Eurhythdia are also acknowledged. BS is a member of the Institut Universitaire de France.

Author contributions

A. W., H. D., B. S. and M. L. designed the study and wrote the paper. A. W. and M. L. constructed the mathematical model, adjusted it to experimental data and carried out the numerical simulations.

Clock genes profiles in the fed and fasted states Fig. S4 shows a complete set of temporal profiles for gene expression and NAD+ levels in the fed state and in the fasting state, compared to the normal diet state. Typically, the fed state is characterized by a marked decrease of the amplitude of the oscillations, whereas the fasting state is associated with a strong phase shift and an increase of the amplitude.

Clock genes profiles in the fed and pharmacologically treated states Fig. S5 shows a complete set of temporal profiles for gene expression and NAD+ levels in the normal diet state, in the fed state, and in the fed state when agonist treated. In all figures, the optimal timing is shown in the middle column. To emphasize that a mistimed treatment is more harmful than none, two other timings, advanced and delayed by 8 hours are also shown for comparison.

As can be seen, almost all profiles can be restored simultaneously to their physiological shape with the optimal agonist delivery timing. The only two exceptions are Ror, whose amplitude becomes significantly higher than in normal diet conditions, but with the correct phase, and NAD+ which does not recover its small peak at ZT5 . Note that Ror overexpression would not necessarily be detrimental given the small dependance of Bmal1 transcription on ROR in our model. The physiological role of the small NAD+ peak has not yet been clarified, and therefore the consequences of its absence are difficult to predict. What is more certain is that the spectacular rescue of the big circadian NAD+ peak at ZT13 would likely be highly beneficial given the many physiological roles of the NAD+-dependent SIRT1 deacetylase, probably impaired in the fed state.

Quality of the expression profile rescue as a function of administration time

To determine how critical the choice of the administration time is, we computed for different administration times the residual error between the Bmal1, Rev-Erb, and Nampt profiles generated by the agonist-treated clock model in the fed condition and the original experimental profiles in the normal chow condition, all other characteristics of the agonist pulse being as given in Table S3K. How the residual error changes with administration time is shown in Fig. S6.

The optimal rescue was obtained when the agonist was administered at ZT 13.7, with a residual error of 9.8. Comparing this value to the residual error of 5.3 obtained between the profiles generated by the clock model in the normal chow condition and the experimental profiles in the same condition, we see that the physiological profiles are relatively well recovered, given that residual errors can range up to 840 for poorly chosen administration times.

The optimal time window is here relatively short, since the residual error remains below twice the minimum value between ZT13.3 and ZT14.1, corresponding to a duration of approximately 50 minutes. However, more detailed studies using realistic agonist pulse profiles are needed before definitive conclusions can be made. Rev-Erb agonist is delivered in the fed state at three different times of the day. Relates to Fig. 6.

Inventory of supplemental Information

Figure S6

Residual error between some time profiles generated by the agonist-treated clock model and the target experimental profiles. Relates to Fig. 6.

Table S1 List of variables of the mathematical model. Relates to Fig. 1.

Table S2 List of the differential equations and mathematical relations defining the mathematical model. Relates to Fig. 1.

Table S3 List of the kinetic constants obtained by parameter identification. For the convenience of the reader, Table S3 is split into Tables S3A-S3L.

Supplemental Experimental Procedures

Supplemental Results