A Mathematical Model to Characterize the Role of Light Adaptation in Mammalian Circadian Clock

In response to a light stimulus, the mammalian circadian clock first dramatically increases the expression of Per1 mRNA, and then drops to a baseline even when light persists. This phenomenon is known as light adaptation, which has been experimentally proven to be related to the CRTC1-SIK1 pathway in suprachiasmatic nucleus (SCN). However, the role of this light adaptation in the circadian rhythm remains to be elucidated. To reveal the in-depth function of light adaptation and the underlying dynamics, we proposed a mathematical model for the CRTC1-SIK1 network and coupled it to a mammalian circadian model. The simulation result proved that the light adaptation is achieved by the self-inhibition of the CRTC1/CREB complex. Also, consistently with experimental observations, this adaptation mechanism can limit the phase response to short-term light stimulus, and it also restricts the rate of the phase shift in a jet lag protocol to avoid overly rapid re-entrainment. More importantly, this light adaptation is predicted to prevent the singularity behavior in the cell population, which represents the abolishment of circadian rhythmicity due to desynchronization of oscillating cells. Furthermore, it has been shown to provide refractoriness to successive stimuli with short gap. Therefore, we concluded that the light adaptation generated by the CRTC1-SIK1 pathway in the SCN provides a robust mechanism, allowing the circadian system to maintain homeostasis in the presence of light perturbations. These results not only give new insights into the dynamics of light adaptation from a computational perspective but also lead us to formulate hypotheses about the related physiological significance.


Supplementary Material 1 Kinetic equations for the CRTC1-SIK1 pathway and for its coupling to the circadian clock
The molecular processes can be described as the following ordinary differential equations. Here the symbols indicate the concentrations of the elements in the model.
The ordinary differential equations for the CRTC1-SIK1 pathway are listed as follows: For the mammalian circadian clock, we listed all the 21 differential equations, eq(12)-eq(32), based on the model which is developed by Henry P. Mirsky et al (Mirsky et al., 2009, supporting information). Experimental results have shown that the nuclear CREB/CRTC1 complex can up-regulate the expressions of Per1 and Per2 mRNA by binding to the CRE elements in their promoters (Jagannath et al., 2013). Therefore, to couple the CRTC1/SIK1 pathway to the mammalian circadian clock, we added a term which is dependent on the concentration of nuclear CREB/CRTC1 complex (CC) in the equations for the time evolutions of Per1 mRNA and Per2 mRNA .
The symbols of variables in the model are listed in Table S1. Table S2 and Table  S3 exhibit the information of parameters (definitions, parameter symbols and parameter values) for the CRTC1-SIK1 model and mammalian circadian clock, respectively.
The main purpose of the CRTC1-SIK1 model is to qualitatively analyze the dynamical mechanism of light adaptation. Therefore, the parameter values for the CRTC1-SIK1 model are chosen based on the qualitative change of experimental observations: (1) The self-inhibition of SIK1 should take effect. (2) The output of the CRTC1-SIK1 module doesn't not disrupt the oscillation of circadian clock.

Parameter analysis
To exhibit the generality of the conclusions in this study, we randomly generated 10000 sets of parameters for the CRTC1-SIK1 module. The range of the parameters are listed in Table S4. The maximal value of each parameter is determined by 100 times of the value in Table S2. For each parameter set, we evaluated the magnitude of light adaptation by the ratio of the stabilized level of CC to the maximal value upon prolonged light exposure.
The magnitude of light adaptation = the stabilized level of CC the maximal value of CC Therefore, smaller ratio implies stronger light adaptation. The simulation results show that 3728 out of 10000 parameter sets yield ratios smaller than 0.8, which implies the occurrence of light adaptation. Among these parameter sets, if Sik1 mRNA is knocked down, the CRTC1-SIK1 model exhibits weakened light adaptation for 3149 parameter sets ( Figure S1). The ratios of light adaptation are weakened more than 20% for 2882 parameter sets. Therefore, this further confirm the conclusion that, the negative feedback loop of SIK1 is necessary for the generation of sufficient light adaptation.
After that, by excluding the sets those destroy the self-oscillation of circadian clock from the 2882 parameter sets, we obtained 2287 parameter sets to test the phase shift response to short-term light pulse. As the circadian clock with CRTC1-SIK1 module exhibits stable oscillation, a 30-min light pulse is applied at 6h after the peak of Per1 mRNA. For most parameter sets, the phases of Per1 mRNA are delayed. By measuring the phase shifts (the time interval between the first peaks of Per1 mRNA after the light pulse and without light pulse) in WT and Sik1 knockdown model, we can observe that the mutant condition exhibits larger phase shifts for 2285 out of 2287 parameter sets ( Figure S2A). The phase shift is increased by 20% in Sik1 knockdown model for 1680 parameter sets. We also found that when the amplitude of the circadian oscillation is small, the phase shifts of WT and Sik1 knockdown model are much larger. Figure S2B-C exhibits an example for this situation, where the difference of phase shifts between WT and Sik1 knockdown model is -3.38h. Thus, the conclusion that the light adaptation generated by CRTC1-SIK1 module prevents the excessive phase shift still holds true for most parameter sets. We also noticed that the phase shift in WT is significantly larger than that in Sik1 mutant model for two parameter sets. Therefore, we checked the magnitude of light adaptation for these parameter sets and found that both WT and Sik1 mutant model exhibit significant light adaptation ( Figure S2D-E). Thus, some other factor may also contribute to the phase shift The 2287 random parameter sets also provided a further test of the result that the light adaptation can reduce the possibility of singularity. We choose the light stimulus at 6.1h after the peak of Per1 mRNA for 7 hours, because it is more likely to generate singularity (Figure 7). The result showed that, 41.8% (957 out of 2287) parameter sets can yield MPDs larger than 5 hours for 100 Sik1 knockdown cells. By contrast, only 23.2% (531 out of 2287) parameter sets generate MPDs larger than 5 hours in 100 WT cells. The situation is similar for the variances of the phases. For 18.45% (422 out of 2287) parameter sets, the differences of the variances are smaller than -5 (the variance in WT minus that in Sik1 knockdown model). However, for only 4.63% (106 out of 2287) parameter sets, the differences of the variances are larger than 5 ( Figure S4).
The above results of the parameter disturbance imply that the regulatory structure of CRTC1-SIK1 network (negative feedback loop) is the essential reason for the light adaptation and preventing the excessive changes of circadian genes.

Initial conditions
Initial conditions for the simulations in this article are listed as follows. This initial condition corresponds to CT0h, which is determined by assuming that Per1 mRNA peaks at CT8.5h according to previous experimental results. CRTC1PC=0

Figure Legend
Figure S1. The differences of the magnitude of light adaptation between WT and Sik1 knockdown model for 3728 parameter sets. 10000 parameter sets are randomly generated for the CRTC1-SIK1 module. 3728 out of 10000 parameter sets yield the magnitude of light adaptation smaller than 0.8, which implies the occurrence of light adaptation. For each parameter set, we used a vertical bar to represent the difference of the magnitude of light adaptation between WT and Sik1 knockdown model (the magnitude of light adaptation in WT minus that in Sik1 knockdown model). It's obvious that the light adaptation is weakened in mutant model for most parameter sets. The maximal synthesis rate of Sik1 mRNA dependent on CREB/CRTC1 complex V1sik1 is decreased by 90% to mimic the knockdown of Sik1.  Table S4, the initial conditions are listed in the section of Initial conditions.   Figure S2A are used to test the the result that the light adaptation can reduce the risk of singularity. The light stimulus is applied at 6.1h after the peak of Per1 mRNA for 7 hours. For each parameter set, we calculated the variances of the phases in 100 slightly