Temperature rhythms are systemic cues that efficiently entrain individual clocks in cultured cells and tissue explants (Morf and Schibler, 2013). A possible mechanism involves the post-transcriptional regulation of the cytoplasmic abundance of Clock mRNA by cold-inducible RNA-binding protein (CIRBP) (Morf et al., 2012) and the transcriptional induction of the nuclear Per2 mRNA by heat shock factor 1 (HSF1) (Saini et al., 2012). Experimental evidence shows that cold and heat exposure can activate HSF1 (Cullen and Sarge, 1997) and upregulate the nuclear Per2 transcript (Chappuis et al., 2013; Fischl et al., 2020), so we focus on the HSF1-associated mechanism in our model.

Without loss of generality, temperature rhythms are modeled as square wave functions. Such an approximation enables us to capture the essential oscillation characteristics such as the phase, period, and peak-to-trough amplitude (Eq. 1). Thermosensing and thermoregulation are closely interconnected (Xiao and Xu, 2021). Since temperature regulates nearly all the physiological processes in living organisms, thermosensors are ubiquitous and indispensable in sensing, transducing, and processing external and internal temperature fluctuations (Xiao and Xu, 2021). Thermosensitive transient receptor potential channels (thermo-TRPs) are molecular thermosensors shown to be present in the plasma membranes of both central and peripheral cells and participate in detecting temperature changes (Wang and Siemens, 2015; Kashio, 2021), thus playing a critical role in the entrainment of peripheral clocks (Bromberg et al., 2013; Poletini et al., 2015). Temperature is one of the gating forces to the thermo-TRPs activation (Matta and Ahern, 2007; Yao et al., 2011; Xiao and Xu, 2021), which means that they can be fully activated only by temperature changes that either fall within a specific temperature range or exceed a particular temperature threshold. This implies that both magnitude and changes in temperature are detected. Such temperature sensing process is expected to be essential for temperature entrainment and is assumed to be the first step of the system’s response to a temperature input. Using an indirect response modeling (IRM) approach (Foteinou et al., 2009; Foteinou et al., 2011), we hypothesize two intermediate signals, denoted as TS ₁ and TS ₂ , to represent the essential components involved in this cascade. Eqs 1, 2 describe their kinetics.

By using the average to represent, for simplicity, the magnitude of the input temperature oscillation, Eq. 2 encapsulates the temperature-induced thermo-TRPs activation upon detecting both temperature magnitude and change as well as their inactivation as a result of the homeostatic regulation. This transient activation of thermo-TRPs in the membrane in response to the temperature rhythm might enable the influx of cation ions and the change of membrane fluidity, which then initiate other downstream intracellular events (Vriens et al., 2014; Wang and Siemens, 2015; Guihur et al., 2022). Eq. 3 generalizes the surge (denoted by a Hill coefficient) and the vanishment of the effector of the thermosensor’s activation.

HSF1 is hypothesized to be involved in the mammalian temperature entrainment through the HSF1-induced heat shock response (HSR) pathway (Buhr et al., 2010). When there is no stress in the canonical HSR pathway, heat shock proteins (HSPs) sequester HSF1 monomers and suppress HSF1 transcriptional activity. Upon heat stress, the denaturation of specific thermolabile proteins leads to the dissociation of HSP chaperones from inactive HSF1/HSP complexes. Then, hyperphosphorylation, trimerization, and translocation of HSF1 occur sequentially (Rieger et al., 2005; Reinke et al., 2008; Scheff et al., 2015; Mazaira et al., 2018; Dudziuk et al., 2019; Masser et al., 2020; Guihur et al., 2022), inducing transcriptional activity of HSF1. The active HSF1 trans-activators stimulate the production of HSPs by binding the heat shock elements (HSEs) on the promotor. The resulting excessive HSPs in turn engage in the termination of HSF1 activation with the same sequestration strategy (Gomez-Pastor et al., 2018; Dudziuk et al., 2019). Independent of this molecular chaperone displacement theory, HSF1 can also be activated through different mechanisms, such as temperature-induced intrinsic structural response and other inter-molecular interactions (Hentze et al., 2016; Gomez-Pastor et al., 2018; Mazaira et al., 2018). However, exactly how HSF1 is activated by temperature variations within the physiological range and how HSR participates in the entrainment process remains unclear.

Therefore, by considering the contribution of thermosensors to thermoregulation and focusing on the HSF1-induced excessive HSPs that play an essential role in recovering HSF1 to its inactivated state, we assume that the HSF1 responds to a temperature input (TS ₂ ) as is transduced by the thermosensing process while being inhibited by the excessive HSPs. The mechanisms are expressed via Eq. 4, which delineates the transduction of the temperature signal leading to the HSPs-mediated activation and the homeostatic inactivation of HSF1. Besides this canonical pathway, HSPs are also found to promote the noncanonical dissociation of active trimer HSF1 from the HSEs on target DNA (Kijima et al., 2018). This is described by HSPs inhibiting their HSF1-induced production in Eq. 5 along with a 1^st-order protein degradation of HSPs.

Moreover, studies suggest that the HSF1 activity, such as its nuclear level, phosphorylation degree, and HSE-binding, exhibits circadian rhythm (Reinke et al., 2008; Hirota et al., 2016). In contrast, the Hsf1 mRNA and total HSF1 protein levels remain constant throughout the day (Buhr et al., 2010; Hirota et al., 2016). As such, the transcription dynamics of Hsf1 mRNA are not included. t e m p e r a t u r e   T , °C = 37 + ∆ T 2 ,   0 ≤ Z T < 12 37 − ∆ T 2 ,   12 ≤ Z T < 24 (1) d T S 1 d t = v a c t 0 , t s 1   T a v g + v a c t 1 , t s 1 T − T a v g − v i n a , t s 1   T S 1 (2) d T S 2 d t = v b , t s 2   T S 1 n k b , t s 2 n + T S 1 n − v d , t s 2   T S 2 (3) d a c t H S F 1 d t = v a c t , h s f 1 H S F 1 t o t − a c t H S F 1 i n d H S P 1 + k T   T S 2 K T + T S 2 − v i n a , h s f 1 a c t H S F 1 (4) d i n d H S P d t = v b , h s p a c t H S F 1 / i n d H S P   − v d , h s p i n d H S P (5)

The active HSF1 transduces the circadian oscillations of the temperature signal to the peripheral clock genes by inducing Per2 transcription (Saini et al., 2012). In addition, an active HSF1:CLOCK/BMAL1 interaction on the Per2 promotor was observed after a heat shock pulse and speculated to be involved in the temperature resetting process (Tamaru et al., 2011). Thus, This interaction is hypothesized to stabilize the HSF1:HSE binding and be competitive with the HSP noncanonical function. These mechanisms are included in the HSF1 entraining term in Eq. 6, which models the induction effect of the active HSF1 trans-activator on the Per/Cry expression as an indirect response.

The intrinsic dynamics of the peripheral clock gene network (Eqs 6–12) are modeled based on previous works (Becker-Weimann et al., 2004; Geier et al., 2005; Mavroudis et al., 2012; Mavroudis et al., 2014; Li and Androulakis, 2021). The network consists of transcriptional translational feedback loops, incorporating intertwined positive and negative feedback loops (Becker-Weimann et al., 2004; Geier et al., 2005). Through a positive feedback loop, the nuclear PER/CRY protein indirectly activates the transcription of Bmal1 mRNA. The translation to cytoplasmic BMAL1 protein and its translocation to the nucleus (nucBMAL1) lead to an increase in the production of CLOCK/BMAL1 heterodimer. On the other hand, the nuclear PER/CRY inhibits the stimulation of the transcription of Per/Cry by CLOCK/BMAL1. d P e r / C r y m R N A d t = v 1 b   C L O C K / B M A L 1 + c k 1 b   1 + n u c P E R / C R Y k 1 i p + C L O C K / B M A L 1 + c   1 + k h s f 1   a c t H S F 1 K h s f 1 + a c t H S F 1 + k h s f 1 , c i i n d H S P C L O C K / B M A L 1 − k 1 d   P e r / C r y m R N A (6) d P E R / C R Y d t = k 2 b   P e r / C r y m R N A q − k 2 d   P E R / C R Y − k 2 t   P E R / C R Y + k 3 t   n u c P E R / C R Y (7) d n u c P E R / C R Y d t = k 2 t P E R / C R Y − k 3 t   n u c P E R / C R Y − k 3 d   n u c P E R / C R Y (8) d B m a l 1 m R N A d t = v 4 b   n u c P E R / C R Y r k 4 b r + n u c P E R / C R Y r − k 4 d   B m a l 1 m R N A (9) d B M A L 1 d t = k 5 b   B m a l 1 m R N A − k 5 d   B M A L 1 − k 5 t   B M A L 1 + k 6 t   n u c B M A L 1 (10) d n u c B M A L 1 d t = k 5 t   B M A L 1 − k 6 t   n u c B M A L 1 − k 6 d   n u c B M A L 1 + k 7 a   C L O C K / B M A L 1 − k 6 a   n u c B M A L 1 (11) d C L O C K / B M A L 1 d t = k 6 a   n u c B M A L 1 − k 7 a   C L O C K / B M A L 1 − k 7 d   C L O C K / B M A L 1 (12)

Fifteen newly introduced parameters associated with the temperature signaling cascade and entrainment process were estimated, while the remaining parameters were set to the values determined earlier (Mavroudis et al., 2012). The parameter estimation aimed to satisfy the following criteria at the single-cell level: 1) the concentrations of TS ₁ , TS _2, and activated HSF1 remain positive; 2) when subject to a 12-h warm/12-h cold (W12/C12 37 ± 1.5   ℃ ) temperature oscillation, the clock genes must be entrained to a 24-h period; 3) the phase of cytoplasmic PER/CRY fits experimental data (Saini et al., 2012). In this study, the phase of a component is defined as the periodically stable difference between its peaking time and the onset of the warm phase of the temperature cycle (ZT0). The period is calculated by averaging the time differences between successive peaks during the component’s steady oscillation stage. The estimation also accounts for the fact that the T S signaling cascade would leamplify of the signal and, furthermore, that the level of (activation of) H S F 1 is likely similar to that of the PCGs (Kaneko et al., 2020). The nominal parameter values of the model are summarized in Table 1.

A population of cells and individuals is simulated in the in silico experiments in this study. The cell population is used to understand the ensemble (average) behavior of peripheral cell oscillators, representative of tissue behavior (Yamazaki et al., 2002; Liu et al., 2007; Abraham et al., 2010), while members of the ensemble of individuals are utilized to approximate personalized responses (intra- and inter-individual variability). We generate 1,000 cells for the cell population by sampling the parameters associated with the temperature signaling cascade and peripheral clock gene network. To simulate inter-individual variability, we assume subjects to have individualized temperature sensing ability, which is attained by sampling the associated parameters ( v a c t 0 , t s 1 , v a c t 1 , t s 1 and k b , t s 2 ). Sampling is accomplished using the Sobol method (Rao and Androulakis, 2017; Scherholz et al., 2020).

To explore how the circadian rhythm of the temperature signal affects the state of clock genes at the peripheral tissue level, we evaluate the synchronicity of the cell population. In particular, the synchronization, defined as R _syn,j , is assessed by quantifying cell oscillators’ deviations from mean levels (in terms of the clock gene component j ) and is calculated by dividing the variance of the mean field by the variance of each oscillator (Mavroudis et al., 2012). R s y n , j = y ¯ j 2 − y ¯ j 2 1 N ∑ i = 1 N y j , i 2 − y j , i 2 (13) where y ¯ j = 1 N ∑ i = 1 N y j , i (14) y ¯ j = 1 t T ∑ t = t 1 t T y ¯ j t (15)

In Eq. 13, y j , i is the time-course vector output generated by the model equations, where the indexes j and i represent the component and the cell (N total), respectively. y ¯ j is the time-course vector of averages over the population of N cells (Eq. 14) and ∎ is the time average (timespan from t 1 to t T ) (Eq. 15). A minimum value of 0 of R _syn,j indicates an entirely desynchronized state and a maximum value of 1 demonstrates a full synchronization. In this study, R _syn,j is calculated for the actHSF1 and/or Per/Cry mRNA components for a timespan of 1200 h, except for the validation and circadian disruption experiments in which it is continuously calculated for every 720 h and 24 h, respectively.

The circadian patterns of core body temperature are controlled by the mammalian master clock in the SCN, which is entrained by light/dark cycles (Edery, 2010; Morf and Schibler, 2013; Coiffard et al., 2021). Changes in the phase and period of the light/dark cycles may affect the circadian body temperature rhythm’s phase, period, and amplitude. This is manifested clearly in shift work, a type of circadian disruption, and has been recorded in many in vivo experiments (Knauth et al., 1978; Reinberg et al., 1980; Knauth et al., 1981; Reinberg et al., 1988; Logan and McClung, 2019). Since such mechanisms are absent in cell cultures, we hypothesize that imposed temperature rhythms can play the same zeitgeber role in cultured cells, offering the opportunity to mimic circadian disruption in vitro (Abraham et al., 2010; Stowie et al., 2019).

By applying persistent perturbations to the imposed temperature rhythm, permanent or alternating “temperature shift” are simulated in our study to examine the long-term effect of circadian disruption on peripheral cell synchronization and the inter-individual variability of responses to circadian disruption. As mentioned earlier, the temperature pattern is simulated as a series of Heaviside step functions (Eq. 1). In other words, during the warm period, the temperature is set at its maximum, whereas during the cold period, the temperature is set at its minimum. Combining different patterns results in various zeitgeber schedules. Thus, an alternating shift schedule (ASS) consists of a combination of successive inverted patterns (hereinafter termed “normal” and “reversed”, respectively): n o r m a l   T   T n = T max ,   0 ≤ Z T < 12   T min ,   12 ≤ Z T < 24 (16) r e v e r s e d   T   T r = T min ,   0 ≤ Z T < 12   T max ,   12 ≤ Z T < 24 (17)

The following temperature schedules are examined: 1) 5-day reversed (R) and 2-day normal (N) shift, 2) 3-day reversed and 4-day normal shift, 3) 6-day reversed and 1-day normal shift every 7 days (essentially simulating, in shorthand, 7–5:2/−3:4/−6:1 (R:N) ASS); 4) 20-day reversed, and 8-day normal shift, 5) 12-day reversed and 16-day normal shift every 28 days (in shorthand, a 28–5:2/−3:4 (R:N) ASS). Schedules (i)-(iii) and (iv)-(v) hold the same rotation periods, respectively, whereas (i) and (iv), (ii), and (v) have the same shift windows, respectively. A detailed depiction of T n , T r , and ASSs is provided in Figure A1.

The modeling and computational analyses were implemented in MATLAB R2020b, and the codes are available in GitHub (https://github.com/IPAndroulakis/Temperature-induced-CR).

To simulate the peripheral entrainment subject to an imposed temperature rhythm, we constructed the model with two parts: (a) temperature signaling cascade and (b) peripheral clock gene network (Figure 1; see details in Materials and Methods). In the temperature signaling cascade, we consider a generalized thermosensor-transductor dynamic model and the downstream effector represented by an HSF1-mediated HSR pathway. The peripheral clock gene network is built upon a gene regulatory network model (Jong, 2002) as our earlier works (Mavroudis et al., 2012; Mavroudis et al., 2014; Mavroudis et al., 2015), while temperature rhythms are modeled as earlier described, with a 12 h warm (inactive) period between zeitgeber time (ZT) 0 h–12 h and the cold (active) period between ZT 12 h–24 h (nocturnal animals).

We hypothesize that temperature rhythms drive the phase and period of the peripheral cell-autonomous oscillators. Consequently, all the clock components’ rhythms are expected to be entrained and maintain stable temperature-phase relations (Abraham et al., 2018). Figure 2A reflects the temporal dynamics of representative components (temperature signal (TS ₂ ), HSF1 activity, clock genes’ mRNA, and their cytoplasmic proteins) when cells are entrained to a nominal W12/C12 37 ± 1.5   ℃ temperature rhythm. Figures 2B,C show that the model reproduces the cytoplasmic PER/CRY protein phase and the antiphase relationship between Per/Cry and Bmal1 mRNA as observed experimentally (Brown et al., 2002; Saini et al., 2012).

Our model further captures the entraining and resetting effect of temperature rhythm on the synchronization of a population of cells. Figure 3A depicts the ensemble dynamics of Per/Cry mRNA of a cell population subject to a switch between introduction and removal of temperature rhythm. In the absence of temperature rhythm t < 10   d , the system experiences a constant temperature T = 37   ℃ . The system is then exposed to temperature oscillations and returns to the isothermal condition at t = 65   d . We observe the onset and stabilization of robust oscillations induced by temperature rhythms, which dampen once the temperature is fixed again. The system was entrained to both normal and reversed temperature profiles (Eq. 16). Depending on the temperature profile, the ensembles are entrained to the appropriate phases (Figure 3B), whereas in both cases, the ensembles are entrained to the 24-h temperature period (Figure 3C). Finally, once the temperature rhythms are imposed, the system reaches a high level of synchronization (Figure 3D), and once the zeitgeber is removed, the ensemble of peripheral cells quickly relaxes to the dispersed intrinsic phase and period distributions.

Interestingly, the qualitative differences in the dynamics of synchronization and desynchronization are apparent: the system requires a substantially shorter time to synchronize than desynchronize. This model prediction can be explained by Per2’s property of immediate early gene (Saini et al., 2012), suggesting that during the period following receiving an input signal, Per/Cry in each cell oscillator functions as an immediate early regulator. However, once the driving force is removed, Per/Cry functions solely as a core clock gene, leading the network to relax more gradually via its intrinsic mutual interactions between components and interlocking feedback loops.

Since zeitgeber’s strength is critical for entrainment and synchronization, we evaluated the synchronization response to variations in the temperature amplitude. By varying the peak-to-trough values, i.e., amplitude ( Δ T = 1 − 5   ℃ ), we simulated the increase in zeitgeber’s strength, which leads to an amplification of the a c t H S F 1 ensemble average (Figure 4A) and consequently the Per/Cry mRNA ensemble average (Figure 4B). However, the former is primarily the result of the temperature-induced amplitude amplification of the single-cell a c t H S F 1 oscillation, whereas the latter results from the increased synchronization of cells leading to robust Per/Cry mRNA oscillations. The increase in ensemble average amplitude as a result of entraining the cells is also indicated by the tightening of phase and period distribution (Figures 4C,D). Finally, Figure 4E quantifies the enhancement in peripheral synchronization as a function of the temperature oscillation amplitude. Our prediction that HSF1 activity does not exhibit a circadian rhythm when there is no fluctuation in temperature are consistent with experimental evidence suggesting that HSEs are not transactivated at a constant temperature (Reinke et al., 2008; Tamaru et al., 2011). The model also captures the observations (Refinetti, 2020) that low-amplitude temperature rhythm can partially synchronize the cells (Refinetti, 2020), and cellular oscillators are synchronized more efficiently with increasing temperature amplitudes. However, the complete elimination of temperature oscillations is approximately equivalent to a loss of the entraining signal, as confirmed by the prediction for k h s f 1 = 0 . Furthermore, our results in Figure 4 indicate that an inherent/intrinsic robust oscillation enables the peripheral clock gene network to be capable of gradually/progressively (with more resistance) adapting to external rhythmic forces. In contrast, the temperature signaling cascade, lacking its own inherent rhythm in the absence of input signals, exhibits an immediate and facile adaptation to the same external forces. This difference can also be observed in Figures 4C–E that the actHSF1’s period and phase adapt with high synchrony to all temperature rhythms with varying amplitudes, which is denoted by the narrowly distributed individual phases and periods and the synchronization index ( R s y n ) of 1.

Similarly, by altering the zeitgeber’s average ( 32 − 40   ℃ with fixed amplitude Δ T = 3   ℃ ), which is rationalized by the experimental evidence showing the core body temperature in humans can maintain its circadian rhythm even in plasmodium parasite-induced fever (Lell et al., 2000), we found that increasing the average of temperature enhanced the average of the single-cell a c t H S F 1 oscillation and thus the a c t H S F 1 ensemble average oscillation (Supplementary Figure S1A). However, this did not translate to an apparent change in peripheral synchronization (Supplementary Figure S1B), which is consistent with experimental findings suggesting that elevated temperatures do not lead to loss of peripheral circadian rhythms and only slightly impact the rhythms’ current characteristics (Saini et al., 2012). Combined with the previous results, it indicates that the variation (amplitude), as opposed to the magnitude (average) of temperature oscillations, is more critical for the phase and period entrainment and the synchronization of peripheral cellular oscillators (Saini et al., 2012).

To further elucidate how the period of temperature oscillations influences the synchronization of peripheral clocks, we analyzed the effect of period mismatch between the ensemble average of peripheral cell oscillators (intrinsic period ∼   23.8   h ) and zeitgeber on ensemble entrainment. Following the approach proposed by Schmal et al., 2015, we characterized the ensemble synchronization level and average entrainment phase on the zeitgeber period - zeitgeber amplitude parameter plane (Supplementary Figure S2). As expected, the Arnold tongue indicates that period mismatch leads to desynchronization that can be partially overcome by increasing temperature amplitude; in turn, decreasing period mismatch can compensate for the weak synchronization power of lower temperature amplitudes (Supplementary Figure S2A). These are also reflected by the peripheral clocks being able to adapt to a broader range of period length when subjected to increasing amplitudes of temperature rhythms, along with more pronounced divergences in the ensemble average entrainment phase and its contrast with the zeitgeber’s phase (defined as the onset of the cold period) (Supplementary Figures S2B, C). In particular, our results suggest that the ensemble average entrainment phase is exquisitely sensitive to temperature rhythms, as evidenced by a phase advance (relative to temperature oscillation) in cells entrained to shorter periods (Supplementary Figure S2C). This finding may reflect the concept of overcompensation and is consistent with the experimental result that a modest change in zeitgeber’s period can trigger a significant phase shift, leading to a substantial difference between the phases of the external time cue and internal clock (Duffy et al., 2001).

The previous results established relations between the synchronization behavior and dynamics of peripheral clocks and unperturbed temperature oscillations. In this subsection, we will investigate and analyze the response of a system exposed to alternating temperature schedules (i.e., perturbed temperature oscillations), which we describe as “alternating shift schedules (ASSs).”

Our first in silico experiments mimic the equivalence to “shift work,” where the system is repeatedly exposed to one zeitgeber schedule followed by a different schedule, each lasting for a certain period of time. Figure 5 shows the dynamics of ensemble phases and synchronization level of peripheral clocks subject to different ASSs (Figure A1). Two parameters define the ASS: a) the total period of the minimum non-repeating unit of schedule and b) the internal breakdown. For example, “7–5:2” implies that we consider a 7-day period during which the system is exposed to one oscillatory pattern for 5 days, and then the pattern is inverted for 2 days. Then, the 7-day cycle repeats in the schedule. Similarly, “7–4:3” implies a 7-day period with 4 days in one pattern and 3 days in the inverted. On the other hand, “28–5:2” means that we consider a 28-day period during which the system is exposed to one temperature pattern for 20 days and then an inverted pattern for 8 days (5:2 is the ratio of days in each pattern). Equivalently, “28–3:4” implies a 28-day non-repeating unit, with 12 days in one pattern and 16 days in the inverted. Specifically, in the simulations that follow, the system is first exposed to a “normal” temperature pattern ( T n , Eq. 16) and then subjected to a schedule that alternates between T n  and  T r (Eq. 17) at t = 50 d .

Figure 5A indicates that the ensemble average of peripheral clocks in a system following a 7–5:2 (R:N) schedule eventually adopts oscillatory characteristics resembling the dominating temperature pattern with a phase advance of −12 h while maintaining robust synchronization in the long run. The manifestation is quite different under the 7–4:3 (R:N) schedule (Figure 5B), where the system never manages to attain a clear phase and achieve synchronization. Thus, this schedule induces significant desynchrony. Circadian misalignment and desynchronization of cellular clocks are major contributing factors to chronic disease prevalence. Our observations are consistent with the experimental finding showing the three-shift workers are more vulnerable (probably by having a more dysregulated immune system (Castanon-Cervantes et al., 2010)) than the five-shift workers, resulting in a higher prevalence of common infections among members of the first category (Mohren et al., 2002).

We then considered slowly rotating schedules (i.e., longer time on each pattern before a change happens), generating more extended periods for the minimum non-repeating units. The peripheral ensemble average phase (Per/Cry mRNA peaking time) fluctuates significantly upon exposure to a 28–5:2 (R:N) schedule (Figure 5C) and effectively loses rhythmicity under the 28–3:4 (R:N) schedule (Figure 5D). Compared with the alternating schedules with a higher frequency, slowly rotating schedules allow the system to achieve a higher synchronization level during the adaptation process while leading to an internal destabilization, despite allowing the system to achieve a higher maximal R _syn during the adaptation process. The destabilization can potentially bring about an allostatic accumulation of the consequent internal adverse effects of the shift schedules (Burgess, 2007).

Our early in silico experiments identified temperature amplitude and period as critical factors that can be manipulated to affect the circadian alignment of the system. To analyze the impact of amplitude further, we examined it in conjunction with ASSs. In Figure 6, the dynamics of ensemble average phases and R _syn of peripheral clocks are reevaluated under the 7–5:2 and 7–3:4 schedules as a function of ∆ T in the T r pattern by examining the cases of high ( ∆ T r = 5   ℃ ), nominal ( Δ T r = 3   o C ) or low ( ∆ T r = 1   ℃ ) amplitude. We observe that high T r amplitude strengthens synchronization, while low T r amplitude reduces it for the 7–5:2 ASS (Figure 6A). We predict an increase of R _syn from ∼   0.5 − 0.6 (partial synchronization) under nominal ∆ T r to ∼   0.8 − 0.9 (nearly complete synchronization) under high ∆ T r . In contrast, a decrease to 0 under low ∆ T r is observed, meaning a loss of the entrainment of Per/Cry mRNA occurs when the T r amplitude is decreased substantially. Since the system is inclined to adapt to the characteristics of the most dominant temperature pattern under the 7–5:2 schedule, an increase in the amplitude of dominant zeitgeber pattern improves the system’s adaptation, whereas an amplitude decrease weakens it. However, the response is reversed under the 7–3:4 ASS (Figure 6B), where a low-amplitude T r rhythm enhances the entrainment and synchronization of peripheral clocks. This is likely because the 7–3:4 schedule tends to desynchronize the system. Therefore, reducing the impact of any pattern in the perturbed zeitgeber schedule diminishes the loss of synchronization.

Recent evidence (Rajaratnam et al., 2013; Rao and Androulakis, 2019; Scherholz et al., 2020) suggests that adaptation to zeitgeber characteristics strongly depends on an individual’s ability to sense and be influenced by them. Therefore, we hypothesize that individual peripheral synchronization is driven by, among others, differences in personalized temperature sensing as indicated by v a c t 0 , t s 1 , v a c t 1 , t s 1 and k b , t s 2 for computational simplicity (see details in Materials and Methods). These parameters are theoretically designed to represent an individual’s temperature-induced activation and subsequent effector response, reflecting the sensitivities of temperature sensing, processing, and transduction that are physiologically crucial for the temperature entrainment event. To test this hypothesis, we sampled these parameters and determined the peripheral synchronization of Per/Cry mRNA as a marker for an individual’s internal coherency/circadian adaptation. For the same values of v a c t 0 , t s 1 , v a c t 1 , t s 1 and k b , t s 2 , we estimated the synchronization levels under a “nominal” temperature schedule (normal T, Eq. 16) and three ASSs: 7–6:1 (R:N), 7–5:2, and 7–3:4 (R:N). Maladaptation to alternating zeitgeber has potential adverse effects (McGowan and Coogan, 2013). To assess more accurately how personalized ASSs impact synchronization, we performed a multi-linear regression to determine the likely impact of each individualized factor on the entrainment of a population of peripheral cell clocks: R s y n = β 0 + β 1 v a c t 0 , t s 1 + β 2 v a c t 1 , t s 1 + β 3 k b , t s 2 + ϵ . The results are summarized in Table 2 and indicate that lower sensing sensitivities to temperature magnitude ( v a c t 0 , t s 1 ) generally grant individuals abilities to maintain a higher mean R _syn,ASS . In contrast, higher sensing sensitivities to temperature variation ( v a c t 1 , t s 1 ) only contribute for the 7–3:4 (R:N) ASS (also see Supplementary Figure S3).

We then evaluated the individual’s ability to remain entrained under an ASS as a function of the robustness of its circadian rhythms under the nominal temperature schedule by plotting the R s y n , A S S vs. R s y n , n o m (Figure 7). Interestingly, it is observed that the more robust the synchronization before exposure to ASSs, the higher the adaptation to ASSs. Moreover, confirming our earlier observations, the less vulnerable the schedule (7–6:1), the higher the synchronization, while the most susceptible ASS (7–3:4) exhibits a sharp curve of the relationship between the personalized synchronization levels before and after the ASS. Since a phase-maintained, or said, minimal repeated re-entrainment behavior was hypothesized to serve as a long-term mechanism for a subject to prevent adverse consequences of being forced to persistently readapt to alternating patterns of shift work (McGowan and Coogan, 2013), we postulate that a tolerant individual would have a more dynamically stable and higher level of internal synchronization, which is defined by a higher mean R s y n , A S S over the entire course of adaptation. Based on this, these results qualitatively capture experimental observations noting that subjects tolerant to long-term shift work are likely to have a more robust homeostatic circadian rhythm (Reinberg et al., 1980).