Entrainment of peripheral clock genes by cortisol

Panteleimon D. Mavroudis, Jeremy D. Scheff, Steve E. Calvano, Stephen F. Lowry, Ioannis P. Androulakis


Circadian rhythmicity in mammals is primarily driven by the suprachiasmatic nucleus (SCN), often called the central pacemaker, which converts the photic information of light and dark cycles into neuronal and hormonal signals in the periphery of the body. Cells of peripheral tissues respond to these centrally mediated cues by adjusting their molecular function to optimize organism performance. Numerous systemic cues orchestrate peripheral rhythmicity, such as feeding, body temperature, the autonomic nervous system, and hormones. We propose a semimechanistic model for the entrainment of peripheral clock genes by cortisol as a representative entrainer of peripheral cells. This model demonstrates the importance of entrainer's characteristics in terms of the synchronization and entrainment of peripheral clock genes, and predicts the loss of intercellular synchrony when cortisol moves out of its homeostatic amplitude and frequency range, as has been observed clinically in chronic stress and cancer. The model also predicts a dynamic regime of entrainment, when cortisol has a slightly decreased amplitude rhythm, where individual clock genes remain relatively synchronized among themselves but are phase shifted in relation to the entrainer. The model illustrates how the loss of communication between the SCN and peripheral tissues could result in desynchronization of peripheral clocks.

  • clock synchronization
  • biological oscillations
  • circadian rhythms
  • mathematical model

exogenous daily signals (light, feeding, and temperature) give rise to endogenous circadian rhythms in virtually all physiological functions. Many circadian oscillations are ultimately driven by the activity of the suprachiasmatic nucleus (SCN) (61), which responds to light by emitting signals that control body temperature, endocrine hormone secretion, and activity (16). The genes responsible for circadian timekeeping in the SCN include Per1–3 and Cry1–2, which, in collaboration with the CLOCK/BMAL1 heterocomplex, form an orchestrated feedback loop that is kept in sync by light-induced signaling and leads to an overall circadian rhythm. The signal transduction of circadian rhythms from the SCN to peripheral tissues consists of centrally mediated signals that synchronize peripheral free-running oscillators. In fact, the molecular machinery responsible for oscillations in the SCN is also generally found in peripheral cells, allowing them to maintain some oscillatory characteristics even in the absence of a central entraining signal.

Independent of the central pacemaker, peripheral tissues exhibit circadian rhythmicity through the oscillations of clock genes (7). These peripheral circadian oscillations have profound physiological importance in regulating the circadian functions of numerous tissues such as the heart, liver, and blood (25, 39, 58). In vitro, peripheral cells initially retain circadian rhythmicity (9, 74), but eventually coherent rhythmicity in the ensemble average disappears. Experiments in rat-1 fibroblasts revealed that the observed dampening of clock gene expression in tissue explants is not the result of lost rhythmicity in individual cells; rather, circadian rhythms are lost as the individual oscillators fall out of sync in the absence of an entrainer, resulting in a blunted ensemble average (55, 72). Therefore, the characteristics of endogenous entraining signals and their mechanisms of action on peripheral tissues are of critical importance in understanding the dynamics of circadian rhythms.

Although there are many endogenous signals that may be responsible for the circadian entrainment of peripheral oscillators (66), glucocorticoids (cortisol in humans) are particularly interesting as putative entrainers. Cortisol secretion from the adrenal glands is regulated by the hypothalamic-pituitary-adrenal (HPA) axis, and it follows a circadian pattern with a zenith in the early morning and a nadir late at night. Circadian rhythms of glucocorticoids regulate expression of peripheral genes through gene/receptor binding to a specific DNA sequence at the promoter region of target genes named the glucocorticoid response element (GRE) (57). Particularly for clock genes, experiments in human peripheral blood mononuclear cells (PBMCs) and fibroblasts (7, 17) established that glucocorticoids, after binding to glucocorticoid receptors, induce the expression of the clock gene Per1 through a GRE in the Per1 genome sequence (73). Further experiments in mesenchymal stem cells revealed the presence of GRE in Per2 and E4bp4 clock gene locus that were continuously occupied with glucocorticoid receptor upon treatment of cells with the synthetic glucocorticoid dexamethasone (65). In addition, genomic deletion of the GRE in Per2 resulted not only in the failure of glucocorticoids to stimulate Per2 response, but also in dampened expression of other clock genes (i.e., Rev-Erbα), illustrating the downstream regulation of the circadian clock network through glucocorticoid receptor binding to GRE. Collectively, these experiments confirm that glucocorticoids, through glucocorticoid receptor signaling, directly stimulate the synchronized rhythmicity of Per1, Per2, and E4bp4 in peripheral cells, further regulating the peripheral clock gene network. Therefore, although there may be multiple systemic circadian signals that entrain peripheral clock genes, one of our underlying hypotheses is that we can use cortisol as a representative entrainer of peripheral cells to explore the dynamics of clock gene synchronization and entrainment by systemic cues. Cortisol is also intriguing as a circadian entrainer given its central role in the inflammatory response and the recent observation that acute changes in both cortisol and clock gene expression occur in peripheral blood leukocytes in response to endotoxemia (38). The disruption of circadian rhythms in cortisol is associated with fatigue, weight loss, insomnia, coronary heart disease, and tumor progression (29, 31, 47, 64). There has also been interest in cortisol circadian rhythmicity as a predictor of breast cancer survivor (63).

The importance of deciphering clock gene dynamics due to their role in regulating the circadian function of numerous tissues, such as heart, liver, and blood, as well as the complexity of clock gene network and entrainment characteristics, motivates the need for mathematical modeling of the peripheral clock network. Mathematical approaches can be of great help in understanding the underlying dynamics and also predicting clinical outcomes and intervention strategies. Several mathematical models have been proposed to investigate and describe the dynamics of clock genes (4, 5, 11, 12, 30, 34, 37, 46, 51). These models, while varying considerably in their underlying assumptions, their degree of complexity, and their method of implementation, all converge to the inclusion of a negative feedback loop that represents the Per/Cry genes and the CLOCK/BMAL1 heterodimer. In line with the experimental evidence described above, we propose a mathematical model of peripheral clock genes that incorporates cortisol as a systemic entrainer. This computational representation linking central and peripheral oscillators is leveraged to study the entraining properties of a central circadian signal on a population of peripheral cells by integrating models of circadian cortisol production (19), glucocorticoid pharmacodynamics (43, 57), and peripheral clock genes (11). To account for heterogeneity between individual peripheral cells, the model is formulated as a system of stochastic differential equations (SDEs).

We observe that cortisol rhythmicity induces peripheral clock gene entrainment and synchronization in an amplitude and frequency-dependent manner. While homeostatic entraining rhythms stimulate a homogeneous circadian pattern to the population of peripheral cells, the loss of circadian amplitude provokes a desynchronization among the population of cell phases. This biological shift from synchronization to desynchronization progresses through a dynamical state where the individual cells retain a relative phase coherency but are phase-shifted relative to the entrainer. Concerning entrainer frequency, peripheral cells remain synchronized only for cortisol frequencies relatively close to the individual cell frequencies. In addition, we observe that even when cells are totally entrained by cortisol, synchronization varies during the day pointing its lowest values when cells are near their nadir or zenith level.



Cortisol production and signal transduction.

Peripheral circadian clocks are entrained by systemic cues, such as the circadian release of cortisol. The circadian production of cortisol is modeled based on the “two rates” model (19) where a zero order production term (RF) switches between two different values representing a high production rate of cortisol in the morning and a low production rate the rest of the day (Eq. 1). We have previously used this model due to its simplicity while still retaining the ability to fit experimental data (62). Due to the physiological importance of glucocorticoids, there has been much interest in the pharmacokinetics and pharmacodynamics of glucocorticoids to understand their behavior in response to endogenous glucocorticoid production by the HPA axis or exogenous treatment. Here, we apply a glucocorticoid pharmacodynamic model which incorporates the steps of the glucocorticoid signal transduction pathway (Eqs. 25) (57), culminating in transcriptional regulation of clock genes. To ultimately regulate the transcription of glucocorticoid responsive genes, cortisol diffuses in the cytoplasm of the target cells and binds with free glucocorticoid receptors (R) forming a glucocorticoid-receptor complex (FR), which eventually translocates to the nucleus [FR(N)] where it dimerizes and binds to GRE. This binding leads to decreased levels of glucocorticoid receptor mRNA (mRNAR, Eq. 2), through a process known as homologous downregulation (56). Then, part of the receptor is recycled back to the cytoplasm and is used anew to bind to glucocorticoids. dFdt=RF+kin,Fenkout,F.FRF={kin,RF10tF1<mod(t,24)<tF2tF2<mod(t,24)<tF1 1 dmRNARdt=ksyn_Rm(1FR(N)kacy7IC50_Rm+FR(N)kacy7)kdgr_RmmRNAR2 dRdt=ksyn_RmRNAR+RfkreFR(N)konFRkdgr_RR3 dFRdt=konFRkTFR4 dFR(N)dt=kTFRkreFR(N)5

It is also known that glucocorticoid receptor binding to GRE is inhibited by the action of the CLOCK/BMAL1(y7) heterodimer that acetylizes glucocorticoid receptor nuclear complex [FR(N)] in its hinge region and blocks its forward action (53). This is modeled by incorporating an inhibitory term based on FR(N) and y7 in the equation mRNAR (Eq. 2). kac parameter that represents the acetylation rate has been chosen adaptively so as not to change the expression levels of the model compartments that were fitted to experimental data. Other parameters were the same as the original models (57, 62) (Table 1).

View this table:
Table 1.

List of parameters used in the model

Chronic emotional or physical stress can manifest itself through flattened cortisol rhythmicity (21, 47). Whether this blunted cortisol rhythm is followed by increased (hypercortisolemia) or decreased (hypocortisolemia) levels of cortisol is highly dependent on the various pathologies (21, 50, 64). Particularly, in a study investigating cortisol rhythms in severely burned adults and children, stressed patients exhibited a blunted cortisol rhythmicity due to higher values of cortisol nadir levels (40). In line with this experimental evidence, we explored the effects of the loss of circadian amplitude in cortisol by increasing cortisol's nadir values (Fmin) and holding cortisol's zenith (Fmax) constant. To do this, we increased the lower value of the zero order production term of cortisol, RF (Eq. 1), while holding the higher value constant. In addition, we explored the behavior of our model as cortisol loses its homeostatic frequency pattern while holding its amplitude constant. We further did this by changing the modulus term [mod(t,24)] of Eq. 1. This provides a consistent way to evaluate the effects of decreased cortisol rhythmicity, in line with experimental evidence.

Peripheral clock genes regulation by cortisol.

We model the interaction between peripheral circadian clocks and cortisol by considering the transcriptional effect of activated glucocorticoid receptor on components of the clock gene network. This peripheral network incorporates a positive and a negative feedback interaction through which peripheral clock genes exhibit oscillatory activity (11). In particular, the periodic expression of clock genes is driven by Per and Cry inhibiting the activity of the CLOCK/BMAL1 dimer (negative feedback) and stimulating Bmal1 gene transcription (positive feedback). Through the negative feedback loop, the heterocomplex CLOCK/BMAL1 activates the transcription of period (Per) and cryptochrome (Cry) genes (y1) upon binding to the E-box promoter region. After the expression of PER/CRY proteins (y2) in the cytoplasm, they translocate to the nucleus (y3) where they inhibit their own transcription by shutting off the transcriptional activity of the CLOCK/BMAL1 heterocomplex (y7). Through the positive feedback loop the nuclear compartment of PER/CRY protein (y3) activates indirectly Bmal1 mRNA (y4) transcription, which after its translation to BMAL1 protein (y5) and its translocation to the nucleus (y6) increases the expression of CLOCK/BMAL1 heterodimer. In reality, PER/CRY proteins regulate Bmal1 transcription through inhibiting the CLOCK/BMAL1-mediated expression of Ror (α, β, and γ) and Rev-Erb (α and β) genes (26, 51, 59). While both simpler and more complex models of the circadian network exist, this system represents a good intermediate level of complexity. This network of interactions, without the entraining effect of cortisol, produces the rhythmic expression of clock genes with a period of 23.4 h. Similar periods have been shown to be adopted by peripheral clock genes oscillators in human fibroblasts (24.5 h with standard deviation of 45 min among individuals) (15) and in human PBMCs (23.91–24.45 h) (13). The activated nuclear glucocorticoid receptor complex directly regulates the expression of clock genes by binding to GRE. Similar to the glucocorticoid pharmacodynamic signal transduction pathway described above, CLOCK/BMAL (y7) antagonizes this binding to GRE, which, based on experimental data (65, 73), is present in the promoter region of Per1 and Per2. As such, glucocorticoids, inhibited by CLOCK/BMAL, entrain the local network of clock genes by regulating the transcription of Per1–2. This circadian regulation then propagates through to the other model variables, imposing a constant phase on the peripheral oscillators based on the circadian rhythm of cortisol. The behavior of these core circadian components of the peripheral cells is described by Eqs. 612 (Table 1). dy1dt=v1b(y7+c)k1b(1+(y3k1i)p)k1dy1+kcFR(N)kacy76 dy2dt=k2by1qk2dy2k2ty2+k3ty37 dy3dt=k2ty2k3ty3k3dy38 dy4dt=v4by3rk4br+y3rk4dy49 dy5dt=k5by4k5dy5k5ty5+k6ty610 dy6dt=k5ty5k6ty6k6dy6+k7ay7k6ay611 dy7dt=k6ay6k7ay7k7dy712

The entraining effect of cortisol nuclear complex that takes into account receptor's acetylation by CLOCK/BMAL1 (y7) is modeled similar to other pharmacodynamic models simulating GRE binding (76) as an additive term [kc*FR(N)/kacy7] on the equation of Per/Cry mRNA (y1), Eq. 6. An additive term in the equation for Per/Cry mRNA was also used by Geier et al. (34) to simulate the light entrainment of clock genes in the SCN, given that light differentially regulates the individual mRNAs that comprise the lumped variable (Per/Cry mRNA, y1). A similar phenomenon is hypothesized to occur with clock gene responses to glucocorticoids, as only Per1 and Per2 have active GREs (65, 73). Although in reality Bmal1 transcription (positive feedback) is indirectly regulated by PER/CRY proteins through the antagonistic effects of ROR (α, β, and γ) and REV-ERB (α and β) transcription factors on the retinoic acid-related orphan receptor response element (RORE) in the promoter region of Bmal1 gene (26, 51, 59), it is not the scope of our work to fully describe the components of this network but, rather, to use a fairly simple representation based on which we can investigate the dynamics of clock gene entrainment by an entraining signal (cortisol).

Thus, we leverage this relatively simple mathematical model of the central pacemaker as a representative peripheral framework to test and explore the dynamics of entrainment of peripheral clocks by a centrally regulated entrainer. The network diagram of the integrated model (Eqs. 112) is displayed in Fig. 1.

Fig. 1.

Schematic figure of the model. Two components formulate the framework of the model. The 1st is a pharmacodynamic compartment through which cortisol F diffuses to cytoplasm, binds the glucocorticoid receptor (R), forming the complex FR, translocates to the nucleus FR(N), and regulates the translation of mRNAR and Per/Cry mRNA, and the 2nd includes the clock gene regulatory positive and negative feedback loops (y1y7).


SDE formulation.

All biological processes [e.g., protein folding (27), gene expression (44), protein-protein interaction (10)] proceed with some level of uncertainty. Relevant to the transcription of clock genes, the binding of a transcription factor such as CLOCK/BMAL1 to a small number of discrete regions to modulate transcription is a highly stochastic process (44). Furthermore, cortisol secretion as well as its signaling pathways contain stochasticity (22, 45), which has been modeled from both the perspectives of hormone release (14) and downstream transcriptional regulatory effects (2, 3, 49). Our model incorporates stochastic expressions for cortisol and clock genes. This uncertainty leads the individual cells, in absence of an entraining signal, to adopt stochastic phases and periods and thus fall out of sync.

The system of ordinary differential equations defined in Eqs. 112 was translated to chemical Langevin stochastic differential equations (CLEs) (35). Generally, this formulation assumes that molecules of the reacting species interact through a number of chemical reactions, each one having a probability of occurrence in an infinitesimal time interval proportional to the rate constant of the reaction. This probability is referred as the propensity of the reaction (Table 2). The CLE recurrence formula used in this work is defined as: xi(t+τ)=xi(t)+j=1Mvjiaj(xi(t))τ+j=1Mvijaj(xi(t))τNj(0,1)13 Where xi (t) are the molecules of the i reacting species (Eq. 112) at time t, τ is the fixed time interval length, j is the different type of the total M = 31 chemical reactions (Table 2), vji is the change in the number of the i molecules caused by the j reaction, and aj is the propensity of the j reaction. Nj are statistically independent, normal (0,1) random variables. The time interval length (τ) has been chosen adaptively so that the propensity of reactions in each time step is slightly changed. To translate the determinist model into stochastic, we needed to know the volume of the cell. Similar to Forger and Peskin (30), this volume was chosen so that the peak of PER/CRY protein molecule count (y2, Eq. 7) is nearly 5,000 molecules, in accordance with experimental evidence from liver cells.

View this table:
Table 2.

Elementary reactions and reaction propensities involved in the model

Quantification of synchronicity, entrainment, and periodicity.

To study how circadian rhythms in cortisol affect the state of clock genes, we evaluate metrics of both synchronicity (i.e., how similar are different cells) and entrainment (i.e., how similar cell periods and phases are to the entrainer). To assess synchronization of clock genes, we measure the deviation of Per/Cry mRNA levels from their mean value. Therefore, in our analysis we incorporated the degree of synchronization for the variable yk, Rsyn,k (33, 71) (Eq. 14), which represents the ratio of the variance of the mean field over the mean variance of each oscillator. Rsyn,k=y¯k2y¯k21Nj=1N(yk,j2yk,j2)14

In Eq. 14, yk,j is the time-course vector output generated by Eqs. 612, where the first index k represents the variable, and the second index j represents the cell (N total). Overbar variables (ȳk) are time-course vectors of averages over the population of N cells (Eq. 15) and brackets represent time averages (Eq. 16). y¯k=1Nj=1Nyk,j15 y¯k=1Tt=t1Ty¯k(t)16

Rsyn,k is calculated for the Per/Cry mRNA variable (k = 1) for a time span of 500 h except for the in silico experiments evaluating circadian variation of synchronization where it was calculated for a period of 2 h. Rsyn,k has a minimum value of zero when the cells are completely desynchronized and a maximum value of 1 when cells are fully synchronized.

In addition to intercellular synchronization, we evaluate the level of entrainment of peripheral cells by cortisol's rhythmic pattern. For this reason we quantify the degree of phase locking of Per/Cry mRNA. Accordingly, we calculate the entrainment index “ρ1,” which is an entropy-based metric that takes into consideration the individual cell period distribution (69). The entrainment index for any variable k of the model is defined as: ρk=1SkSk,max17 Where Sk=l=1NPk,l · ln(Pk,l) is the entropy of the discrete period distribution for the variable k, l is the number of total N bins, Pk,l is the normalized occupancy of the lth bin for the kth variable, and Sk,max = ln(N).

Lastly, we evaluated the phase coherence of the population of cells as the cortisol amplitude was varied by calculating the standard deviation of the cell phases again for the variable Per/Cry mRNA (σΦPer/Cry mRNA). This was accomplished by using the fast Fourier transformation (FFT) (24) and then finding the phase angles (ΦPer/Cry mRNA) over the same time span that were used to calculate Rsyn,1 and ρ1(500 h). Phase angles were found relative to 0° of the trigonometric cycle, which can be translated to a constant external time point in the time domain. Because the phase of cortisol is held constant throughout all simulations, these phase angles can be also considered to be relative to the phase of cortisol. Incorporating this frequency domain metric facilitated not only a check on the validity of the two aforementioned metrics, but also the calculation of the phase shift between the entrainer and the population. This procedure was repeated for every amplitude of cortisol tested. The standard deviation of cell phases is expected to change in the opposite direction of Rsyn,k and ρk, taking a value of zero when cells are perfectly entrained and increasing as the cells become desynchronized. We also use the FFT to find the frequency and period of each cell's oscillations. In the results shown below, all of the aforementioned metrics are applied after time t ∼800 h.


A population of 1,000 cells is used in all of in silico experiments performed in this work. Initially, to understand how the stochastic nature of the system impacts the population characteristics independent of entrainment, we investigated the distribution of single cell phases and periods in absence of cortisol. In Fig. 2A we examine single cell phases as the coupling strength (kc) is set to zero. Phases of single cells adopt a uniform distribution as they maintain values through the entire regime from 0 to 2π. Relative to individual cell periods, the stochastic dynamics induce a normally distributed pattern with a mean period of 23.4 h (Fig. 2B).

Fig. 2.

Distribution of single cell phases and periods when no entrainer is present. A: unentrained single cell phases adopt a uniform distribution possessing values through the entire regime from 0 to 2π. B: individual cell periods adopt a normally distributed pattern with mean period equal to 23.4 h.

In Fig. 3, we explore the synchronization that cortisol imposes on the ensemble average of the population. From time t = 0 h until time t = 200 h, cells are independent from cortisol regulation as kc is set to zero. This causes the average to have a very low rhythmicity due to the incoherent rhythms in the single cell level. At time t = 200 h, kc is set to its nonzero value, imposing entrainment by cortisol on the peripheral oscillators. Cortisol affects the ensemble of cells in two ways. First, it stimulates a robust pattern with large amplitude, and second, it imposes its frequency on the peripheral clocks so that the population average has a 24 h period. At time t = 1,000 h, cortisol's effect is again removed by setting kc to zero and the ensemble progressively loses its rhythmicity. Thus, Fig. 3 shows both the synchronization and desynchronization that occur as an entrainer is added and then removed from the system.

Fig. 3.

Cortisol entrainment to Per/Cry mRNA (y1) compartment (1,000 cells). Ensemble average profile of Per/Cry mRNA (y1) before and after the presence of cortisol. Cortisol entrainment (200 < t < 1,000 h) results in a robust expression signal at the population level in contrast with desynchronized states (t < 200 and t > 1,000 h) where the population signal is weaker.

Amplitude-Dependent Synchronization

The impact of cortisol's amplitude on the synchronization of peripheral genes is depicted in Fig. 4. As the amplitude of cortisol increases, moving rightward on the x-axis, the cell synchronization increases as denoted by the gradually increasing Rsyn,1 and ρ1 metrics and the decreasing standard deviation of cell phases. In particular, our model reveals three qualitative regimes of synchronization. In the first regime (I), individual cells are nearly fully desynchronized shown by the small values of the Rsyn,1 and ρ1 entrainment metrics and the high standard deviation of cell phases. As cortisol's amplitude increases to reach its homeostatic value, individual cells gradually become entrained. The second regime (II) is an intermediate state in which some of the individual cells begin to fall in sync. In the third regime (III), individual cell maintain a high degree of synchronization and the standard deviation of cell phases is small. To explore the single cell dynamics when changing the entrainer's amplitude, we further examine single cell phases (ΦPer/Cry mRNA) and periods for different cortisol amplitudes. Figure 5 shows the distribution of cell phases (ΦPer/Cry mRNA) as the amplitude of cortisol decreases. We perform eight in silico experiments for different ratios of cortisol amplitudes. At the homeostatic level of cortisol, where amplitude (amp = 1 − Fmin/Fmax) is 0.85 (Fig. 5A), the population of cells approximates a normal distribution. As cortisol gradually loses its rhythmicity, although the distribution of cell phases retains its normally distributed pattern, the population mean drifts to higher values (Fig. 5, B–D). As the phase of cortisol remains constant, a change in the mean peripheral circadian phase denotes a phase shift relative to the entrainer. At the same time, the standard deviation slowly increases. In Fig. 5, E and F, the population reaches the second regime (II) of synchronization and some cells gradually adopt different phases. Further decreasing cortisol's amplitude the standard deviation increases, culminating in a uniform distribution when cortisol amplitude has become almost constant (regime I; Fig. 5, G and H).

Fig. 4.

Cortisol's amplitude dependent synchronization of peripheral cells (1,000 cells). In small cortisol amplitudes (I) individual cells are desynchronized as it is denoted by the small values of Rsyn,1, ρ1, and the high standard deviation of cell phases (σΦPer/Cry mRNA). As cortisol amplitude approaches its homeostatic value, individual cells pass through an intermediate regime of synchronization (II), where some of the cells are gradually becoming synchronized. Regime III corresponds to the entrained state of the population where the cells are nearly fully synchronized as it is denoted by the high values of Rsyn,1 and ρ1 metrics and the low value of σΦPer/Cry mRNA.

Fig. 5.

Distribution of individual cell phases (ΦPer/Cry mRNA) for 8 cortisol amplitudes (1,000 cells). While in large cortisol amplitudes (amp = 1 − Fmin/Fmax) distribution of single cell phases adopt a normally distributed pattern, as the amplitude of cortisol decreases, the normal distribution becomes gradually uniform. From A to H cortisol amplitude is decreasing. μ and σ are the mean and standard deviation of population phases, respectively.

In Fig. 6 we evaluate single cell periods for different cortisol amplitudes. When cortisol's amplitude results in the third regime of synchronization (III, Fig. 4), the single cell periods are narrowly distributed with mean values equal to the circadian period of the entrainer (Fig. 6, A–D). Further decrease of cortisol's amplitude induces the emergence of a second distribution (Fig. 6E) that indicates that cells are beginning to fall out of sync since they adopt periods different from that of the entrainer. This state is gradually reached by the majority of the cells (Fig. 6, F and G) until it becomes the sole distribution of the population (Fig. 6H) when all of the cells become unentrained.

Fig. 6.

Distribution of individual cell periods for 8 cortisol amplitudes (1,000 cells). In large cortisol amplitudes (amp = 1 − Fmin/Fmax), distribution of cell periods adopt a narrow distribution centered around the circadian period of the entrainer (24 h). As the amplitude decreases, a 2nd distribution of cell periods surges, which ultimately becomes the solely distribution of the population. From A to H cortisol amplitude is decreasing. μ and σ are the mean and standard deviation of population periods, respectively.

Frequency-Dependent Synchronization

In addition to assessing entrainment properties as the amplitude of the entrainer changes, we examine the synchronization of the cell population upon changing the period of the entrainer (Fig. 7). High levels of synchronization are observed only for a certain range of cortisol periods near the period of the free-running peripheral circadian oscillator (23.4 h).

Fig. 7.

Cortisol's frequency dependent synchronization of peripheral cells (1,000 cells). Synchronization as calculated with Rsyn,1 and ρ1 metrics is observed only for entrainer periods relatively close to the individual cell period (23.4 h).

In Figs. 8 and 9 we evaluate the impact of changing cortisol's period on the single cell level. In particular, Fig. 8 shows individual cell phases for six values of entrainer's period. When cortisol's period is very different from that of a single cell (Figs. 8, A and B), cell phases adopt a uniform distribution with high standard deviation. As cortisol's period approaches the period of individual cells, their individual phases gradually converge and the standard deviation decreases (Fig. 8C) until the point where the distribution becomes normally shaped (Fig. 8D). For greater period values, individual cells fall out of sync and the population again adopts a uniform distribution (Fig. 8, E and F).

Fig. 8.

Distribution of individual cell's phases for several cortisol periods (1,000 cells). As entrainer period remain highly different than that of individual cells (23.4 h), the cell phases adopt a uniform distribution. For cortisol periods close to that of individual cells, there is a gradual concentration of phases under a normal distribution. From A to F cortisol period is increasing. μ and σ are the mean and standard deviation of population periods, respectively.

Fig. 9.

Distribution of individual cell periods for several cortisol periods (1,000 cells). As cortisol period approaches or departs from the period of individual cells (23.4 h) we see the rise of a 2nd distribution denoting that cells are gradually becoming synchronized or desynchronized respectively. From A to F cortisol period is increasing. μ and σ are the mean and standard deviation of population periods, respectively.

We further examine individual cell periods (Fig. 9) upon changing cortisol's period. In response to setting cortisol's period to 10 h, individual cell periods are concentrated in a normal distribution (Fig. 9A). As the entrainer period approaches the regime of high synchronization (Fig. 7 regime of high values of Rsyn,1 and ρ1), we observe again the rise of a second distribution (Fig. 9B) that gradually becomes the dominant state of the population (Fig. 9C). This narrow distribution increases in prominence with higher entraining periods (Fig. 9D). Finally, as cortisol's period moves away of the synchronization regime, it induces the appearance of a new distribution (Fig. 9E). As we further increase the period of cortisol, the majority of the cells go to this new distribution until it eventually becomes the sole distribution of the population for a cortisol period equal to 39 h (Fig. 9F).

Circadian Variation of Clock Gene Synchronization

Lastly, we investigated variation in clock gene synchronization throughout a 24 h period. Figure 10A shows single cell profiles of Per/Cry mRNA in response to homeostatic cortisol rhythms that cause peripheral cells to retain a high level of entrainment. For the 24 h period, we calculate the Rsyn,1 synchronization metric for successive time windows of 2 h and observe that synchronization varies in a circadian manner, reaching lower values when Per/Cry mRNA approaches its nadir or zenith levels (Rsyn,1 = 0.64 or 0.13, respectively). Focusing on these regions of low synchronization, Fig. 10, B and C, reveals that single cells have uncorrelated expression profiles with high variances in cell patterns. In contrast, in regions of high synchronization (Fig. 10C), cells adopt coherent profiles with highly correlated dynamics, leading to high values of Rsyn,1.

Fig. 10.

Circadian variation of clock gene synchronization throughout a 24 h period. A: individual cell Per/Cry mRNA expression for homeostatic cortisol rhythms. The synchronization metric (Rsyn,1) has been calculated for consecutive time windows of 2 h and has been placed over each time interval. Different colors denote different windows where the metric has been calculated. B: representation of a small number of cells for the 6–8 h window that indicates the highly uncorrelated profiles of single cells. C: a small number of cell profiles for the regime of maximum desynchronization (12–14 h). D: a small number of cells for a regime of high synchronization (18–20 h). The x- and y-axes of the inset plots are same to these of the main figure A.


The synchronization of endogenous rhythms to external commanding zeitgebers is crucial for the maintenance of biological fitness. In this paper we studied the underlying dynamics of the entrainment of peripheral clocks by a single entrainer (cortisol), focusing on the dependency of cell synchronization on the characteristics of the entraining rhythmic cortisol pattern.

Introducing an external rhythmic force such as cortisol circadian secretion to an oscillating variable such as Per/Cry mRNA (at t = 200 h in Fig. 3) produces a system equivalent to a forced oscillator (68). The effect of cortisol on the peripheral clock is dependent on two characteristics: the frequency and the amplitude of the entrainer. Experiments in chronically jet-lagged and SCN-ablated mice indicate that the loss of glucocorticoid circadian frequency was associated with accelerated tumor growth (28, 29). Also, in humans, cortisol circadian frequency has been found to be disrupted in breast and ovarian cancer patients (70). Heterogeneity in cortisol amplitudes arise due to HPA dysregulation in chronic stress (54). The effective amplitude sensed by the peripheral circadian module in our model is largely dependent on the choice of the parameter kc, which is multiplied with FR(N) in Eq. 6, thus representing the coupling strength between central and peripheral circadian systems (12, 37). In the simulations performed here, this parameter was fixed to a value that was chosen relative to the following qualitative limitations. First, in response to normal cortisol rhythms, the body is at its homeostatic state and peripheral clock gene responses are entrained to cortisol's frequency. Second, when cortisol has completely lost its amplitude variability (FminFmax), the individual oscillators lose synchronization and individual cells adopt random phases and periods. This latter modeling constraint is supported by experiments on SCN-ablated animals (60, 77) that show the loss of circadian rhythmicity in cortisol along with desynchronization in peripheral cells. These limits define a small range of parameters (0.005 < kc < 0.013) inside of which we chose the kc value (kc = 0.009). However, the qualitative results of the model would not be different if we had chosen any other value within that range.

The absence of an entrainer induced a dispersed distribution of single cell phases (Fig. 2A) and a normally distributed pattern in cell periods (Fig. 2B). In particular, unentrained oscillators (t < 200 h and t > 1,000 h, Fig. 3) have a blunted ensemble average similar to the experimentally observed flattening of clock gene rhythmicity when cells are in in vitro unentrained conditions (45), while entrained states lead to a robust rhythm. As was theoretically hypothesized by Balsalobre et al. (9) and experimentally tested by Nagoshi et al. (55) and Welsh et al. (72), this flattening of population-level rhythmicity we observed is an outcome of phase drifting between robust oscillators rather than the loss of individual cellular rhythms. In the in silico experiment of Fig. 4, our model predicts the desynchronization of peripheral cells when cortisol's amplitude fall below a critical value (regime I and II). Given the relationship between circadian rhythms and recovery from critical illness (18, 21), effects of perturbations to cortisol's circadian rhythmicity are particularly interesting because, in addition to its role as a circadian entrainer described here, cortisol plays a critical role as an anti-inflammatory hormone. Furthermore, experimental evidence from severely burned adults indicates that the amplitude of cortisol in stressed patients is diminished, corresponding to our second regime of synchronization (amplitude = 1 − Fmin/Fmax ∼ 0.38) (41). Recently, human endotoxemia experiments showed that there is both an acute increase in cortisol and a time-of-day-independent resetting of clock gene expression in response to endotoxin treatment (38), which further points to an interaction between cortisol, circadian rhythmicity, and inflammation. Godin and Buchman (36) hypothesized that systemic inflammation might lead to uncoupling between biological oscillators and that this loss of interorgan communication may be important in the progression and recovery of critically ill patients. In this context, the loss of central diurnal signals, leading to the desynchronization of peripheral clocks manifested as greater signal regularity in the ensemble level (Fig. 3), can be viewed as a loss of communication between the SCN and peripheral tissues. Thus, a chronically stressed patient might not retain the beneficial properties of circadian rhythms in peripheral tissues (32).

Given the presence of stochastic noise in our system and that there is no intercellular coupling among the peripheral cells in our model, it is not surprising that cells fall out of phase as the rhythmic entraining input decreases. While it is known that the entrainment of SCN neurons depends on intercellular coupling, cells in peripheral tissues display different entrainment dynamics that can generally be explained by their lack of strong coupling (1, 48, 75). As a result, the degree and rapidity with which the individual cells fall out of phase depend upon the amount of stochastic noise, which in our case is determined by the formulation of the CLE algorithm and the reaction propensities, and the amplitude of the entraining signal. However, Fig. 4 reveals an unintuitive dynamic of peripheral cell entrainment, where individual cells maintain relative synchronicity as they drift away from their cortisol-entrained phases. Evidence supporting this finding comes from feeding experiments (78), where temporally restricted feeding induces not only a phase shift in the expression of clock genes but also flattened cortisol rhythmicity. Furthermore, it has been observed that PBMC of cancer patients, under exposure to a blunted systemic cortisol rhythm, maintain a phase-shifted profile (6). Our model predicts this phase-shifting regime when cortisol has slightly lost its amplitude rhythmicity (Fig. 5, A–D, μ value). In this regime, individual cells are normally distributed and the standard deviation is low, denoting the high level of synchronization in the population of cells. For further decreases of cortisol's amplitude, individual cells become desynchronized and adopt phases through the entire range from 0 to 2π (Fig. 5E), culminating in a uniform distribution where cells are fully desynchronized (Fig. 5, G and H).

In our model, we consider the transcriptional effects of glucocorticoids on peripheral clock genes using a fairly simple mathematical model that incorporates negative and positive feedback loops mediated by the PER/CRY protein complex. Additionally, our model takes into account important nontranscriptional effects of glucocorticoids. Recent experiments of Nader et al.(53) and Charmandari et al. (20) show a direct linkage between clock gene rhythmicity and glucocorticoid effects both in HeLa and HCT116 cell lines as well as in PBMC. Specifically, the CLOCK/BMAL1 transcription factor acetylizes the glucocorticoid-receptor complex at a multiple lysine cluster in its hinge region, producing inefficient binding of the complex to GRE in a time-dependent manner. These two signals (CLOCK/BMAL1 production and glucocorticoid/receptor binding to GRE) are in rhythmic synchrony, resulting in decreased tissue sensitivity to cortisol in the morning when cortisol levels are high and increased sensitivity at night when cortisol reaches its nadir. Consequently, a phase shift in CLOCK/BMAL1 clock gene rhythmicity may indirectly affect the impact of glucocorticoids in peripheral tissues through altering the tissue's sensitivity to glucocorticoids.

Our simulations show that even when the cells are phase shifted relative to the entrainer, their individual period remains equal to that of the entrainer (Fig. 6, A–D). In particular, cortisol imposes its period for all the amplitudes of the third regime (III) of Fig. 4. In contrast with the corresponding individual phases (Fig. 5, A–D), the deviation of cell periods for the third regime of entrainment is very small. This is an outcome of entrainment dynamics. Generally, a population of cells (oscillators) that is synchronized by an entrainer retains two characteristics. First, all individual cells adopt the period of the entrainer (Fig. 6, A–D), and second, the individual cell phases are phase-locked relative to the entrainer. This is also true for our model in the regime of synchronization (Figs. 5, A–D, and 6, A–D). However, the stochastic nature of cell dynamics leads different cellular oscillators to adopt slightly different periods and phases (Fig. 2B). Consequently, while the individual cell periods can only slightly deviate from the circadian period, there is a wider regime of individual cell phases that retain the phase locking characteristic. For the case where cortisol has nearly no rhythmicity (Figs. 5H and 6H), individual cell phases are uniformly distributed (Fig. 6H), while individual cell periods adopt a normal distribution due to the absence of an entraining rhythm (Fig. 5H). At this point, it is interesting to note that our model predicts that the presence of even a very small amplitude-entraining signal induces a reduction in the mean period of the population relative to unentrained cells (Fig. 2B vs. Fig. 5H). Changes in population periodicity relative to entrainer concentration has also been observed in synthetic oscillators (67).

As cortisol's amplitude decreases, cells gradually fall out of sync. In particular, we observe the increasing prominence of a second distribution of individual cell periods (Fig. 6E) that further increases (Fig. 6, F and G) and finally becomes the only period distribution of the population (Fig. 6H). This dynamic denotes the “movement” of the population from the regime of synchronization to the regime of desynchronization. This transition happens in a gradual manner. Similar to what was observed for cell phases (Fig. 5E), where we see a gradual dispersion of cell phases that corresponds to the second regime of synchronization (Fig. 4, II), for cell periods we see a transition of cell periods from a regime where they are distributed narrowly around the circadian period of cortisol to a regime where they are distributed with a different mean period.

It is well established that disruption of cortisol circadian periodicity is associated with detrimental outcomes (28). Whether the disruption of cortisol rhythm is a cause or a consequence is an ongoing research problem. In our model we approached this issue by examining the consequence of a disrupted central systemic entrainment on the expression of peripheral circadian clocks.

Figure 7 shows that high levels of synchronization are achieved only for cortisol periods close to the individual cell period (23.4 h). In particular, although individual cell phases adopt a nearly uniform distribution for cortisol periods far from the individual cell frequency (Fig. 8, A and B), as cortisol's period approaches that of individual cells, we see a gradual concentration of cell phases around a mean value and a concurrent reduction of standard deviation (Fig. 8C) until the regime of entrainment (Fig. 8D). The dynamics of our model lie in qualitative accordance with the experimental results of Mondragon-Palomino et al. (52). In their experiment, they used a synthetic genetic oscillator with arabinose as the entraining signal and observed that for entrainer periods near the free synthetic oscillator period, oscillators' phases adopt a dense distribution, whereas for periods away from this regime the synthetic oscillators adopt a “double trough” uniform distribution. In the same experiment they also observed the formation of a bimodal distribution similar to what we observed (Fig. 9, B and E) as the entrainer period approaches and departs from values that induce the entraining regime.

Both Figs. 6 and 9, in accordance with experimental evidence (52), reveal the existence of two distinct dynamic states. Disruption of cortisol rhythms, either by flattening its amplitude or by changing its circadian period, leads the peripheral clock genes to a new nonhomeostatic steady state where individual cells have lost their entraining characteristics respective to both their individual cell periods and individual cell phases. However, our model predicts that even for homeostatic cortisol amplitudes and frequencies, individual cells retain a dynamically varying level of synchronization throughout the day.

There is a circadian time structure associated with many diseases. For example, patients suffering from cardiovascular diseases are at a higher risk of cardiac death (30%) in the morning (23, 39). Peripheral clocks may play a crucial role in this context since their rhythmic output directly regulate circadian generation of heart beat rhythms (42). In addition, the experiments of Balsalobre et al. (8) showed that liver's clock genes profiles after injection of dexamethasone produces different phase-shifting behavior in the population depending on the time that dexamethasone is administered, further illustrating the dynamic nature of peripheral clock gene behavior. Our modeling effort revealed that the level of synchronization, even at states of high entrainment (homeostatic cortisol amplitude and frequency), is dynamic. When calculating the Rsyn,1 metric in consecutive 2 h time windows for a period of 24 h, we see that synchronization fluctuates, with minimum values at Per/Cry mRNA's nadir and zenith (Fig. 10A). In these regions, single cell profiles follow uncorrelated dynamics leading the denominator of the Rsyn,1 formula (mean variance of each oscillator) to increase (Fig. 10, B and C). In contrast, at other times, cells adopt highly coherent profiles (Fig. 10D). This more local analysis reveals that, although in macroscale cells are highly entrained by cortisol's rhythmic pattern, there are certain periods of the day where cell oscillators are desynchronized (nadir and zenith levels; Fig. 10, B and C). This further implies that there are certain times throughout the day when the host may be more vulnerable to external stimuli leading to desynchronization.

In summary, our model is able to predict and evaluate a number of experimental observations, including the necessity of a significant systemic rhythm for the synchronization of the peripheral clocks (Figs. 3 and 4) and that even mild changes in the amplitude pattern of the entrainer can have indirect effects through altering the orchestration of local networks (e.g., acetylization induced by CLOCK/BMAL1) (Figs. 5 and 6). Furthermore, upon changing either the amplitude or frequency of the entrainer, the transition from the unentrained to entrained state and vice versa happens in a gradual manner, as can be observed in the bimodal distributions of Figs. 6 and 9. Our results show that dysregulated cortisol rhythmicity can induce desynchronization of peripheral clock genes, uncoupling them from the systemic circadian network.


P. D. Mavroudis and I. P. Androulakis acknowledge support from National Institute of General Medical Sciences (NIGMS) Grant GM-082974. P. D. Mavroudis, J. D. Scheff, S. E. Calvano, and S. F. Lowry are supported, in part, by NIGMS Grant GM-34695.


No conflicts of interest, financial or otherwise, are declared by the author(s).


Author contributions: P.D.M. and I.P.A. conception and design of research; P.D.M. performed experiments; P.D.M. analyzed data; P.D.M., J.D.S., and I.P.A. interpreted results of experiments; P.D.M. prepared figures; P.D.M. and J.D.S. drafted manuscript; P.D.M., J.D.S., and I.P.A. edited and revised manuscript; P.D.M., J.D.S., S.E.C., S.F.L., and I.P.A. approved final version of manuscript.


  1. 1.
  2. 2.
  3. 3.
  4. 4.
  5. 5.
  6. 6.
  7. 7.
  8. 8.
  9. 9.
  10. 10.
  11. 11.
  12. 12.
  13. 13.
  14. 14.
  15. 15.
  16. 16.
  17. 17.
  18. 18.
  19. 19.
  20. 20.
  21. 21.
  22. 22.
  23. 23.
  24. 24.
  25. 25.
  26. 26.
  27. 27.
  28. 28.
  29. 29.
  30. 30.
  31. 31.
  32. 32.
  33. 33.
  34. 34.
  35. 35.
  36. 36.
  37. 37.
  38. 38.
  39. 39.
  40. 40.
  41. 41.
  42. 42.
  43. 43.
  44. 44.
  45. 45.
  46. 46.
  47. 47.
  48. 48.
  49. 49.
  50. 50.
  51. 51.
  52. 52.
  53. 53.
  54. 54.
  55. 55.
  56. 56.
  57. 57.
  58. 58.
  59. 59.
  60. 60.
  61. 61.
  62. 62.
  63. 63.
  64. 64.
  65. 65.
  66. 66.
  67. 67.
  68. 68.
  69. 69.
  70. 70.
  71. 71.
  72. 72.
  73. 73.
  74. 74.
  75. 75.
  76. 76.
  77. 77.
  78. 78.
View Abstract