Modeling clock-related metabolic syndrome due to conflicting light and food cues

Most organisms possess a light- and food- entrainable circadian clock system enabling their adaptation to daily environmental changes in sunlight and food availability. The mammalian circadian system is composed of multiple clocks throughout the body. These local clocks are entrained by nutrient, neural, endocrine and temperature cues and drive diverse physiological functions including metabolism. In particular, the clock of the pancreatic β cell rhythmically regulates the transcription of genes involved in glucose-stimulated insulin secretion. Perturbations of this fine-tuned oscillatory network increase the susceptibility to diseases. Besides chronic jet lag and shift work, common perturbations are ill-timed eating patterns which can lead to metabolic troubles (such as hypoinsulinemia). We have built a mathematical model describing the clock-dependent pancreatic regulation of glucose homeostasis in rodents. After calibrating the model using experimental data, we have investigated the effect of restricting food access to the normal rest phase. Our simulations show that the conflict between the light-dark cycle and the feeding-fasting cycle creates a differential phase shift in the expression of core clock genes (consistent with experimental observations). Our model further predicts that this induces a non-concomitance between nutrient cues and clock-controlled cues driving metabolic outputs which results in hypoinsulinemia, hyperglycemia as well as in a loss of food anticipation.


Model building
Our purpose was to obtain a mathematical model as simple as possible, but able to reproduce experimental phenotypes. To this end, we have built a model composed of 13 ordinary differential equations describing the time evolution of the concentration of clock components as well as of elements of the glucose-insulin module. The variables of the model are defined in Table S1. A few additional equations are used to form the external forcing functions representing food intake (meal) and SCN-driven cues. Concentration of exocytosis factors mRNA [EXO] Concentration of exocytosis factors protein Table S1: Definition of the variables of the mathematical model.
For the sake of simplicity, we have made several approximations. The different isoforms of each clock component were lumped together: Per 1,2,3 are merged into one variable, Per, and Cry1,2 are merged into Cry. Furthermore, we do not distinguish the cytoplasmic and nuclear forms of each protein. In addition, we neglect some core clock components. First of all, the activator CLOCK was not explicitely modeled because its expression is generally considered to be constitutive [10]. We thus suppose that it is the protein BMAL1 that is driving the dynamics of the positive arm of the clock: this enabled us to reduce the number of equations by considering that BMAL1 is playing the role of the complex CLOCK-BMAL1. We also do not incorporate the nuclear receptor ROR. Indeed, the latter seems to have a less important role than the nuclear receptor REV-ERB as it does not oscillate in every tissue [19]. Despite its simplicity, the model successfully reproduces experimental time series and phenotypes.
The transcription functions are detailed in Table S2 and reflect the molecular regulation of the core clock genes transcription: the activator complex CLOCK-BMAL1 (here represented by the protein BMAL1) first binds on the promotor of clock-controlled genes. This establishes a transcriptionally active state. Then, the repressors PER and CRY arrive which sets the beginning of a repressive state [7]. The transcription functions (Table S2) based on Hill functions were first successfully used in [14] and then in [18]. In section 2, we show how they can be derived using the quasi-steady state assumption. It should be noted that a basal gene transcription is also allowed which is reflected in the fold change parameters; g X (that represents the ratio of gene expression (e.g. transcription rate) in the presence and absence of transcription factors [1]). and f Rev,SCN in Eqs. S1 and S8). As the activation of clock gene expression by systemic cues generally occurs at sites that are distinct from those of CLOCK-BMAL1 regulation (see [16]), the function describing systemic depend-transcription is multipled by the one describing the clock-dependent transcription.

Non-shifting SCN cues
As the maximal neural activity occurs peaks during the light phase [5], we model this SCN output by a Fourier series with the appropriate phase: The parameter values are period = 24h, a 0 =0.2726, a 1 =-0.16, b 1 =0.14, a 2 =-0.01, b 2 =-0.01 and a 3 =0.25.

Nutrient cues
In our model, food intake is modelled by a periodic forcing function which qualitatively reproduces the number of pellets that animals consume per hour when food access is ad libitum. The feeding activity pattern corresponds to the alternation between starvation phases at the beginning of the light phase and feeding phases at the beginning of the dark phase [11]. To model this specific pattern (see Figure 1A in the main text), we choose a flexible forcing function, meal, which can reproduce a wide range of waveforms: The variables x and y used in the upper expressions are solutions of the ordinary differential equations: The parameter values are cfood =3.29, period = 24h, steep= 0.5, dur =3h and pha=17h for nighttime feeding and 5h for daytime feeding.
We use this forcing function because it is very flexible and this allows the comparison between different food-access paradigms such as ad libitum food access (ad lib) versus food access restricted to the night (tRF). Panda et al have shown that for tRF, the food intake function is more narrow and has a higher amplitude compared to the ad lib situation [4]. Our simulations show that the expression of some clock genes as well as insulin levels are boosted for a steeper food intake function and this is in agreement with Panda' results [4] ( Figure S6).

Equations for the clock-controlled exocytosis
It has been shown that BMAL1 and CLOCK bind to enhancers of genes coding for diverses factors involved in insulin exocytosis ( [12]). These elements are involved in the budding, transport, tethering and fusion of vesicles participating in insulin secretion. The phases of these cycling transcripts are variable: some of them peak 8h after Bmal1 mRNA while others peak at the transition between light and dark. The earlier elements include for example Vps13a,c,d which affect trafficking [12]. The latter elements include transcripts of factors such as Rhoa and Rhob that activate the protein kinase C (PKC) which contribute to initiate vesicle fusion [12]. Here, for the sake of simplicity, we only model the dynamics (mRNA and protein) of one generic element involved in exocytosis that we call Exo: we assume that its transcription is driven by the positive element BMAL1 and repressed by the negative arm of the clock (Eq. S15). We use a transcription function (f Exo,clock ) that is similar to those used for the core clock (see Table S2). The mRNA Exo is then translated into a protein, EXO (Eq. S16). While we did not constrain the parameters describing Exo mRNA dynamics, we did for those of EXO protein because we wanted its phase to be such that the EXO protein and glucose cues act simultanously to induce optimal insulin secretion under normal physiological conditions.
It should be mentioned that the β cell clock likely regulates insulin exocytosis in multiple ways. While in our model, the coupling of the clock to metabolism is mainly achieved via BMAL1 in agreement recent experimental results [12]), it is known that some exocytosisrelated genes are also regulated by REV-ERBα. Indeed, Vieira et al 2012 showed that the expression of some exocytosis-related factors is reduced after a knockdown of Rev-Erbα [17] but modeling these specific regulations in more detail is beyond the scope of this study. Instead, our view was to keep our model as simple as possible (one generic clock-controlled exocytosis factor), focusing on the important point that the maximum accumulation of our generic exocytosis factor has to coincide with the peak of glucose.

Equations for Insulin and Glucose
We model insulin dynamics as the balance between its secretion and its clearance (Eq. S17). Insulin secretion is modelled as a sigmoidal function of glucose levels which is consistent with experimental data (see references in [15]). We also use a sigmoidal function for the effect of the protein EXO on insulin secretion to account for the complexity of the exocytosis process. We consider that both the levels of glucose and exocytosis proteins (EXO) have to be high to promote a high insulin secretion: this corresponds to an AND gate and therefore, we multiply the functions describing the effect of glucose and EXO on insulin secretion by each other (see f Ins,glucose and f Ins,EXO in Table S2).
We model blood glucose dynamics as a balance between its production and its uptake (Eq. S18). Glucose production is modelled as the sum of a basal production rate and a meal-dependent production rate. We consider that the basal production rate incorporates inter alia the glucagon-dependent rise of blood glucose. Glucose uptake is written as a michaelian function of insulin level.

Derivation of transcription functions
In order to clarify how we obtain the transcription functions of Table S2, we here show as an example how we derive the transcription function used to describe the dynamics of Per mRNA. The transcription of Per is controlled by the clock as well as by systemic cues. If the two cues act independently, the transcription function can be written as T P er = f P er,clock · f P er,nutrient (S19) We thus have dP er dt = vmP · f P er,clock · f P er,nutrient − dmP · P er (S20) The dynamics of Per transcription can be described as follows: first, the binding of the activator (here BMAL1 ) on the promotor of the Per gene establishes a transcriptionally active state and leads to the synthesis of Per mRNA. Then, the repressor complex PER-CRY (PC ) binds to the E-box, which sets the beginning of the repressive state: the transcription of the gene Per is then inhibited. This set of reactions can be written as: where E0 corresponds to the unbound form of the E-box, E1 to the E-box with BMAL1 bound and E2 to the E-box with BMAL1 and PER-CRY bound. The parameters k 1 , k −1 , k 3 and k −3 are binding and unbinding rates while the parameters k 0 and k 2 are transcription rates. A conservation relation can also be written: If we translate the above reactions into ordinary differential equations by using mass action kinetics, we obtain: Taking into account the cooperative binding to the E-box leads to where the cooperativity is characterised by the Hill coefficients haP and hiP .
3 Parameter values

Parameter estimation
To our knowledge, the kinetic parameters have not been measured in pancreas. Here there are all estimated by parameter estimation methods. We used the Hooke and Jeeves pattern search optimisation method to minimise an objective function defined as the weighted sum of the square of the distance between in silico data points and experimental data from mice pancreas [11]. In these experiments, mice were kept in 12h light-dark conditions (DD), fed ad libitum and sacrified every 4 hours at ZT0 or ZT12 (with ZT0 corresponding the start of the light period and ZT12 to the start of the dark period). The data we used correspond to clock gene expression and metabolic time profiles. We first normalised the data (that is, divided them by the mean expression of each gene) and then interpolated them by Fourier series to facilitate the adjustment of the model to experimental data. The comparison between in silico and experimental time profiles is illustrated in Figure 3 in the main text. In [11], data originate from at least three experiments and are reported as mean +/-SEM. SEM is the standard error of mean which characterises the precision of the estimated mean of the population. Contrary to the standard deviation, it does not represent the variability within the data. To avoid misleading interpretation, SEM is not shown Figure 1 and 3 of the main text.
Parameter values are given in Tables S3-S5. Note that most transcription rates on the one hand and most degradation rates on the other hand lie in the same range except those of Rev-Erbα (vmR, dmR). When strongly reducing the value of the synthesis and the degradation rates of Rev-Erbα while keeping other parameters at their default value, the robustness of the oscillations during daytime feeding is lost. A possible explanation could be that as the profile of Rev-Erbα is far from a sinusoidal one, larger parameters values are needed to describe it. However, we cannot exclude that there are other sets of parameter values for which vmR and vmR are much lower and which would also reproduce the data.
We should also point out that given the simplified nature of the model, a single parameter typically encompass various processes (RNA maturation, post-translational modifications, transport).

Sensitivity to parameter values
One might wonder if the "gene-specific phase shift" property depends on a specific set of parameter values or if it is general for the described model. To answer this question, we investigate how the phase shift of Bmal1 and Per mRNAs is influenced by the variations of a few relevant parameters. Variations in these parameter values modify either the strength of the input of the zeitgebers (light or/and food) or the strength of the different core clock feedback loops. Figure S2 shows that when inverting the feeding schedule: • the higher the strength of food intake, the more Bmal1 mRNA shifts ( Figure S2A-B).
• the higher the strength of light cues (transmitted through SCN-driven), the less Bmal1 mRNA shifts ( Figure S2C-D).
• the higher the activation of Per by BMAL1, the less Per and Bmal1 mRNAs shift . This is probably because in this case, more information about the non-shifting light cues is transmitted to Per mRNA through Bmal1 mRNA (and Rev-Erb mRNA) ( Figure S2E-F).
• the higher the inhibition threshold of Rev-Erb by PER-CRY, the less Bmal1 shifts. This is probably because in this case, less information about the shifting food cues is transmitted to Rev-Erb mRNA and then to Bmal1 mRNA through Per mRNA ( Figure S2G-H).
• the higher the inhibition threshold of Bmal1 by REV-ERB, the more Bmal1 mRNA shifts. This is probably because in this case, less information about non-shifting light cues is transmitted to Bmal1 mRNA through Rev-Erb mRNA ( Figure S2I-J).
This analysis thus shows how light and food information are transmitted to the positive arm (here represented by Bmal1 ) and the negative arm (here represented by Per ) of the core clock. In particular, it indicates that the gene-specific phase shift property does not depend on a specific set of parameter values but can instead easily be obtained by modulating a few key parameters that control how feeding and light cues are transmitted to the different components of the core clock.    Right Y-axis: The phase shift difference between Bmal1 and Per mRNAs (which corresponds to the phase difference Bmal -Per in nighttime feeding (NF) minus the phase difference Bmal -Per in the shifted feeding regime ) is also represented as a function of the strength of SCN-driven cues (red curve). (B) As in A in the case of a 2h delay in food intake. (C) Insulin levels compared for nighttime feeding (blue curve) and a 2h advance in food intake (red curve). (D) As in C in the case of a 2h delay in food intake. The alternation of white and grey bars symbolises the LD cycles. Figure S4: Effect of a shift from nighttime feeding (NF, blue curves) to daytime feeding (DF). Comparison between the assumed physiological situation (red curves) and the case where peripheral clocks would uncouple from the SCN (green curves). The alternation of white and grey bars symbolises the LD cycles.  Figure S5: Resynchronisation time to the new feeding schedule (shift from nighttime feeding to daytime after two days) for different architectures, namely "No central clock" condition vs "No clock at all" condition). (A-B) Bmal1 expression and Insulin levels when the eating schedule is shifted from nighttime feeding to daytime (red curve) in the "No central clock" condition, to be compared with the nighttime feeding (blue curve). (C-D) Bmal1 expression and Insulin levels when the eating schedule is shifted from nighttime feeding to daytime (red curve) in the "No clock at all" condition, to be compared with the nighttime feeding (blue curve). The alternation of white and grey bars symbolises the LD cycles. (E) Time required to resynchronise to the new feeding schedule ( when shifting from nighttime feeding to daytime) compared for the "No central clock" condition (green bar) and the "no clock at all" condition (red bar). Figure S6: Comparison of clock gene expression and levels of metabolic actors for two different food-access paradigms: ad libitum food access (light blue curves) versus food access restricted to the night (dark blue curves). The alternation of white and grey bars symbolises the LD cycles.