Complementary phase responses via functional differentiation of dual negative feedback loops

Multiple feedback loops are often found in gene regulations for various cellular functions. In mammalian circadian clocks, oscillations of Period1 (Per1) and Period2 (Per2) expression are caused by interacting negative feedback loops (NFLs) whose protein products with similar molecular functions repress each other. However, Per1 expression peaks earlier than Per2 in the pacemaker tissue, raising the question of whether the peak time difference reflects their different dynamical functions. Here, we address this question by analyzing phase responses of the circadian clock caused by light-induced transcription of both Per1 and Per2 mRNAs. Through mathematical analyses of dual NFLs, we show that phase advance is mainly driven by light inputs to the repressor with an earlier expression peak as Per1, whereas phase delay is driven by the other repressor with a later peak as Per2. Due to the complementary contributions to phase responses, the ratio of light-induced transcription rates between Per1 and Per2 determines the magnitude and direction of phase shifts at each time of day. Specifically, stronger Per1 light induction than Per2 results in a phase response curve (PRC) with a larger phase advance zone than delay zone as observed in rats and hamsters, whereas stronger Per2 induction causes a larger delay zone as observed in mice. Furthermore, the ratio of light-induced transcription rates required for entrainment is determined by the relation between the circadian and light-dark periods. Namely, if the autonomous period of a circadian clock is longer than the light-dark period, a larger light-induced transcription rate of Per1 than Per2 is required for entrainment, and vice versa. In short, the time difference between Per1 and Per2 expression peaks can differentiate their dynamical functions. The resultant complementary contributions to phase responses can determine entrainability of the circadian clock to the light-dark cycle.


Introduction
Complex gene regulatory networks are responsible for diverse cellular functions, such as transcriptional switches, adaptation, noise filtering, and genetic oscillations [1][2][3]. These gene regulatory networks often include multiple feedback regulations. Naturally, redundancy in multiple feedback regulations confers robustness on the system, but understanding whether and how redundant feedbacks acquire different functions is less resolved. Here we address this question with the interacting negative feedback loops (NFLs) in the mammalian circadian clock system.
Almost all known organisms possess circadian clocks that can set a subjective time for an individual in a constant environment and regulate behavioral and physiological rhythms. A fundamental characteristic of these circadian clocks is their ability to entrain to zeitgebers, including light-dark (LD) cycles. Entrainment properties of circadian clocks have been studied by measuring responses to brief light signals under conditions of constant darkness. In mammals including human, a light signal administered at subjective dawn advances the clock, whereas one at subjective dusk or night delays the clock [4,5]. Such advance and delay of the circadian clock are termed phase shifts. Furthermore, plotting phase shifts as a function of the time of light administration results in a phase response curve (PRC), which predicts how the circadian clock entrains to LD cycles [4,6,7].
Zeitgebers, including light signals, shift the phase of the circadian clock by affecting its molecular machinery. Negative feedback regulations of circadian clock genes with time delays generate oscillations in their gene expression with a nearly 24-hour period. In mammals, the circadian clock genes Period1 (Per1) and Period2 (Per2) encode the transcriptional repressors of their own promoters ( Fig 1A). Transcription of Per1 and Per2 mRNAs is induced by the CLOCK-BMAL1 complex through its binding to the E-box element (CACGTG) in their promoters. Each of the PER1 and PER2 proteins forms a complex with co-repressor CRYPTO-CHROME1 or 2 (CRYs), then represses the transcriptional activity of CLOCK-BMAL1, closing an NFL. Also, since the molecular structures and functions of PER1 and PER2 proteins are similar [8,9], the repression of CLOCK-BMAL1 activity by PER1 and PER2 results in mutual repression between them (Fig 1A and 1B). The dual NFLs of Per1 and Per2 confer the regularity of gene expression rhythms, increasing the robustness of the circadian clock system [10,11]. Interestingly, one of the notable differences between Per1 and Per2 is the peak time of gene expression. In mice and rats, Per1 expression peaks at subjective midday, and is followed by Per2 expression, with about a 4-hour delay in the central pacemaker tissue, the suprachiasmatic nucleus (SCN) as indicated by in-situ hybridization [12][13][14] and q-PCR (e.g. ref. [15] and S1B Fig in ref. [16]) for the two Per genes. Moreover, in human U2OS cells, 2-hour difference between the expression peaks of PER1 and PER2 mRNAs was also observed [17]. Previous work has identified positive feedback of Per2 as a possible mechanism for the 4-hour peak time difference [18]. In addition, a separate study indicated that a functional noncanonical E-box (CATGTG) in the Per2 promoter also regulated the peak time of Per2 expression [15]. However, it remains unclear whether the peak time difference between Per1 and Per2 reflects their different functions in entrainment to LD cycles, as well as the determination of period and amplitude of oscillation.
In addition to sustained rhythm generation, Per1 and Per2 also play key roles in phase responses of the mammalian circadian clock to light signals. Light signals induce the transcription of both Per1 and Per2 mRNAs in the SCN (Fig 1A) [12,[19][20][21][22][23]. The resultant elevated levels of PER1 and PER2 proteins contribute to phase shifts of the clock and entrainment to LD cycles [7,24,25]. In addition to the difference in the peak times of Per1 and Per2 expression described above, experimental studies have shown that the magnitude of mRNA induction by light signals may also differ between them. For example, induction of Per1 mRNA by a short light pulse was stronger than that of Per2 in rat and hamster SCNs [12,26]. Yet, the biological significance of different induction levels between Per1 and Per2 by light signals for entrainment has also not been addressed.
Motivated by these experimental observations, here we use mathematical modeling and simulations to analyze the phase responses of two interacting NFLs to light signals. We show that a time difference between the expression peaks of the two repressors in the NFLs, as observed in Per1 and Per2, leads to functional differentiation in their phase responses to light signals. Specifically, the repressor with an earlier peak time (i.e. Per1) mainly contributes to phase advance, whereas the repressor with a later peak time (i.e. Per2) contributes to phase delay. Due to these complementary contributions, the ratio of light-induced transcription rates between Per1 and Per2 determines the proportion of phase advance and delay zones in a PRC, and thereby determines the entrainability of a circadian clock to a 24-hour LD cycle. Thus, the dual NFLs can possess different dynamical roles according to the peak time difference, resulting in the tunability of phase responses to external input signals.

Time delay model for two interacting NFLs
To analyze phase responses of oscillations to light signals, we consider two interacting NFLs composed of transcriptional repressor genes P1 and P2, corresponding to mammalian Per1 and Per2, respectively ( Fig 1B). Their protein products bind to each promoter and repress transcription in both. In general, transcription, translation, and transport of these protein products between the cytoplasm and nucleus require a certain time to be completed. Thus, the NFLs inherently include time delays, and such delayed negative feedbacks are known to produce self-sustained oscillations [1,27,28]. To describe the delayed feedbacks, we adopt a set of differential equations with time delay parameters [29][30][31]. The model includes P1 and P2 gene products, namely the levels of mRNAs m i available for translation in the cytoplasm, and those of proteins p i in the nucleus, as variables (i = 1,2). The time evolution of these mRNA and protein levels are described as: where β is the maximum light-independent transcription rate, α is the degradation rate of mRNA, v is the translation rate, and μ is the degradation rate of protein. We assume the linear degradation of both mRNAs and proteins. The first term of Eq (1a) represents the regulation of transcription independent of light signals, such as through E-box. We describe transcriptional repression by the repressor proteins with a Hill function of the coefficient n and the dissociation constant K of the repressor proteins to the promoters. To obtain this Hill function, we assume the cooperative binding of repressor proteins where binding of a protein to one of the binding sites in the promoter makes the protein of same type more likely to bind to neighboring binding sites [29,32]. However, the competition for the promoter would not affect results presented below, because we introduce the peak time difference between m 1 and m 2 later, and resultant separate expression of proteins eventually relax the competition. τ i in Eq (1a) represents time required for the processes to produce matured mRNAs m i available for translation, such as splicing, modification, and transport from the nucleus to the cytoplasm ( Fig 1C). Hence, the first term in Eq (1a), describing the rate of increase of the matured mRNAs at time t, is determined by the protein levels at the onset of transcription t − τ i , p 1 (t − τ i ) and p 2 (t − τ i ). Note that the timing of transcriptional repression by repressor proteins is the same for both P1 and P2 genes. For example, if the levels of P1 and P2 proteins become high enough to repress transcription at time t r , this effect of repression is reflected at time t = t r + τ 1 in m 1 , and at t = t r + τ 2 in m 2 . Similarly, T i in Eq (1b) represents the time required to translate mRNA and transport the products into the nucleus (Fig 1C). In the mouse SCN, the time difference between the expression peaks of a Per mRNA and the corresponding PER protein was about 4-6 hours [33]. Based on this data, we set T 1 = T 2 = 5 h unless mentioned otherwise in this study.
With appropriate values of time delays and reaction parameters, the model generates a stable limit cycle (Fig 1D and 1E). Although the regularity of autonomous rhythms in Per1 or Per2 deficient SCN was impaired [10,11], for simplicity, we choose a parameter regime where each single NFL can sustain rhythms even in the absence of the other NFL. A previous study indicated that a positive feedback loop (PFL) of Per2 was considered as a possible mechanism underlying the 4-hour peak time difference between Per1 and Per2 [18]. However, for simplicity, Eq (1) does not include a PFL. Instead, the peak time difference between m 1 and m 2 is controlled by the difference in time delays Δτ = τ 2τ 1 in transcription through E-box in Eq (1a) (Fig 1D and 1E). In fact, another experimental study showed that the peak time of Per2 expression was also regulated by a functional non-canonical E-box in the Per2 promoter in mice [15]. Note that Δτ and the peak time difference between m 1 and m 2 are equivalent in Eq (1) (S1A and S1B Fig). In addition, Δτ also reflects the difference in phases measured by the first Fourie components of m 1 and m 2 (S1A and S1C Fig), validating Δτ as a measure for phase difference. To focus on just the effect of the peak time difference between P1 and P2, we first assume that the maximum light-independent transcription rate, translation rate, degradation rate, and dissociation constants are the same between the two NFLs in Eq (1). We examine the effects of differences in the values of reaction parameters between P1 and P2 later in the result section.
The function γ i in Eq (1a) describes the transcription of mRNA induced by light signals. We assume the following rectangular function for this light-induced transcription [34]: ( where t l is the time at which a light signal is administered and T d is the duration of the light signal (Fig 1F and 1G). In this paper, t l = 0 indicates light administration starting at the time of an m 1 trough (Fig 1F and 1G). Since the expression levels of Per1 mRNA start to increase at dawn in the mouse SCN [16,35], time around a m 1 trough roughly corresponds to subjective dawn, and time 12-hour after corresponds to subjective dusk. A light signal induces the transcription of the repressor P i at a rate � i . � t i in Eq (1a) represents the time delay in mRNA production induced by light signals. We assume that � t 1 < � t 2 , because the levels of Per1 mRNA are elevated quicker than that of Per2 after light administration in the mouse SCN [9,36].
To compute the phase shift caused by the light signals, we simulate the time evolution of mRNAs and proteins with a light pulse at time t l (perturbed) and those without a light pulse (unperturbed). We simulate 50 cycles after the administration of the light signal to stabilize the trajectory (Fig 1F and 1G). Then, we calculate the peak time difference Δϕ of P1 mRNAs by subtracting the peak time of the perturbed trajectory from that of the unperturbed one. A positive value of Δϕ indicates phase advance due to the light signal (Fig 1F), whereas a negative value of Δϕ indicates phase delay ( Fig 1G). Thus, we define Δϕ as the phase shift induced by the light signal. In later sections, we may denote the phase shift as a function of the lightinduced transcription rates of P1 and P2, as Δϕ = Δϕ(� 1 , � 2 ).
In this study, we analyze the dependence of phase responses to light signals on the time delay parameters τ i in Eq (1a). We examine phase shifts within the ranges of the time delays where the period of autonomous oscillation T p nears 24 hours. Changes in the time delay parameters may also affect T p . Therefore, to better compare the magnitude of phase shifts with different parameter values, we scale the duration of a light signal T d in Eq (2) with T p as T d ! (T p /24)T d [34]. The details of numerical simulations and values of parameters are described in S1 Text and S1 File.
To quantify the shape of a PRC, we measure the unsigned areas of its advance (Δϕ � 0) and delay (Δϕ < 0) zones, which we refer to A and D (A � 0, D � 0), respectively. Then, we compute the fraction of A to the total area: If R is close to one (zero), the typical response of the NFLs to light signals is phase advance (delay). R should be a function of the light-induced transcription rates � i in Eq (2), R = R(� 1 , � 2 ).

Dependence of the period and amplitude on the two NFLs
We start the analysis by examining the difference in dynamical functions between two NFLs in the regulation of autonomous oscillations. In Fig 1E, we introduce a 4-hour peak time difference between m 1 and m 2 with Δτ = τ 2τ 1 = 4 h based on experimental data on mammalian Per1 and Per2. With this peak time difference, we observe that the time at which m 1 starts to decrease is close to the time when p 1 surpasses the value of the dissociation constant K in Eq (1a), and the decrease in m 2 follows 4-hour later ( Fig 1E). If p 2 /p 1 < 1 due to the peak time difference between the two proteins, the light-independent transcription rate in Eq (1a) can be approximated as β/(1 + (p 1 /K) n (1+(p 2 /p 1 ) n )) � β/(1+(p 1 /K) n ) where we omit time delay parameters for notational simplicity. In addition, when p 1 = K with p 2 /p 1 < 1, the light-independent transcription rate decreases to nearly the half of its maximum β/2. We confirm this approximation by observing the timeseries of the light-independent transcription rate and the levels of the two proteins with Δτ = 4 h (S2A- S2C Fig). Hence, the passage time at which p 1 surpasses K determines the onset of a repressed state where the transcription rate is less than β/2. We change this passage time of p 1 with respect to K by shifting the value of the time delay in protein translation T 1 in Eq (1b) with the condition τ 1 + T 1 < τ 2 + T 2 (S2D- We then examine the increase of mRNA. We notice that when p 2 decreases to p 2 = K, the light-independent transcription rate increases to nearly the half of its maximum β/2 (S2A- S2C  Fig). At the time p 2 = K, p 1 /p 2 < 1 in the presence of peak time difference (Fig 1E). Then, the light-independent transcription rate can be approximated as β/(1 + (p 2 /K) n ((p 1 /p 2 ) n + 1)) � β/ (1+(p 2 /K) n ) = β/2 (S2A- S2C Fig). Therefore, the time at which p 2 becomes smaller than K sets the onset of induction state where the transcription rate is larger than β/2. This onset of induction state is delayed by the increase in T 1 because it extends the time interval of p 2 increase as described above (S2F, S2G and S2N Fig). Thus, the increase in T 1 lengthens the autonomous period as well as amplitude (S2L and S2N Fig). The increase in the time delay T 2 for P2 translation in Eq (1b) also lengthens the autonomous period, as T 2 delays the time for p 2 to become smaller than K (S2H-S2K and S2O Fig). In summary, in the presence of peak time difference, the repressor with an earlier peak time determines the onset of transcriptional repression, whereas the other repressor with the later peak determines the onset of mRNA transcription.

Complementary contributions of two interacting NFLs to PRC
Next, to study contributions of each NFL to phase responses, we draw PRCs with the induction of either P1 or P2 mRNA by light signals (Fig 1F and 1G). For simplicity, we assume that time delays in light-induced transcription are the same as those for light-independent transcription We first analyze the two identical NFLs Δτ = τ 2τ 1 = 0 (Fig 2A). In this case, the waveforms of P1 and P2 mRNAs and proteins are the same. Correspondingly, PRCs obtained only by P1 induction (� 1 = 0.5, � 2 = 0 in Eq (2)) are identical to those by P2 induction (� 1 = 0, � 2 = 0.5) (Fig  2A). The PRC for simultaneous P1 and P2 inductions (� 1 = � 2 = 0.5) is approximately the sum of the phase shifts caused by each NFL, indicating the additivity of the two PRCs, Δϕ(� 1, � 2 ) � Δϕ(� 1 ,0)+Δϕ(0,� 2 ). This additivity of PRCs probably originates from the additivity of phase sensitivity or infinitesimal PRCs in phase reduction theory [37][38][39].
Then, we introduce a peak time difference between P1 and P2 by increasing τ 2 . As the difference in time delays Δτ becomes large, the difference in the contributions of P1 and P2 to PRCs increases (Fig 2B-2D). Notably, the induction of P1 mRNA, which peaks earlier, contributes to phase advance rather than phase delay. The area of the advance zone A in the PRC, i.e. the area above zero phase shift, increases with Δτ ( Fig 2D). In contrast, the induction of P2 mRNA contributes to phase delay as R decreases with Δτ ( Fig 2D). The PRC obtained by simultaneous induction of both P1 and P2 mRNAs is nearly the sum of each contribution, retaining the additivity of PRC (Fig 2B and 2C). Thus, the time difference between expression peaks of the two repressors separates their contributions to PRCs such that each NFL complements the other.
To quantify the magnitude of complementarity in the phase responses, we define an index: where R(�, 0) is the area fraction of advance zone in a PRC obtained only by P1 induction, with rate � defined by Eq (3), and R(0, �) is the area fraction obtained only by P2 induction. c becomes large with Δτ as the contributions of two NFLs to phase shifts become complementary ( Fig 2D). Interestingly, the index c increases until Δτ � 4, then it slowly decreases with Δτ ( Fig 2D), indicating the existence of optimal Δτ for complementarity. Subsequently, we examine the mechanism for this functional differentiation in phase responses (Figs 2G-2J and S3). P1 protein determines the time at which the light-independent  (1) and (2)  transcription rate decreases to its half maximum as described above. If a light signal further induces transcription of P1 mRNA when mRNA levels are increasing, P1 protein levels exceed the dissociation constant K in Eq (1a) earlier, due to increased translation by the excess mRNA ( Fig 2G). Therefore, the light-independent transcription rate decreases to β/2 earlier (S3A Fig), which results in lower accumulation of both P1 and P2 proteins ( Fig 2G). Consequently, the degradation of these proteins to levels lower than K takes less time, advancing the phase of oscillation. In contrast, light induction of P1 mRNA at the time interval where its levels are decreasing does not affect the light-independent transcription, because P2 protein levels are still abundant (Figs 2I and S3C). At this time interval, conversely, the light induction of P2 mRNA delays the time at which the light-independent transcription rate increases to β/2 by prolonging the transcriptional repression by excess P2 protein (Figs 2J and S3D). On the other hand, induction of P2 mRNA when its levels are increasing only weakly advances the phase of oscillation (Figs 2H and S3B). This is because P1 protein is already abundant, masking the effect of P2 induction. Note that these observed phase shifts are correlated with the transient change in the amplitude of P2 protein, like the covariation of the amplitude and period shown in the previous section.
The above analysis suggests that differences in time at which p 1 and p 2 become larger or smaller than the dissociation constant K are key to the complementary phase responses. These time differences can be influenced by not only Δτ but also time delays in translation T i . Indeed, if T 2 = T 1 − Δτ, the passage time at which p 2 becomes smaller than K is the same as the corresponding passage time of p 1 (S4A Fig). In this case, (τ 2 + T 2 )/(τ 1 + T 1 ) = 1 and the index c is almost zero (S4D Fig). As we increase T 2 , this passage time of p 2 with respect to K becomes later than that of p 1 (S4B and S4C Fig) and c grows (S4D Fig). Similarly, as T 1 becomes larger than T 2 , and closer to T 1 = T 2 + Δτ, the value of index c decreases. Thus, the complementary phase responses require the condition (τ 2 + T 2 )/(τ 1 + T 1 ) > 1. The experimentally observed 4-hour time difference between expression peaks of Per1 and Per2 mRNAs with similar time delays in translation T 1 � T 2 is a way to satisfy this condition.
Next, we examine the dependence of these complementary phase responses on the parameters involved in light induction of P1 and P2 mRNAs. We confirm that the complementary contributions of dual NFLs to PRCs are preserved with different values of time delays in lightinduced transcription � t i in Eq (1a), supporting the robustness of the results (S5A- S5C Fig). In contrast, we find that the index c decreases with an increase in light-induced transcription rate � (Fig 2E). A longer duration of the light signal T d also decreases c ( Fig 2F). As the values of these parameters increase, P1 and P2 inductions also cause substantial phase delay and advance, respectively (S5D- S5I Fig). A strong light induction of P1 mRNA at its trough results in p 1 larger than both the dissociation constant K and p 2 , delaying the phase of oscillation (S5J Fig). Similarly, by a strong light induction of P2 mRNA at its trough timing, p 2 exceeds K before p 1 , advancing the phase (S5K Fig). However, as long as the PRCs remain continuous and of type-1 in simulations, c never reaches zero (Figs 2E and 2F and S5D-S5I), indicating that this complementarity is an inherent property of the two interacting NFLs with the time difference between the expression peaks of the repressors.

Dependence of complementary phase responses on reaction parameters
In the previous section, we have studied the phase responses of the interacting NFLs as described in Eq (1), in which only time delays differ between the two loops. Since Per1 and Per2 gene products have highly similar protein sequences [9], their other reaction parameter values (e.g. translation and degradation rates, and dissociation constants to promoters) are expected to be similar. In fact, both PER1 and PER2 proteins are incorporated into the repressor complex with CLOCK-BMAL1 and CRY1/2 [40], confirming their similar dissociation constants. In addition, the degradation of both PER1 and PER2 is regulated by Casein kinase 1ε/δ [41]. Hence, it is worth examining whether the similar reaction parameters of P1 and P2 facilitate their complementary contributions to PRCs more than different values would. To investigate this possibility, we extend Eq (1) to incorporate differences in parameter values between the two NFLs (S1 Text). Then, we introduce the ratios of each reaction parameter between the two loops by nondimensionalizing the extended version of Eq (1) (S1 Text). Subsequently, we study the dependence of the complementarity index c on these reaction parameter ratios, keeping Δτ = 4 h (S6A-S6E Fig).
We find that the complementarity index c is maximized if the ratios of the reaction parameters are close to one (i.e. the reaction parameters are similar between the two loops) (S6A- S6E  Fig). If the production rate of the P1 protein is above or its degradation rate is below that of P2 protein, P1 induction by light signals tends to cause both phase advance and delay. In such cases, the peak value of the P1 protein becomes higher than that of P2 protein (S6F Fig). Consequently, p 1 levels are reduced to less than K near the same time as p 2 levels are. Hence, the passage time of p 1 with respect to K shifts to later than that of p 2 by the light-induced elevation of p 1 levels, delaying the time at which the light-independent transcription rate increases to β/ 2. In other words, complementary contributions of phase responses by a peak time difference are more likely to occur if the amplitudes of mRNAs and proteins are similar between two NFLs. Additionally, Hill coefficients in both NFLs should be large for such complementary contributions (S6G and S6H Fig).
If we constrain the reaction parameters in the two NFLs to be identical, as in Eq (1), the nondimensional model has only two reaction parameters -the ratio of protein degradation rate to mRNA degradation rate μ/α and the ratio of effective protein production rate to the dissociation constant vβ/(α 2 K) (S1 Text). If the ratio vβ/(α 2 K) is close to 1 in S6I-S6K Fig, the system converges to a steady state and no oscillation occurs. As this ratio becomes large, a stable limit cycle emerges via a Hopf bifurcation, and light induction of either P1 or P2 mRNA results in a continuous type-1 PRC in simulations. However, if the ratio is further increased, the PRC shifts to type-0. Therefore, we calculate the complementarity index c within the region where the PRC shape is type-1. c depends on vβ/(α 2 K) nonmonotonically and the peak is located near the Hopf bifurcation point (S6I- S6K Fig). As vβ/(α 2 K) becomes larger, P2 induction by light signals starts to cause phase advance. If a light induction of P2 mRNA occurs at around the trough of p 2 with a high value of vβ/(α 2 K), p 2 levels are more likely to surpass K due to the light induction. The elevated p 2 levels above K at such timing reduce the subsequent peak value of p 2 , leading to the earlier relief of transcriptional repression.
The complementary index c also depends on the ratio of degradation rate of protein to that of mRNA, μ/α. As μ/α becomes large, c peaks at a larger value of vβ/(α 2 K) and its peak value increases (S6I- S6K Fig). Rapid protein degradation prevents the accumulation of P1 and P2 proteins after the light induction. Hence, after the light induction of P1 mRNA at its decrease, P1 protein levels remain lower than P2 protein levels. Similarly, P2 protein levels stay lower than K after the light induction of P2 mRNA at its trough. In this way, fast protein degradation facilitates the functional differentiation of P1 and P2 in phase responses. A previous experimental study reported that the half-life of PER2 protein fused with a luciferase reporter (PER2::LUC) in the mouse SCN slices was about 1.9 hours [42]. The data for the half-lives of Per1 and Per2 mRNAs in the SCN is currently not available. However, in NIH3T3 cells, the half-life of Per2 mRNA was 0.9 hours and in embryonic stem cells, it was 2.9 hours [30,43,44]. Using these available data, we estimate μ/α = 0.47~1. 52

Complementary phase responses in dual NFLs including multiple states of mRNA and protein
So far, we have described the dynamics of repressors in dual NFLs using delay differential equations for the ease of controlling their peak time difference. Another description of time delays in NFLs is to model different states of mRNA and protein explicitly as the Goodwintype models [1,27,28,[45][46][47][48][49][50]. To examine whether complementary phase responses occur regardless of the description of time delays, here we analyze models that include multiple states of mRNAs and proteins as separate variables in ordinary differential equations (ODEs; S1 Text and S7 and S8 Figs).
First, to show how the inclusion of multiple states of mRNA affects period, amplitude and phase responses, we consider a single NFL with one repressor (S7A Fig). The model describes multiple states of mRNA (m 11 , m 12 , . . ., m 1u ) and proteins (p 11 ,p 12 , . . ., p 1r ) with ODEs. m 11 is the levels of transcribed nascent mRNA in nucleus and m 1u is the levels of matured mRNA at cytoplasm available for translation. Similarly, p 11 is the levels of translated nascent repressor protein in cytoplasm and p 1r is those of functional protein in nucleus. For simplicity, we assume linear chains of state transition from m 1i − 1 to m 1i , and p 1i − 1 to p 1i with time constants η and λ, respectively (S7A Fig). The details of the model are described in S1 Text. As the time constant η for the state transition of mRNA increases, the period decreases, whereas the amplitude increases (S7B and S7C Fig). In contrast, both the period and amplitude of oscillation increase with the state number of mRNA u (S7E and S7F Fig). Thus, these two parameters determine the delays in the production of functional mRNA. PRCs of this single NFL include both phase advance and delay zones (S7D and S7G Fig).
Next, we consider two interacting NFLs based on the description described above (S1 Text and S8A Fig). We assume that the time constants of state transition and state numbers are different between P1 and P2 mRNAs. A larger state number together with a smaller time constant of P2 mRNA than those of P1 mRNA generates peak time difference between their functional mRNAs with nearly 4 hours (S8B- S8E Fig). We observe complementary phase responses of the dual NFLs to light signals in this description as well (S8F and S8G Fig): light-induced transcription of P1 mRNA mainly causes phase advance, whereas that of P2 mRNA causes phase delay. Therefore, we conclude that complementary phase responses to light signals occur in the presence of the peak time difference, regardless of the description of time delays in dual NFLs.

Diversity of PRC shapes resulting from different ratios of light-induced transcription rates
It is known that PRC shapes differ among mammalian species. For example, PRCs of mice include a larger delay zone than advance zone, whereas those of rats and hamsters include a larger advance zone [4,51]. How does this diversity of PRC shapes arise? The complementary contributions of Per1 and Per2 to the PRC described above could be the reason for the different mammalian PRCs. To verify this prediction, we simulate phase shifts at various values of light-induced transcription rate of P1, � 1 in Eqs (1) and (2), with a fixed total rate � t = � 1 + � 2 . Thus, that of P2 is � 2 = � t − � 1 . Hereafter, we use Eqs (1) and (2) with identical values of the reaction parameters between the two NFLs. We quantify the difference in PRC shapes by the area fraction of advance zone R(� 1 ,� t -� 1 ) defined in Eq (3).
If the peak time difference between the two NFLs is small, PRC shapes by different ratios � 1 /� t are almost the same (Fig 3A and 3C). In contrast, if the peak time difference is large, the PRC shape strongly depends on � 1 /� t (Fig 3B and 3C). When � 1 /� t is small, the PRC includes a large delay zone, i.e. a small advance zone, resulting in small R. As � 1 /� t increases, the area of advance zone gradually expands and R reaches close to one. As mentioned in the introduction, the magnitude of light-induced transcription of Per1 mRNA in rat and hamster SCNs has been found to be greater than that of Per2 mRNA [12,26]. Hence, PRCs with a larger advance zone than delay zone in these animals [4,51] coincide with the prediction of the simulation. In summary, our theoretical results indicate that different PRC shapes can be generated depending on the ratio of light-induced transcription rates between two repressors in dual NFLs in the presence of time difference in their expression peaks.

Entrainment to a 24-hour LD cycle
As described in the previous section, the ratio of light-induced transcription rates between the two NFLs in the presence of the time difference between expression peaks of the two repressors Δτ controls the magnitude and direction of phase shifts at each subjective time. The magnitude of phase shifts by light signals determines whether a circadian clock with a certain autonomous period can entrain to an LD cycle. Hence, we identify the ratio of light-induced transcription rates that enables a circadian clock with the peak time difference Δτ and autonomous period T p to entrain to a 12:12 LD cycle (Fig 4). We control T p in this analysis by changing the time delays in translation T 1 and T 2 in Eq (1b) (S1 Text). As in the previous section, we fix the total rate as � t = � 1 + � 2 , and change the value of � 1 to examine entrainment. Fig 4A-4C shows the range of entrainment in the T p and � 1 /� t parameter space. For the identical NFLs Δτ = 0, the range of entrainment is symmetric about � 1 /� t = 0.5 (Fig 4A). To maximize the range of entrainment along the T p axis, the ratio of P1 induction should be either close to one or close to zero. This suggests that if the period difference is large, only one NFL should be light-responsive to attain entrainment with a fixed total rate � t .
With the peak time difference Δτ > 0, the range of entrainment becomes asymmetric about � 1 /� t = 0.5 (Fig 4B and 4C). If T p is shorter than the period of the LD cycle, the value of � 1 /� t needs to be small for entrainment to occur (Fig 4B and 4C). In other words, stronger light induction of P2 than P1 is necessary to delay the clock (Fig 4D and 4F). Conversely, if T p is longer than the period of the LD cycle, the value of � 1 /� t needs to be greater than 0.5 for such slower clocks to entrain to the LD cycle by increasing their speed (Fig 4B-4G). Thus, the repressor that should be induced by light signals more strongly than the other is determined based on the complementary contributions of the two NFLs and the autonomous period. In addition, the range of entrainment with peak time difference Δτ > 0 partly encompasses that with Δτ = 0, as shown in the bottom left (T p < 24 and � 1 /� t < 0.5) and top right (T p > 24 and � 1 /� t > 0.5) corners of Fig 4B and 4C. Thus, a light-induced transcription rate of P1 mRNA satisfying � 1 /� t > 0.5 (� 1 /� t < 0.5) in the presence of peak time difference can entrain a circadian clock with a longer (shorter) autonomous period compared to its absence. The ratio � 1 /� t also determines the peak times of P1 and P2 mRNAs in an entrained rhythm (S9 Fig). Taking a human autonomous period (~24.5 hours [52]) as an example, P1 and P2 mRNAs peak earlier as the ratio � 1 /� t increases (S9A-S9D Fig). Since the peak time of P1 and P2 mRNAs in a LD cycle would represent human chronotypes [53], this result suggests that the variations in the ratio of light-induced transcription rates between the two Per genes can be one of the determining factors of human chronotypes.
Taken together, two interacting NFLs have complementary effects on phase shifts and, as a result, the ratio of their light-induced transcription rates determines entrainability to LD cycles.

Discussion
In this study, we examined two interacting NFLs for genetic oscillations to better understand the mammalian circadian clock system. In mammals, Per1 and Per2 consist of interacting NFLs that are responsible for circadian genetic oscillations. Peak time differences between Per1 and Per2 expression have been observed in the central pacemaker tissue, SCN [12][13][14][15][16]. Furthermore, both Per1 and Per2 play key roles in entrainment to LD cycles [7,24,26,54]. The lack of one of the two Per genes resulted in irregularity in autonomous rhythms in SCN, and increased sensitivity of circadian rhythms to genetic background [10,11], suggesting that the dual NFLs of Per genes enhance robustness of the circadian clock. However, besides their role as redundant transcriptional repressors, functional differences between Per1 and Per2 gene products have yet to be revealed. Our results predict the novel and distinct roles of Per1 and Per2 in phase responses to light signals: light induction of Per1 mRNA mainly contributes to phase advance, whereas induction of Per2 mRNA contributes to phase delay.
Our analysis revealed that difference in passage time at which the levels of PER1 and PER2 proteins become higher or lower than the dissociation constant to E-box is key to their functional differentiation in phase responses to light signals. A previous immunohistochemical study observed that the PER1 and PER2 protein expression peaked at CT12 in the mouse SCN [33]. Since the modification such as phosphorylation of PER proteins influences their repressor assembly and activities [49,55], the complementarity of PER1 and PER2 proteins must be reflected by the kinetics of their presence on E-box. In fact, a Chip-seq analysis using mouse liver indicated that PER2 protein stays at the E-boxes of circadian clock gene promoters 1~4-hour later than PER1 [56], suggesting that the timepoint of PER2 protein levels below the dissociation constant is later than that of PER1. Similar experiments using the SCN would be worth to evaluate our conclusion.
We addressed phase responses of circadian rhythms to light signals in a single SCN neuron, since its phase responses is the basis of those at tissue and individual levels. In the SCN, neurons interact with each other via various neurotransmitters, such as vasoactive intestinal peptide and arginine vasopressin [57]. In addition to keeping precise ticking [10], the intercellular coupling may also modulate responsiveness of oscillators to environmental signals [58,59]. Phase responses of coupled circadian oscillators have been quantified by neuronal firing rate in SCN slices [21,60]. Induction of Per1 mRNA in the rat SCN slices by glutamate advanced the circadian phase of the neuronal firing rate [21]. Importantly, the phase advance by glutamate was inhibited in the presence of an antisense oligodeoxynucleotide against Per1 mRNA, but not by antisense oligodeoxynucleotide against Per2 [21]. Thus, these results are consistence with our theoretical prediction in a single neuron.
The complementary contributions of interacting NFLs to phase responses may underlie the diverse PRC shapes observed in different animal species. Even within mammals, the area ratios of advance to delay zones in PRCs differ among species [4]. One way to modify the shape of a PRC is by tuning reaction parameters within NFLs, such as transcription and degradation rates of Per1 and Per2 mRNAs [34], although such changes may also affect the period and amplitude of oscillation. Another possibility could be to gate light input signals to the two Per genes in the SCN depending on time of day, to modulate the ratio of advance and delay zones [61,62]. However, such gating would require an elaborate clock-dependent mechanism. On the other hand, our current study revealed a new mechanism to create diverse PRC shapes only by changing the ratio of light-induced transcription rates between Per1 and Per2 in the presence of peak time difference in the SCN.
The PRC shape has been considered important for stable entrainment to LD cycles [4,7,63]. With a phase oscillator model, a previous theoretical study analyzed an optimal PRC shape for stable entrainment to LD cycles, when a circadian clock system includes two light-sensitive reactions [64]. Interestingly, if there is a phase difference between two variables influenced by these light-sensitive reactions, to maximize entrainability, each of these should contribute to phase responses in a complementary way, with one variable mainly contributing to phase advance, whereas the other contributes to phase delay [64], like P1 and P2 in the current study. However, the abstract phase oscillator model cannot address questions whether and how such PRCs are realized by gene regulations in the circadian clock. By modeling gene regulatory network, our current study offers a possible mechanism for complementary phase responses with the two light-responding clock genes, Per1 and Per2. Mammalian circadian clocks also include other redundant NFLs such as Cry1 and Cry2, and Rev-erbα and Rev-erbβ [40]. A future modeling approach similar to the current one might reveal functional differentiation of these redundant feedback loops. In addition, we notice that dual NFLs with a 4-hour difference between expression peaks of the two repressors generate an interval in a PRC where phase shifts are close to zero, termed a dead zone [4,34], but those with a 1-hour difference do not (Fig 3A and 3B). Thus, the dual NFLs of Per1 and Per2 would be responsible for not only diverse PRC shapes but also characteristics of type-1 PRC observed in mammals.
In rats and hamsters, the area of the advance zone in a PRC for behavioral rhythms is larger than that of the delay zone [4,51]. The current theoretical results predict that such larger advance zones result from the greater light-induced transcription of Per1 than Per2. Consistent with this prediction, Per1 expression in their SCNs was found to be strongly induced by short light signals at subjective night, whereas Per2 expression was weaker than Per1 [12,26]. Furthermore, our model predicts that if the period of an LD cycle is shorter than the autonomous period of animals, stronger induction of Per1 than Per2 mRNAs is required to entrain to the LD cycle, and vice versa. Experimental results for rats and hamsters are consistent with this prediction (S1 Table) [26,54], suggesting that the ratio of light-induced transcription rates between Per1 and Per2 would be constrained by the difference between the autonomous and the LD periods in these animals.
In mice, administration of light signals that caused phase delays in behavioral rhythms at subjective night resulted in weak induction of Per1 mRNA but strong induction of Per2 mRNA in the dorsal region of the SCN (S1 Table) [23,25]. The current model suggests that such strong Per2 induction in the dorsal SCN may underlie the larger delay zone in mouse PRCs [4,7]. Conversely, when a light signal caused phase advance in behavioral rhythms at subjective late night, the induction levels of Per1 mRNA were higher than those of Per2 mRNA in the dorsal region [23,25]. Thus, the correlations between relative induction levels of Per genes in the dorsal region of the mouse SCN and the direction of phase shifts are consistent with current theory. In the ventral region of the mouse SCN, however, Per1 mRNA was induced by light signals that caused a phase advance [23], phase delays [19,65] and even no phase shift [23] in behavioral rhythms (S1 Table), probably reflecting a major role of this region in the reception of light signals from retina. In summary, with few exceptions, experimental observations from mammals have consistently indicated that the induction of Per1 mRNA contributes to phase advance, whereas that of Per2 mRNA contributes to phase delay.
Although autonomous rhythms in SCN cells are compromised in Per1 or Per2 deficient mice, their behavioral rhythms are almost indistinguishable from WT [10,11]. A previous experiment measured the PRCs of mice lacking either functional Per1 or Per2 [7]. These mutant mice possessed only a single NFL, but this remaining NFL could cause both phase advance and delay, according to our theoretical results. Therefore, an experiment to verify the current model predictions would require the inhibition of light inputs into one of the two Per genes without any effect on the circadian expression of both genes. Light signals lead to the activation of cAMP response element binding protein (CREB), which induces transcription of Per1 and Per2 after binding to the CRE element in their promoter regions. Since circadian expression of the two Per genes depends not on CRE but on the E-box element in their own promoters, isolation of CRE mutant mice of Per1 or Per2 should satisfy the above requirement. Alternatively, the complementary contributions of Per1 and Per2 to phase responses could be examined by specific induction of one of the two Per genes. Although such chemical inducers are not available at present, the current study should provide motivation for their screening for the treatment of circadian rhythm disorders and adjustment for different chronotypes. Previous genome-wide association studies associated morningness with single nucleotide polymorphisms near human PER1 and PER2 loci [66,67], confirming the relevance of these two PER genes to human chronotypes. The current study predicts that a chemical compound that specifically induces PER1 mRNA may be able to cause phase advances only, regardless of administration time, and could be useful for the treatment of delayed sleep phase syndrome [68,69]. Conversely, a compound that targets PER2 specific induction might be used to delay the clock at any zeitgeber time, which would be useful for the treatment of advanced sleep phase syndrome [69,70]. Thus, as our simulations suggested, the administration of such compounds may adjust the phase of entrainment to be more desirable for daily life.
Although we focused on the interlocked NFLs of the two Per genes, the mammalian circadian clocks include other positive and negative feedback loops [40,71]. Importantly, the expression of clock genes containing E-box in their promoters, such as Per1 and Rev-erbs, tends to peak at earlier time of a day, whereas the expression of those containing ROR element in their promoters, such as Bmal1, peaks at later time [16,40,72]. The different phases of these clock gene expression would reflect their specific functions important for sustaining robust circadian oscillation [31,48,50] and entrainment. How these multiple feedback loops modulate the complementary phase responses by the two Per genes is an important question that should be addressed in future study. Since functional differentiation between PER1 and PER2 proteins is predicted to depend on the similar repression activities of both proteins, it is not necessary to revise our conclusion as long as the other feedback loops do not affect the relative repression activities.
In conclusion, we revealed the complementary contributions of two interacting NFLs to phase responses and entrainment. Our results suggest that the peak time differences between transcription factors in the multiple feedback loops lead to their functional differentiation. Such functional differentiation may be key to understanding temporal dynamics in complex gene regulatory networks.