## Abstract

Bioluminescence rhythms from cellular reporters have become the most common method used to quantify oscillations in circadian gene expression. These experimental systems can reveal phase and amplitude change resulting from circadian disturbances, and can be used in conjunction with mathematical models to lend further insight into the mechanistic basis of clock amplitude regulation. However, bioluminescence experiments track the mean output from thousands of noisy, uncoupled oscillators, obscuring the direct effect of a given stimulus on the genetic regulatory network. In many cases, it is unclear whether changes in amplitude are due to individual changes in gene expression level or to a change in coherence of the population. Although such systems can be modeled using explicit stochastic simulations, these models are computationally cumbersome and limit analytical insight into the mechanisms of amplitude change. We therefore develop theoretical and computational tools to approximate the mean expression level in large populations of noninteracting oscillators, and further define computationally efficient amplitude response calculations to describe phase-dependent amplitude change. At the single-cell level, a mechanistic nonlinear ordinary differential equation model is used to calculate the transient response of each cell to a perturbation, whereas population-level dynamics are captured by coupling this detailed model to a phase density function. Our analysis reveals that amplitude changes mediated at either the individual-cell or the population level can be distinguished in tissue-level bioluminescence data without the need for single-cell measurements. We demonstrate the effectiveness of the method by modeling experimental bioluminescence profiles of light-sensitive fibroblasts, reconciling the conclusions of two seemingly contradictory studies. This modeling framework allows a direct comparison between in vitro bioluminescence experiments and in silico ordinary differential equation models, and will lead to a better quantitative understanding of the factors that affect clock amplitude.

## Introduction

In mammals, circadian rhythms are endogenous oscillations in gene transcription responsible for coordinating daily changes in physiology. Although the suprachiasmatic nucleus (SCN) in the brain serves as the body’s master pacemaker, cells found in peripheral tissues also oscillate in a circadian manner (

1

). These peripheral clocks process systemic and SCN-mediated entraining cues, buffering against rhythmic changes in energy availability to maintain metabolic homeostasis (2

). The amplitude of circadian transcription is a relevant factor, and has been shown to play a critical role in phase resetting and entrainment (3

, 4

). Recent studies have further highlighted the importance of high peripheral clock amplitudes in maintaining metabolic health. Mice lacking an intact clock have been shown to develop metabolic disease (5

), while low-amplitude clock oscillations, whether caused by diet (6

) or age (7

), have also been tied to metabolic disorders. An understanding of how clock amplitudes are regulated is therefore a topic of ongoing research with potential therapeutic applications (8

). Unlike the network of oscillators in the SCN, in which intercellular coupling maintains robust amplitudes even in the absence of external cues, peripheral oscillators are thought to lack a direct mechanism to spontaneously synchronize (9

). As a result, populations of peripheral clocks are likely synchronized by common external cues, with stochastic effects and cell heterogeneity driving entrained populations gradually toward desynchrony.The development of immortalized peripheral oscillator cell lines with bioluminescent reporters has allowed high-throughput analysis of the responses of circadian rhythms to genetic and pharmacological manipulation (

10

, 11

). In addition to providing experimental tractability, these systems provide detailed information on the amplitude of oscillations in gene transcription, a measure often lacking from earlier experiments using wheel-running activity. As a result, these cell lines have proven useful in studying core clock connectivity and stoichiometry (12

). However, because in vitro experiments typically measure entire cultures of cells, data collected at the population-level can obscure the response of the clock at the scale of the gene-regulatory network. Even when individual cells can be recorded, stochastic noise hinders accurate amplitude determination. As a result, two studies using similar perturbations to understand the mechanism of light-induced amplitude reduction reached differing conclusions, in which either single-cell amplitude reduction or population-level desynchrony was identified as the dominant factor (13

, 14

).Mathematical models have long been used to understand the results of circadian experiments (

15

, 16

), aided by definitions and computational techniques designed to match modeling predictions to experimental data. One such definition is the response function, a general technique that maps a change in an output variable to a temporary change in parameters (17

). For instance, the phase response curve (PRC) has been used to characterize the entrainment behavior of both experimental and mathematical systems (18

, 19

, 20

), and in analyzing the synchrony of populations of oscillators (21

). Accurate and efficient numerical routines for finding infinitesimal PRCs have therefore been developed (22

, 23

, 24

). In addition to changes in phase, phase-dependent changes in amplitude are an important factor in understanding circadian rhythms. Amplitude response curves (ARCs) were first used in the clock literature to comprehend the effects of light pulses on simple phase-amplitude models, and were useful in predicting phase singularity behavior (25

). ARCs have also been used to characterize perturbations to groups of oscillators through desynchrony (14

, 26

), at the level of the single oscillator (27

, 28

), and even in studying the effect of entrainment phase on SCN rhythm amplitude (29

). However, previous definitions of the ARC are inconsistent between studies and do not simultaneously consider amplitude effects at the single-cell and population level.The choice of modeling framework dictates the type of amplitude response that can be predicted. Ordinary differential equation (ODE) models of gene regulation are capable of describing the amplitude and phase-resetting behavior of single cells, but fail to capture the collective dynamics of a population of oscillators. Likewise, phase-only models correctly capture the change in synchrony of a population, but do not capture fluctuations in amplitude of individual oscillators. Explicit stochastic simulations of populations of cells are capable of realistically capturing both single-cell and population-level effects, and have been successfully used to understand the response of coupled oscillators to external VIP perturbation (

30

). However, these methods are computationally expensive, and cannot be used for an analytical understanding of amplitude response.In this study, we describe approaches to quantify amplitude change in a population of noninteracting oscillators. By exploiting the independence of each oscillator, we derive computationally efficient methods to approximate the mean dynamics of full stochastic simulations. Additionally, our method allows the calculation of ARCs at both the single-cell and population level, allowing the behavior of the system to be quickly profiled. Specifically, we use ODE models to describe the transient amplitude response at the single-cell level, coupled with a phase probability density function to describe population-level dynamics. After a perturbation with a finite duration, a limit cycle oscillator undergoes a transient change in amplitude. When the mean expression level of a population of oscillators is also considered, an additional change in amplitude is incurred due to the change in synchrony of the population, which persists until synchrony is changed by subsequent perturbations. We therefore observe a separation in timescales between the effects on clock output mediated at the single-cell level and those mediated by population synchrony, allowing the source of an amplitude change to be qualitatively inferred by inspection of bioluminescence data. Understanding the mechanisms and consequences of both types of amplitude regulation will be important in understanding how peripheral amplitudes are maintained, and may lead to the design of pharmacological or behavioral strategies to boost circadian amplitudes.

## Materials and Methods

### Basic definitions

ODE models take the general form

in which

The period is the smallest

The phase variable

$\frac{dx}{dt}=f\left(x\left(t\right),p\right),$

(1)

in which

*x*(*t*) represents the concentrations of the state variables, such as mRNA and protein concentrations;*f*contains information on the production, degradation, and reactivity of the states; and*p*represents the kinetic parameters governing reaction kinetics. Limit cycle models are ODE models in which the solution approaches a steady-state oscillatory trajectory, satisfying:$\underset{t\to \infty}{\text{lim}}\left[x\left(t+T\right)-x\left(t\right)\right]=0.$

(2)

The period is the smallest

*T*> 0 for which Eq. 2 holds. The points on the stable limit cycle are denoted by*x*^{γ}(*θ*), with each point assigned to a value of a phase variable*θ*∈[0,2*π*). For convenience, time in Eq. 1 can be rescaled such that the period is 2*π*:$\tilde{t}=\frac{2\pi}{T}t;\phantom{\rule{1em}{0ex}}\tilde{f}=\frac{T}{2\pi}\phantom{\rule{0.25em}{0ex}}f;\phantom{\rule{1em}{0ex}}\frac{dx}{d\tilde{t}}=\tilde{f}\left(x\left(\tilde{t}\right),p\right).$

(3)

The phase variable

*θ*is therefore defined on the limit cycle as $\theta =\tilde{t}\phantom{\rule{0.25em}{0ex}}\text{mod}\phantom{\rule{0.25em}{0ex}}2\pi $, with*θ*= 0 assigned to a unique and identifiable point.### Perturbations to limit cycle systems

In this study we restrict our analysis to temporary perturbations capable of entraining an oscillatory system, excluding permanent parameter changes arising, for instance, from genetic knockout experiments. The simplest entraining perturbation involves adding or removing components to a limit cycle system, resulting in a perturbed trajectory

This trajectory evolves according to Eq. 1, eventually returning to

In addition to perturbations directly to the state of the system, oscillators can also be perturbed by temporary changes to the parameters. For a parameter pulse Δ

The pulse trajectory ${x}_{\tilde{d}}\left(\tilde{t}\right)$ is then integrated from $\tilde{t}=-\tilde{d}\to 0$, at which point the pulse is removed. A temporary parameter pulse is therefore equivalent to a state perturbation, but with the perturbation at

A schematic depicting a perturbed and reference trajectory exemplar is shown in Fig. 1

*x*(*t*). Here the initial conditions are determined by the strength of the perturbation and the phase at which it is applied:$x\left(0\right):={x}^{\gamma}\left({\theta}_{0}\right)+\Delta x\left(0\right).$

(4)

This trajectory evolves according to Eq. 1, eventually returning to

*x*^{γ}. It is useful to express this trajectory by the deviation from the limit cycle:$\Delta x\left(\tilde{t}\right)=x\left(\tilde{t}\right)-{x}^{\gamma}\left(\tilde{t}+{\theta}_{0}\right).$

(5)

In addition to perturbations directly to the state of the system, oscillators can also be perturbed by temporary changes to the parameters. For a parameter pulse Δ

*p*of duration $\tilde{d}$ (in radians) that ends at*θ*_{0}, the oscillator deviates from the limit cycle trajectory according to the following:${x}_{\tilde{d}}\left(-\tilde{d}\right)={x}^{\gamma}\left({\theta}_{0}-\tilde{d}\right),$

(6)

$\frac{d{x}_{\tilde{d}}}{d\tilde{t}}=\tilde{f}\left({x}_{\tilde{d}}\left(\tilde{t}\right),p+\Delta p\right).$

(7)

The pulse trajectory ${x}_{\tilde{d}}\left(\tilde{t}\right)$ is then integrated from $\tilde{t}=-\tilde{d}\to 0$, at which point the pulse is removed. A temporary parameter pulse is therefore equivalent to a state perturbation, but with the perturbation at

*t*= 0 defined by$\Delta x\left(0\right)={x}_{\tilde{d}}\left(0\right)-{x}^{\gamma}\left({\theta}_{0}\right).$

(8)

A schematic depicting a perturbed and reference trajectory exemplar is shown in Fig. 1

*A*.The response of the system to very small perturbations is more efficiently calculated using ODE sensitivity analysis (

can be calculated directly by integrating

with

31

). The sensitivity matrix,$S\left(t\right)=\frac{dx\left(t\right)}{dx\left(0\right)}=\underset{\Delta x\left(0\right)\to 0}{\text{lim}}\frac{\Delta x\left(t\right)}{\Delta x\left(0\right)},$

(9)

can be calculated directly by integrating

$\frac{d}{dt}S\left(t\right)=\frac{df\left(x\left(t\right),p\right)}{dx}\phantom{\rule{0.25em}{0ex}}S\left(t\right)$

(10)

with

*S*(0) =*I*.### Phase response curves

Quantifying phase changes after a perturbation has been well studied and is particularly relevant for models of circadian rhythms (

This result follows from the fact that in the limit of an infinitely short and small parameter pulse,

22

, 24

). A phase response curve (PRC) maps the change in phase resulting from the same perturbation applied at each initial phase. Infinitesimal PRCs, the derivative of the phase change with respect to the perturbation, can be defined for state and parameter-impulse perturbations (24

). Methods for efficiently calculating these quantities using ODE sensitivity analysis have been developed (24

), with the important result that the parameter- and state-impulse PRCs can be related by the following Jacobian matrix:$\frac{d}{d\tilde{t}}\frac{d\theta}{dp}=\frac{d\theta}{dx}\frac{d\tilde{f}}{dp}.$

(11)

This result follows from the fact that in the limit of an infinitely short and small parameter pulse,

$\Delta x\left(0\right)\to \frac{d\tilde{f}}{dp}\phantom{\rule{0.25em}{0ex}}\tilde{d}\phantom{\rule{0.25em}{0ex}}\Delta p.$

### Phase-diffusion model

Large populations of oscillators are typically described using phase-only models (

The shape of $p\left(\theta ,\tilde{t}\right)$ changes as the cells advance in time. Stochastic effects cause the population to gradually desynchronize as slight cycle-to-cycle deviations are propagated throughout the population (

Due to the rescaling of

The solution of Eqs. 13–16 is well characterized, with $p\left(\theta ,\tilde{t}\right)$ evolving in time as the convolution of the initial conditions with a wrapped normal distribution with mean $\tilde{t}$ and standard deviation $\sqrt{2d\tilde{t}}$ (

in which the wrapped normal distribution () is defined as

Because the convolution of two normal distributions is also a normal distribution, it is efficient when possible to describe

32

), in which the state of each oscillator is represented only by its phase, *θ*. The synchrony of the population can be modeled using a probability density function $p\left(\theta ,\tilde{t}\right)$ that describes the probability of finding an oscillator at each phase (33

). The usefulness of probability density functions in describing the phase and amplitude responses of populations of circadian cells has been previously shown (14

). As with all probability density functions,${\int}_{0}^{2\pi}p\left(\theta ,\tilde{t}\right)\phantom{\rule{0.25em}{0ex}}d\theta =1.$

(12)

The shape of $p\left(\theta ,\tilde{t}\right)$ changes as the cells advance in time. Stochastic effects cause the population to gradually desynchronize as slight cycle-to-cycle deviations are propagated throughout the population (

34

). For an infinite population of oscillators, these effects are well described by a Fokker-Plank equation (35

):$\frac{\partial p}{\partial \tilde{t}}=\frac{\partial p}{\partial \theta}+d\frac{{\partial}^{2}p}{\partial {\theta}^{2}}.$

(13)

Due to the rescaling of

*t*, the mean period of the population is 2*π*. Here, the*∂p*/*∂θ*term, analogous to convection, describes the mean oscillatory period, whereas the*∂*^{2}*p*/*∂θ*^{2}term describes the diffusion of phases across [0,2*π*). The phase diffusivity parameter*d*(in units of inverse radians) describes the speed with which the population desynchronizes and can be fit to experimental data (36

). Equation 13 has periodic boundary conditions, with initial condition *ϕ*(*θ*) as the phase population at*t*= 0:$\text{BCs}:\phantom{\rule{1em}{0ex}}p\left(0,\tilde{t}\right)=p\left(2\pi ,\tilde{t}\right),$

(14)

$\frac{\partial p}{\partial \theta}\left(0,\tilde{t}\right)=\frac{\partial p}{\partial \theta}\left(2\pi ,\tilde{t}\right),$

(15)

$\text{IC}:\phantom{\rule{1em}{0ex}}p\left(\theta ,0\right)=\varphi \left(\theta \right).$

(16)

The solution of Eqs. 13–16 is well characterized, with $p\left(\theta ,\tilde{t}\right)$ evolving in time as the convolution of the initial conditions with a wrapped normal distribution with mean $\tilde{t}$ and standard deviation $\sqrt{2d\tilde{t}}$ (

37

):$p\left(\theta ,\tilde{t}\right)=\varphi \left(\theta \right)\ast \mathcal{W}\mathcal{N}\left(\theta ;\tilde{t},\sqrt{2\phantom{\rule{0.25em}{0ex}}d\phantom{\rule{0.25em}{0ex}}\tilde{t}}\right),$

(17)

in which the wrapped normal distribution () is defined as

$\mathcal{W}\mathcal{N}\left(\theta ;\mu ,\sigma \right)=\frac{1}{\sigma \sqrt{2\pi}}\sum _{k=-\infty}^{\infty}\text{exp}\left[\frac{-{\left(\theta -\mu +2\pi k\right)}^{2}}{2{\sigma}^{2}}\right].$

(18)

Because the convolution of two normal distributions is also a normal distribution, it is efficient when possible to describe

*ϕ*(*θ*) as a normal distribution with mean*μ*_{0}and standard deviation*σ*_{0}, such that $p\left(\theta ,\tilde{t}\right)$ can be found analytically through$p\left(\theta ,\tilde{t}\right)=\mathcal{W}\mathcal{N}\left(\theta ;{\mu}_{0}+\tilde{t},\sqrt{{\sigma}_{0}^{2}+2{d}^{2}{\tilde{t}}^{2}}\right).$

(19)

### Numerical simulations

ODE models were simulated in the software PYTHON (Python Software Foundation, https://www.python.org/psf/), using the computer algebra package CASADI (https://github.com/casadi/casadi/wiki) (

39

). Stochastic simulations were performed using the STOCHKIT2 package (http://sourceforge.net/projects/stochkit/) (40

). The codes used to generate the figures are included in the Supporting Material.## Results and Discussion

A signal to a population of limit cycle oscillators can affect amplitude in two ways:

- 1.Individual cells are perturbed from their limit cycles, and exhibit transient dynamics before settling back to steady-state amplitudes; and
- 2.The phases of the population are changed, resulting in a permanent change in population synchrony.

### Definition of an amplitude metric at the single-cell level

After a perturbation to a single cell, the path that the perturbed trajectory

The amplitude change metric is defined as

This amplitude metric was chosen over alternatives, such as peak-trough distance, for its analytical tractability. The integrand in Eq. 21, abbreviated by

*x*(*t*) takes in returning to the limit cycle will have a different amplitude than the unperturbed trajectory. Although such an amplitude change will only have a finite duration, it plays an important role when perturbations are repeatedly received by the clock, such as a peripheral oscillator entrained to daily metabolic stimuli. To define an amplitude change metric for such a case, we compare a perturbed trajectory*x*(*t*) to a phase-shifted limit cycle reference*y*(*t*), for which*x*(*t*) →*y*(*t*) for sufficiently long times. Because*x*(*t*) approaches the reference as*t*→ ∞, the means of both trajectories are equal and can be calculated by$\mu :={\int}_{0}^{2\pi}\frac{{x}^{\gamma}\left(\theta \right)}{2\pi}\phantom{\rule{0.25em}{0ex}}d\theta .$

(20)

The amplitude change metric is defined as

$\begin{array}{c}\Delta A\left(x\left(t\right),y\left(t\right)\right):={\int}_{0}^{\infty}{\left(x\left(t\right)-\mu \right)}^{2}-{\left(y\left(t\right)-\mu \right)}^{2}\phantom{\rule{0.25em}{0ex}}dt\\ ={\int}_{0}^{\infty}h\left(t\right)\phantom{\rule{0.25em}{0ex}}dt.\end{array}$

(21)

This amplitude metric was chosen over alternatives, such as peak-trough distance, for its analytical tractability. The integrand in Eq. 21, abbreviated by

*h*(*t*), compares the variances of the reference and perturbed trajectories. When*h*(*t*) > 0, the trajectory is further from the mean than the reference, and similarly when*h*(*t*) < 0, the trajectory is closer to the mean than the reference. Thus the overall amplitude change can be calculated by integrating*h*(*t*) until the two trajectories converge, returning an amplitude value for each state variable. Amplitude change for a two-dimensional oscillator is easy to visualize graphically. In Fig. 1,*B*–*D*, the same state perturbation Δ*x*(0) applied at two different phases results in opposite amplitude changes, depending on whether the perturbation shifts the trajectory to the interior of the limit cycle (reduced amplitudes) or to the outside the limit cycle (increased amplitudes). Trajectories for the second state variable, $Y\left(\tilde{t}\right)$, and the corresponding integrand for the amplitude change equation, $h\left(\tilde{t}\right)$, demonstrate how this transient change is quantified.### Single-cell amplitude response curves

Similar to the PRC, we denote the phase-dependent amplitude change after a perturbation as an amplitude response curve (ARC), using the metric presented in Eq. 21. In Fig. 2, we calculate the PRC and ARC for a state perturbation of three different strengths. The weakest perturbation results in Type-1 (weak) resetting, in which the PRC is continuous, whereas the strongest perturbation results in type 0 (strong) resetting, in which the PRC is discontinuous. This transition occurs once the state perturbation is strong enough to push the trajectory over the unstable fixed point located at the middle of the limit cycle. For a perturbation of intermediate strength, at a critical phase the trajectory will be pushed close to the fixed point and take a long time to recover the steady-state amplitude. The behavior at this singularity point is indicated by a sharp dip in the ARC for this perturbation strength.

Infinitesimal versions of response curves are often more general and easier to compute than those that track a specific perturbation strength (

Although this quantity could be calculating by using a very small Δ

Simplifying the integrand in Eq. 22,

Taking the limit of this integrand as Δ

Note that in Eq. 29,

Here, the first term of the integrand $\left(S-\tilde{f}\dot{\theta}\right)$ tracks the distance from the perturbed trajectory to the limit cycle, which decays to zero as

and may be calculated from the state-impulse version with the following relationship:

As with Eq. 11, this equivalency reflects the fact that for a pulse of infinitely short duration, a parameter change is equivalent to changing the state of the system along the direction specified by the Jacobian. Convergence between the methods in Eqs. 30–32 and finite-difference approaches is shown in Fig. S1 in the Supporting Material, which demonstrates that the numerically efficient differential ARCs remain representative even for moderate-strength perturbations.

17

). We therefore derive an expression for the infinitesimal ARC, defined as$\frac{dA}{dx}:=\underset{\Delta x\left(0\right)\to 0}{\text{lim}}\frac{\Delta A\left(x\left(\tilde{t}\right),\phantom{\rule{0.25em}{0ex}}{x}^{\gamma}\left(\tilde{t}+{\theta}_{0}+\Delta \theta \right)\right)}{\Delta x\left(0\right)}={\int}_{0}^{\infty}\underset{\Delta x\left(0\right)\to 0}{\text{lim}}\frac{h\left(t\right)}{\Delta x\left(0\right)}\phantom{\rule{0.25em}{0ex}}dt.$

(22)

Although this quantity could be calculating by using a very small Δ

*x*(0), it is more accurate and efficient to derive a direct method to calculate Eq. 22 using the ODE sensitivities defined in Eq. 10. For simplicity, we define ${t}_{\theta}=\tilde{t}+{\theta}_{0}$, which allows the time variable for the perturbation to vary from 0 → ∞ while tracking the appropriate phase on the limit cycle. Because Δ*θ*→ 0 as Δ*x*(0) → 0, we Taylor-expand the limit cycle trajectory around*t*_{θ}:${x}^{\gamma}\left({t}_{\theta}+\Delta \theta \right)={x}^{\gamma}\left({t}_{\theta}\right)+\frac{d{x}^{\gamma}\left({t}_{\theta}\right)}{d\theta}\Delta \theta +O\left(\Delta {\theta}^{2}\right)$

(23)

$={x}^{\gamma}\left({t}_{\theta}\right)+\tilde{f}\left({x}^{\gamma}\left({t}_{\theta}\right)\right)\Delta \theta +O\left(\Delta {\theta}^{2}\right).$

(24)

Simplifying the integrand in Eq. 22,

$\frac{h\left(t\right)}{\Delta x\left(0\right)}=\frac{1}{\Delta x\left(0\right)}\left[{\left({x}^{\gamma}\left({t}_{\theta}\right)+\Delta x\left(\tilde{t}\right)-\mu \right)}^{2}-{\left({x}^{\gamma}\left({t}_{\theta}+\Delta \theta \right)-\mu \right)}^{2}\right]$

(25)

$=\frac{1}{\Delta x\left(0\right)}\left[{\left({x}^{\gamma}\left({t}_{\theta}\right)+\Delta x\left(\tilde{t}\right)-\mu \right)}^{2}-{\left({x}^{\gamma}\left({t}_{\theta}\right)+\tilde{f}\left({x}^{\gamma}\left({t}_{\theta}\right)\right)\Delta \theta -\mu \right)}^{2}\right]$

(26)

$=\frac{1}{\Delta x\left(0\right)}\left[\left(\Delta x\left(\tilde{t}\right)-\tilde{f}\left({x}^{\gamma}\left({t}_{\theta}\right)\right)\Delta \theta \right)\left(\Delta x\left(\tilde{t}\right)+\tilde{f}\left({x}^{\gamma}\left({t}_{\theta}\right)\right)\Delta \theta +2\left({x}^{\gamma}\left({t}_{\theta}\right)-\mu \right)\right)\right].$

(27)

Taking the limit of this integrand as Δ

*x*(0) → 0 cancels several differential terms, and allows the remainder to be substituted with quantities that can be calculated using ODE sensitivity analysis:$\underset{\Delta x\left(0\right)\to 0}{\text{lim}}\frac{h\left(t\right)}{\Delta x\left(0\right)}=2\left[\right(\underset{\Delta x\left(0\right)\to 0}{\text{lim}}\frac{\Delta x\left(\tilde{t}\right)}{\Delta x\left(0\right)})-\tilde{f}\left({x}^{\gamma}\left({t}_{\theta}\right)\right)(\underset{\Delta x\left(0\right)\to 0}{\text{lim}}\frac{\Delta \theta}{\Delta x\left(0\right)})]\left({x}^{\gamma}\left({t}_{\theta}\right)-\mu \right)$

(28)

$=2\left(S\left(\tilde{t}\right)-\tilde{f}\left({x}^{\gamma}\left({t}_{\theta}\right)\right)\frac{d\theta}{dx}\right)\left({x}^{\gamma}\left({t}_{\theta}\right)-\mu \right).$

(29)

Note that in Eq. 29,

*dθ*/*dx*represents the derivative of Δ*θ*with respect to the perturbation, and is therefore a scalar quantity. The infinitesimal state-impulse ARC may therefore be calculated directly from the ODE sensitivity matrix and the PRC, allowing the amplitude change for an infinitesimal perturbation to be calculated exactly and efficiently:$\frac{dA}{dx}={\int}_{0}^{\infty}2\left(S\left(\tilde{t}\right)-\tilde{f}\left({x}^{\gamma}\left({t}_{\theta}\right)\right)\frac{d\theta}{dx}\right)\left({x}^{\gamma}\left({t}_{\theta}\right)-\mu \right)\phantom{\rule{0.25em}{0ex}}dt.$

(30)

Here, the first term of the integrand $\left(S-\tilde{f}\dot{\theta}\right)$ tracks the distance from the perturbed trajectory to the limit cycle, which decays to zero as

*t*→ ∞. The second term (*x*^{γ}–*μ*) weights this distance by whether (or not) the deviation occurs above or below the oscillatory mean, yielding negative amplitude changes when the trajectory is perturbed closer to the mean. Just as with the parameter-impulse PRC, the infinitesimal parameter-impulse ARC is defined as$\frac{d}{d\tilde{t}}\frac{dA}{dp}:=\underset{\tilde{d},\phantom{\rule{0.25em}{0ex}}\Delta p\to 0}{\text{lim}}\frac{\Delta A}{\tilde{d}\phantom{\rule{0.25em}{0ex}}\Delta p},$

(31)

and may be calculated from the state-impulse version with the following relationship:

$\frac{d}{d\tilde{t}}\frac{dA}{dp}=\frac{dA}{dx}\frac{d\tilde{f}}{dp}.$

(32)

As with Eq. 11, this equivalency reflects the fact that for a pulse of infinitely short duration, a parameter change is equivalent to changing the state of the system along the direction specified by the Jacobian. Convergence between the methods in Eqs. 30–32 and finite-difference approaches is shown in Fig. S1 in the Supporting Material, which demonstrates that the numerically efficient differential ARCs remain representative even for moderate-strength perturbations.

### Population-level response curves

We next briefly describe how to calculate the PRC and ARC for a phase-density model. These approaches are commonly used in understanding phase models (

However, it is easier to estimate the mean and standard deviation of the perturbed population $\stackrel{\u02c6}{p}\left(\theta ,\tilde{t}\right)$ using directional statistics (). A population defined on the unit circle can be described by a complex variable $z=\rho {e}^{i\overline{\theta}}$, where $\overline{\theta}$ is the mean phase and

In Eq. 35, we avoid calculating the perturbed population explicitly by instead integrating over the new phases at the prior population density function. Population-level amplitude and phase responses can therefore be calculated by

Population-level PRCs and ARCs can be tabulated by solving Eqs. 34–37 for populations $p\left(\theta ,\tilde{t}\right)$ with different mean phases. It is important to note that the ARC at the population level strongly depends on the slope of the PRC, as has been shown previously (

14

, 33

), and are presented here to match previous definitions for single cells. A phase transition curve, g(*θ*) =*θ*+ Δ*θ*, maps the phase of an oscillator before perturbation to its phase after the perturbation. Because individual oscillators are neither created nor destroyed during the perturbation, it is possible—yet numerically difficult—to directly calculate the phase probability distribution after perturbation using the standard change of variables relation:$\stackrel{\u02c6}{p}\left(\theta ,\tilde{t}\right)\phantom{\rule{0.25em}{0ex}}dg\left(\theta \right)=p\left(\theta ,\tilde{t}\right)\phantom{\rule{0.25em}{0ex}}d\theta .$

(33)

However, it is easier to estimate the mean and standard deviation of the perturbed population $\stackrel{\u02c6}{p}\left(\theta ,\tilde{t}\right)$ using directional statistics (). A population defined on the unit circle can be described by a complex variable $z=\rho {e}^{i\overline{\theta}}$, where $\overline{\theta}$ is the mean phase and

*ρ*, the synchronization index, is related to the standard deviation of the population. For*ρ*= 1, the population is clustered about one mean phase, whereas for*ρ*= 0 the population is evenly balanced across the unit circle. Complex variables for the population before and after perturbation can be calculated via$z:={\int}_{0}^{2\pi}{e}^{i\theta}p\left(\theta ,\tilde{t}\right)\phantom{\rule{0.25em}{0ex}}d\theta ,$

(34)

$\stackrel{\u02c6}{z}:={\int}_{0}^{2\pi}{e}^{i\theta}\stackrel{\u02c6}{p}\left(\theta ,\tilde{t}\right)\phantom{\rule{0.25em}{0ex}}d\theta ={\int}_{0}^{2\pi}{e}^{ig\left(\theta \right)}p\left(\theta ,\tilde{t}\right)\phantom{\rule{0.25em}{0ex}}d\theta .$

(35)

In Eq. 35, we avoid calculating the perturbed population explicitly by instead integrating over the new phases at the prior population density function. Population-level amplitude and phase responses can therefore be calculated by

$\Delta \overline{\theta}=\angle z-\angle \stackrel{\u02c6}{z},$

(36)

$\Delta \rho =\left|z\right|-\left|\stackrel{\u02c6}{z}\right|.$

(37)

Population-level PRCs and ARCs can be tabulated by solving Eqs. 34–37 for populations $p\left(\theta ,\tilde{t}\right)$ with different mean phases. It is important to note that the ARC at the population level strongly depends on the slope of the PRC, as has been shown previously (

14

).### Population-level mean expression profiles

To efficiently capture the population-level effects of bioluminescence experiments, we couple the detailed single-cell ODE model to a phase-density model. Previous work has used limit cycle models to estimate population-level parameters, such as desynchronization rate, for phase-only models (

36

). In this study, we use an ODE model to calculate the response to a perturbation at each phase, and subsequently take the weighted average of these responses according to the phases of the cells in the population.Assuming each oscillator in the population follows the dynamics described by

*x*^{γ}(*θ*), the unperturbed mean population-level expression, $\overline{x}\left(\tilde{t}\right)$, can be found by taking the weighted average of the expression level over our current population:$\overline{x}\left(\tilde{t}\right)={\int}_{0}^{2\pi}{x}^{\gamma}\left(\theta \right)p\left(\theta ,\tilde{t}\right)\phantom{\rule{0.25em}{0ex}}d\theta .$

(38)

Phase-diffusion models can explain why the gradual damping from experimental population-level data closely resembles an exponentially damped sinusoid, a result that has been shown experimentally (

Due to the smoothing effect of phase diffusion, higher-frequency sinusoidal components of the limit cycle are damped faster than lower-frequency components, resulting in exponentially decaying sinusoids even for limit cycles that are not sinusoidal in sufficiently disperse populations. The analytical result in Eq. 42 allows us to easily estimate the phase diffusivity parameter in Eq. 13 by simply fitting an exponentially damped sinusoid to detrended bioluminescence data.

9

) and computationally (36

). Our method similarly demonstrates exponential decay: for the idealized system *x*^{γ}(*θ*) = cos(*θ*) starting from a synchronized state,$\overline{x}\left(\tilde{t}\right)=\frac{1}{\sqrt{4\pi d\tilde{t}}}{\int}_{-\infty}^{\infty}\text{cos}\left(x\right)\text{exp}\left(-\frac{{\left(x-\tilde{t}\right)}^{2}}{4d\tilde{t}}\right)\phantom{\rule{0.25em}{0ex}}dx$

(39)

$=\mathrm{\Re}\left[\frac{1}{\sqrt{4\pi d\tilde{t}}}{\int}_{-\infty}^{\infty}\text{exp}\left(ix\right)\text{exp}\left(-\frac{{\left(x-\tilde{t}\right)}^{2}}{4d\tilde{t}}\right)\phantom{\rule{0.25em}{0ex}}dx\right]$

(40)

$=\mathrm{\Re}\left[{e}^{\left(i-d\right)\tilde{t}}\right]$

(41)

$={e}^{-d\tilde{t}}\text{cos}\left(\tilde{t}\right).$

(42)

Due to the smoothing effect of phase diffusion, higher-frequency sinusoidal components of the limit cycle are damped faster than lower-frequency components, resulting in exponentially decaying sinusoids even for limit cycles that are not sinusoidal in sufficiently disperse populations. The analytical result in Eq. 42 allows us to easily estimate the phase diffusivity parameter in Eq. 13 by simply fitting an exponentially damped sinusoid to detrended bioluminescence data.

Using our continuous approximation to population-level dynamics, we demonstrate an example of an unperturbed trajectory using Model 1 (adapted from Novák and Tyson (

41

); see the Supporting Material). In Fig. 3, an initially jagged population density smoothes and widens over time as cells desynchronize, similar to the effects of diffusion. Although each cell’s expression level follows the limit cycle, the population amplitude gradually damps with time as cells with diverse phases are averaged together.Next we describe how to calculate population-level mean expression after a perturbation. The perturbed trajectory $\stackrel{\u02c6}{x}\left(\tilde{t}\right)$ can be decomposed into contributions from two sources:

First, for long times after the perturbation, each individual oscillator will have returned to the limit cycle ${x}^{\gamma}\left(\theta \right)$, but with a new population density $\stackrel{\u02c6}{p}\left(\theta ,\tilde{t}\right)$. This steady-state perturbed trajectory ${\stackrel{\u02c6}{x}}_{ss}\left(\tilde{t}\right)$ can be found by

For very long times after perturbation, the new phase probability density could be approximated using the initial mean and standard deviation found through Eq. 35, because jagged profiles in the population density will eventually smooth to a normal distribution. However, for shorter times after perturbation, $\stackrel{\u02c6}{p}\left(\theta ,\tilde{t}\right)$ must be calculated numerically.

${\stackrel{\u02c6}{x}}_{ss}\left(\tilde{t}\right)={\int}_{0}^{2\pi}{x}^{\gamma}\left(\theta \right)\stackrel{\u02c6}{p}\left(\theta ,\tilde{t}\right)\phantom{\rule{0.25em}{0ex}}d\theta .$

(43)

For very long times after perturbation, the new phase probability density could be approximated using the initial mean and standard deviation found through Eq. 35, because jagged profiles in the population density will eventually smooth to a normal distribution. However, for shorter times after perturbation, $\stackrel{\u02c6}{p}\left(\theta ,\tilde{t}\right)$ must be calculated numerically.

The second contribution to the perturbed population trajectory comes from deviations in the limit cycle oscillations in each cell. We calculate the population-level effect of these deviations by averaging over the deviations that occur at each phase. We define the deviation trajectory $\delta x\left(\theta ,\tilde{t}\right)$ for each phase as the distance between the perturbed trajectory and the phase-adjusted reference:

Because the perturbed trajectory ultimately converges with the phase-adjusted reference, deviations will converge to zero. Because the phase change, Δ

Here the first term is equivalent to ${\stackrel{\u02c6}{x}}_{ss}\left(\tilde{t}\right)$, the steady-state perturbed trajectory, whereas the second term decays to zero as individual oscillators return to their steady-state amplitude. The accuracy of the continuous phase-diffusion model was tested by explicitly simulating a perturbation to a population of 225 uncoupled oscillators using a stochastic simulation algorithm (

$\delta x\left({\theta}_{0},\tilde{t}\right):=x\left(\tilde{t}\right)-{x}^{\gamma}\left(\tilde{t}+{\theta}_{0}+\Delta \theta \right),$

(44)

$\therefore \phantom{\rule{0.25em}{0ex}}\underset{\tilde{t}\to \infty}{\text{lim}}\delta x\left(\theta ,\tilde{t}\right)=0.$

(45)

Because the perturbed trajectory ultimately converges with the phase-adjusted reference, deviations will converge to zero. Because the phase change, Δ

*θ*, associated with a perturbation at each phase is likely not known before calculating the perturbed trajectory, it is difficult to tabulate deviation trajectories associated with each final phase,*θ*_{0}+ Δ*θ*. It is therefore more straightforward to find the average effect of single-cell perturbations at the population level by weighting the deviations by the phase density function before perturbation. The population-level response to a perturbation, $\stackrel{\u02c6}{x}\left(\tilde{t}\right)$, is therefore defined as$\stackrel{\u02c6}{x}\left(\tilde{t}\right)={\int}_{0}^{2\pi}{x}^{\gamma}\left(\theta \right)\stackrel{\u02c6}{p}\left(\theta ,\tilde{t}\right)+\delta x\left(\theta ,\tilde{t}\right)p\left(\theta ,\tilde{t}\right)\phantom{\rule{0.25em}{0ex}}d\theta .$

(46)

Here the first term is equivalent to ${\stackrel{\u02c6}{x}}_{ss}\left(\tilde{t}\right)$, the steady-state perturbed trajectory, whereas the second term decays to zero as individual oscillators return to their steady-state amplitude. The accuracy of the continuous phase-diffusion model was tested by explicitly simulating a perturbation to a population of 225 uncoupled oscillators using a stochastic simulation algorithm (

40

, 42

). Good agreement between the stochastic population and continuous approximation was shown in phase shift, single-cell level amplitude change, and alteration of population synchrony (see Fig. S2, Movie S1, and Movie S2 in the Supporting Material). These results indicate that the use of a continuous probability function is justified in the case of cultured cellular reporter systems, which may contain up to 10^{5}individual cells per culture (9

). In addition, the method was verified by similar simulations using the OREGONATOR (43

), a model containing only mass-action terms, and a detailed mechanistic model of circadian rhythms (44

) at a variety of phase timings (see Fig. S3). The continuous approximation is therefore able to capture the population-level dynamics of a wide variety of limit cycle models and perturbations.We next demonstrate the usefulness of single-cell and population level ARCs in predicting the mean response of a population. Using a detailed single-cell model (Model 2), we calculate the response of the system, specifically the average CRY2 protein expression level, to a temporary increase in

*Per*transcription rate. In Fig. 4*A*, we plot the resulting single-cell and population-level response curves. In this case, the population-level phase change is a slightly smoothed version of the single-cell PRC, because the population has an averaging effect on incoming perturbations with each cell receiving the input at a slightly different internal phase. This smoothing follows from the consistent definition of phase between the single-cell and population level. In the ARC plot, however, the shape of the single-cell and population-based ARC are different, because they describe different types of changes (finite versus sustained) in the output trajectories.Using these response curves, we demonstrate characteristic features of amplitude change mediated at the single-cell and population levels. In Fig. 4

*B*, the perturbation at*θ*_{1}does not cause a significant change in phase or population synchrony, yet strongly reduces single-cell amplitudes. The resulting population-level trajectory appears damped for the first cycle, until the oscillators within the population return to the limit cycle and normal amplitudes are restored. In Fig. 4*C*, the perturbation at*θ*_{2}reduces the population synchrony while increasing single-cell amplitudes, yielding amplitudes that are transiently increased before the population settles to a lower-amplitude, desynchronized trajectory.These examples demonstrate the qualitative differences between amplitude change mediated at the single-cell and population levels. Their characteristic features allow us to distinguish between these two sources of amplitude change without explicitly recording single-cell amplitudes: amplitude change induced at the population level will be sustained, but may be masked in the short term by changes to single-cell level amplitudes. Similarly, a change induced at the single-cell level will be evident only at short timescales. These results underscore the importance in considering both single-cell and population-mediated amplitude change in predicting the effects of daily stimuli on clock amplitudes.

### Application to experimental data

We conclude by using the proposed modeling framework to reconcile two experimental studies that measured the effect of light on the circadian amplitude of photosensitive fibroblast cells (

13

, 14

). Both studies sought to determine the main factor of amplitude reduction after a transient light pulse, in which either desynchrony alone (14

) or a combination of single-cell amplitude reduction and desynchrony (13

) was identified as the dominant factor.Data on the response of melanopsin-responsive and control NIH3T3 mouse fibroblasts cells is shown in Fig. 5

*B*(14

). The experimental data is denoised using a discrete wavelet transform (45

). In the experiment, a light pulse is given to desynchronize a colony of cells, resulting in drastically reduced amplitudes. A second (longer) light pulse subsequently resynchronizes the population, resulting in increased amplitudes. To capture this phenomenon in silico, we used a recent model of the core circadian feedback circuit (44

). To find parameters for the phase density function, we assumed an initial phase population with *σ*= 0.1 and calculated a phase diffusivity parameter*d*= 0.104 from the experimental control trajectory.The next step in capturing the experimental data was finding a perturbation capable of desynchronizing the system. To find such a perturbation, we calculated population-level ARCs at several different knockdown strengths for each parameter in the model. We selected the degradation rate of

*Per*mRNA from the feasible parameters due to PER’s known induction by CREB after photoperturbation (46

). Although not an exact match of the experimental system, reduction in the degradation rate yields a similar response to increasing transcription. The parameter was reduced by 28.5% during the light pulses, a value fit to maximize light-induced desynchrony. The output variable, PER2-luciferase in the experimental system, was chosen to be represented by *Cry1*mRNA—an E-box activated gene that is buffered from the direct effects of the parameter perturbation. The single-cell and population-level response curves for this parameter choice are shown in Fig. 5*A*, which demonstrates the strong desynchronization that occurs at*θ*≈*π*/4. Because the reporter used in the experimental system likely has a phase lag from the corresponding mRNA, we did not try to match the initial phase of the simulation to experiment. Instead, the initial phase of simulation was chosen such that the first pulse occurred when the system was at the phase that corresponded with the minimum of the population-level ARC.Population-level trajectories were found by solving for the resulting limit cycle trajectories at every phase, and subsequently finding the weighted average of these trajectories according the phase density (as described in Eqs. 38 and 46). The simulated control and perturbed trajectories, shown in Fig. 5

*B*, closely match the experimental results, in which the model captures both the phase shifts and amplitude modulation of the light pulses. The phase probability density function for the light-sensitive model trajectory is shown in Fig. 5*C*for several representative time points. Changes to the synchronization of the population from each perturbation are readily apparent, with both perturbations inducing a bimodal distribution in the phase density. Although experimental data on individual cell phases was unavailable for this data set, bimodal distributions in SCN neuron firing after a phase shift have previously been seen experimentally (47

). As higher-frequency features, these bimodal distributions dissipate as the phases diffuse. Amplitude change in this case is mediated at both the population- and single-cell levels, as indicated by their representative ARCs. The contributions from each source are summarized in Fig. S4.Several quantitative differences between the model and experiment do appear. Most importantly, the acute amplitude change after the light pulse is not correctly captured by the single-cell model. From the experimental system, it appears the light pulse temporarily reduces PER2-luciferase amplitudes immediately after perturbation, separate from the overall change in population synchrony. This result is most obvious after the second light pulse, where rhythms require some time before their maximum amplitude is reached. Such a delay is indicative of perturbations to the limit cycle system, suggesting a contribution from single cells in determining the population amplitude. Light-induced amplitude suppression in fibroblasts is therefore likely mediated at the single-cell level at short timescales, and at the population level for longer timescales. To improve the fit of the model to the data, the response curves could be fit to match experimentally predicted values. Infinitesimal PRCs and ARCs calculated at the single-cell level are numerically efficient to evaluate, and could be more readily incorporated into a parameter estimation algorithm than explicit stochastic simulations of a population of oscillators.

## Conclusion

In this study, we have described what we believe to be new tools for understanding amplitude change in independent cells based on analytically tractable methods for simulating a large population of oscillators. These tools synthesize single-cell level ODE sensitivity analysis metrics with population-level measures of synchrony, allowing the quantification of clock amplitude at multiple levels of biological organization. Using these tools, we have demonstrated how ODE models can be directly compared to experimental bioluminescence profiles, which, as of this writing, serve as a standard system to investigate clock perturbations. More generally, we have demonstrated the differences between acute and prolonged changes in amplitude after temporary perturbations, and characterized the mechanisms that give rise to each.

Although our examples have focused on light-mediated perturbations, these approaches could also be applied to pharmacological perturbations. Some difficulty arises in obtaining such data in cultured cells, however, as entraining perturbations would require that pharmacological agents be introduced for only a finite duration. Because medium or temperature changes are often enough to resynchronize cultured cells, a transient application is difficult to achieve experimentally. The search for clock-enhancing molecules has therefore tended to focus on constant drug concentrations: for instance, dose-dependent period or amplitude change after inhibition of CK1

*δ*or similar targets (48

). For in vivo systems, however, pharmacokinetics dictates a finite duration of action for both naturally secreted hormones and pharmacological therapies. More effective treatments might therefore be designed by explicitly accounting for such a transient response, perturbing peripheral clocks at the proper phase for inducing resynchronization and an increase in single-cell level amplitudes.Maintaining robust rhythmicity in peripheral oscillators has been linked to protection against metabolic disease. It is therefore interesting that liver cells have not developed a stronger mechanism of intercellular coupling, such as in the SCN, to maintain robust amplitudes in the absence of external cues. Perhaps the ability to quickly reentrain to a phase shift, which is often slowed by coupling (

4

), has historically been more advantageous than protection against highly variable food intake. Regardless, the necessities of modern society often dictate irregular circadian-metabolic regimes that damp peripheral clock oscillations. For such instances, we hope that the mathematical theory presented here will permit the development of optimized treatment regimes.We thank Yongqiang Wang for the helpful discussions, and gratefully acknowledge the three anonymous reviewers whose comments/suggestions helped improve this work.

This work was supported by the National Institutes of Health/National Institute of General Medical Sciences under award No. 1R01GM096873-01 and by the Institute for Collaborative Biotechnologies through grant No. W911NF-09-0001 from the U.S. Army Research Office.

## Supporting Material

- Document S1. Three models and four figures

- Movie S1. Stochastic Simulation of a Two-State Limit Cycle Oscillator

- Movie S2. Perturbation to a Population of Stochastic Oscillators

## References

- Physiological significance of a peripheral tissue circadian clock.
*Proc. Natl. Acad. Sci. USA.*2008; 105: 15172-15177 - System-driven and oscillator-dependent circadian transcription in mice with a conditionally active liver clock.
*PLoS Biol.*2007; 5: e34 - The amplitude of circadian oscillations: temperature dependence, latitudinal clines, and the photoperiodic time measurement.
*J. Biol. Rhythms.*1991; 6: 299-313 - Coupling governs entrainment range of circadian clocks.
*Mol. Syst. Biol.*2010; 6: 438 - Disruption of the clock components CLOCK and BMAL1 leads to hypoinsulinemia and diabetes.
*Nature.*2010; 466: 627-631 - Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet.
*Cell Metab.*2012; 15: 848-860 - SIRT1 mediates central circadian control in the SCN by a mechanism that decays with aging.
*Cell.*2013; 153: 1448-1460 - Spatiotemporal separation of PER and CRY posttranslational regulation in the mammalian circadian clock.
*Proc. Natl. Acad. Sci. USA.*2014; 111: 2040-2045 - Bioluminescence imaging of individual fibroblasts reveals persistent, independently phased circadian rhythms of clock gene expression.
*Curr. Biol.*2004; 14: 2289-2295 - High-throughput chemical screen identifies a novel potent modulator of cellular circadian rhythms and reveals CKI
*α*as a clock regulatory kinase.*PLoS Biol.*2010; 8: e1000559 - Cell type-specific functions of period genes revealed by novel adipocyte and hepatocyte circadian clock models.
*PLoS Genet.*2014; 10: e1004244 - Network features of the mammalian circadian clock.
*PLoS Biol.*2009; 7: e52 - Reciprocity between phase shifts and amplitude changes in the mammalian circadian clock.
*Proc. Natl. Acad. Sci. USA.*2007; 104: 20356-20361 - Melanopsin-dependent photo-perturbation reveals desynchronization underlying the singularity of mammalian circadian clocks.
*Nat. Cell Biol.*2007; 9: 1327-1334 - Toward a detailed computational model for the mammalian circadian clock.
*Proc. Natl. Acad. Sci. USA.*2003; 100: 7051-7056 - Modeling feedback loops of the mammalian circadian oscillator.
*Biophys. J.*2004; 87: 3023-3034 - Design principles underlying circadian clocks.
*J. R. Soc. Interface.*2004; 1: 119-130 - A Functional analysis of circadian pacemakers in nocturnal rodents.
*J. Comp. Physiol. A Neuroethol. Sens. Neural Behav. Physiol.*1976; 106: 253-266 - Oscillator model reduction preserving the phase response: application to the circadian clock.
*Biophys. J.*2008; 95: 1658-1673 - Robust entrainment of circadian oscillators requires specific phase response curves.
*Biophys. J.*2011; 100: 2557-2565 - Molecular mechanisms that regulate the coupled period of the mammalian circadian clock.
*Biophys. J.*2014; 106: 2071-2081 - Sensitivity analysis of oscillatory systems.
*Appl. Math. Model.*1984; 8: 328-340 - Isochron-based phase response analysis of circadian rhythms.
*Biophys. J.*2006; 91: 2131-2141 - Sensitivity measures for oscillating systems: application to mammalian circadian gene network.
*IEEE Trans. Automat. Contr.*2008; 53: 177-188 - Refinement of a limit cycle oscillator model of the effects of light on the human circadian pacemaker.
*J. Theor. Biol.*1998; 192: 455-465 - Modeling circadian rhythm generation in the suprachiasmatic nucleus with locally coupled self-sustained oscillators: phase shifts and phase response curves.
*J. Biol. Rhythms.*1999; 14: 460-468 - Cardiac atrial circadian rhythms in PERIOD2:LUCIFERASE and per1:luc mice: amplitude and phase responses to glucocorticoid signaling and medium treatment.
*PLoS ONE.*2012; 7: e47692 - Phase-amplitude response functions for transient-state stimuli.
*J. Math. Neurosci.*2013; 3: 13 - Amplitude of the SCN clock enhanced by the behavioral activity rhythm.
*PLoS ONE.*2012; 7: e39693 - A neuropeptide speeds circadian entrainment by reducing intercellular synchrony.
*Proc. Natl. Acad. Sci. USA.*2013; 110: E4355-E4361 - Sensitivity analysis in chemical kinetics.
*Annu. Rev. Phys. Chem.*1983; 34: 419-461 - Collective synchronization in populations of globally coupled phase oscillators with drifting frequencies.
*Phys. Rev. E Stat. Nonlin. Soft Matter Phys.*2006; 73: 011104 - Chemical Oscillations, Waves, and Turbulence.Springer, Berlin, Germany1984
- Robustness of the noise-induced phase synchronization in a general class of limit cycle oscillators.
*Phys. Rev. Lett.*2004; 93: 204103 - A theoretical analysis of neuronal variability.
*Biophys. J.*1965; 5: 173-194 - Dynamical signatures of cellular fluctuations and oscillator stability in peripheral circadian clocks.
*Mol. Syst. Biol.*2007; 3: 93 - Stochastic Models, Information Theory, and Lie Groups. Vol. 1. Birkhäuser, Boston, MA2009
- Directional Statistics.John Wiley, New York2009
Andersson, J. 2013. A General-Purpose Software Framework for Dynamic Optimization. Ph.D. thesis. Faculty of Engineering, KU Leuven, Leuven, Belgium.

- STOCHKIT2: software for discrete stochastic simulation of biochemical systems with events.
*Bioinformatics.*2011; 27: 2457-2458 - Design principles of biochemical oscillators.
*Nat. Rev. Mol. Cell Biol.*2008; 9: 981-991 - Exact stochastic simulation of coupled chemical reactions.
*J. Phys. Chem.*1977; 84: 2340-2361 - Oscillations in chemical systems. IV. Limit cycle behavior in a model of a real chemical reaction.
*J. Chem. Phys.*1974; 60: 1877 - Identification of small molecule activators of cryptochrome.
*Science.*2012; 337: 1094-1097 - Wavelet-based time series analysis of circadian rhythms.
*J. Biol. Rhythms.*2011; 26: 454-463 - Ca
^{2+}/cAMP response element-binding protein (CREB)-dependent activation of Per1 is required for light-induced signaling in the suprachiasmatic nucleus circadian clock.*J. Biol. Chem.*2003; 278: 718-723 - Phase resetting of the mammalian circadian clock relies on a rapid shift of a small population of pacemaker neurons.
*PLoS ONE.*2011; 6: e25437 - Small molecule modifiers of circadian clocks.
*Cell. Mol. Life Sci.*2013; 70: 2985-2998

## Article Info

### Publication History

Editor: Leah Edelstein-Keshet.

Accepted:
October 1,
2014

Received:
June 30,
2014

### Identification

### Copyright

© 2014 Biophysical Society. Published by Elsevier Inc.

### User License

Elsevier user license | How you can reuse

Elsevier's open access license policy

Elsevier user license

## Permitted

### For non-commercial purposes:

- Read, print & download
- Text & data mine
- Translate the article

## Not Permitted

- Reuse portions or extracts from the article in other works
- Redistribute or republish the final article
- Sell or re-use for commercial purposes

Elsevier's open access license policy