Light entrainment of the SCN circadian clock and implications for personalized alterations of corticosterone rhythms in shift work and jet lag

The suprachiasmatic nucleus (SCN) functions as the central pacemaker aligning physiological and behavioral oscillations to day/night (activity/inactivity) transitions. The light signal entrains the molecular clock of the photo-sensitive ventrolateral (VL) core of the SCN which in turn entrains the dorsomedial (DM) shell via the neurotransmitter vasoactive intestinal polypeptide (VIP). The shell converts the VIP rhythmic signals to circadian oscillations of arginine vasopressin (AVP), which eventually act as a neurotransmitter signal entraining the hypothalamic–pituitary–adrenal (HPA) axis, leading to robust circadian secretion of glucocorticoids. In this work, we discuss a semi-mechanistic mathematical model that reflects the essential hierarchical structure of the photic signal transduction from the SCN to the HPA axis. By incorporating the interactions across the core, the shell, and the HPA axis, we investigate how these coupled systems synchronize leading to robust circadian oscillations. Our model predicts the existence of personalized synchronization strategies that enable the maintenance of homeostatic rhythms while allowing for differential responses to transient and permanent light schedule changes. We simulated different behavioral situations leading to perturbed rhythmicity, performed a detailed computational analysis of the dynamic response of the system under varying light schedules, and determined that (1) significant interindividual diversity and flexibility characterize adaptation to varying light schedules; (2) an individual’s tolerances to jet lag and alternating shift work are positively correlated, while the tolerances to jet lag and transient shift work are negatively correlated, which indicates trade-offs in an individual’s ability to maintain physiological rhythmicity; (3) weak light sensitivity leads to the reduction of circadian flexibility, implying that light therapy can be a potential approach to address shift work and jet lag related disorders. Finally, we developed a map of the impact of the synchronization within the SCN and between the SCN and the HPA axis as it relates to the emergence of circadian flexibility.

www.nature.com/scientificreports/ The SCN can be divided into two distinct sub-regions based on anatomical differences and neurochemical content: a ventrolateral core and a dorsomedial shell region surrounding the core. Neurons of the core contain several neurotransmitters: vasoactive intestinal polypeptide (VIP), calretinin, neurotensin (NT), and gastrinreleasing peptide (GRP), with VIP the most prevalent neurotransmitter. VIP is released rhythmically from the core and induces Per1 and Per2 expression in neurons for both the core and the shell by acting at the receptor VPAC2 5,6 , accomplishing the synchronizing function of the SCN core. The shell is characterized by neurons containing arginine vasopressin (AVP) colocalized with gamma aminobutyric acid (GABA) 7 . Directly below the SCN, the optic chiasm is located which projects the photic input to the SCN. In the core, the non-oscillatory calbindin (CalB)-containing cells are retinorecipient, thereby process light-input directly from the optic chiasm, subsequently conveying photic information to the oscillatory cells in the core. The core in turn transduces the entraining information to the shell through "master synchronizer" neurotransmitter VIP 8,9 . The core projects densely to the shell, while the shell projects sparsely back to the core, enabling photic information transduction indirectly to the shell through the core [10][11][12] .
Entrained by the central pacemaker (SCN), clock genes are also expressed rhythmically in peripheral tissues, primarily mediated via the rhythmic regulation of glucocorticoids (GCs) secreted by the hypothalamic-pituitary-adrenal (HPA) axis 13 . The hypothalamic paraventricular nucleus (PVN) releases a corticotrophin-releasing hormone (CRH) which activates the release of adrenocorticotropic hormone (ACTH) from the pituitary gland and in turn induces the secretion of GCs (cortisol in humans, corticosterone in rodents) from the adrenal glands 14 . Subsequently, GCs suppress upstream regulatory pathways in the PVN and the pituitary gland mediated through glucocorticoid receptors, thus forming a negative feedback loop resulting in the robust oscillation of GCs. This process is synchronized to photic cues (i.e., entrained by light) by the SCN-derived neurotransmitters. AVP exerts negative control over basal plasma corticosterone (CORT) concentrations by inhibiting the secretion of CRH macroscopically, resulting in a precise systemic circadian pattern with GCs peaking at the beginning of the active phase (day for humans, night for rodents) 15,16 . Eventually, the rhythmically produced HPA hormones play a central role in the response to environmental variations, enabling the host to actively re-establish homeostasis 17 . When the HPA axis becomes distorted upon exposure to irregular photic signals, allostatic load accumulation results in GC phase shifts and amplitude dampening which may potentially result in numerous pathologies 18,19 .
The importance of understanding the dynamics of photic entrainment of circadian rhythms motivates the need for characterizing the endogenous pathways, interactions between the SCN and the HPA axis, and allostatic state of the system in adverse environments. Considering the difficulty in measuring fluctuations of the physiological state of the host, mathematical modeling approaches can be particularly helpful in generating the experimentally verifiable hypotheses and deciphering the underlying dynamics to predict clinical outcomes and intervention strategies.
In this paper, we proposed a mathematical model incorporating the entraining effects of photic Zeitgebers on the core, the shell, and eventually the HPA axis. Our model reproduces a series of experimental observations including (1) the phase relation under 12 h/12 h light/dark entrainment; (2) temporal signal transient across the core, the shell, and the HPA axis; (3) a variety in individual's ability to adapt to perturbed schedules such as jet lag and shift work. We predict that the individualized circadian adaption tolerances to jet lag and shift work are highly distinct among the population. Far from being coincidental, we found that an increased rhythmic flexibility results from a higher level of light sensitivity and AVP coupling strength. Thus, our work provides insights into the individualized allostatic adaptation strategies and functional trade-offs by investigating complex physiological entrainment architectures between the SCN and the HPA axis.

Results
General model structure. To simulate synchronization within the SCN and between the SCN and the HPA axis, our model is constructed in two levels. At the SCN level, we use a single-cell gene regulatory network (GRN) model for neurons in both the core and the shell to describe the intracellular circadian dynamics of clock genes and proteins. The coupling within the SCN is accomplished by VIP neurotransmitters released from the core. At the HPA axis level, we consider the corticosterone dynamics which are regulated by the AVP signal released from the shell. The overall schematic of the model is shown in Fig. 1.

Intracellular oscillators in the SCN.
The entrainment dynamics of neurons in the SCN are based on an interlocked circadian oscillator model in mammals [20][21][22] . Neurons in the SCN core and shell share the same transcriptional and translational feedback structure which is described by seven ordinary differential equations. The intracellular network incorporates a positive and a negative feedback loop through which the SCN neurons exhibit self-sustained oscillations ( Fig. 1; details described in "Materials and methods").
Synchronization of the core and the shell. To entrain the biological clock in the SCN to the photic Zeitgeber, light/dark cycles are modeled as a simple step function ("Materials and methods"), which is active for 12 h during the light period between Zeitgeber time (ZT) 0−12 h and inactive during the dark period between ZT 12−24 h . The circadian release of neuropeptides mediates intercellular synchronization of oscillators in the SCN. We assume that the coupling between the core and the shell is accomplished by the neurotransmitter VIP upon PER/CRY protein activity as suggested by previous studies 23,24 . The releasing of VIP activates Per/Cry transcription in both the SCN core and shell, therefore functions as an entrainer for the non-photosensitive shell.
Corticosterone regulation by AVP. The  www.nature.com/scientificreports/ the skeletal structure of the Goodwin oscillator, previously proposed by Sriram et al. 25 and later modified by Rao et al. 26 to encompass the synchronization of the HPA axis activities by light/dark cycle via the SCN in both nocturnal and diurnal species (this study focuses on nocturnal species). In addition, we account for the synchronization effect of AVP on the HPA axis. Experiments demonstrated that AVP micro-infusions to rodents' hypothalamus led to notable downregulation of adrenal corticosterone 27 . Moreover, shell-mediated secretion of AVP is crucial for the maintenance of PVN neuronal firing 28,29 , indicating that the autonomous oscillations of the HPA axis are entrained by the SCN, where AVP plays a crucial role. Therefore, we hypothesize that AVP, as a neurotransmitter output from the SCN, entrains the circadian profile of corticosterone by regulating the degradation of the PVN's output, CRH.
Model calibration. Twelve newly introduced parameters associated with light entraining signal, coupling effects of neurotransmitters VIP and AVP, and dynamic rates of VIP and AVP were estimated independently (denoted by '*' in Table 1), while the remaining parameters were set to earlier determined values 21,26,30 , but slightly adjusted by scaling factors. The main attempt was to limit the parameter adjustment to avoid potential overfitting. This parameter estimation methodology allows to modify the system properties such as phase, period, and amplitude while maintaining the overall relation among the existing parameters. The parameter estimation aimed at satisfying the following criteria: (i) in constant-dark condition, the three compartments must be synchronized (i.e., the periods of different compartments are identical); (ii) when exposed to a 12 h /12 h light/dark entrainer (Zeitgeber), the system must be entrained (i.e., the periods of different compartments equal to 24 h ); (iii) the intrinsic periods of the SCN core and shell should be in qualitative agreement with the experimental observations where the intrinsic period of the shell is less than the core, and both periods should not exceed 24 h ; and (iv) the model should capture the correct phase relation among key oscillatory components: Per/Cry mRNA (core and shell), Bmal1 mRNA (core and shell), VIP, AVP, and CORT. To account for coupling effects across light, the core, the shell, and the HPA axis, we use three coupling parameters v l , v c2 and v coe which represents light sensitivity of the core, the coupling strength of VIP, and the coupling strength of AVP, respectively (see "Materials and methods"). Since these three parameters couple the three independent compartments (the core, the shell, and the HPA axis), we define intrinsic periods of the three compartments as their periods in the absence of coupling (i.e., v l , v c2 , v coe = 0 ). The estimation assigned the core, the shell, and the HPA axis intrinsic periods of 23.8 h , 21.2 h , 24.2 h , respectively. The 8% shorter period of the shell compared to the core qualitatively agrees with experimental findings 31 . Besides, the model reproduced a light-induced phase delay between the core and the shell of about 3 h as observed experimentally 11 . The nominal parameter values and  www.nature.com/scientificreports/ descriptions are included in Table 1. We assume that the nominal parameter set represents the optimum calibration state for our model. Figure 2a depicts the dynamics of seven representative components (clock genes' mRNA, neurotransmitters, and corticosterone) as the model gets entrained by a rhythmic light/dark stimulus ( L12/D12 ). Entrainment drives the SCN and the HPA axis autonomous rhythms to adopt period of the Zeitgeber (light), and as a result, all the components' phases track the light/dark cycles (Fig. 2b). Specifically, the almost anti-phasic relation between Bmal1 and Per/Cry mRNA is reproduced by the model, driven by the activated form of BMAL1 20 . The peaking time of the clock genes' mRNA in the shell exhibits a time delay compared to the core since the expression of PER/CRY in the core drives the secretion of VIP, which ultimately entrains the shell. Experimental evidence suggests Per1 mRNA expression in shell lags that of the core by about an hour, following a light pulse 11 . This feature is qualitatively captured by our model (albeit predicting a slightly greater delay of about 3 h ), suggesting the sequential information transduction within the SCN from the core to the shell. Similarly, VIP and AVP exhibit a phase delay as VIP peaks earlier during the subjective day, while AVP peaks around dusk. In the HPA axis, corticosterone secretion reaches a maximum at the onset of the dark period, which is the beginning of the active phase for nocturnal species.
Since light is the key entrainer, we evaluated the model's response to light by determining the phase response curve (PRC). A 3 h light stimulus nearly 12 times higher than the rhythmic light intensity was introduced at different subjective time, after the system had acclimated to constant darkness (DD). Since in DD there is no light/ dark Zeitgeber time (ZT), we defined the point where corticosterone peaks as circadian time twelve (CT12), following experiments setting the subjective time to twelve at the beginning of the nocturnal animal's active phase. Figure 3 depicts the PRC by measuring the phase changes between the perturbed and unperturbed corticosterone profiles. The resulting PRC is of type I, where both phase advance and phase delay can be produced continuously depending upon the timing of the photic perturbation. Phase delay is observed when the light is introduced at the descending phase of CORT (between CT 11 and CT 22). In experimental observations, a typical PRC exhibits phase delay in the early subjective night and phase advance in the late subjective night, with little phase-shifting occurring during the subjective day. Hence, the subjective day portion of the PRC is often referred to as the "dead zone" 32 . In our model, light stimulus causes phase delay when introduced from CT 11 to CT 22, while a slight phase advance is predicted when light pulses are introduced during the early subjective day and late subjective night. Therefore, our model predicts a large delay-to-advance ratio, which is consistent with experimental observations 33 .
Model synchronization dynamics. An essential and innovative feature of the model is the existence of coupling links across different independent compartments. Understanding the synchronization mechanisms of the system as well as the role of the three coupling parameters is critical to revealing the underlying physiological relations. In our work, we focus on VIP, AVP, and CORT as the key outputs to represent the three compartments: the core, the shell, and the HPA axis, respectively. It is worth mentioning that VIP and AVP are not only outputs of their corresponding compartments but also entrainers to the subsequent compartment. The computational "actograms" of VIP, AVP, and CORT show the synchronization procedure of our system (Fig. 4). As shown in the beginning stage of the actograms, the SCN core and shell, as well as the HPA axis, oscillate independently www.nature.com/scientificreports/ with their intrinsic periods with the absence of the three coupling parameters (i.e., v l = 0, v c2 = 0, v coe = 0 ). As the uncoupled system gets to steady-state, VIP coupling strength is firstly restored to its nominal value, resulting in the synchronization of AVP indicated by the parallel phase patterns of VIP and AVP. Subsequently, AVP coupling strength is restored to its nominal value after the new steady-state has been reached, resulting in the synchronization of CORT. Finally, with the re-establishment of light sensitivity, VIP, AVP, and CORT get entrained by the light/dark cycle, indicated by the vertical phase patterns with a 24 h period oscillations. These results clearly illustrate the synchronization dynamics of the system, demonstrating the regulation and entrainment of the inner compartments by the external Zeitgeber. Next, we investigated the entrainment range differences (i.e., the period mismatch between the subsystems and the Zeitgeber that leads to entrainment) among the three compartments, and more importantly, how the values of the coupling parameters affect the entrainment. To evaluate the entrainment range of the system, we adopted the protocol suggested by Schmal et al. 34 , whereby the entrainment status is characterized on the "Zeitgeber period -Zeitgeber photoperiods" parameter plane. Within the entrainment region, Zeitgeber enforces its period in the circadian compartments so that the compartments get entrained. The resulting entrainment zones for AVP, VIP, and CORT are shaped in the form of skewed "Arnold onions", Fig. 5a. For each Arnold onion, the bottom and upper points are affected by the free-running periods in constant-dark (DD) and constantlight (LL), respectively. The tips (photoperiod χ = 0 or 1 ) of the onion pointing to Zeitgeber period-axis are given by periods of the system under DD ( τ DD ) and LL ( τ LL ), respectively. We determined that the Arnold onions are skewed toward the left since the system shortens its period with increasing light intensity. The widest entrainment range is found near the equinoctial photoperiod (i.e., Zeitgeber photoperiod χ = 0.5 ). As intuitively expected, the entrainment range diminishes moving from the core to the shell and eventually to the periphery (i.e., VIP > AVP > CORT ). Area difference between the Arnold onions of two neighboring compartments  Table 1. All components oscillate with a period of 24 h once the system gets entrained. Bmal1 mRNA oscillates almost anti-phasic to Per/Cry mRNA. mRNAs in the shell oscillate with a phase delay of 3 h compared to mRNAs in the core. White/dark patterns denote light/dark settings. (b) Experimentally determined phase relation between the representative components (left) and model predictions (right) are consistent 63 www.nature.com/scientificreports/ represents the entraining ability of one over the other, where a smaller area difference corresponds to a better entraining ability. We determined the area differences between VIP and AVP, AVP and CORT as a function of the corresponding coupling strength values. As shown in Fig. 5b, as the coupling strength of VIP (i.e., v c2 ) increases, the entraining ability of VIP over the shell increases, while as the coupling strength of AVP (i.e., v coe ) increases, the entraining ability of AVP over the HPA axis firstly increase and then decreases. The non-monotonic profile of the coupling effect might be due to the nonlinear nature of the circadian oscillators. More details related to the effects of coupling strengths will be discussed later.
In-silico population and their response to perturbed light schedules. As the primary outcome of the HPA axis, the CORT circadian phenotype should be maintained within strict physiological bounds to support downstream homeostatic mechanisms. However, significant individual variability exists within these physiological bounds, and this difference potentially contributes to personalized synchronization strategies that allow for differential responses to transient and permanent light schedule changes. We hypothesize that individual synchronization is driven by, among others, the inter-individual differences in coupling parameters within the SCN and between the SCN and the HPA axis as reflected in parameters v l , v c2 and v coe . These parameters reflect physiologically crucial light sensitivity and the coupling strengths of VIP and AVP. The interplay between these values results in synchronization heterogeneity in a compensatory manner to maintain CORT circadian rhythms within reasonable bounds. A virtual population was created using Sobol sampling 35,36 where the distribution of three coupling parameters was identified such that the circadian profiles of simulated CORT lie within ± 20% of the nominal CORT profile (Fig. 6a). This sampling approach helps to determine the subspace of three parameters ( v l , v c2 and v coe ) that corresponds to homeostatic CORT rhythms. We find that v l and v c2 have relatively wide acceptable distribution ranges while the range of v coe values is narrow (Fig. 6b). This selection result is consistent with our observation in Fig. 5b, since the increase of v l and v c2 leads to a stronger entraining effect for the SCN core and shell, while for v coe , only a narrow range of its value induces strong coupling effects of the HPA axis. Jet lag and shift work increase the incidence of chronic pathologies, such as diabetes, depression, and cancer, driven by the temporal misalignment of the endogenous circadian oscillations with the external Zeitgebers. Previous studies 26,37 have shown pronounced individual differences in jet lag and shift work tolerance, while the correlation between jet lag tolerance and shift work tolerance is less studied. To better understand the interindividual differences in response to perturbed light schedules and the underlying trade-offs between the ability to adapt to jet lag and shift work, we evaluated the capacity of the system to respond to perturbations in light schedules reflecting various forms of jet lag and shift work.
To characterize the response of the model to jet lag, we first conducted a jet lag simulation under the nominal parameter set. A 6 h advanced light phase shift (simulating eastward traveling) is introduced once the host has been acclimated to a regular light schedule. The change in light schedule triggered a gradual shift of the peaking time of VIP, AVP, and CORT towards the new light regimen (Fig. 7). Eventually, the system reaches a new steadystate characterized by stable 24 h oscillation with a 6 h phase advance. The phase shift is most rapid at the early stages of the transitions and then decreases exponentially (Fig. 7b). Our model predictions are consistent with experimental observations 38 . We consider a component to be resynchronized when its phase shift is larger than 5.8 h . By testing the resynchronization time of 5 components (Per/Cry mRNA in the core, VIP, Per/Cry mRNA in the shell, AVP, CORT) among all the selected individuals, we observed an increase in resynchronization time as the entraining signal transitions from upstream to downstream, resulting in a minimum resynchronization www.nature.com/scientificreports/ time of Per/Cry mRNA in the core and a maximum resynchronization time of CORT in the HPA axis (Fig. 8). Analogous findings were reported by Kiessling et al. 39 , where a strong heterogeneity in resynchronization time under jet lag schedule was found between different organs and tissues following an entraining signal transduction order from the SCN to the adrenal, and then to the peripheral cells. Our model qualitatively captured this adaptation heterogeneity between different compartments due to the hierarchical structure of the physiological oscillatory system. Adrenal glucocorticoids play a key role in circadian re-entrainment under jet lag and shift work. Manipulation of the phase-shifting of adrenal glucocorticoid circadian rhythms regulates the speed of behavioral re-entrainment. Therefore, we focused our characterization on the CORT when investigating inter-individual responses to jet lag and shift work. Figure 9a records the resynchronization time of CORT for all the selected individuals under a 12 h jet lag perturbation. As light sensitivity and APV coupling strength ( v l and v coe ) increase, the time to resynchronization decreases, indicating that individuals who are more sensitive to light and AVP bear a better jet lag adaptation ability. The inter-individual diversity in jet lag resynchronization time is distinctive, with more than a threefold difference among the simulated population. Another intriguing feature is the direction of jet lag. Experimental studies have shown that eastward (phase advance) travel is typically worse than westward (phase delay) travel for human beings, due to the fact that the average intrinsic period of the circadian clock in humans is longer than 24 h 40 . Similarly, the simulated system shows a stronger tendency to exhibit phase advances than phase delays since its intrinsic period is shorter than 24 h . To quantify the ability of each individual to exhibit phase delay in jet lag, we recorded the CORT shifting direction for each individual under 23 shifting schedules (+1 h , +2 h , …, +12 h , −1 h , −2 h , …, −11 h; " + " means phase advance schedule while " − " means phase delay schedule). We define the phase delay index as the number of the schedules under which CORT exhibits phase delay shifting (rightward in actograms). Our results show that the maximum value of phase delay index is 10, which is smaller than the overall phase delay schedule numbers (which is 11), meaning all the virtual individuals have a stronger tendency to shift advancingly. The level of this tendency varies among the population and is mainly affected by the light sensitivity and AVP coupling strength (Fig. 9b). Moreover, an individual's ability to exhibit phase delay shift is positively related to its resynchronization rate under jet lag, suggesting the existence of the inner correlation between the resynchronization direction and adaptation rate.
To understand both long-term and short-term effects of shift work, we simulated two types of shift work schedules: transient shift work and permanent alternating shift work. The former denotes a temporary perturbation 26 in which an 8-day transient inversion in the light/dark schedule was introduced. The maximum phase difference between perturbed and unperturbed CORT oscillations was measured as a gauge of tolerance to short term shift work (a larger phase difference represents a lower level of tolerance). Once again, our model predicts a parametric dependence of the transient shift work tolerance (Fig. 10). Specifically, a higher level of light sensitivity and AVP coupling strength leads to a lower level of transient shift work tolerance. Combining these with the jet lag findings, our model predicts that those individuals whose circadian clock adapts more efficiently to jet lag are more susceptible to transient shift work. This trade-off was also indicated in our earlier work 26 .
Following this, we sought to determine the effects of alternating shift work patterns. We examined the system response to rotating light/dark cycles with 2 days inversion every 7 days (essentially simulating a 7-day schedule: 5 days normal and 2 days reversed shift). Following exposure to alternating shift work schedule, a distinct www.nature.com/scientificreports/ inter-individual variation was observed (Fig. 11a,b). Based on whether an individual can maintain their period within 24 h ± 0.05 h or not, the population was assigned in one of the two groups. Figures 11a,b represent the actograms upon exposure to the fast-rotating schedule for individuals who failed (group A) and succeeded (group B) to maintain a 24 h oscillation. Individuals in group A appear to manifest a "free running" status given that their period is smaller than 24 h . Individuals in group B show a greater propensity to maintain the normal light/ dark cycle. While individuals in group B adapt, their phases are slightly shifted compared to their initial phases. The phase-locked behavior in alternating shift work of group B was also reported in experimental observations and was hypothesized to serve as a mechanism in the long term to prevent detrimental consequences of oscillations having to repeatedly resynchronize to alternating patterns of shift work 41 . Therefore, we hypothesize that individuals belonging to group B can tolerate the fast-rotating schedule better than individuals in group A. The parametric space distributions of the two groups are shown in Fig. 11c  We define the regulatory plasticity of individuals from group B to be higher than those from group A, since individuals from group B have a shorter resynchronization time and can adapt to the external shifting directions more easily. This plasticity also accounts for the trade-offs between the tolerance to transient shift work and the alternating shift work, since a lower tolerance of the former and a higher tolerance of the latter both result from the emergence of a higher regulatory plasticity.

Discussion
Entrainment of the internal physiological clock by the external environment is essential to maintaining homeostasis. As the central effector that controls organismal behavior, the HPA axis conveys entrainment signals from the master pacemaker (SCN) to peripheral cells. When light/dark schedules become irregular, the SCN output to the HPA axis becomes disturbed, resulting in glucocorticoids rhythm, phase and period disruptions which have been linked to numerous pathologies 18 . The complexity of the circadian timing in our body requires a systematic approach to investigate a multitude of interactions. While previous mathematical models have incorporated the entraining effects of light on the HPA axis 21,22,26,36 , to our knowledge, this study is one of the first attempts to There is a monotonically increasing balance between v l and v coe that as v l increases, v coe increases in a relatively narrow range. The main, left, and above views show the existence of thresholds for v l and v c2 such that they generate a robust coupled system that matches the experimental rhythms of CORT. www.nature.com/scientificreports/ construct an integrated multi-compartment model combing the heterogeneous light entrainment of the SCN core and shell, and the subsequent entrainment of the HPA activity through efferent output signals from the SCN. Our model explored the temporal dynamic and light responses for the three compartments: the SCN core and shell, and the HPA axis. We observed the existence of temporal information transfer across these compartments under regular (Fig. 2) and perturbed light/dark schedules (Figs. 7 and 8). Under a regular schedule, the components' peaking phases follow a spatially and temporally ordered pattern, as entrainment signals transduce from the environment (light input) through the core to the shell, and eventually the HPA axis. Between the core and the shell, there is a time lag for the light-induced Per expression to spread since neurons in the core project to the shell, but not the reverse 12 suggesting that temporal information is transmitted from the core to the shell 42 . In agreement with jet lag experiments 39 , our model predicts considerable transient desynchronization of the system, where the adaptation time-delay varies moving from the SCN core to the shell and eventually the HPA axis. The simulated alternating circadian phase shift also causes desynchronization among different oscillators (Fig. 11a,b). These misalignments among different compartments are likely to cause detrimental effects related to jet lag and shift work.
We hypothesize that the differences in coupling parameters values (defined by v l , v c2 and v coe ) define the inter-individual regulatory plasticity. Our results indicate the existence of a balance between the three regulatory mechanisms to match the desired corticosterone levels, Fig. 6b. While AVP coupling strength lies within tight bounds, as indicated by the narrow distribution of v coe values, both light sensitivity and VIP coupling strength show high variability among individuals, with the establishment of their threshold values as indicated in Fig. 6b ( v l , > 0.1 , v c2 > 0.8 ). Such high diversity in light sensitivity has also been experimentally observed with more than a 50-fold difference across individuals 43 . In addition to stimulating the visual system, light elicits responses www.nature.com/scientificreports/ related to physiological and neurobehavioral functions, known as the non-image forming (NIF) effects of the light 44 . Light sensitivity has a critical effect on NIF functions, such as the circadian entrainment, and shifts the timing of circadian rhythms. Previous experiments have found a potential consistency in the tolerance to disrupted schedules, as both male and young groups were reported to have an improved ability to cope with shift work and jet lag compared to female and aging groups 45,46 . Our model suggests that the physiological cause of this consistency may be related to the coupling strengths. A lower light sensitivity has been reported among female and aging groups 47,48 . On the other hand, earlier studies have shown that chronic jet lag increases mortality in aged mice, and aging groups have difficulty recovering from simulated jet lag 49,50 , while social jet lag is more detrimental to female in their cognitive ability compared to male 51 , and female may struggle more with the effects of jet lag or night shifts than male 52 . Our model predictions show that individuals with lower light sensitivity are more vulnerable to shift work and take longer to recover from jet lag, implying that the attenuation of circadian entrainment in jet lag and shift work among aging or female groups may be caused by their lower light sensitivity. In other words, inter-individual differences in light sensitivity may explain differential vulnerability to circadian disruption and subsequent impact on human health. Based on this, our model suggests that light therapy could be used to improve circadian entrainment under perturbed schedule among vulnerable population. Apart from light sensitivity, VIP coupling strength is also essential in maintaining intracellular and behavioral rhythmicity in the SCN. Studies have shown that mice bearing a null mutation of the VIP receptors cannot maintain normal circadian rhythms. Moreover, mice lacking VIP in the SCN have abnormal circadian activity and reduced synchrony among neurons 53 . Although VIP is also expressed through out other parts of the mammalian brain such as cortex, retina, and superior colliculus, a recent study by Mazuski et al. shows that VIP from the SCN is crucial in sustaining daily rhythms synchronization. Therefore, VIP has been hypothesized to be necessary for the intra-SCN connectivity and generation of the SCN output to hypothalamic areas 54 . In accord with the implications of these findings, our model predicts that keeping VIP coupling strength above a minimum value inside the SCN is crucial for maintaining the rhythmicity throughout the system. The insufficiency of VIP coupling strength can be due to a lack of VIP neurons or VIP receptors. Furthermore, our model also predicts that the role of VIP coupling strength plays in perturbed circadian schedules is negligible since there is almost no variation in response of the system as v c2 changes (Figs. 9, 10 and 11). Earlier work has shown that AVP released by neurons of the SCN strongly affects the depressive disorder among patients. Once over-expressed or overreleased, AVP may contribute to hyper-anxiety and disturbed rhythmicity 55 . Our model incorporates the effect of AVP levels and consistent with the experiments, our results suggest a narrow range for AVP coupling strength to maintain normal behavior homeostasis. More importantly, we predict that, within this range, AVP coupling strength increases as light sensitivity increases. Both light sensitivity and AVP coupling strength significantly affect an individual's flexibility when exposed to perturbed schedules. We explored the system response under non-rhythmic constant darkness condition (Fig. 3). We tested the PRC as a probe for understanding the circadian regulatory mechanism of the system. As a type I PRC, the system displays relatively small phase shifts with a continuous transition between phase delays and advances. The phase-resetting indicated by the PRC curve compensates for the fact that the free-running period of the system is not equal to 24 h . In other words, light can reset the clock to equalize the period of the entrained oscillator to the period of its cycle. With the verified entrainment ability of the system, we further investigated the difference in the ability to get entrained by measuring the area of Arnold onions for three compartments (Fig. 5). The left-skewed Arnold onions indicate that constant light induces a decrease in the intrinsic period of the system. This prediction does not abide by Aschoff 's first rule in which higher light intensity is associated with longer circadian periods in nocturnal organisms. This disobedience may be due to the missing of one pathway in our model where constant light lengthens the nocturnal circadian period by elevating levels of mPER2 protein and www.nature.com/scientificreports/ constitutively enhance the phase-delaying limb of the molecular oscillator 56 . Our model also predicts a coupling parametric dependence of entrainment ranges. Specifically, the entrainment range of VIP and AVP gets increased as their corresponding coupling sensitivity v l and v c2 gets larger, while the entrainment range of CORT firstly increases then decreases as v coe increases with the maximal achieved at v coe ≈ 0.8 . The parametric dependencies shown in Fig. 5b accounts for the parameter sets distribution after selection shown in Fig. 6b, since the prerequisite for individuals to fit in the nominal CORT phenotype is to establish a well coupled and entrained multi-compartment system. Finally, we explored inter-individual adaptation mechanisms in perturbed schedules. Studies show that during long-term shift work, complete adjustment is seldom achieved 38 . By categorizing the simulated population into a susceptible group (group A) and a resistant group (group B), our results also indicate the lack of entertainment for both group A and B in alternating shift work schedules. Individuals from group A oscillate with an unstable period smaller than 24 h , while although individuals from group B establish a stable 24 h period, they also fail to completely track the external rotating schedule. By tracking the previous steady light/dark schedule, group B establishes a mechanism in the long term to prevent the maladaptive consequences of rhythms having to constantly re-entrain to a rotating pattern of shift work. Thus, though both groups A and B failed to be entrained by the rotating schedule, resistant individuals exhibit a less pronounced de-synchronization. Because the susceptible individuals are also less prone to exhibit phase delay shifts in jet lag (Fig. 9b), we hypothesize that their higher susceptibility to alternating shift work mainly results from their stronger propensity to shift in the same direction as their intrinsic period (i.e., phase advance). This tendency to align their rhythm to their free-running pattern also contributes to a lesser phase shift in transient shift work of group A. Overall, our model predicts that individuals with higher light sensitivity tend to have enhanced ability to cope with jet lag and alternating www.nature.com/scientificreports/ shift work, but not transient shift work. The mathematical model of the SCN and the HPA axis that we used in this study has several limitations. For instance, the model failed to reproduce the "dead zone" part of the phase response curve due to the fact that the system didn't include the phase-dependent sensitivity (gating) effects 30 . Moreover, it is well understood that apart from its classic circadian oscillation, the HPA axis also exhibits prominent ultradian pattern of discrete pulsatile release of glucocorticoids. Apart from that, glucocorticoids have also been shown to have feedback effects on central circadian clocks, adding complexity to the circadian dynamics of glucocorticoids. Nevertheless, studies have shown that many essential properties of glucocorticoid rhythms can be explained using circadian limit cycle oscillators of the HPA axis 25,26,36,[57][58][59] . Therefore, with future modifications of the model that gradually incorporate the mentioned properties, we expect that the key features of our simulation results might be preserved.
In summary, our model describes the molecular rhythms observed in the SCN and the HPA axis and describes the entrainment pathways at multiple levels. Focusing on the coupling parameters between the three compartments (SCN core, SCN shell, and the HPA axis), we investigated the likely drivers of individualized susceptibility under different circadian perturbations, proposing a better interpretation of the relation between the genotypic diversity and reactivity of a population composed of phenotypically similar individuals. Since the model incorporates the essential structure and processes of the circadian system, future work includes the multi-cell scale models describing the synchronization of neurons in the SCN; the introduction of other Zeitgebers such as sleep and metabolism; the inclusion of more parameters to account for inter-individual variety.

Materials and methods
SCN compartment. In the clock network of both the core and shell, the nuclear translocated PER/CRY proteins (y 3 and y 11 ) downregulate their synthesis by stimulating Bmal1 gene transcription (positive feedback) and inhibiting the activity of the CLOCK/BMAL1 heterocomplex (negative feedback). The heterocomplex CLOCK/BMAL1 (y 7 and y 15 ) functions as an activator for the transcription of Per/Cry genes (y 1 and y 8 ). In the positive feedback loop, after the transcription of Per/Cry mRNA, the expressed PER/CRY proteins (y 2 and y 10 ) translocate to the nucleus (y 3 and y 11 ). The PER/CRY proteins in the nucleus activate Bmal1 mRNA (y 4 and y 12 ) transcription, which further increases the expression of CLOCK/BMAL1 heterodimer (y 7 and y 15 ) after the translation to BMAl1 protein (y 5 and y 13 ) and the protein translocation to the nucleus (y 6 and y 14 ). In the negative feedback loop, PER/CRY proteins in the nucleus shut off the transcriptional activity of the CLOCK/BMAL1 heterocomplex. The positive and negative feedback structures essentially result in the autonomous oscillators. The behavior of the clock components in the core and the shell is described by Eqs. (2)(3)(4)(5)(6)(7)(8) and (10)(11)(12)(13)(14)(15)(16), respectively.
Light/dark cycles are modeled as a L12/D12 step function (Eq. 1). We introduce the light effect by an additional saturation of photoreceptors term in Eq. (2) encompassing the light-induced transcriptional activation of Per/Cry mRNA in the light-sensitive core. This additive term accounts for the independence of photic-induced transcription from CLOCK/BMAL1 driven transcription 60 . The Michaelis-Menten equation is used to account for the saturation of photoreceptors 61 .
To ensure that the neurotransmitter is released quickly after Per/Cry mRNA transcription, we assume the neurotransmitters are released upon cytosolic PER/CRY protein activity as suggested by other modeling works 23,24 . Correspondingly, VIP is induced by PER/CRY protein in the core (Eq. 9), while AVP is induced by PER/CRY protein in the shell (Eq. 17). The coupling between neurons in the core and the shell is accomplished by the neurotransmitter VIP, which leads to Per/Cry mRNA transcription in both the core and the shell (Eq. 2 and 10). In the VIP coupling terms, Hill coefficients are used to control the steepness of the Per/Cry promoter feedback loop.

Data availability
Simulations and analysis were performed in the MATLAB_R2020b (The MathWorks, http:// www. mathw orks. com). Additional data and materials generated in this work is available in a public GitHub repository at https:// github. com/ IPAnd roula kis/ Light-SCN-HPA.