A Phenomenological Mouse Circadian Pacemaker Model.

Mathematical models have been used extensively in chronobiology to explore characteristics of biological clocks. In particular, for human circadian studies, the Kronauer model has been modified multiple times to describe rhythm production and responses to sensory input. This phenomenological model comprises a single set of parameters which can simulate circadian responses in humans under a variety of environmental conditions. However, corresponding models for nocturnal rodents commonly used in circadian rhythm studies are not available and may require new parameter values for different species and even strains. Moreover, due to a considerable variation in experimental data collected from mice of the same strain, within and across laboratories, a range of valid parameters is essential. This study develops a Kronauer-like model for mice by re-fitting relevant parameters to published phase response curve and period data using total least squares. Local parameter sensitivity analysis and parameter distributions determine the parameter ranges that give a near-identical model and data distribution of periods. However, the model required further parameter adjustments to match characteristics of other mouse strains, implying that the model itself detects changes in the core processes of rhythm generation and control. The model is a useful tool to understand and interpret future mouse circadian clock experiments.

The suprachiasmatic nucleus (SCN) is the site of the central clock in mammals and is responsible for regulating the body’s circadian cycles and influencing the response to various stimuli. Many important rhythms such as the sleep–wake cycle and body temperature cycle are circadian, and are thus controlled by the SCN. It has long been of great interest within the field of circadian biology to construct a working model of the SCN. However, due to the inherent coupled complexity of biological systems, determining the connections between genes, protein networks, neurons, and neuronal circuits and their relation to behavior has been a formidable task in both theoretical and experimental biology.

A phenomenological model of the human circadian pacemaker, originally by Kronauer (1990), has been developed, studied, and iterated upon over the last few decades (Phillips et al., 2011; St Hilaire et al., 2007a; Forger et al., 1999; Jewett et al., 1999; Kronauer et al., 1999; Klerman et al., 1996). Additional phenomena such as non-photic stimuli (St Hilaire et al., 2007a, 2007b) have been incorporated, while modifications to the model and changes to the parameter values have improved its accuracy in producing human phase response curve (PRC) and resetting data.

Other examples of general phenomenological models for mammals have been developed to understand the dynamics of the circadian clock, such as those by Daan and Berde (1978), Kawato and Suzuki (1980), and Schmal et al. (2020). However, to our knowledge, there does not exist a correspondence to the model by Kronauer (1990) that describes the circadian behavior in mice.

There are many mechanistic models that are used to decipher the molecular interactions underlying the circadian clock. These include the models of Herzel (Becker-Weimann et al., 2004a, 2004b), Forger (Forger and Peskin, 2003; Kim and Forger, 2012; Dewoskin et al., 2014), and others (Podkolodnaya et al., 2017; Mirsky et al., 2009; Lema et al., 2000). In addition, there are single-cell models of Relógio et al. (2011) and Shiju and Sriram (2017), and the multi-cell models of Vasalou and Henson (2010), To et al. (2007), and DeWoskin et al. (2015). However, these models are computationally expensive and difficult to corroborate with behavior. A simpler, phenomenological model is much more convenient for analysis and comparison with the available mice behavioral data for photic and non-photic inputs.

Given the similarities of the structure of the circadian clock between mice and humans, we reason that using the Van der Pol oscillator equations of the Kronauer model, originally developed for the human pacemaker, should also be suitable to model the mouse pacemaker. However, the parameter values used in the human model are not applicable to mice and do not result in good agreement with mouse experimental data. Furthermore, because the original human pacemaker parameters were fit based on an optimization approach that minimized the perpendicular distance to the data, these parameters are only suitable when describing the average process and fail to describe specific samples in the population.

Many authors have measured the free-running period of the wild-type C57BL/6 mouse in constant conditions and have reported widely varying results. It has been discussed in several of these studies that physical characteristics, such as age, may play an important role. In particular, it has been noted that the free-running period of the C57BL/6J mouse for wheel-running activity under constant dark conditions appears to lengthen with age (Mayeda et al., 1997; Possidente et al., 1995; Valentinuzzi et al., 1997). Other important factors that have been considered are the effects of access versus restriction to a running wheel (Capri et al., 2019); different types of white light (Alves-Simoes et al., 2016); particular wavelengths of light (Hofstetter et al., 2005); sub-strains of C57BL/6 mice, particularly that of the 6J and 6N varieties (Capri et al., 2019; Ebihara et al., 1978); and perhaps even the presence of light-induced retinal damage (Gonzalez, 2018).

In this article, we develop a fully functional phenomenological model of the mouse circadian clock. We focus on one of the most commonly used mouse strains: the wild-type C57BL/6 mouse. We use the existing formulation of the model (Kronauer, 1990), but re-fit the parameters to experimental PRC data. We demonstrate the robustness of the fit mouse parameters by validating against other PRCs available in the literature and comparing the results with that of the human model. Then, because of the large variation in the free-running period within and across labs, we perform a local sensitivity analysis to determine the most relevant parameters that affect the simulated period in constant conditions. With these relevant parameters, we deduce the range of potential values by parameter variation and fitting the simulated free-running period to the experimentally observed free-running periods in the literature. Finally, we demonstrate some applications of the mouse model: The model is used to deduce the maximal phase response to photic stimuli and the model is adapted to describe a mutant mouse.

Methods

Mathematical Model

To model the circadian pacemaker for mice, we begin with Kronauer’s oscillator model (Kronauer, 1990), derived from the Van der Pol oscillator, the proto-typical limit-cycle oscillator. The equations describing the evolution of the state of the circadian pacemaker, known as Process P, are as follows:

In the absence of light , Equations 1 and 2 generate a limit-cycle oscillation in and with a period of . The parameter defines the stiffness of the oscillator, and the coefficients in Equation 1 have been selected so that the amplitude of the oscillation is approximately 1. The correction factor forces the constant darkness (DD) period of the oscillator to be .

Process P is driven by light, through Process L, according to the following equations:

The variable is interpreted as the proportion of activated photoreceptors; is the proportion of ready photoreceptors. The rate at which ready photoreceptors are converted to their active state is a function of light level given by Equation 6. The photoreceptors generate photic drive to the clock that is proportional to their activation rate with proportionality constant . This drive is modulated by the state of clock as described by Equation 5. The parameter determines the circadian sensitivity modulation of the photic drive to light. The rate is the rate at which the photoreceptors are converted back to their ready state. This rate is considered to be independent of light since Walch et al. (2015) found that light-dependent values did not significantly enhance fits of multi-electrode array data. The constant is the conversion factor from minutes to hours for the rates and .

The parameter determines the photic drive strength. In previous works, it was suggested that for , as the strength of increased, the average observed model period decreased (Jewett et al., 1999). Hence, this requirement in the model agreed with Aschoff’s fourth rule, which states that the free-running period of diurnal organisms shortens as the light intensity is increased. However, in this study, it is determined that also satisfies Aschoff’s third rule, which states the opposite effect for nocturnal organisms: Increasing the light intensity increases the free-running period. To confirm this, our simulations showed that increasing the drive strength also increased the average model period under constant light (LL).

MATLAB code for simulating the model and generating a PRC is available on ModelDB (McDougal et al., 2017) at http://modeldb.yale.edu/267250.

Fitting Data Sources

The model parameters were determined by fitting to data sets from Pendergast et al. (2010) and Vajtay et al. (2017) simultaneously. The experimental protocols for these data sets are described below.

Pendergast et al. (2010) set up singly housed wild-type C57BL/6J mice in a cage where they free-ran for 6 days in DD. The authors used linear regression to determine the onset of activity at 12 h circadian time (CT 12) on the following day and, at the appropriately calculated CT, the cages were transferred to another light-tight box and given light of 150 lx (measured at the cage bottom), where the mice ran for 15 min. Then, the mice were returned to DD and free-ran for a week.

Vajtay et al. (2017) entrained wild-type C57BL/6J mice (7-8 males) for 2 weeks under LD12:12 and then released them into DD. Every 14 days, the lights were turned on for 3 h at appropriate CTs. In this experiment, the authors note that the light was turned on such that the midpoint of the light exposure correlated with the appropriate CT of the mouse. The light intensities were in the range of 250 to 400 lx at bedding level. Each mouse had a total of 8 randomly timed light exposures. The phase shifts were then averaged for each 3-h bin, giving 5 to 9 measurements for each data point.

Validation Data Sources

For model validation, we did not alter our fit parameters and, following the respective light protocols and experimental setups, compared the simulation results against data from Vitaterna et al. (2006) and Comas et al. (2006).

Vitaterna et al. (2006) entrained C57BL/6J mice to a LD12:12 cycle for 1 week and subsequently released them into DD. After 3 weeks of free-running in DD, a 6-h light pulse was given at an appropriate CT, calculated from the period and onset of activity (CT 12). The data were collected for an additional 10 days in DD after the light pulse.

In Comas et al.’s (2006) experiments, C57BL/6J and C57BL/6J-OlaHsd mice were entrained for 2 weeks in LD12:12 so that the initial phases were all the same. The mice were then released into DD and given a pulse of light of 100 lx every 11 days.

As a further validation of the model parameters, and to demonstrate the sensitivity of the parameters to specific strains of mice, we also compared our model with the data for a non-C57BL/6J mouse. Ouk et al. (2019) used mice from a CBA × C57BL/6J background. The mice were entrained to an LD12:12 cycle for 7 to 10 days, followed by 10 to 14 days in DD. Then, at appropriate CTs, the mice were exposed to 15 min of white light at 300 lx. The data were collected for at least 10 to 14 days after the pulse to determine the phase shift calculations.

Fitting Procedure

The model fitting was performed in MATLAB using global parallel optimization, with the method of total least squares. The simulation was set up as follows. First, the “animal” is entrained to a 400 lx LD12:12 cycle, using an arbitrary initial condition on for 50 days. Then, using the results for as an initial condition, we simulated for 11 cycles in DD, where each cycle length corresponded to the intrinsic period of the mouse, . Ignoring the first 7 days, the periodic minima were measured. Finally, using the output from in DD as the initial condition, we ran the simulation again, but with a light pulse at the appropriate CT. Again, ignoring the first 7 days, the minimum of is calculated. The phase difference is defined as the difference between the initial before the pulse and after, with negative values corresponding to phase delay and positive values corresponding to phase advance.

The experiments did not report any error in the independent (time) variable; however, there is inherent uncertainty in determining the precise CTs for the light stimulus, which are calculated through linear regression. For this reason, it is more sensible to fit data according to the perpendicular distance between the data and the model curve. Thus, our objective function is defined to be the squared sum of the minimum Euclidean distance between each data point, , and its (closest) point on the curve generated by the model, :

in which denotes the parameters of the model. We construct 2 loops for the minimization. The outer loop, in parallel, minimizes the Euclidean distance between the data and model, as given by Equation 7. The inner loop finds the corresponding model point that is the closest to the data point . We also incorporated simultaneous fitting for the free-running period of the mouse under LL: The periods were calculated and fit against the provided average experimental period of 25.3 h (Pendergast et al., 2010).

Parameters and were fixed in accordance with previous model parameters as and . For the remaining parameters, we varied . We chose specific bounds on these parameters as larger intervals yielded unstable results in the simulations. Moreover, to mitigate any further instability in the resultant phase difference, the computational code ignored high error parameters by specifying a maximum allowed mean squared error between the model and data points. The parameters were bounded as , , , and .

We modified parameter , which is determined as a correction to fix the measured period of the oscillator to its intrinsic period . Due to the variable nature of biological parameters, such as that of the intrinsic period , it is unnecessary to compute through the Poincaré-Lindstedt method (Strogatz, 1994) for any choice of intrinsic period and oscillator stiffness to sufficiently high order. Instead, we chose to programmatically fit the parameter so that the measured period is sufficiently close to that of the chosen intrinsic period. The objective function was defined to be the mean squared sum of the error between the intrinsic period and calculated periods.

Results

For our nominal parameters, it was determined that for with an intrinsic period of , provided the best fit for the observed model period in DD to match the intrinsic period.

The optimal parameters determined by the fit are provided in Table 1.

Table 1.

A comparison of model parameters between the human model (St Hilaire et al., 2007b) and our mouse model.

Parameter	Human	C57BL/6J	Units	Description
α_0	0.1	1.8	min–1	Maximum response rate to light
β	0.007	0.005	min–1	Decay rate of light input
k	0.55	0.20	1	Photic drive strength
b	0.4	0.59	1	Modulation constant
G	37	52	1	Scaling factor
p	0.5	0.64	1	Light expansion exponent
τ_x	24.2	23.6	h	Intrinsic period

Note, the intrinsic period is not a fit parameter of the model equations (Equation 2)—It is determined directly from data.

The model parameters were fit against data from Pendergast et al. (2010) and Vajtay et al. (2017). The comparison between the optimal fit parameters and the experimental data is shown in Figure 1.

Figure 1.

Fit model and experimental data as provided in Pendergast et al. (2010) and Vajtay et al. (2017). The phase difference was calculated as the average of the difference between the minimum of during a constant darkness simulation prior to the light pulse, and the minimum of after the light pulse, ignoring the first 10 days. Abbreviation: CT = circadian time.

The model was validated against data from Vitaterna et al. (2006) with the mouse model parameters (Table 1) providing a very good match with the experiment, as shown in Figure 2. The experimental data involved measuring the phase shift in 61 wild-type mice in response to 6-h light pulses at appropriate CTs. For a summary of the experimental protocols, see “Validation Data Sources” section. Data from a representative PRC for a 6-h light pulse in Nakamura et al. (2016) are included as well for comprehensiveness. In Figure 2, we also plot and compare our model PRC with the PRC generated using the human model parameters (Table 1). To be consistent, we have used identical prior photoperiods, LD12:12, for both the mouse and human simulations. One can clearly see that these parameters do not yield an appreciable phase response for the identical experimental protocol: The human model results in a maximum phase delay of 1.76 h and a maximum advance of 1.45 h, significantly lacking in both the advance and delay zones present in the data.

Figure 2.

Validation against data for a light pulse of 6 h (Vitaterna et al., 2006; Nakamura et al., 2016) using the optimal mouse model parameters, compared with the human model parameters (Table 1). As in Figure 1, the phase difference was calculated as the average of the difference of minima before and after the light pulse. See Table 1 for parameter values. Abbreviation: CT = circadian time.

We have also compared our model simulations with that of Comas et al. (2006) for light pulses of duration 1, 3, 4, 6, 9, and 12 h, shown in Figure 3. We emphasize that the parameters (Table 1) were completely unaltered. From the figures below, it is clear the model provides an excellent match with the data, capturing the slight advance region in the early CTs and transition into the dead zone for pulses of 1 to 6 h, as well as the lack of a dead zone in the 9- and 12-h pulses.

Figure 3.

Validation with experimental data (Comas et al., 2006) for light pulses of duration (a) 1 h, 9 h; (b) 3 h, 6 h; and (c) 4 h, 12h at 100 lx. The mouse model demonstrates an excellent qualitative agreement with the data, clearly showing the lack of dead zones for the 9- and 12-h pulses as well as the gradual shift from slight advance to the dead zone to the delay zone in the other pulses. Abbreviation: CT = circadian time.

Figure 4 demonstrates another validation of the model. With the C57BL/6J mouse model parameters (Table 1), comparison of the generated model curves to the experimental data presented in Ouk et al. (2019) results in a poor match, particularly near the region of phase delay near CT 15. This is expected since the mice used in Ouk et al. (2019) are not wild-type C57BL/6J mice, but rather mice on a CBA × C57BL/6J background. There are several experimental discrepancies in circadian photosensitivity results between the CBA and C57BL/6 mouse, which are discussed by Yoshimura et al. (1994). The maximum absolute error between the data and model curves in Figure 1 is 1.34 and 1.10 h, and the average estimated error is 0.463 and 0.486 h for Pendergast et al. (2010) and Vajtay et al. (2017), respectively. Comparatively, in Figure 4, the maximum absolute error between the data and model curves is 2.47 h, while the average estimated error is 0.818 h. Due to this, we cannot expect our mouse model, fit for C57BL/6 mice, to provide a close comparison with these experimental data.

Figure 4.

Comparison with experimental data (Ouk et al., 2019): The model parameters determined from the fit in Figure 1 are valid for a C57BL/6J mouse—The mice used in Ouk et al. (2019) are on a CBA × C57BL/6J background. As a result, the model does not have good agreement with the data especially in the phase delay region near CT 15. Abbreviation: CT = circadian time.

Local Sensitivity Analysis

It is well known that there is a great degree of intrinsic variability in the wild-type mouse free-running periods under both DD and LL. Some reasons for the intervariability between mice have been determined, such as sourcing from different labs (Capri et al., 2019), exposure to different wavelengths of light (Hofstetter et al., 2005) or even different types of white light (Alves-Simoes et al., 2016), and measuring the periods at different times in the mouse’s life cycle (Possidente et al., 1995). Even with these factors controlled, there is still significant variability between mice within the same environment. For these reasons, the parameters provided in this study should also include a range of possible numerical values.

First, we performed a local sensitivity analysis of the mouse model using a central finite difference method with a one-at-a-time (OAT) approach (Saltelli et al., 2008) to deduce the most significant parameters. For a parameter , the oscillator period’s sensitivity to changes in can be estimated as

where is defined as the standard unit vector of the ith dimension of the space containing .

Using Equation 8, we perturbed each parameter and determined the corresponding normalized sensitivities of the oscillator period at nominal for the lighting conditions DD, LD12:12, and LL, of which the latter two (conducted at 150 lx) are listed in Table 2. Comparing these protocols, we can see that the sensitivity of the circadian period is extremely low in LD. Meanwhile, the perturbations have a significantly more noticeable effect on the period in LL, particularly for parameters . The sensitivity of the period in DD was also determined; however, it was negligible.

Table 2.

Period sensitivities to a parameter perturbation of under LL and LD12:12 lighting conditions.

Parameter	Nominal Value	Units	s_LL(10^−3)	s_LD(10^−6)
α_0	1.8	min–1	3.8	27
β	0.005	min–1	97	485
k	0.20	1	–31	26
b	0.59	1	103	706
G	52	1	101	585
p	0.64	1	–10	–317

Abbreviation: LL = constant light. Period sensitivities to constant darkness are not presented here as it was determined to be effectively negligible: Sensitivity values were on the order of .

We plot the model periods under LL, LD, and DD in Figure 5 for each perturbed parameter, up to variations of 25%. As reflected in Table 2, the period in LD and DD does not change appreciably under perturbations. The most relevant parameters that affect the period are ; the decay rate, ; the proportionality constant for the photic drive; and ; the modulation constant.

Figure 5.

Model period under different lighting conditions: (a) LL, (b) LD12:12, and (c) DD by perturbing each model parameter by with an OAT approach. The period under (a) LL is sensitive to the parameters, decreasing with and increasing with . The period under (b) LD and (c) DD is insensitive to all parameters. Abbreviations: LL = constant light; DD = constant darkness; OAT = one-at-a-time.

The sensitivity analysis conducted suggests that the periods in LL are particularly sensitive to changes in the parameters . In addition, the sensitivity analysis also illustrates that the model period in DD is insensitive to parameter perturbation. This is expected because of the parameter , which was modified so that the model period is equal to the intrinsic period .

We also demonstrate the effects of varying and on an example PRC from Comas et al. (2006). We did not include the other parameters because ’s effect is nearly identical to , while the effects of the other parameters were less dramatic. In Figure 6, we compare the PRCs from the optimal parameters in Table 1 with the PRCs generated assuming the parameters and were varied, one at a time, by . It can be seen that increasing parameter increases the peak-to-peak amplitude of the PRC. Increasing , on the contrary, increases the effect of phase delay and widens the range of the dead zone (see “Discussion” section). The 1.25b PRC appears to be a closer match to the experimental data with an average weighted PRC error of 1.08 h as opposed to the optimal with an error of 1.1 h; however, the period of the oscillator in LL was determined to be 25.9 h, which is significantly greater than the mean period of 25.3 h that the original parameters were fit to. While Comas et al. (2006) did not report the free-running periods in their experiments, it is very unlikely the mean period in their study would be greater than 2 standard deviations away from the general mean mouse period under LL (see Figure 7).

Figure 6.

Comparison between a 3-h pulse phase response curve from Comas et al. (2006) with the optimal parameters and perturbations of the parameters . Increasing increases the effect of phase delay and widens the range of the dead zone. Increasing increases the effect of phase delay and phase advance. Abbreviation: CT = circadian time.

Figure 7.

Period determined by model simulation under DD (blue histogram) and LL (yellow histogram) assuming a normal distribution for parameters . Black lines represent a period distribution assuming mean and standard deviations of m, and , under DD and LL, respectively. Data are plotted as mean ± standard deviation in a swarm chart, where each row represents a different source of data. Filled-in circles represent DD data, while empty circles represent LL data. Rows that include multiple data points are for different experimental conditions. “Ralph lab” refers to the (unpublished) running wheel data collected from the Ralph Chronobiology Lab (Department of Psychology, University of Toronto). Abbreviations: DD = constant darkness; LL = constant light. Color version of the figure is available online.

Parameter Variability

To determine the ranges for the key parameters, we re-fit with available data in the literature. We collected a wide range of studies, sourced from 12 different labs that reported the DD and/or LL periods of the C57BL/6 mouse. We calculated the total mean and standard deviation by combining the reported standard deviations and means from each data set , using the computational formula for the sample variance:

where is the mean of the ith data set, is the number of mice in the ith data set, is the standard deviation of the ith data set, is the total number of mice from all sets, and is the total mean of the population. The total number of mice exposed to DD from all of the collected data sources was . The mean and standard deviation were calculated to be and , respectively. The total number of mice exposed to LL from all of the data sources was . The mean and standard deviation were calculated to be and , respectively.

The parameter variability was determined as follows. First, the parameters were re-fit so that the observed period in LL was equal to the average experimental period in LL, . Then, assuming that follow a normal distribution with a coefficient of variation, the model was simulated for  days in DD, followed by  days in LL. The corresponding periods in DD and LL were calculated, and their mean and standard deviations determined. Then, the coefficients of variation were varied until the mean and standard deviations of the observed period were sufficiently close to the mean and standard deviations of the data. Data were sourced from a wide range of sources: Vitaterna et al. (1999), Capri et al. (2019), Alves-Simoes et al. (2016), Schwartz and Zimmerman (1990), Kiessling et al. (2017), Pendergast et al. (2010), Ebihara et al. (1978), Ali et al. (2019), Hofstetter et al. (2005), Loh et al. (2013), Valentinuzzi et al. (1997), and from the Ralph Lab at the University of Toronto.

The best coefficients of variation were determined to be for , respectively. In Figure 7, we plot a histogram of periods, assuming the parameters follow a normal distribution with their respective coefficients of variation. The experimental data are plotted as a swarm chart—Each line in the vertical represents data sourced from a specific lab. The DD periods are organized by mouse age, where the lowest point in the swarm chart has a period of with a mouse age of 19 to 22 months. The ages of the mice in Vitaterna et al. (1999) are unknown. The black curves represent a normal distribution of periods assuming the experimental data’s total mean and standard deviation, , and , under DD and LL, respectively. Overall, the standard deviation of the period is greater under LL than DD; however, we suggest this may be due to light being the most important zeitgeber, and so differing protocols of photic stimuli (e.g., different wavelengths of light, types of light) used play a significant role in affecting the measured periods between labs. This argument is strengthened by the fact that for mice from each lab, the standard deviation of the periods under LL is comparable with the standard deviation of the periods under DD. Overall, one can see from Figure 7 that the determined coefficients of variation yield an excellent agreement between the model distribution (histogram) and the distribution determined from the experimental data.

In Table 3, we present the mean and standard deviation for the periods under DD and LL, along with the values of the mean and standard deviation for the model parameters, .

Table 3.

Mean period (h) in constant darkness and constant light and their standard deviations (h) as calculated from the data sources in Figure 7.

	Mean	Value	Standard deviation	Value
Data	μ_DD	23.68	σ_DD	0.2718
	μ_LL	24.96	σ_LL	0.3884
Model parameter	G	48	σ_G	4.8
	b	0.54	σ_b	0.054
	β	0.005	σ_β	0.0015
	τ_x	23.68	σ_τ_x	0.2718

Relevant mean parameter values for and intrinsic period are presented, as based on a fit to the mean periods . Parameter standard deviations were derived by varying the parameter means until the simulated period standard deviation matched that of the data.

Applications

Maximal Phase Response

One use of the mouse model is to search for interesting protocols that result in large or maximal phase shifts. Previous empirical studies have investigated some aspects of inducing this behavior (Comas et al., 2006) but were restricted by the feasibility of large-scale experimentation in deducing the optimal protocol. Thus, we employed our mouse model, in silico, to design and report the optimal protocols of how to achieve maximal peak-to-peak amplitude, advance, and delay.

To determine the ideal light stimulus to induce the maximal peak-to-peak amplitude, advance, and delay, we fixed the light intensity at each chosen lux level (Table 4) and allowed for variation of light ON–OFF switching times, for pulse experiments. Our numerical simulations consisted of minimizing the negative of the absolute value of the difference between the maximum phase shift and minimum phase shift, the negative of the absolute value of the maximum of the phase advance, and the minimum of the phase delay for each response behavior, respectively.

Table 4.

Largest peak-to-peak amplitude of a phase response curve (∆ (h)), largest phase advance (+∆ (h)), and largest phase delay (–∆ (h)) as calculated through an optimization procedure of switching times for light stimuli for varying light intensities (lx).

Intensity (lx)	1 pulse			2 pulses
	Δ ∆ (h)	+Δ (h)	−Δ (h)	Δ ∆ (h)	+Δ (h)	−Δ (h)
100	5.60	2.19	3.49	7.34	3.30	4.19
200	5.74	2.25	3.57	7.60	3.42	4.33
300	5.81	2.28	3.61	7.71	3.47	4.39
400	5.84	2.29	3.63	7.77	3.50	4.42
500	5.87	2.30	3.64	7.82	3.53	4.45
600	5.88	2.31	3.65	7.85	3.54	4.46
700	5.90	2.32	3.66	7.87	3.55	4.47
800	5.91	2.32	3.67	7.89	3.56	4.48
900	5.92	2.33	3.67	7.91	3.57	4.49
1000	5.92	2.33	3.68	7.92	3.57	4.50

The optimal protocol for inducing the maximal amplitude, phase advance, and phase delay was found to be a 2-pulse light stimulus, with the ON–OFF switching times determined to be [0, 7.7, 22.3, 24], [0, 9.2, 21.9, 24], and [0, 7.1, 23.3, 24], respectively. For a single pulse, the maximal amplitude, advance, and delay can be induced by durations of 8.6, 10.5, and 7.8 h, respectively.

Our convention for the ON–OFF switching times is that light is ON at 0, OFF at the subsequent time, and this alternates as ON/OFF until the end of the stimulus where the light must be switched off so the simulation can be run in DD.

We restrict the maximal possible length of time of the stimulus to be 24 h. Within this time, we are able to provide any kind of stimulus (e.g., 2 pulses, 5 pulses). Note that this linear representation is not meant to be reflective of the circular clock but instead is a mechanism of representing multiple pulses of specific lengths and separations, which start from determined appropriate CTs in PRC experiments. As a concrete example, let us assume that the experimentalist constructing PRCs has determined the appropriate CTs. Then, say at CT 0, they implement our protocol, . This involves a pulse of length  h, followed by darkness of length  h, and finally another pulse of  h, after which the light is shut off, which will allow the animal to free-run in DD.

For a number of pulses , the simulations converged to a 2-pulse protocol. The single pulse consistently performed worse in inducing maximal phase shifts. The best protocol for inducing maximal amplitude, advance, and delay was determined to be a 2-pulse light stimulus with ON–OFF switching times of [0, 7.7, 22.3, 24], [0, 9.2, 21.9, 24], and [0, 7.1, 23.3, 24], respectively. We repeated this for a range of light intensities, from 100 to 1000 lx. However, past 400 lx, the simulations were purely for expository purposes. Typical mouse circadian experiments do not use very high light intensities (Eckel-Mahan and Sassone-Corsi, 2015). The maximal shifts are provided in Table 4, where the best shifts possible due to a single pulse of light versus 2 pulses of light are compared. For a single pulse, the switching times for maximal amplitude, advance, and delay were determined to be [0, 8.6], [0, 10.5], and [0, 7.8], respectively. Note that our convention for ON–OFF switching times always starts at 0; for experimental purposes, the 0 would be the location of every appropriate CT that the light pulse is administered. Furthermore, we emphasize that the points and are not points on the clock, but reflect how long the stimulus experiment should be: For example, in , the represents the start of the pulse and is on until have elapsed, followed by of DD, and then finally a short pulse of , after which the light is shut off and the mouse free-runs in DD.

It has been reported in previous studies on C57BL/6J-OlaHsd mice (Comas et al., 2006) and on humans (Dewan et al., 2011) that increasing duration is more effective in inducing phase changes of the circadian clock than that of increasing light intensity. We can see a similar result in our mouse model (Figure 8): the effect of light intensity on inducing phase shifts does cause some noticeable change, but the duration, and more specifically the pattern of light stimulus, is more effective.

Figure 8.

Phase response curves for varying light intensities using the optimal protocols (given above each panel) for (a) maximum advance and for (b) maximum delay. The optimal protocols for inducing maximal advance and maximal delay are the sets of ON–OFF switching times [0, 9.2, 21.9, 24] and [0, 7.1, 23.3, 24], respectively. Abbreviation: CT = circadian time.

The mouse study (Comas et al., 2006) only examined experimentally different durations of 100 lx single pulses and concluded the optimal duration for inducing the maximal amplitude is 9 h. Our model confirms that around 8.5 to 9 h indeed result in the maximal amplitude for a single pulse; however, it instead suggests that the 2-pulse protocol [0, 7.1, 23.3, 24] is more effective overall.

Comparison With a Mutant Mouse

We emphasize that the fit model parameters (Table 1) are specifically for the C57BL/6 mouse—Other mice or mutant strains that alter or possibly influence the circadian clock would not correlate well with the model parameters. One example of this is the casein kinase 1 epsilon deficient (CK1 ) mouse, as reported in Etchegaray et al. (2009). In their study, males with a floxed allele and the Protamine-Cre transgene were backcrossed to the C57BL/6J background for at least 10 generations. Figure 9 shows the difference between the original parameters fit for the C57BL/6J mouse versus the CK1 mouse. It is clear that the original parameters do not yield an appreciable fit to the PRC data. While the new parameters seem to perform worse near the phase delay region, we note that the error bars are significantly larger. The parameters determined were , , , , and . Given the lack of advance and delay regions compared with the C57BL/6J mouse, it is unsurprising that the parameters , responsible for the photic drive, were significantly lower for the CK1 mouse. Parameter , responsible for the circadian modulation, is larger—This parameter predominantly controls the size of the dead and delay zones (see “Discussion” section). Since the dead zone encompasses most of the PRC, it is reasonable to expect a larger . The parameter is typically responsible for the effects of duration versus intensity of light. The effects of the difference between the 2 values here are not significant, but the minimizer determined to be the optimal value. It is clear from this example that mutations or strain differences require the identification of specific parameter changes that accommodate the inborn alterations to reproduce and predict outcomes.

Figure 9.

Phase response curve comparison between experimental data from Etchegaray et al. (2009) for the CK1 mouse with our re-fit model parameters and the original optimal parameters for the C57BL/6 mouse. Abbreviation: CT = circadian time.

Discussion

Our results have revealed that there is a clear difference between the parameters responsible for photic input in humans compared with the C57BL/6 mouse. This discrepancy in parameter values between mice and humans can be explained as follows. Experiments involving nocturnal rodents, especially mice, use low-intensity light because of their susceptibility to retinal damage at high light intensities compared with humans (Gonzalez, 2018). Nocturnal rodents are also significantly more sensitive to low-intensity light than humans, avoiding dim lights and being able to condition to low-intensity light signaling (Eckmier et al., 2016). While it has been shown that human circadian systems also undergo phase shifts to lower light intensities (e.g., 180 lx; Boivin et al., 1996), the phases only shift with a mean of ∼∼ . For similar light intensities, the mouse circadian pacemaker phase shifts by significantly greater amounts. This is apparent from the experimental PRCs in Figures 1 and 2, where there are substantial phase shifts to light intensities between 150 and 400 lx. Moreover, it can be seen from human PRCs (Boivin et al., 1996; St Hilaire et al., 2012; Khalsa et al., 2003) that similar significant phase shift magnitudes occur only with markedly greater intensities between and orders of magnitude. Thus, we expect that due to the heightened sensitivity of the mouse to low-intensity light, the parameters and , responsible for effects of photic drive, should also be higher. As discussed in the previous analysis, is responsible for the effects of pulse duration versus light intensity. A larger causes a variation in low light intensities to result in greater variation in phase response. For larger intensities, the difference between phase response as a function of lx is less pronounced, as can be seen in Figure 8: as the intensity is increased, the effect on the phase shift is less and less noticeable.

We admit that the lux unit can be problematic because it is based on the perceived brightness correlated to the human visual system (Peirson et al., 2018), and thus lends to potential issues in experimental data (Alves-Simoes et al., 2016; Hofstetter et al., 2005) as discussed prior. However, the use of lux has been accepted as a standard in the field due to the wide availability of and easy access to lux meters. Furthermore, above saturation, the difference between rhodopsin and melanopsin spectral sensitivities should not affect the interpretation of experiments. Regardless, our model takes into account the effects of different types of artificial light and provides a range of suitable parameters using parameter variation that reliably reproduces the characteristic data and behavior.

We note that the parameter k was slightly lower in mice than in humans. Previous authors studying human PRCs determined that the parameter was required to shorten the period of the pacemaker. However, in our study, we determined that with for nocturnal organisms, light lengthens the period of the pacemaker in accordance with Aschoff’s rules. Specifically, increasing the intensity of light increases the free-running period in LL. For our set of parameters with a fixed light intensity, as increases, so too does the free-running period under LL. The suitable value of that matches the experimental period, however, turned out to be k ∼∼ .

The circadian modulation, communicated through parameter , should be larger than that of the human model to correlate with the stronger modulation of sensitivity in mice during subjective day—the so-called dead zone, where phase shifts are rather negligible. During subjective day, when is positive, a larger value of causes the drive strength to be reduced further. Moreover, during subjective night, both and are negative, forcing the drive strength to be greater. This modulation also impacts the location of maximal sensitivity to photic stimuli, which occurs between CT 15 and CT 17. Examining Figure 2 indicates that this larger choice of indeed fits both the dead zone and the location of maximal delay more accurately than that of the human model parameters.

We have presented a mouse model of the original human model developed by Kronauer (1990). The model parameters were fit to available mouse data from the literature. The model was validated against other sets of experimental data and compared against the human model. For the same protocol, the mouse parameters yielded a significantly better match with the data than that of the human parameters. Sensitivity analysis conducted revealed that the measured period is only affected by perturbations under LL lighting conditions. Due to the reported variability of periods across and within labs, the most relevant parameters affecting the period were varied and compared with data from a wide range of sources to deduce a range of values for the key parameters. Using the mouse model to search for the optimal protocol for inducing maximal amplitude, advance, and delay determined that a 2-pulse light stimulus of switching times [0, 7.7, 22.3, 24], [0, 9.2, 21.9, 24], and [0, 7.1, 23.3, 24], respectively, was best. To further validate this model, we propose that future studies attempt other light protocols, especially the 2-pulse protocol as suggested by the optimization. The application of the mouse model to the CK1 mutant revealed that the parameters are indeed mouse-strain specific. The mouse model may be adjusted accordingly based on general physiology when attempting to study other strains.