The Synergistic Role of Light-Feeding Phase Relations on Entraining Robust Circadian Rhythms in the Periphery.

The feeding and fasting cycles are strong behavioral signals that entrain biological rhythms of the periphery. The feeding rhythms synchronize the activities of the metabolic organs, such as liver, synergistically with the light/dark cycle primarily entraining the suprachiasmatic nucleus. The likely phase misalignment between the feeding rhythms and the light/dark cycles appears to induce circadian disruptions leading to multiple physiological abnormalities motivating the need to investigate the mechanisms behind joint light-feeding circadian entrainment of peripheral tissues. To address this question, we propose a semimechanistic mathematical model describing the circadian dynamics of peripheral clock genes in human hepatocyte under the control of metabolic and light rhythmic signals. The model takes the synergistically acting light/dark cycles and feeding rhythms as inputs and incorporates the activity of sirtuin 1, a cellular energy sensor and a metabolic enzyme activated by nicotinamide adenine dinucleotide. The clock gene dynamics was simulated under various light-feeding phase relations and intensities, to explore the feeding entrainment mechanism as well as the convolution of light and feeding signals in the periphery. Our model predicts that the peripheral clock genes in hepatocyte can be completely entrained to the feeding rhythms, independent of the light/dark cycle. Furthermore, it predicts that light-feeding phase relationship is a critical factor in robust circadian oscillations.

Introduction

In mammalian organisms, the circadian machinery (ie, the biological time-keeping mechanisms which maintain 24-hour oscillations in cellular and physiological processes) is organized as a hierarchical, interconnected network of clocks.[1] The peripheral clocks are distributed across all cells and need to be synchronized within and across tissues and organs. The master, or central, clock in the suprachiasmatic nucleus (SCN) is mainly linked to light/dark cycle and synchronizes the network of peripheral clocks by controlling timing of feeding and activity cycles.[2,3] Rhythmic hormonal and metabolic signaling establishes the homeostatic phase relations among the various clocks, whereas these rhythmic signals play a major role in immune[4-6] and metabolic functions[7] conferring adaptive advantages by means of anticipatory control mechanisms.[8] These robust circadian patterns arise as a result of interactions between endocrine, immune, autonomic, and central nervous systems[9] regulated in a precise temporal manner by the central and peripheral circadian clocks[10,11] coordinating the diurnal variation of important physiological mediators.

Aside from the light/dark cycle, food intake is also a strong zeitgeber, ie, a cue given by the environment to the host to reset the internal body clock. In organs governing energy homeostasis, such as the liver, feeding rhythms are a more potent zeitgeber compared with the light/dark cycle, metabolically synchronizing peripheral clocks independent of the SCN.[12,13] Therefore, it is logical to suggest that there is an interconnection between the circadian clock machinery and metabolism, and indeed, disruption of one seems to affect the other. Several studies have shown evidence that disruption of circadian rhythms can facilitate metabolic syndrome; shift work and sleep deprivation dampen circadian rhythms and are linked to obesity,[14] whereas diabetes and cardiovascular diseases are linked with disruption of daily patterns in food intake in both children and adults.[15-17] To understand the cause and progression of such diseases, it is imperative to investigate the mechanisms of interplay between circadian clocks and metabolism.

The circadian clocks and metabolism seem to bidirectionally influence each other. In mammals, metabolic activities are under the regulation of daily feeding rhythms as well as the peripheral clock machinery. In turn, the feeding rhythms influence the circadian rhythms of key clock components via enzymatic reactions and transcriptional regulations.[18-20] Previous literature has found evidence for the interaction between circadian clocks and metabolism in relation to the well-established core circadian clock machinery. At a molecular level, clock mutant mice show a decrease in metabolic rate,[21] and liver-specific Bmal1 knockout mice exhibit an interruption in glucose homeostasis.[22] However, time-restricted feeding studies in mice show that metabolic cues influence circadian rhythmicity. Interestingly, restricted feeding schedule restored oscillations of some peripheral clock genes (PCGs) in clock-deficient mouse livers,[23] suggesting that behavioral adjustments can restore the disrupted circadian rhythms.

The above observations lead to the realization that environmental or behavioral cues can influence specific, low-level targets. To explore this hypothesis further in the context of circadian rhythms and metabolism, we first need to answer the question of how feeding rhythms entrain the PCGs in peripheral tissues, independently or synergistically with entraining the SCN. The focus on liver is driven by the major role it plays in maintaining the energy balance and the fact that its activity is regulated by nutrient intake. Once the entrainment mechanism is established, the synergistic or antagonistic role of competing light/dark cycle and feeding rhythms in the periphery should be explored to understand in depth the implications of light-feeding alignment and misalignment on the proper functioning of PCGs. In this work, we tackle these 2 questions by mathematically modeling the entrainment of the periphery by the light/dark cycle and the feeding/fasting cycle. Previous works have mathematically modeled the entrainment of peripheral clocks by either light/dark cycle[24,25] or feeding/fasting cycle,[26] but the convoluting effects of the 2 zeitgebers have not been studied mathematically. Our model consists of a central compartment which takes in the environmental and behavioral cues as inputs and a peripheral compartment which represents a human hepatocyte as it is more potent to the feeding rhythms rather than the SCN. The peripheral compartment will encompass the core clock machinery that oscillates autonomously and human sirtuin 1 (SIRT1)-mediated enzymatic reactions that influence the rhythms of PCGs.

Although tissue-specific rhythms are observed throughout the body, the main driver for the core circadian clock in most tissues, including liver and kidney, is the positive and negative feedback loops between 2 protein complexes, CLOCK/BMAL1 heterodimer and PER/CRY complex.[27-30] The PER/CRY complex inhibits the CLOCK/BMAL1-mediated transcription of PER and CRY proteins while stimulating the expression of BMAL1. Nutrient sensors such as adenosine monophosphate–activated protein kinase (AMPK), SIRT1, and poly(ADP-ribose) polymerase 1 (PARP1) exhibit circadian behavior and interact with the key molecules of the core circadian clocks, whereas also playing key roles in gluconeogenesis[31,32] and other metabolic activities. Among these nutrient sensors, we chose SIRT1 as a key candidate for bridging the circadian clocks and metabolism because it forms a complex network of regulations with components of the core clock machinery. Previous works on modeling circadian rhythms in the liver have recognized the importance of SIRT1 activity and incorporated it in their work.[26,33]

SIRT1 is a class III histone deacetylase, a homolog of Sir2 (silence information regulator 2) in yeast.[19] Its activity takes place in the nucleus, modulating lipid/protein/carbohydrate metabolism and enhancing mitochondrial activity,[34] in addition to modulating a variety of biological activities such as oncogenesis and aging.[35] SIRT1 has a clear role in metabolism as its enzymatic activity requires binding of nicotinamide adenine dinucleotide (NAD+) into its catalytic site along with the substrate, effectively acting as a sensor for the energy state of the cell. Modeling the activity of SIRT1 would require an accurate portrayal of the NAD+ level in the cell as it directly activates SIRT1. In addition to synthesizing NAD+ from amino acids, cells can also recover NAD+ from the NAD+ salvage pathway.[36] In the NAD+ salvage pathway, nicotinamide (NAM) is released from NAD+, NAM is converted to nicotinamide mononucleotide (NMN), and NMN is converted back to NAD+, completing the cycle.

The activity of SIRT1 has complex relationships with the NAD+ salvage pathway and the peripheral clocks. First, SIRT1 activity requires binding of NAD+, thus NAD+ has an activating effect on SIRT1. Second, SIRT1 activity is inhibited by NAM, the precursor of the rate-limiting step. Third, expression of nicotinamide mononucleotide adenylyltransferase (NMNAT), the rate-limiting enzyme of the NAD+ salvage cycle, is under the regulation of SIRT1. SIRT1 is also a cotranscription factor for the Nampt gene along with PCGs. In summary, SIRT1 functions as a nutrient sensor, being under the influence of the energy state of the cells, activated by NAD+. It is also under the effect of the circadian rhythmicity presented by NAD+, NAM, and nicotinamide phosphoribosyltransferase (NAMPT).

Based on the aforementioned observations, we propose a mathematical model that describes the interactions among SIRT1, NAD+ salvage cycle, and PCGs in human hepatocyte, to study the effects of light and feeding cycles on circadian dynamics in the periphery. We validate the model by comparing predictions under aligned and misaligned light/dark and feeding/fasting cycles of differing intensities, and comparing them to known behaviors under ad libitum and time-restricted feeding conditions. We further perturb the model by introducing delayed feedings of different feeding durations and explore the dynamics of PCGs under conflicting zeitgebers. Our model successfully entrains the PCGs of the hepatocyte to the feeding/fasting cycle away from the central compartment, agreeing with experimental studies on mammals. Our model predictions characterize the specific phase relations between light and feeding cycles that yield the strongest entrainment of PCGs. Furthermore, it is suggested that maintaining robust oscillatory patterns for SIRT1 is critical to drive robust oscillations of PCGs.

Materials and Methods

We propose a mathematical model to study the dynamics of the interactions of the circadian and metabolic components under the entrainment of 2 environmental cues: photoperiod (light/dark cycle) and nutrition (feeding cycle). The model was built upon our previous works,[24,37] introducing novel elements describing the dynamics of cellular energetics, PCGs, and cortisol level. The model is essentially composed of 2 compartments: a central compartment receiving the light/dark cycle and feeding information as inputs and processing them via modulations of the hypothalamic-pituitary-adrenal (HPA) axis and a peripheral compartment representing a human hepatocyte expressing the convoluted effects of signaling molecules emanating from the central compartment. Cortisol, entrained by the light/dark cycle, and NAD+, entrained by feeding, exert synergistic influences on the PCGs and metabolic enzymes in the periphery. A schematic representation displaying the network of signals transmitted from the 2 environmental cues to the downstream PCGs is shown in Figure 1.

Figure 1.

Schematic representation of the light- and feeding-entrained model. The HPA axis is entrained to the light/dark cycle, and NAD+ availability is entrained to food availability. Cortisol and SIRT1 in the periphery entrain the peripheral clock genes. HPA indicates hypothalamic-pituitary-adrenal; NAD+ nicotinamide adenine dinucleotide; SIRT, sirtuin 1.

Central compartment

The central compartment receives photic and metabolic cues as inputs and processes them. The signal transduction in this compartment governs the daily oscillations of the signaling molecules that modulate downstream activities in the periphery. The light/dark cycle governs cortisol release through the HPA axis, consisting of negative feedback between corticotropin-releasing hormone (CRH), adrenocorticotropin (ACTH), and cortisol.[24,37] The feeding rhythms transmit through a series of transit compartment signals eventually regulating NAD+ oscillations in the periphery. Cortisol and NAD+ will each modulate the intertwined network of the core clock machinery and metabolic enzymes in the periphery. We now consider these elements in greater detail.

In our model, the environmental cues are represented as step functions, similar to earlier studies.[24,25,38] Although the environmental and behavioral cues do not appear and disappear abruptly in nature, most time-restricted feeding experiments involving mammals provide illumination and food in an on/off manner.[39-41] Therefore, representing the light/dark cycle and feeding rhythms as step functions suffice both for recreating laboratory conditions in silico and for simplicity of the model. Equation (1) models the feeding rhythms and equation (4) models the light/dark cycle. Under the baseline case when feeding pattern is synchronized to the light schedule, a constant food supply is available between the hours of 6 am and 6 pm, and a constant supply of light source is also present between the hours of 6 am and 6 pm. Equation (1) was later modified for testing various feeding durations and delays relative to the light schedule. Equations (2) and (3) send the feeding signal through 2 transit compartments each with a delay of 3 hours (τ) to reflect the observation that NAD+ level peaks 5 to 6 hours after the beginning of the active period in rat liver,[42] suggesting a delay in feeding signal modifying the redox reactions between NAD+ and NADH, modifying the cellular NAD+/NADH ratio.[43]

The light signal regulates the release of cortisol from the central compartment to the periphery, driven by the self-sustained oscillations of CRH and ACTH in the HPA axis.[24,44] The entrainment of these components by the light/dark cycle is described by equations (4) to (7). These equations represent a Goodwin oscillator, modified to include Michaelis-Menten kinetics in the terms that represent synthesis and degradation of each component to avoid the use of unrealistically high Hill coefficients.[45] The production of CRH in the hypothalamus in equation (5) leads to the secretion of ACTH in the anterior lobe of the pituitary gland in equation (6). Equation (7) describes ACTH acting on the adrenal cortex, which produces the glucocorticoid hormones (mainly cortisol in humans). Then, cortisol negatively regulates CRH and ACTH via receptor-mediated activities, creating a negative feedback loop and maintaining sustained oscillations in the absence of environmental cues. Equations (8) to (11) describe the signal transduction pathway that forms the cortisol-receptor complex. These equations were derived from a corticosteroid pharmacodynamic model.[46] Receptor messenger RNA (mRNA) dynamics is described by a zeroth-order production and first-order degradation terms in equation (8). The indirect response model accounts for the experimentally observed downregulation of receptor mRNA upon methylprednisolone (MPL) treatment. The receptor protein dynamics is modeled in equation (9), mediated by the receptor mRNA, degradation, binding rate to the ligand cortisol, and transfer in and out of the nucleus. Equations (10) and (11) describe the receptor binding to cortisol and the resulting cortisol-receptor complex translocating to the nucleus to inhibit CRH and ACTH secretion. Finally, equation (12) describes the translocation of pro-inflammatory cytokines to the central compartment, which has an HPA-activating effect represented as stimulation of ACTH and cortisol using indirect response.[25]

Light is sensed by the hypothalamus in the form of photic signals communicated via the retinal ganglion cells in the eye.[47] The hypothalamic SCN integrates the photic input and regulates the physiological circadian rhythms of the cortisol and the periphery.[11,48] The SCN mediates secretion of light-induced arginine vasopressin (AVP).[49] Given that decreased secretion of AVP leads to increased cortisol levels,[50] stimulating both CRH and ACTH release,[51] we ultimately regulate the production of cortisol by light-induced degradation of CRH (equation (5)). We modified the HPA axis to reflect the signals received by the SCN from feeding. Ventromedial arcuate nucleus (vmARC) is often associated with satiety because metabolically active hormones such as ghrelin, leptin, glucocorticoid, insulin, and their receptors are expressed at high levels in it.[52-55] These hormones, along with glucose, modulate the electrical activity of vmARC.[56-58] The vmARC forms a complex with the subepidermal layer of the median eminence and have reciprocal connections with the SCN.[59] Therefore, we hypothesized that cellular energy state change caused by feeding in the periphery is transmitted to the SCN via electric signaling through the vmARC. Furthermore, studies showed that corticosterone rhythms change to show a food-anticipatory peak upon a change of feeding schedule.[60] Motivated by the above observations, we created a transit compartment variable, EntF, lumping the signaling activity from the periphery to the SCN via vmARC. Cortisol production responds to this signal in the HPA axis using Michaelis-Menten kinetics in equation (7). Parameter k in equation (7) represents the coupling strength of the feeding signal to the SCN.

Peripheral compartment

The peripheral compartment represents a human hepatocyte and encompasses a network of transcription, translation, and enzymatic reactions involving the PCGs and metabolic enzymes in the liver. This compartment considers the dynamics of signaling molecules entrained to the environmental and behavioral cues as inputs while also feeding information back to the central compartment to modulate the HPA axis activity. In our model, the oscillation of NAD+ is primarily governed by the feeding rhythms while also interacting with clock components, whereas the cortisol rhythms are under the influence of both light and feeding cycles. These molecules then activate a set of downstream reactions involving the peripheral clocks and metabolic enzymes, which form a complex, bidirectional network with each other.

The cortisol dynamics from the central compartment dominates the dynamics of cortisol in the peripheral compartment. Cortisol is present ubiquitously in all tissues of human body, as it performs a broad spectrum of physiological functions essential for life, including regulation of metabolism, immune response, and maintaining salt/water balance. However, the cellular response to cortisol varies in magnitude and specificity of action, both across different tissues and within the same tissue.[61] For example, kidney cells in the distal tubules have low cortisol sensitivity, although they do contain functional cortisol receptors and genes that contain glucocorticoid response elements because cortisol is metabolized by the enzyme 11β-hydroxysteroid dehydrogenase (11β-HSD).[62,63] In the liver, large amounts of glucocorticoid receptors (GRs) are present[64] and are involved in the regulation of metabolism including gluconeogenesis and lipogenesis. Cortisol controls the circadian clock machinery across different types of tissues. Dexamethasone, a glucocorticoid hormone analog, transiently changes the circadian gene expression phase in liver, kidney, and heart.[65] The diffusion of cortisol to the cytoplasm in hepatocyte is modeled as a transit compartment in equation (13). Once in the cytoplasm, cortisol can bind the mineralocorticoid receptor (MR) and the GR and activate them. The dynamics of the cortisol-receptor interactions has been previously modeled by our group.[37] In equations (14) and (17), we assume that cortisol activates either the MR or GR by phosphorylation. The production and degradation of the receptors are described with Michaelis-Menten kinetics. In these equations, MRT and GRT represent the total receptor concentrations, where MR and GR represent the activated form of the receptors that are already phosphorylated. Subsequent terms represent the binding of the receptor with cortisol and the recycled receptors. The dynamics of the receptor-ligand complex is described by equations (15) and (18), controlled by the binding rate and importation rate into the nucleus. Once the cortisol-receptor complex is translocated into the nucleus, the dynamics follow equations (16) and (19), controlled by the importation rate and recycle rate back to the cytoplasm.

Once the cortisol-MR complex and cortisol-GR complex translocate to the nucleus, they exert influence on the dynamics of pro-inflammatory cytokines and the PCGs. Equations describing the circadian rhythmicity of pro-inflammatory cytokine and its receptor were adapted from previous works by our group.[25,37] Although the role of pro-inflammatory cytokines was modeled for an immune subsystem previously, the influence of cytokines on cortisol dynamics was still implemented for modeling the periphery representing a hepatocyte, as shown in equation (7). The cytokines were adapted based on the justification that liver serves as an immunological organ.[66] In particular, human hepatocytes express pro-inflammatory cytokines such as interleukin 8, tumor necrosis factor α (TNF-α), and growth-related oncogene (GRO)-α, GRO-β, and GRO-γ[67] on bacterial infection, showing that hepatocytes both initiate and amplify inflammatory responses. Furthermore, lipopolysaccharide-induced endotoxic shock has been shown to exaggerate the abolition or alteration of circadian rhythms in the liver of animals subjected to chronic jet lag,[68] suggesting that mediation of PCGs and cortisol by pro-inflammatory cytokines should be included in the model for a liver subsystem.

It has been well established that cortisol functions as a critical driver for secretion of circadian cytokines, such as IFN-γ, interleukin 1, and TNF-α.[69] Furthermore, studies suggest a GR-mediated cytokine inhibition as experimental evidence shows that treatment with GR antagonist reduces expression of cytokines.[70-73] The inhibitory activity of the GR and cortisol complex is expressed in equation (20). In the same equation, the indirect inhibition of BMAL1 is modeled to reflect the diurnal rhythms of PCG-mediated pro-inflammatory cytokines. Participation of BMAL1 in the cytokine rhythms is supported by Bmal1 mRNA experiments which found that BMAL1-deficient mice myeloid cells showed exacerbated immune responses and increased secretion of cytokines under challenges with endotoxin or bacteria.[74-76] After translation, cytokines bind to cytokine receptor and form a complex represented as PR in Equation (24), which feeds back to the cytokine mRNA due to autocrine effects.[77] We simulated the effect of cytokine/cytokine receptor complex with an indirect response in equation (20). The translation of cytokine mRNA is described by equation (21). The translated cytokine protein stimulates the HPA axis as previously mentioned and indirectly induces the translation of Per/Cry mRNA, which will be described later. The cortisol-mediated upregulation of cytokine receptors was modeled in equation (22) as indirect stimulation of the receptor mRNA. Equation (23) describes the translation of cytokine receptor, and equation (24) represents the dynamics of the cytokine/cytokine receptor complex.

Although the light/dark cycle entrains the PCGs via the cortisol-receptor complex, the feeding rhythms entrain the PCGs via NAD+ and SIRT1. The mathematical representation of feeding entrainment on the PCGs involves the rhythmic response of NAD+ and SIRT1 to feeding and fasting states, which react with multiple components in the circadian machinery. The dynamics of NAD+ in the periphery is mainly driven by 2 reactions: (1) reduction to NADH and (2) the NAD+ salvage pathway. It is suggested that while feeding, continuous glycolytic throughput would result in the accumulation of NADH and therefore limit the NAD+ concentration.[78] In contrast, fasting would decrease the glycolytic NAD+ demand, resulting in higher NADH oxidation and higher NAD+ concentration.[78] Therefore, in equation (25), which describes the changing concentration of NAD+ in the periphery due to food availability and the NAD+ salvage cycle, feeding signal (through the transit compartments in equations (2) and (3)) is the driving force for the inhibition term. The transit compartments for feeding were added to the model based on the experimental observations that NAD+ level peaks 5 to 6 hours after the beginning of the active period in rat liver,[42] suggesting a delay in feeding signal modifying the redox reactions between NAD+ and NADH. The regeneration of NAD+ from NADH by oxidation reaction is modeled using Michaelis-Menten kinetics. In equation (25), nad represents the combined concentration of NAD+ and NADH together; therefore, the quantity (nad-NAD) represents the concentration of NADH available to be oxidized, and the rate of NAD+ regeneration is dependent on it. Aside from NADH oxidation, the NAD+ salvage pathway is another way cells generate NAD+,[79] and the reactions are modeled using Michaelis-Menten kinetics in equations (25) to (27). In the salvage pathway, NAM is released from NAD+ in adenosine diphosphate-ribose transfer reactions, represented as a degradation term in equation (25) describing the change in NAD+ concentration, and represented as a synthesis term in equation (26) describing the change in NAM concentration. Nicotinamide is converted to NMN by rate-limiting enzyme NAMPT catalyzation. Because the conversion of NAM to NMN is dependent on the availability of NAMPT, the synthesis term for NMN is multiplied by the NAMPT concentration in equation (27). Finally, NMN is converted back to NAD+ by an enzyme called NMNAT, completing the cycle of the NAD+ salvage pathway.

NAD+ acts the representative agent that communicates the energy state of the periphery to the SCN in our model. As mentioned previously, the electric signaling activity of vmARC is regulated by metabolically active hormones such as ghrelin, leptin, glucocorticoid, insulin, and their receptors, along with glucose.[52-58] Furthermore, tracer experiments using fluorophore-conjugated cholera toxin B show that the SCN is one of the sites that is involved in reciprocal communication with the vmARC.[59] Therefore, we hypothesize that information about food availability in the periphery is transmitted to the SCN via electric signaling through the vmARC and lump the hormonal and electric signaling activity into 1 transit compartment in equation (28). This transit compartment, termed EntF, stimulates cortisol secretion in the HPA axis in the central compartment, as shown in equation (7).

NAD+ is also the activating agent of SIRT1. It is well established that NAD+ activates SIRT1 by directly binding and altering the conformation of the catalytic site to allow binding of substrates, as the crystal structure of SIRT1 in the active conformation has shown.[80,81] Furthermore, SIRT1 is an attractive candidate to study in the context of metabolic diseases both experimentally and in silico because multiple inhibitors and activators, such as indoles/indole derivatives[82] and resveratrol,[83] were studied as possible therapeutics for cancer, obesity, and aging-related heart diseases.[84] The activity of NAD+ and SIRT1 has been adapted for modeling circadian clocks in the liver by other groups[26,33] because they play a major role in communicating the feeding rhythm to the clock genes. The stimulating effect of NAD+ on SIRT1 is modeled in equation (29) using Michaelis-Menten kinetics. The constant sirtT represents the total concentration of SIRT1 protein, both in the inactive conformation and the NAD+-bound active conformation, and SIRT1 represents the activated SIRT1 enzyme. The production rate of active SIRT1 is dependent on the concentration of inactive SIRT1 protein availability as well as the activator (NAD+) concentration. We model the dynamics of active SIRT1 in this manner because experimental evidence suggests that the total SIRT1 protein concentration stays relatively constant, whereas the enzymatic activity level oscillates with a circadian rhythm.[19]

So far, the environmental cues entrain the periphery through 2 different channels. The light/dark cycle entrains the cortisol rhythms in the central compartment through the HPA axis, which drives the dynamics of cortisol-receptor complex in the periphery. The feeding cycle entrains the dynamics of NAD+ which communicates the cellular energy state to the SCN and cortisol in the HPA axis and also the SIRT1 protein. The cortisol-receptor complex and SIRT1 are the 2 main entrainers of PCGs that will produce convoluted downstream effects, and their entrainment dynamics will be described.

The entrainment dynamics of the PCGs by cortisol is based on our previous work.[24] A network of transcriptional and translational feedback loops are incorporated into equations (30) to (36) with important modifications to model the circadian rhythmicity of the PCGs, ultimately resulting in autonomous oscillations of the PCGs. Although these equations were originally developed for an immune subsystem such as macrophages or neutrophils, the basic structure behind the core clock machinery involving the PER/CRY and CLOCK/BMAL1 complexes is similar across different tissues in the body[27-30]; however, the rhythmic patterns of clock components may exhibit varying properties depending on the tissue type.[85] Therefore, these equations are successful in describing the PCG behavior in the liver. The negative and positive feedback loops consisting of the 2 protein complexes have been implemented in mathematical modeling of a hepatocyte in another work,[33] further justifying the utilization of this mechanism in modeling a hepatocyte. In equation (30), the entrainment of Per and Cry mRNAs by cortisol is taken into account based on the consideration that cortisol and GR complex binds to the promoter region of the Per1 and Per2 genes.[86] The transcription of the Per and Cry genes is stimulated when CLOCK/BMAL1 complex binds to an Ebox enhancer. In addition, translated PER/CRY protein translocates to the nucleus and inhibits its own translational activity, forming a negative feedback loop. The exponent p to the PER/CRY protein in the nucleus in equation (30) is a Hill function coefficient and is used to describe the switch-like behavior of the translational activities. The role of pro-inflammatory cytokines is also incorporated into this equation via indirect response, based on studies that show that Per1, Cry1, and Cry2 expressions are induced upon treatment with pro-inflammatory cytokines such as IL-6 or TNF-α.[87-89] The dynamics of PER/CRY protein is represented by equation (31) and is composed of the complex formation rate, degradation rate, and nuclear import and export rates. The PER/CRY protein dynamics inside the nucleus is governed by equation (32), which is also described by the nuclear import/export rate as well as the degradation rate. When inside the nucleus, PER/CRY protein also suppresses the activation of REV-ERBα, which is promoted by CLOCK/BMAL1 complex.[90] In turn, REV-ERBα negatively regulates transcription of Bmal1. Therefore, PER/CRY indirectly induces Bmal1 transcription, and this relationship was modeled as a Hill equation with a coefficient r in equation (33). The dynamics of the BMAL1 protein is dependent on translation, degradation, and nuclear import/export rates, as shown in equation (34). Inside the nucleus, BMAL1(N) binds to the CLOCK protein, which is assumed to be in an access of BMAL1(N) in this model. The dynamics of BMAL1(N) and the CLOCK/BMAL1 complex is shown in equations (35) and (36).

The feeding rhythms entrain the PCGs via NAD+ and SIRT1, but the interaction between PCGs and NAD+-activated SIRT1 is bidirectional. SIRT1 exerts influence on the rhythms of some clock components and is also under the regulation of them at the same time. The activated SIRT1 protein binds to the CLOCK/BMAL1 complex, one of the key components that drive the core circadian machinery in the periphery. Once it forms the CLOCK/BMAL1/SIRT1 complex, it promotes the expression of NAMPT.[91] As previously mentioned, NAMPT controls the rate of regeneration of NAD+, the activator of SIRT1, from the NAD+ salvage cycle. Therefore, SIRT1 self-regulates its activation through interacting with core peripheral clock components. The dynamics of CLOCK/BMAL1/SIRT1 complex is modeled in equation (37). In this equation, k8 represents the association rate of CLOCK/BMAL1 complex and SIRT1, k8 represents the dissociation rate of CLOCK/BMAL1/SIRT1 complex back to CLOCK/BMAL1 and SIRT1, and k9 represents the degradation rate of this complex. The regeneration of CLOCK/BMAL1 from CLOCK/BMAL1/SIRT1 complex is present in equation (36). The activity of NAMPT is modeled in equation (38) using first-order production (k10) and degradation (k10) terms. By binding to the CLOCK/BMAL1 complex, SIRT1 also takes away the availability of CLOCK/BMAL1 complex that promotes the expression of PER/CRY. Therefore, accumulation of activated SIRT1 enzyme will inhibit the translation of Per/Cry mRNA. SIRT1 affects the rhythm of PER/CRY with its deacetylating activity as well. In equation (32), SIRT1 indirectly stimulates the nucPER/CRY degradation. This modification is introduced because SIRT1 is known to facilitate degradation of PER2 protein by deacetylation of PER2.[18,19] Supporting experimental evidence also shows abnormally high PER2 accumulation in SIRT1-deficient mouse embryonic fibroblasts.[92]

In summary, our model uses SIRT1 as the key molecule that connects the cellular energy state resulting from feeding state to the PCGs, whose rhythms are now entrained by signals from both light/dark and feeding/fasting cycles. The dynamics of NAD+ (a SIRT1 activator) is modeled to reflect the change in redox reaction from feeding as well as the NAD+ salvage cycle that takes input from the circadian dynamics from CLOCK/BMAL1 heterodimer. The interaction between NAD+ and cortisol is also captured to reflect adjustments to the cortisol dynamics upon change in feeding schedule, introducing a modification to the HPA axis:

Parameter Estimation and Sensitivity Analysis

Our model aims to study the convoluted effects of light/dark cycle and feeding rhythms on the PCGs in human hepatocyte by mathematical modeling and naturally involves numerous equations with parameters representing signal transduction and reaction rates. However, quantified level or activity for the model components is not always available throughout the day, and constraining the parameters becomes an extremely difficult task. Often, extrapolations are made from diverse animal studies across different species, adding more uncertainties to parameter estimation. This issue is common in systems biology and the field has recognized that qualitatively capturing biological phenomena while making accurate and relevant predictions is more valuable than constraining every parameter in the model.[93] To that end, where new components or reactions were introduced in our model, parameters were estimated to capture the qualitative characteristics of key components including phase relations and amplitude changes at varying environmental conditions. In this section, we first explain what specific qualitative features of physiology were considered in our estimated parameters, listed in Table 1, and describe how the sensitivity analysis was performed.

Table 1.

Parameter values and descriptions.

No.	Parameter	Value	Units	Description/reference
1	kp 1	0.7965	µM h−1	Rate constant of CRH production[24]
2	Kp 1	1.0577	µM	Dissociation constant for CRH production[24]
3	Vd 1	0.5084	µM h−1	Rate of CRH enzymatic degradation[24]
4	Kd 1	1.9627	µM	Michaelis constant of CRH enzymatic degradation[24]
5	kfp	0.15	µM−1	Efficiency of P on ACTH and F stimulation/estimated[25]
6	kp 2	0.6857	µM h−1	Rate of ACTH production[24]
7	Kp 2	1.0577	µM	Dissociation constant for ACTH production[24]
8	Vd 2	0.5129	µM h−1	Rate of ACTH enzymatic degradation[24]
9	Kd 2	0.3069	µM	Michaelis constant of ACTH enzymatic degradation[24]
10	kp 3	1.0302	µM h−1	Rate of F central production/estimated[25]
11	k n	1.2		Coupling constant of cortisol to cellular energy state
12	Vd 3	0.3618	µM h−1	Rate of F central enzymatic degradation[24]
13	Kd 3	0.4695	µM	Michaelis constant of F central enzymatic degradation[24]
14	ksynRm	2.9	fmol g−1 h−1	Synthesis rate or glucocorticoid receptor mRNA[46]
15	IC 50Rm	26.2	nmol L−1 mg protein−1	Concentration of FR(N) at which mRNA, R synthesis drops to its half[46]
	R 0	540.7	nmol L−1 mg protein−1	Baseline value of free cytosolic glucocorticoid receptor[46]
	Rm 0	25.8	fmol g−1	Baseline value of glucocorticoid receptor mRNA[46]
16	kdgrRm	ksynRm/Rm0		Degradation rate of glucocorticoid receptor mRNA[46]
17	ksynR	(R0/Rm0)*kdgrR		Synthesis rate of free cytosolic receptor[46]
18	rf	0.49		Fraction of cortisol recycled[46]
19	kre	0.57	h−1	Rate of receptor recycling from nucleus to cytoplasm[46]
20	kon	0.00329	L nmol−1 h−1	Second-order rate constant of glucocorticoid receptor binding[46]
21	kdgrR	0.0572	h−1	Degradation rate of cytosolic glucocorticoid receptor[46]
22	kt	0.63	h−1	Rate of receptor translocation to the nucleus[46]
23	τ	0.25	h	Delay in cortisol production following ACTH stimulation
24	kMR	0.34	nM h−1	Base transcription rate of MR[37]
25	kF,MR	1.1011	1	Maximum extent of Fperiphery-mediated activation of MR[37]
26	KF,MR	0.5	nM	Michaelis constant for Fperiphery-mediated activation of MR[37]
27	MRT	1.45	nM	Total MR concentration[37]
28	KMR	0.21	nM	Michaelis constant for MR production[37]
29	kMR,deg	0.70	nM.h−1	Degradation rate for MR[37]
30	KMR,deg	1.65	nM	Michaelis constant for degradation of MR[37]
31	kb,MR	0.00329	nM−1 h−1	Degradation rate for cortisol/mineralocorticoid receptor binding[46]
32	kon,MR	1	L nmol−1 h−1	Second order rate constant of mineralocorticoid and receptor binding[37]
33	kt,MR	1	h−1	Rate of mineralocorticoid receptor translocation to the nucleus[46]
34	kre,MR	1	h−1	Rate of mineralocorticoid receptor recycling from nucleus to cytoplasm[37]
35	kGR	1.18	nM h−1	Base transcription rate of GR[37]
36	kF,GR	15	1	Maximum extent of Fperiphery-mediated activation of GR[37]
37	KF,GR	30	nM	Michaelis constant for Fperiphery-mediated activation of GR[37]
38	GRT	1.81	nM	Total GR concentration[37]
39	KGR	0.74	nM	Michaelis constant for GR production[37]
40	kGR,deg	1.52	nM h−1	Degradation rate for GR[37]
41	KGR,deg	1.05	nM	Michaelis constant for degradation of GR[37]
42	kb,GR	0.00329	nM−1 h−1	Degradation rate for cortisol/glucocorticoid receptor binding[46]
43	kon,GR	1	L nmol−1 h−1	Second-order rate constant of mineralocorticoid and receptor binding[46]
44	kt,GR	1	h−1	Rate of mineralocorticoid receptor translocation to the nucleus[46]
45	kre,GR	1	h−1	Rate of mineralocorticoid receptor recycling from nucleus to cytoplasm[46]
46	kmRNARpin	0.61	µM h−1	Base transcription rate of mRNARP[37]
47	kfr 2	0.8	1	Maximum extent of FMR(N)-mediated transcription of mRNARP[37]
48	Kfr 2	0.5	µM	Michaelis constant for FMR(N)-mediated transcription of mRNARP[37]
49	kpc	0.3	1	Maximum extent of BMAL1-mediated suppression of mRNAP and mRNATLR4 estimated
50	Kpc	25	µM	Michaelis constant for BMAL1-mediated suppression of mRNAP and mRNATLR4 estimated
51	kmRNARpout	0.19	h−1	Degradation rate of mRNARP[37]
52	kinP	0.29	h−1	Translation rate of P[37]
53	koutP	1.06	h−1	Degradation rate of P[37]
54	kmRNARpin	0.61	µM h−1	Base transcription rate of mRNARP[37]
55	kfr 2	0.8	1	Maximum extent of FMR(N)-mediated transcription of mRNARP[37]
56	kmRNARpout	0.19	h−1	Degradation rate of mRNARP[37]
57	kinRp	1.11	h−1	Translation rate of RP[37]
58	kd	0.14	µM−1 h−1	PR binding rate constant[37]
59	koutRp	0.26	h−1	Degradation rate of RP[37]
60	koutPRp	1.3	h−1	Dissociation rate of PR[37]
61	v 1b	9	nM h−1	Maximal rate of Per/Cry transcription[94]
62	k 1b	1	nM	Michaelis constant of Per/Cry transcription[94]
63	k 1i	0.56	nM	Inhibition constant of Per/Cry transcription[94]
64	C	0.01	nM	Concentration of constitutive activator[94]
65	P	8		Hill coefficient of inhibition of Per/Cry transcription[94]
66	kf	1.2	nM−1	Efficiency of P on transcription of Per/Cry/estimated[25]
67	kc	0.009	nM h−1	Coupling strength/estimated
68	k 1d	0.12	h−1	Degradation rate of Per/Cry mRNA[94]
69	k 2b	0.3	nM−1 h−1	Complex formation rate of Per/Cry mRNA[94]
70	Q	2		No. of PER/CRY complex forming subunits[94]
71	k 2d	0.05	h−1	Degradation rate of cytoplasmatic PER/CRY[94]
72	k 2t	0.24	h−1	Nuclear import rate of the PER/CRY complex
73	k 3t	0.02	h−1	Nuclear export rate of PER/CRY complex
74	k 3d	0.02	h−1	Degradation rate of the nuclear PER/CRY complex
75	v 4b	3.6	nM h−1	Maximal rate of Bmal1 transcription[94]
76	k 4b	2.16	nM	Michaelis constant of Bmal1 transcription[94]
77	R	3		Hill coefficient of activation of Bmal1 transcription[94]
78	k 4d	0.75	h−1	Degradation rate of Bmal1 mRNA[94]
79	k 5b	0.24	h−1	Translation rate of BMAL1[94]
80	k 5d	0.06	h−1	Degradation rate of cytoplasmatic BMAL1[94]
81	k 5t	0.45	h−1	Nuclear import rate of BMAL1[94]
82	k 6t	0.06	h−1	Nuclear export rate of BMAL1[94]
83	k 6d	0.12	h−1	Degradation rate of nuclear BMAL1[94]
84	k 6a	1	h−1	Activation rate of nuclear CLOCK/BMAL1
85	k 7a	0.1	h−1	Deactivation rate of CLOCK/BMAL1
86	k 7d	0.5	h−1	Degradation rate of CLOCK/BMAL1
87	k m 8a	10	h−1	Association rate of CLOCK/BMAL1/SIRT1
88	k m 8d	20	h−1	Dissociation rate of CLOCK/BMAL1/SIRT1
89	τf	3	h	Delay between feeding and NAD+ reduction to NADH
90	k m 1	5	1	Maximum extent of NADH converting to NAD+
91	k m 1	2	µM	Michaelis constant for NADH converting to NAD+
92	k m 2	40	1	Maximum extent of NMN converting to NAD+
93	k m 2	1	µM	Michaelis constant for NMN converting to NAD+
94	k m 3	5	1	Maximum extent of NAD+ converting to NADH upon feeding
95	k m 3	1	µM	Michaelis constant NAD+ converting to NADH upon feeding
96	k m 4	20	1	Maximum extent of NAD+ converting to NAM
97	k m 4	20	µM	Michaelis constant for NAD+ converting to NAM
98	k m 5	40	1	Maximum extent of NAM converting to NMN, mediated by NAMPT
99	k m 5	5	µM	Michaelis constant for of NAM converting to NMN, mediated by NAMPT
100	k m 11	5	1	Maximum extent of entrainment of cortisol by feeding
101	k m 11	2	µM	Michaelis constant for entrainment of cortisol by feeding
102	k m 12	1	µM h−1	Rate of cortisol-entraining signal degradation
103	k m 6	5	1	Maximum extent of SIRT1 activation mediated by NAD
104	k m 6	1	µM	Michaelis constant of SIRT1 activation mediated by NAD
105	sirtT	5	µM	Sum of active and inactive SIRT1
106	k m 7	2	1	Maximum extent of SIRT1 degradation
107	k m 7	1	µM	Michaelis constant for SIRT1 degradation
108	k m 8a	10	µM−2 h−1	Rate of SIRT1 and CLOCK/BMAL1 complex association
109	k m 8d	20	µM−1 h−1	Rate of SIRT1 and CLOCK/BMAL1 complex dissociation
110	k m 9d	0.1	µM−1 h−1	Rate of CLOCK/BMAL1/SIRT1 complex degradation
111	k m 10a	2	µM h−1	Rate of NAMPT production mediated by CLOCK/BMAL1/SIRT1 complex
112	k m 10d	2	h−1	Rate of NAMPT degradation
113	Nad	1	µM	Sum of NAD+ and NADH

Abbreviations: ACTH, adrenocorticotropin; CRH, corticotropin-releasing hormone; GR, glucocorticoid receptor; mRNA, messenger RNA; NAD+, nicotinamide adenine dinucleotide; NAM, nicotinamide; NAMPT, nicotinamide phosphoribosyltransferase; NMN, nicotinamide mononucleotide; SIRT1, sirtuin 1.

Bold parameters are estimated.

Of particular mention are certain parameters: k in equation (7), the coupling strength of the feeding signal to the SCN, was estimated to capture the behavior of cortisol under conditions where feeding is synchronized to the active (light) period, and feeding is antisynchronized to the active period. When food is available during the active period, there is ample evidence that cortisol level peaks in the morning time for humans.[95] However, when food availability is not synchronized to the light period, the cortisol level peaks slightly before food availability in the anticipation of incoming nutrients,[60] and the oscillation amplitude is expected to decrease.[96] For the equations involving the NAD+ salvage cycle (equations (25)-(27)), the key feature captured in parameter estimation is that NAD+ levels are high during fasting and become lower during feeding.[97] In addition, there is a delay of 5 to 6 hours between feeding start time and NAD+ level decline,[97] and NAMPT tends to peak shortly before the time when NAD+ level peaks.[98] The estimated parameters ensure that these qualities are captured regardless of the presence of conflict between the light/dark cycle and feeding cycle. Estimation of parameters for equation (29) describing the dynamics of SIRT1 captures that SIRT1 level peaks ~3 hours post the first food availability based on mouse liver data.[19] Estimation of parameters describing the interaction between core clock genes and metabolic components (equations (30)-(38)) was performed in 2 steps. First, parameters were optimized for correct phase relations of components relative to light under synchronized schedule. For example, it is well established that PER and CRY proteins peak in the morning time or early active phase.[95,99,100] Once these conditions were met, the parameters were estimated once again at a state in which feeding was antisynchronized to light. The new estimation would yield complete inversion of the metabolic components and PCGs,[60] and decrease in amplitude for PCGs,[96] agreeing with experimental observations. For all estimations, PER/CRY protein was used as the representative sample for PCGs because PER/CRY rhythms are under close influence of both SIRT1 and nuclear cortisol-receptor complex. The finalized values of the estimated parameters and their descriptions are shown in Table 1. The system contains 113 parameters, of which only 35 were adjusted in this study, whereas the remaining were fixed based on prior information.

A sensitivity analysis was performed on the model to gain insight into the PCG dynamics predicted by our model and the influence of the various parameters. The analysis was performed under 2 different entrainment conditions: first with light and feeding signals aligned and next with light and feeding signals antisynchronized. We computed a relative sensitivity coefficient[38] for every parameter in Table 1 using , where p is the tested parameter and y is the measured response. Because we are most interested in studying the dynamics of PCGs where both light and feeding cues exert their influence, we used the profile of PER/CRY to calculate the sensitivity coefficients. PER/CRY level is under the influence of transcriptional activation by the cortisol-receptor complex and degradation by SIRT1, exhibiting convoluted entrainment effects from both light/dark cycle and feeding rhythms. Each parameter was varied by 1%, and relative local sensitivity coefficients for each parameter based on the response variable’s amplitude and phase angles were computed.

Results

Our in silico studies aim to explore the characteristics of the dynamics of the peripheral clock machinery under the control of 2 independent, but strongly interacting, entrainers: light and feeding. The schematics of the model including the HPA axis entrained by light and the periphery representing a hepatocyte entrained by cortisol and feeding are shown in Figure 1. The redox reaction between NAD+ and NADH receives the signal from feeding/fasting cycle, and the resulting NAD+ dynamics stimulates the activation of SIRT1, which serves as a bridge between light and feeding entrainment. SIRT1 then interacts with the PCGs in a bidirectional manner, directly and indirectly regulating the rhythms of PER/CRY via binding and deacetylation reactions while being under the influence of CLOCK/BMAL1 heterocomplex.

In Figure 2, the time profiles of cortisol (Figure 2B), PER/CRY protein (Figure 2C), and SIRT1 (Figure 2D) are shown over a period of 24 hours, under synchronized and antisynchronized light and feeding schedules. Figure 2A shows the light and feeding signals under both conditions. Light signal was on from 6 am to 6 pm for both conditions. The feeding signal was also on from 6 am to 6 pm for the synchronized schedule, whereas the signal was on from 6 pm to 6 am, opposite from the light signal, for the antisynchronized schedule. When feeding is synchronized to light, cortisol and PER/CRY protein peak in the morning time or during the early active period. In contrast, when feeding is antisynchronized to light, PER/CRY phase is completely inverted, whereas cortisol peak shifts to the beginning of the inactive phase. In addition, the amplitudes for cortisol and PER/CRY oscillations are smaller for the antisynchronized schedule, compared with the synchronized schedule. These features are in qualitative agreement with restricted feeding studies performed on mice.[101,102] The biphasic behavior of NAD+ observed in rodent experiments[40,97] was not reproduced by our model; instead, the higher of the 2 peaks which falls during the feeding phase was reproduced. This led to a realistic profile of activated SIRT1 in Figure 2D, peaking 4 to 5 hours after the beginning of the active phase and agreeing with previously published experimental work.[103] The amplitude of SIRT1 does not undergo a change between the synchronized and antisynchronized schedules, responding only in peak time.

Figure 2.

The time profiles of key components throughout the day upon synchronized (black dotted line) and antisynchronized (blue line) light and feeding schedules. Light signal was at 1 from 6 am to 6 pm and at 0 for the rest of the day. (A) Feeding signals synchronized and antisynchronized to light, (B) corresponding cortisol, (C) PER/CRY protein, and (D) SIRT1 profiles. SIRT indicates sirtuin 1.

We then explored how abolishing light rhythms by maintaining light of constant intensity during the light period will affect the PCG dynamics. Feeding was available from 6 am to 6 pm, whereas light was on for the entire 24-hour period (Figure 3A). The resulting time profiles for cortisol, PER/CRY, and SIRT1 were compared against a synchronized schedule for the light and feeding (Figure 3). The results show that under constant light, cortisol (Figure 3B) and PER/CRY (Figure 3C) oscillate albeit with reduced amplitudes compared with light/dark cycles synchronized to feeding. Applying constant light of brighter intensity resulted in greater reduction in amplitude for cortisol and PER/CRY. This observation is in qualitative agreement with experimental observations in mice.[104] The phase angles or cortisol and PER/CRY were not affected by changing intensities of light. However, SIRT1 oscillations remained identical when applying different light schedules and intensities and exhibited no difference in amplitude or phase.

Figure 3.

Time profiles of key components under constant light schedule at different intensities are compared. (A) Feeding signal was at 1 from 6 am to 6 pm and at 0 for the rest of the day. Constant light signals at different intensities, along with a control, in which feeding and light are synchronized. (B) Cortisol, (C) PER/CRY protein, and (D) SIRT1 profiles under the different light intensities. SIRT1 indicates sirtuin 1.

Imposing a constant light signal to the system not only results in a diverged PCG dynamics from the synchronized light/dark cycle but also changes the system response when feeding signals are inverted. The transition time for feeding inversion under constant light conditions was calculated and is shown in Figure 4 along with the transition time under light/dark cycle. The transition time is defined as the number of days to first reach the steady-state phase relations after feeding inversion. Initially, the feeding signal was on from 6 am to 6 pm and then switched to be on from 6 pm to 6 am at midnight. The phase angle at new steady state was first computed, and then each day’s phase angle following feeding inversion was compared with the computed steady-state phase angle to determine when the PER/CRY rhythm first reached the steady-state phase angle. As shown in Figure 4, the phase inversion was achieved more quickly with constant light conditions at both bright (L/L) and dark (D/D) light intensities, compared with the light/dark cycle (L/D).

Figure 4.

Number of days taken to first reach the steady-state phase angle upon feeding inversion at different light schedules. L/D is a 12-hour light, 12-hour dark cycle. D/D is a 24-hour dim light schedule with an intensity of 0.1. L/L is a 24-hour light schedule with an intensity of 1. D/D indicates dark light intensity; L/D, light/dark cycle; L/L, bright light intensity.

In Figure 5, we have emulated the condition of time-restricted feeding (access to food is restricted for specific time intervals during the day without calorie restrictions[13]) and tested the effects of phase relations between light and feeding signals under 2 different feeding durations. In all cases, we assumed that the “caloric contents” of the nutritional cues were equivalent by maintain constant area under the curve (AUC) of the nutritional signal. Either a 12-hour duration feeding signal at an amplitude of 1 (represented by thick lines) or a 6-hour duration feeding signal at an amplitude of 2 (represented by thin dashed lines) was applied with varying start time delays relative to light. The amplitudes for feeding signals with different durations were set such that the AUC of feeding for a given 24-hour period is identical between the 2 feeding durations tested. From the resulting dynamics, the amplitudes and phase angles for cortisol (Figure 5A), PER/CRY (Figure 5B), and SIRT1 (Figure 5C) profiles were shown in compass plots. The amplitudes are represented as the lengths of the arrowheads, and phase angles are represented as the directions of the arrowheads in the plot. The amplitude changes for cortisol were accompanied with changes in phase relations between light and feeding, characterized by smaller amplitudes in the first and fourth quadrants compared with second and third quadrants. However, feeding duration did not affect the amplitude of cortisol oscillations as much as the effects caused by feeding delays. The largest percent difference between amplitudes caused by feeding duration was 14%, observed at the light-feeding phase difference of 20 hours. Compared with this case, the largest percent difference in amplitude due to light-feeding phase difference was 62%, observed between 4 and 16 hours for the 12-hour duration feeding. However, PER/CRY amplitudes were influenced both by feeding duration and feeding delays, characterized by the difference in lengths between thick solid and thin dashed arrowheads as well as the changing amplitude around the plot. For both feeding durations, amplitudes were higher when light-feeding phase difference was small (0-4 hours). Amplitudes decreased as phase difference increased, with a minimum at 16 hours of phase difference. For all light-feeding phase differences, amplitudes were greater for the 6-hour duration feeding compared with 12-hour duration feeding. Furthermore, feeding delays in 6-hour duration schedule resulted in phase-jump behavior, characterized by irregular angles between each thick arrow. For example, in the 12-hour duration feeding case, the difference in phase angle between 8- and 12-hour light-feeding phase difference is about 100°, whereas that between 16- and 20-hour difference is less than 10°. In contrast, SIRT1 phase angles strictly followed the feeding schedules while amplitudes were identical for each feeding duration. SIRT1 amplitude for 6-hour duration was greater than that of the 12-hour duration due to the difference in feeding signal strength.

Figure 5.

(A) Amplitude and phase of cortisol, (B) PER/CRY protein, and (C) SIRT1 at various light-feeding phase relations. Thick lines represent data with 12-hour feeding duration with an amplitude of 1, and thin dashed lines represent data with 6-hour feeding duration with an amplitude of 2. SIRT1 indicates sirtuin 1.

To explore the effects of feeding duration and phase relative to light in more detail, the amplitudes of cortisol and PER/CRY were plotted along an axis of light-feeding phase difference at 30-minute intervals in Figure 6. A positive value of the light-feeding phase difference means that feeding signal was started after the start of light signal. Therefore, feeding was started at ZT0, ZT0.5, ZT1, ZT1.5, etc. Feeding durations of 12 hours at an amplitude of 1, 6 hours at an amplitude of 2, and 2 hours at an amplitude of 6 were tested, keeping the AUC of feeding identical over a 24-hour period. In Figure 6A, we observe that the cortisol amplitude is high when feeding starts shortly after light phase and drops as feeding delay increases. The amplitude is at the maximum value when feeding is started at ZT3 for 12-hour duration feeding, ZT4 for 6-hour duration feeding, and ZT6 for 2-hour duration feeding. The amplitudes reach the minimum values at feeding start time of ZT16 for the 12-hour duration feeding, ZT19 for the 6-hour duration feeding, and ZT23 for the 2-hour duration feeding. The feeding delay time at the minimum amplitude is offset from the time of complete antisynchrony between light and feeding, which is at ZT12. In Figure 6B, the amplitudes for PER/CRY protein are plotted in the identical manner. The influence of feeding duration is apparent from this plot, as the amplitudes for shorter feeding duration are higher than those for longer feeding duration for most of the delay times, from ZT1 to ZT19. The PER/CRY protein oscillates with larger amplitude when feeding starts early morning or a few hours before the light period starts. The PER/CRY amplitude is at maximum when feeding starts 2 hours before light under 12-hour duration feeding schedule. For 6-hour duration feeding, starting feeding at ZT1 gives the maximum amplitude, and starting feeding at ZT3 will give maximum amplitude for 2-hour duration feeding. From then on, amplitudes decline as feeding delay is increased. Unlike the cortisol, there is a range of delay times in which the PER/CRY oscillates at minimum amplitude. For example, a feeding start time between ZT7 and ZT17 will cause PER/CRY to oscillate at the minimum amplitude. Similar ranges exist for 6- and 2-hour duration feedings and are centered around the delays exhibiting minimum correlation coefficients between light and feeding. The slopes leading to and exiting from the minimum ranges are asymmetric and can be characterized by a slow decline of amplitude and quick recovery of amplitude with increasing feeding delays.

Figure 6.

(A) Amplitudes of cortisol and (B) PER/CRY proteins at different feeding delay times relative to light.

Sensitivity coefficients for PER/CRY protein amplitude and phase were computed as described earlier, once with synchronized light and feeding signals and another time with antisynchronized light and feeding signals. In Figure 7, the parameters with 5 largest sensitivity coefficients for each condition based on PER/CRY amplitude were selected, and their sensitivity coefficients under both conditions are shown. Sensitivity coefficients for all parameters can be found in Supplementary Figure 1. When light and feeding signals were synchronized, the most sensitive parameters were k1, K1, K2, V3, and k4. Among these parameters, the first 4 are associated with the Goodwin oscillator in the HPA axis. Last parameter k4 is the Michaelis constant for Bmal1 transcription. When light and feeding signals are antisynchronized, k1, k3, k6, k7, and nad appear to be most sensitive to perturbations. Among these parameters, k1, k3, and nad are associated with the dynamics of NAD+ concentration in the periphery. The other 2 parameters, k6 and k7, describe the dynamics of activated SIRT1. Finally, the amplitude of PER/CRY protein was more sensitive than phase angle when the parameters were perturbed, under both synchronized and antisynchronized light and feeding conditions.

Figure 7.

Sensitivity coefficients for PER/CRY protein under synchronized and antisynchronized light and feeding signals. Sensitivity coefficients were calculated based on PER/CRY amplitude and phase angles.

Discussion

The role of circadian rhythms on metabolic activity has been well established in mammals. Maintaining homeostasis of plasma glucose level throughout the day is extremely important for mammals because insufficient or excessive glucose levels can have detrimental effects on key biological functions such as neuronal activity and balancing of body fluid and electrolytes.[105] Important metabolic pathways that contribute to steady glucose level, such as gluconeogenesis and glycogen metabolism, are linked to proper functioning of the peripheral clock machinery.[22,106] Lipid metabolism has also shown to exhibit circadian activities[107] while also being regulated by BMAL1,[108] a core component of the feedback loop that creates circadian rhythms. Therefore, it is not surprising that circadian disruption can result in metabolic syndrome[16] or that animals with metabolic syndromes exhibit attenuated circadian rhythms in the periphery.[109] Although the metabolic implications from circadian disruption may manifest through multiple unknown mechanisms, there are some distinguishing effects on the clock genes upon changes in meal timing. One of the common and well-established observations is that restricting food access to the animals’ rest period can re-entrain the circadian rhythms in the periphery to be synchronized to the feeding schedule, uncoupling the rhythms of the periphery from the central clock in the SCN. More specifically, when mice with pancreatic adenocarcinoma had 4 hours of access to food during the light period, PER2 and BMAL1 peaks were advanced by 8 hours.[110] In another study, diurnal fed mice showed completely inverted expression of D-site–binding mRNA and protein responsible for some clock output genes.[111] Furthermore, the phase inversion in the periphery is achieved whether the animals are subject to a 12-hour light, 12-hour dark (LD) cycle, or a constant dark (DD) schedule. When feeding time was advanced by 7 hours in mice, the peaks for Per1 and Per2 mRNA, D-site–binding protein, and Cyp7A mRNA in the liver were advanced by 6 to 12 hours.[96]

As shown in Figure 2, our model predicts that cortisol peak will be pushed to the beginning of the feeding time. At the same time, PER/CRY protein rhythms will completely invert upon feeding inversion, re-entraining to the new feeding schedule. The peak was shifted from early active phase to early rest phase, consistent with the literature data. A loss of amplitude is also observed for the PER/CRY rhythms, also consistent with findings that feeding exclusively during rest phase results in reduced oscillations for PER2.[96] The complete re-entrainment of PER/CRY rhythms to feeding is achieved through 2 major channels in our model. First, feeding-entrained NAD+ affects the synthesis of cortisol through a transit compartment, EntF. As a result, cortisol exhibits an altered rhythmic pattern (Figure 2B), peaking near the feeding start time. The peak-delay behavior is reflective of the food-anticipatory rise in glucocorticoid levels in mice subjected to altered feeding schedules.[112] Unlike PER/CRY, cortisol rhythm is not completely shifted by 12 hours because it is strongly entrained by the light signal in the HPA axis. Light negatively regulates the CRH, which promotes ACTH, which in turn stimulates the release of cortisol. The cortisol-receptor complex promotes the expression of PER/CRY with the delayed rhythm, pushing the peak location of PER/CRY back. However, altered cortisol profile alone cannot explain the complete 12-hour phase shift of PER/CRY as the cortisol peak was not shifted by 12 hours. The complete phase inversion is achieved by the re-entrainment of CLOCK/BMAL1 heterocomplex by SIRT1. In our model, SIRT1 binds to CLOCK/BMAL1 to form CLOCK/BMAL1/SIRT1 complex, which serves as a transcription factor for NAMPT. Therefore, SIRT1 depletes the pool of CLOCK/BMAL1 complex and delays the peak time. Because CLOCK/BMAL1 promotes the transcription of Per/Cry genes, the peak time of PER/CRY protein is further delayed and achieves a complete 12-hour inversion.

We hypothesized that imposing a constant light condition will result in attenuated amplitudes for cortisol and PER/CRY rhythms while improving the transition time between feeding rhythm changes. As shown in Figure 3, we confirmed that constant light signals result in reduced amplitudes for cortisol and PER/CRY rhythms. However, SIRT1 amplitude or phase was unaffected by the light intensity because it is tightly controlled by the feeding-entrained NAD+. Although all constant light signals at all intensities result in lower amplitude oscillations for cortisol and PCGs, the intensity of the constant light signal affects their dynamics. Higher light intensity, or bright light, attenuates the oscillation amplitude more than dim light. This prediction agrees with experimental findings in mice.[104] Light signal entrains the cortisol by facilitating the degradation of CRH in equation (5). Brighter light will result in greater degradation of CRH, which yields lower production of ACTH and therefore less cortisol (Figure 3B). Because the PCGs are entrained by cortisol in the periphery, the PER/CRY amplitude is reduced more under bright light in return. Imposing constant light signal also alters the transition time of the PCGs upon feeding inversion. As shown in Figure 4, the PCGs adjust to the new feeding pattern more quickly under constant light condition than a 12-hour light, 12-hour dark cycle. This model behavior is consistent with the experiment that in mice carrying GR null alleles exclusively in the hepatocytes, the PCGs in the liver entrain much faster to the feeding regimen compared with the wild-type mice.[112] Because cortisol and cortisol-receptor interaction are essential in entraining the PCGs by light, we rationalize that removing the light schedule creates a similar environment in silico and observe that faster entrainment to feeding is achieved under constant light conditions. In our model, transition time is faster under constant light because the light/dark cycle acts as a conflicting entrainer and inhibits the transition to the feeding regimen.

In Figure 5, our model predicts that restricting feeding to a few hours during the active period results in higher amplitude oscillations for the peripheral clock machinery. In this figure, 2 different feeding durations (12 and 6 hours) with identical AUC over a given 24-hour period were tested at varying start times relative to light. The PER/CRY amplitudes are the highest when a 6-hour duration feeding is started at ZT0 and ZT4 (Figure 5B). The prediction is in qualitative agreement with animal studies. When a mouse model of diet-induced obesity was subjected to 8 hours of restricted feeding during the active phase, Per2 and Bmal1 mRNA oscillations were enhanced with increased amplitudes, for both high-fat-content and normal diets.[40] The increase in amplitude with shorter feeding duration is most likely due to the increased amplitude in the oscillation of SIRT1. If the feeding duration is shortened, the intensity of feeding signal has to be higher to compensate for the short duration and provide identical AUC of feeding to the system. Because SIRT1 is directly activated by NAD+ which is closely tied to feeding signals, SIRT1 amplitude also increases under shorter feeding duration, as shown in Figure 5C. Then, the SIRT1 dynamics will affect PER/CRY rhythms by facilitating the degradation rate inside the nucleus, eventually decreasing the protein’s self-inhibitory effect on its own transcription rate (equation (30)). In addition to the feeding duration, feeding start time relative to light also affects the amplitude of PER/CRY oscillation. For 12-hour feeding duration, starting feeding between ZT8 and ZT16 resulted in low-amplitude oscillations, whereas feeding start time of ZT12 to ZT20 resulted in low-amplitude oscillations for 6-hour duration feeding. The light-feeding phase differences that yield low-amplitude PER/CRY oscillations also result in lower amplitude oscillations for cortisol (Figure 5A). This phenomenon is due to that light and feeding both exert influence on cortisol, much like PER/CRY, although the light entrainment is much stronger in cortisol secreted from the HPA axis. We also observe from Figure 5 that phase angles progress through the feeding start time in different ways for cortisol, PER/CRY, and SIRT1. The irregularities in phase angles are expected as amplitude death in coupled limit-cycle oscillators is associated with phase-flip behavior,[113] where the phase angle between the 2 oscillators suddenly increases at a threshold difference between the external entrainers. SIRT1’s phase angles over the varying feeding start times are evenly spaced because it is mainly entrained by the feeding signal alone. However, PER/CRY and cortisol phase angles progress through the feeding start times at varying angles because they are coupled with both light and phase signals entering externally.

In Figure 6, the amplitudes for cortisol and PER/CRY through the varying feeding start times are shown in more detail, to identify which light-phase phase relations would yield the highest and lowest amplitude oscillations. Simulation of the model was performed at 30-minute intervals at 3 different feeding durations (12, 6, and 2 hours). The amplitudes of the feeding signals were adjusted to match the AUC of feeding signal over a 24-hour period. From Figure 6A, cortisol oscillation amplitudes exhibit identical trends among the 3 feeding durations where they are initially high when feeding starts soon after the light starting time, decline slowly to reach a minimum, and then recover to the maximum amplitude quickly. However, the specific feeding start time that gives the maximum and minimum amplitudes are different for every feeding duration. For example, under the 12-hour feeding duration condition, the maximum amplitude is reached at feeding start time of ZT3 and the minimum is reached at ZT16. But these are not the times that give maximum and minimum amplitudes for 6- or 2-hour feeding durations. Interestingly, 12-hour delay, or complete inversion between light and feeding, does not result in the most reduced oscillation of cortisol. Such asymmetry is observed for PER/CRY amplitude profiles (Figure 6B) but manifested in a different way. Our model predicts that there is a range of feeding start times where PER/CRY amplitudes are also at a minimum value for each feeding duration. For the 12-hour duration, this minimum range is centered around feeding start time at ZT12. However, the slope of the amplitude profiles before and after the minimum values is asymmetric. The amplitude slowly declines to the minimum and then recovers fast to reach the maximum either at the end of the dark period or at the beginning of the light period. Considering the characteristics of PER/CRY and cortisol together, these results suggest that starting to eat earlier in the day will give more robust circadian rhythms in the periphery than starting to eat later in the day. The model prediction is consistent not only with the common belief that late-night snacks are disadvantageous for the health but also with studies where distributing more calories to breakfast resulted in more weight loss under the same caloric intake.[114] Furthermore, amplitude death in coupled limit-cycle oscillators are associated with phase-flip behaviour,[113] where the phase angle between the 2 oscillators suddenly increases at a threshold difference between the external entrainers. In a symmetric system, where the entrainer strength and the coupling strength between the 2 oscillators are identical, the phase flip is centered around 12-hour delay between the entrainers, provided that the oscillations have a 24-hour period.[115] However, the timing of phase flip moves when either the relative strengths of entrainers or the coupling strengths are asymmetric. Integrating these observations, the offset of the minimum cortisol amplitude and asymmetric PER/CRY profile may suggest that our model is an asymmetric system whose entraining strengths of light and feeding are transmitted to the PCGs with different efficiencies. Because the PCGs eventually entrain to the feeding cycle and the cortisol peak is phase delayed by 7 to 8 hours on a 12-hour feeding inversion, feeding signals have a greater impact on the downstream events than the light entrainment on the PCGs.

The sensitivity analysis reveals that different sets of parameters have highest sensitivity coefficients under synchronized and antisynchronized light and feeding signals. When the 2 external signals are synchronized, parameters associated with the Goodwin oscillator describing the production of cortisol in the HPA axis have the highest sensitivity coefficients, as shown in Figure 7. These parameters include k1 and K1, each describing the CRH production and dissociation rates; K2, which estimates dissociation rate of ACTH; and Vd3, which describes the rate of cortisol degradation. Together, these parameters affect the dynamics of the light-entrained negative feedback loop in the HPA axis (equations (5)-(7)), and it is expected that these parameters are among the most sensitive parameters because the secreted cortisol entrains the downstream PCGs. Parameter k4, the Michaelis constant of Bmal1 transcription, also has a high sensitivity coefficient, likely due to SIRT1-BMAL1 binding interaction as well as its involvement in controlling the dynamics of the pro-inflammatory cytokines. The pro-inflammatory cytokine dynamics and BMAL1 dynamics both directly affect the transcription of Per/Cry, leading to higher sensitivity. When light and feeding are antisynchronized, parameters related to feeding entrainment exhibit higher sensitivity coefficients. This is likely due to that in the antisynchronized state, feeding entrainment decouples the PCGs and cortisol from the Goodwin oscillator from the HPA axis, exerting more influence on the downstream events. Parameters k1, k3, k6, k7, and nad appear to be most sensitive to perturbations. Among these parameters, k1, k3, and nad are associated with the dynamics of NAD+ concentration in the periphery. Parameter k1 is the maximum extent of NADH converting to NAD+, and k3 describes the same for NAD+ converting to NADH on feeding. The combined concentration of NAD+ and NADH is parameter nad. The other 2 parameters, k6 and k7, describe the dynamics of activated SIRT1. They each describe the maximum extent of SIRT1 activation mediated by NAD+ and the maximum extent of SIRT1 degradation. SIRT1 is one of the key components in our model because it delivers the feeding signal to the peripheral clock machinery. Therefore, the sensitivity caused by SIRT1 is due to multiple interactions SIRT1 has with the PCGs. One of the ways in which SIRT1 interacts with the PCGs is by binding to the CLOCK/BMAL1 complex to form the CLOCK/BMAL1/SIRT1 complex, which drives the expression of NAMPT. Meanwhile, the CLOCK/BMAL1 complex is involved in the promotion of PER and CRY expression. In addition, the amplitude of the PCGs appears to be more sensitive than the phase angle of the PCGs, indicating that the phase relations achieved by the 2 entrainers, light and feeding, are robust, whereas the oscillatory strength can change due to perturbations to the parameters.

In summary, our model qualitatively captures the key features of feeding-entrained peripheral clock machinery. Through the simulations, we could relate the robust circadian rhythms under short duration feeding in active phase to the higher amplitude oscillation of SIRT1 due to strong but short feeding signal. Our prediction suggests that controlling the dynamics of SIRT1 may be helpful in restoring or strengthening the oscillations of the peripheral clock machinery, also supported by a study on methylselenocysteine, where restored NAD+ oscillations and enhanced SIRT1 activity resulted in restoration of circadian rhythms of mouse mammary tumor model.[116] Our model also predicts that the phase relation between light and feeding plays an important role in determining the oscillation amplitude for the PCGs, highlighting the need to study the metabolic implications caused by the interplay between these 2 environmental cues. Furthermore, analysis of the model suggests that the reason best feeding time lies earlier in the active phase may be due to the asymmetry between the efficiencies of light and feeding entrainment of the PCGs. Future work could involve studying further downstream metabolic activities such as hepatic gluconeogenesis and lipogenesis to explore the changing dynamics of energy homeostasis in relation to other metabolic genes such as mTOR and AMPK under circadian disruption. Clearly, a number of issues remain to be further examined, including but not limited to, assessing the impact of nutritional composition as well as multiple patterns of metabolic rhythms, such as simulating multiple meals.