Amplitude Metrics for Cellular Circadian Bioluminescence Reporters

      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 (
      • Lamia K.A.
      • Storch K.-F.
      • Weitz C.J.
      Physiological significance of a peripheral tissue circadian clock.
      ). These peripheral clocks process systemic and SCN-mediated entraining cues, buffering against rhythmic changes in energy availability to maintain metabolic homeostasis (
      • Kornmann B.
      • Schaad O.
      • Schibler U.
      • et al.
      System-driven and oscillator-dependent circadian transcription in mice with a conditionally active liver clock.
      ). The amplitude of circadian transcription is a relevant factor, and has been shown to play a critical role in phase resetting and entrainment (
      • Pittendrigh C.S.
      • Kyner W.T.
      • Takamura T.
      The amplitude of circadian oscillations: temperature dependence, latitudinal clines, and the photoperiodic time measurement.
      ,
      • Abraham U.
      • Granada A.E.
      • Herzel H.
      • et al.
      Coupling governs entrainment range of circadian clocks.
      ). 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 (
      • Marcheva B.
      • Ramsey K.M.
      • Bass J.
      • et al.
      Disruption of the clock components CLOCK and BMAL1 leads to hypoinsulinemia and diabetes.
      ), while low-amplitude clock oscillations, whether caused by diet (
      • Hatori M.
      • Vollmers C.
      • Panda S.
      • et al.
      Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet.
      ) or age (
      • Chang H.-C.
      • Guarente L.
      SIRT1 mediates central circadian control in the SCN by a mechanism that decays with aging.
      ), 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 (
      • St John P.C.
      • Hirota T.
      • Doyle 3rd, F.J.
      • et al.
      Spatiotemporal separation of PER and CRY posttranslational regulation in the mammalian circadian clock.
      ). 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 (
      • Welsh D.K.
      • Yoo S.-H.
      • Kay S.A.
      • et al.
      Bioluminescence imaging of individual fibroblasts reveals persistent, independently phased circadian rhythms of clock gene expression.
      ). 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 (
      • Hirota T.
      • Lee J.W.
      • Kay S.A.
      • et al.
      High-throughput chemical screen identifies a novel potent modulator of cellular circadian rhythms and reveals CKIα as a clock regulatory kinase.
      ,
      • Ramanathan C.
      • Xu H.
      • Liu A.C.
      • et al.
      Cell type-specific functions of period genes revealed by novel adipocyte and hepatocyte circadian clock models.
      ). 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 (
      • Baggs J.E.
      • Price T.S.
      • Hogenesch J.B.
      • et al.
      Network features of the mammalian circadian clock.
      ). 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 (
      • Pulivarthy S.R.
      • Tanaka N.
      • Panda S.
      • et al.
      Reciprocity between phase shifts and amplitude changes in the mammalian circadian clock.
      ,
      • Ukai H.
      • Kobayashi T.J.
      • Ueda H.R.
      • et al.
      Melanopsin-dependent photo-perturbation reveals desynchronization underlying the singularity of mammalian circadian clocks.
      ).
      Mathematical models have long been used to understand the results of circadian experiments (
      • Leloup J.-C.
      • Goldbeter A.
      Toward a detailed computational model for the mammalian circadian clock.
      ,
      • Becker-Weimann S.
      • Wolf J.
      • Kramer A.
      • et al.
      Modeling feedback loops of the mammalian circadian oscillator.
      ), 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 (
      • Rand D.A.
      • Shulgin B.V.
      • Millar A.J.
      • et al.
      Design principles underlying circadian clocks.
      ). For instance, the phase response curve (PRC) has been used to characterize the entrainment behavior of both experimental and mathematical systems (
      • Daan S.
      • Pittendrigh C.S.
      A Functional analysis of circadian pacemakers in nocturnal rodents.
      ,
      • Taylor S.R.
      • Doyle 3rd, F.J.
      • Petzold L.R.
      Oscillator model reduction preserving the phase response: application to the circadian clock.
      ,
      • Pfeuty B.
      • Thommen Q.
      • Lefranc M.
      Robust entrainment of circadian oscillators requires specific phase response curves.
      ), and in analyzing the synchrony of populations of oscillators (
      • Kim J.K.
      • Kilpatrick Z.P.
      • Josić K.
      • et al.
      Molecular mechanisms that regulate the coupled period of the mammalian circadian clock.
      ). Accurate and efficient numerical routines for finding infinitesimal PRCs have therefore been developed (
      • Kramer M.A.
      • Rabitz H.
      • Calo J.M.
      Sensitivity analysis of oscillatory systems.
      ,
      • Gunawan R.
      • Doyle 3rd, F.J.
      Isochron-based phase response analysis of circadian rhythms.
      ,
      • Taylor S.R.
      • Gunawan R.
      • Doyle III, F.J.
      • et al.
      Sensitivity measures for oscillating systems: application to mammalian circadian gene network.
      ). 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 (
      • Jewett M.E.
      • Kronauer R.E.
      Refinement of a limit cycle oscillator model of the effects of light on the human circadian pacemaker.
      ). ARCs have also been used to characterize perturbations to groups of oscillators through desynchrony (
      • Ukai H.
      • Kobayashi T.J.
      • Ueda H.R.
      • et al.
      Melanopsin-dependent photo-perturbation reveals desynchronization underlying the singularity of mammalian circadian clocks.
      ,
      • Achermann P.
      • Kunz H.
      Modeling circadian rhythm generation in the suprachiasmatic nucleus with locally coupled self-sustained oscillators: phase shifts and phase response curves.
      ), at the level of the single oscillator (
      • van der Veen D.R.
      • Shao J.
      • Duffield G.E.
      • et al.
      Cardiac atrial circadian rhythms in PERIOD2:LUCIFERASE and per1:luc mice: amplitude and phase responses to glucocorticoid signaling and medium treatment.
      ,
      • Castejón O.
      • Guillamon A.
      • Huguet G.
      Phase-amplitude response functions for transient-state stimuli.
      ), and even in studying the effect of entrainment phase on SCN rhythm amplitude (
      • van Oosterhout F.
      • Lucassen E.A.
      • Meijer J.H.
      • et al.
      Amplitude of the SCN clock enhanced by the behavioral activity rhythm.
      ). 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 (
      • An S.
      • Harang R.
      • Herzog E.D.
      • et al.
      A neuropeptide speeds circadian entrainment by reducing intercellular synchrony.
      ). 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
      dxdt=f(x(t),p),
      (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:
      limt[x(t+T)x(t)]=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π:
      t˜=2πTt;f˜=T2πf;dxdt˜=f˜(x(t˜),p).
      (3)


      The phase variable θ is therefore defined on the limit cycle as θ=t˜mod2π, 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 x(t). Here the initial conditions are determined by the strength of the perturbation and the phase at which it is applied:
      x(0):=xγ(θ0)+Δx(0).
      (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:
      Δx(t˜)=x(t˜)xγ(t˜+θ0).
      (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 d˜ (in radians) that ends at θ0, the oscillator deviates from the limit cycle trajectory according to the following:
      xd˜(d˜)=xγ(θ0d˜),
      (6)


      dxd˜dt˜=f˜(xd˜(t˜),p+Δp).
      (7)


      The pulse trajectory xd˜(t˜) is then integrated from t˜=d˜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
      Δx(0)=xd˜(0)xγ(θ0).
      (8)


      A schematic depicting a perturbed and reference trajectory exemplar is shown in Fig. 1 A.
      Figure thumbnail gr1
      Figure 1Amplitude metrics at the single-cell level measure transient deviations from the limit cycle. (A) Schematic showing trajectories used in the calculation of single-cell phase and amplitude change. The perturbation x(t˜,p) (blue), from the limit cycle solution xγ(t˜,p) (red), ultimately returns to the limit cycle with a phase shift xγ(t˜+Δθ,p) (green). (B–D) The same perturbation applied at two different phases can result in opposite amplitude effects. The state-space representation (B) reveals the path both perturbations take to return to the limit cycle. The time-series representation (C) shows how one perturbation (blue) results in an amplitude decrease, whereas the other (red) results in an amplitude increase. (D) This amplitude change is quantified by integrating h(t˜), the difference in variance from each solution to the limit cycle as defined in Eq. 21, from t = 0 → ∞. (Shaded region) Area under the curve, ΔA, for the perturbations (red and blue). Model adapted from Novák and Tyson (
      • Novák B.
      • Tyson J.J.
      Design principles of biochemical oscillators.
      ). To see this figure in color, go online.
      The response of the system to very small perturbations is more efficiently calculated using ODE sensitivity analysis (
      • Rabitz H.
      • Kramer M.
      • Dacol D.
      Sensitivity analysis in chemical kinetics.
      ). The sensitivity matrix,
      S(t)=dx(t)dx(0)=limΔx(0)0Δx(t)Δx(0),
      (9)


      can be calculated directly by integrating
      ddtS(t)=df(x(t),p)dxS(t)
      (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 (
      • Kramer M.A.
      • Rabitz H.
      • Calo J.M.
      Sensitivity analysis of oscillatory systems.
      ,
      • Taylor S.R.
      • Gunawan R.
      • Doyle III, F.J.
      • et al.
      Sensitivity measures for oscillating systems: application to mammalian circadian gene network.
      ). 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 (
      • Taylor S.R.
      • Gunawan R.
      • Doyle III, F.J.
      • et al.
      Sensitivity measures for oscillating systems: application to mammalian circadian gene network.
      ). Methods for efficiently calculating these quantities using ODE sensitivity analysis have been developed (
      • Taylor S.R.
      • Gunawan R.
      • Doyle III, F.J.
      • et al.
      Sensitivity measures for oscillating systems: application to mammalian circadian gene network.
      ), with the important result that the parameter- and state-impulse PRCs can be related by the following Jacobian matrix:
      ddt˜dθdp=dθdxdf˜dp.
      (11)


      This result follows from the fact that in the limit of an infinitely short and small parameter pulse,
      Δx(0)df˜dpd˜Δp.


      Phase-diffusion model

      Large populations of oscillators are typically described using phase-only models (
      • Rougemont J.
      • Naef F.
      Collective synchronization in populations of globally coupled phase oscillators with drifting frequencies.
      ), 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(θ,t˜) that describes the probability of finding an oscillator at each phase (
      • Kuramoto Y.
      Chemical Oscillations, Waves, and Turbulence.
      ). The usefulness of probability density functions in describing the phase and amplitude responses of populations of circadian cells has been previously shown (
      • Ukai H.
      • Kobayashi T.J.
      • Ueda H.R.
      • et al.
      Melanopsin-dependent photo-perturbation reveals desynchronization underlying the singularity of mammalian circadian clocks.
      ). As with all probability density functions,
      02πp(θ,t˜)dθ=1.
      (12)


      The shape of p(θ,t˜) 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 (
      • Teramae J.N.
      • Tanaka D.
      Robustness of the noise-induced phase synchronization in a general class of limit cycle oscillators.
      ). For an infinite population of oscillators, these effects are well described by a Fokker-Plank equation (
      • Stein R.B.
      A theoretical analysis of neuronal variability.
      ):
      pt˜=pθ+d2pθ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 2p/∂θ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 (
      • Rougemont J.
      • Naef F.
      Dynamical signatures of cellular fluctuations and oscillator stability in peripheral circadian clocks.
      ). Equation 13 has periodic boundary conditions, with initial condition ϕ(θ) as the phase population at t = 0:
      BCs:p(0,t˜)=p(2π,t˜),
      (14)


      pθ(0,t˜)=pθ(2π,t˜),
      (15)


      IC:p(θ,0)=ϕ(θ).
      (16)


      The solution of Eqs. 13–16 is well characterized, with p(θ,t˜) evolving in time as the convolution of the initial conditions with a wrapped normal distribution with mean t˜ and standard deviation 2dt˜ (
      • Chirikjian G.S.
      ):
      p(θ,t˜)=ϕ(θ)WN(θ;t˜,2dt˜),
      (17)


      in which the wrapped normal distribution (
      • Mardia K.V.
      • Jupp P.E.
      Directional Statistics.
      ) is defined as
      WN(θ;μ,σ)=1σ2πk=exp[(θμ+2πk)22σ2].
      (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(θ,t˜) can be found analytically through
      p(θ,t˜)=WN(θ;μ0+t˜,σ02+2d2t˜2).
      (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) (

      Andersson, J. 2013. A General-Purpose Software Framework for Dynamic Optimization. Ph.D. thesis. Faculty of Engineering, KU Leuven, Leuven, Belgium.

      ). Stochastic simulations were performed using the STOCHKIT2 package (http://sourceforge.net/projects/stochkit/) (
      • Sanft K.R.
      • Wu S.
      • Petzold L.R.
      • et al.
      STOCHKIT2: software for discrete stochastic simulation of biochemical systems with events.
      ). 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.
      Therefore, to derive continuous approximations to the dynamics of a large population of noninteracting oscillators, we first demonstrate how ARCs can be found at both the single-cell and population level.

      Definition of an amplitude metric at the single-cell level

      After a perturbation to a single cell, the path that the perturbed trajectory 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
      μ:=02πxγ(θ)2πdθ.
      (20)


      The amplitude change metric is defined as
      ΔA(x(t),y(t)):=0(x(t)μ)2(y(t)μ)2dt=0h(t)dt.
      (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, BD, 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(t˜), and the corresponding integrand for the amplitude change equation, h(t˜), 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.
      Figure thumbnail gr2
      Figure 2Phase and amplitude response curves (PRCs and ARCs) describe how a single oscillator will respond to a perturbation. (A and B) PRCs and ARCS for three perturbations of increasing strengths demonstrate the transition from Type-1 phase resetting for a weak stimulus (blue) to type 0 phase resetting from a strong stimulus (red). (C) The state space representation for the perturbations of varying strengths demonstrate how the trajectories return to the limit cycle. For one particular perturbation (green), the trajectory starts very close to the singularity, and therefore takes a long time to recover. (D) For an oscillator perturbed to the singularity, a strong amplitude reduction ensues. Here, the oscillator with an initial phase θ ≈ 2π/2 is perturbed by Δx = 1.1 at t˜=0. Several cycles are required for normal amplitudes to be restored, corresponding to a dip in the ARC (C, green). To see this figure in color, go online.
      Infinitesimal versions of response curves are often more general and easier to compute than those that track a specific perturbation strength (
      • Rand D.A.
      • Shulgin B.V.
      • Millar A.J.
      • et al.
      Design principles underlying circadian clocks.
      ). We therefore derive an expression for the infinitesimal ARC, defined as
      dAdx:=limΔx(0)0ΔA(x(t˜),xγ(t˜+θ0+Δθ))Δx(0)=0limΔx(0)0h(t)Δx(0)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θ=t˜+θ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γ(tθ+Δθ)=xγ(tθ)+dxγ(tθ)dθΔθ+O(Δθ2)
      (23)


      =xγ(tθ)+f˜(xγ(tθ))Δθ+O(Δθ2).
      (24)


      Simplifying the integrand in Eq. 22,
      h(t)Δx(0)=1Δx(0)[(xγ(tθ)+Δx(t˜)μ)2(xγ(tθ+Δθ)μ)2]
      (25)


      =1Δx(0)[(xγ(tθ)+Δx(t˜)μ)2(xγ(tθ)+f˜(xγ(tθ))Δθμ)2]
      (26)


      =1Δx(0)[(Δx(t˜)f˜(xγ(tθ))Δθ)(Δx(t˜)+f˜(xγ(tθ))Δθ+2(xγ(tθ)μ))].
      (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:
      limΔx(0)0h(t)Δx(0)=2[(limΔx(0)0Δx(t˜)Δx(0))f˜(xγ(tθ))(limΔx(0)0ΔθΔx(0))](xγ(tθ)μ)
      (28)


      =2(S(t˜)f˜(xγ(tθ))dθdx)(xγ(tθ)μ).
      (29)


      Note that in Eq. 29, /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:
      dAdx=02(S(t˜)f˜(xγ(tθ))dθdx)(xγ(tθ)μ)dt.
      (30)


      Here, the first term of the integrand (Sf˜θ˙) 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
      ddt˜dAdp:=limd˜,Δp0ΔAd˜Δp,
      (31)


      and may be calculated from the state-impulse version with the following relationship:
      ddt˜dAdp=dAdxdf˜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 (
      • Ukai H.
      • Kobayashi T.J.
      • Ueda H.R.
      • et al.
      Melanopsin-dependent photo-perturbation reveals desynchronization underlying the singularity of mammalian circadian clocks.
      ,
      • Kuramoto Y.
      Chemical Oscillations, Waves, and Turbulence.
      ), 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:
      pˆ(θ,t˜)dg(θ)=p(θ,t˜)dθ.
      (33)


      However, it is easier to estimate the mean and standard deviation of the perturbed population pˆ(θ,t˜) using directional statistics (
      • Mardia K.V.
      • Jupp P.E.
      Directional Statistics.
      ). A population defined on the unit circle can be described by a complex variable z=ρeiθ¯, where θ¯ 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:=02πeiθp(θ,t˜)dθ,
      (34)


      zˆ:=02πeiθpˆ(θ,t˜)dθ=02πeig(θ)p(θ,t˜)dθ.
      (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
      Δθ¯=zzˆ,
      (36)


      Δρ=|z||zˆ|.
      (37)


      Population-level PRCs and ARCs can be tabulated by solving Eqs. 34–37 for populations p(θ,t˜) 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 (
      • Ukai H.
      • Kobayashi T.J.
      • Ueda H.R.
      • et al.
      Melanopsin-dependent photo-perturbation reveals desynchronization underlying the singularity of mammalian circadian clocks.
      ).

      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 (
      • Rougemont J.
      • Naef F.
      Dynamical signatures of cellular fluctuations and oscillator stability in peripheral circadian clocks.
      ). 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, x¯(t˜), can be found by taking the weighted average of the expression level over our current population:
      x¯(t˜)=02πxγ(θ)p(θ,t˜)dθ.
      (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 (
      • Welsh D.K.
      • Yoo S.-H.
      • Kay S.A.
      • et al.
      Bioluminescence imaging of individual fibroblasts reveals persistent, independently phased circadian rhythms of clock gene expression.
      ) and computationally (
      • Rougemont J.
      • Naef F.
      Dynamical signatures of cellular fluctuations and oscillator stability in peripheral circadian clocks.
      ). Our method similarly demonstrates exponential decay: for the idealized system xγ(θ) = cos(θ) starting from a synchronized state,
      x¯(t˜)=14πdt˜cos(x)exp((xt˜)24dt˜)dx
      (39)


      =[14πdt˜exp(ix)exp((xt˜)24dt˜)dx]
      (40)


      =[e(id)t˜]
      (41)


      =edt˜cos(t˜).
      (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 (
      • Novák B.
      • Tyson J.J.
      Design principles of biochemical oscillators.
      ); 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.
      Figure thumbnail gr3
      Figure 3Synchrony affects population-level amplitude. (A) Stochastic fluctuations cause a population of cells to gradually desynchronize with time. The phase probability density, p(θ,t˜), gradually widens as it advances in phase according to the mean period. (B) The mean amplitude from a population of oscillators is determined by the probability density function p(θ,t˜) and the oscillator’s limit cycle xγ(θ). As time passes, population-level rhythms resemble an exponentially damped sinusoid. To see this figure in color, go online.
      Next we describe how to calculate population-level mean expression after a perturbation. The perturbed trajectory xˆ(t˜) 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γ(θ), but with a new population density pˆ(θ,t˜). This steady-state perturbed trajectory xˆss(t˜) can be found by
      xˆss(t˜)=02πxγ(θ)pˆ(θ,t˜)dθ.
      (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, pˆ(θ,t˜) 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 δx(θ,t˜) for each phase as the distance between the perturbed trajectory and the phase-adjusted reference:
      δx(θ0,t˜):=x(t˜)xγ(t˜+θ0+Δθ),
      (44)


      limt˜δx(θ,t˜)=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, xˆ(t˜), is therefore defined as
      xˆ(t˜)=02πxγ(θ)pˆ(θ,t˜)+δx(θ,t˜)p(θ,t˜)dθ.
      (46)


      Here the first term is equivalent to xˆss(t˜), 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 (
      • Sanft K.R.
      • Wu S.
      • Petzold L.R.
      • et al.
      STOCHKIT2: software for discrete stochastic simulation of biochemical systems with events.
      ,
      • Gillespie D.T.
      Exact stochastic simulation of coupled chemical reactions.
      ). 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 105 individual cells per culture (
      • Welsh D.K.
      • Yoo S.-H.
      • Kay S.A.
      • et al.
      Bioluminescence imaging of individual fibroblasts reveals persistent, independently phased circadian rhythms of clock gene expression.
      ). In addition, the method was verified by similar simulations using the OREGONATOR (
      • Field R.J.
      • Noyes R.M.
      Oscillations in chemical systems. IV. Limit cycle behavior in a model of a real chemical reaction.
      ), a model containing only mass-action terms, and a detailed mechanistic model of circadian rhythms (
      • Hirota T.
      • Lee J.W.
      • Kay S.A.
      • et al.
      Identification of small molecule activators of cryptochrome.
      ) 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.
      Figure thumbnail gr4
      Figure 4Response curves describe different types of amplitude change. A population of cells is simulated using a detailed model of circadian rhythms (Model 2). Each cell in the population is subjected to a 20% increase in the vtp parameter (Per transcription rate) for d = π/2, inducing both a phase and amplitude change. (A) The PRCs are shown for both the single-limit cycle oscillator and population average (left). The population-level ARC is shown together with the single-cell amplitude response, with both curves normalized to σ = 1. (B and C) Population-level rhythms, xˆ(t˜) resulting from perturbations given at the indicated phases in panel A (red and cyan lines, respectively). The population transitions from the unperturbed population, x¯(t˜), to the steady-state perturbed population, xˆss(t˜). (B) A perturbation at this phase yields a transient amplitude reduction resulting from perturbations at the single-cell level. However, there is little change in population synchrony. Population-level rhythms are therefore transiently damped before regaining normal amplitudes. (C) Oscillators are desynchronized, but with a transient increase in limit cycle amplitude. It therefore takes some time before the amplitude reduction from desynchrony is realized. To see this figure in color, go online.
      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 (
      • Pulivarthy S.R.
      • Tanaka N.
      • Panda S.
      • et al.
      Reciprocity between phase shifts and amplitude changes in the mammalian circadian clock.
      ,
      • Ukai H.
      • Kobayashi T.J.
      • Ueda H.R.
      • et al.
      Melanopsin-dependent photo-perturbation reveals desynchronization underlying the singularity of mammalian circadian clocks.
      ). Both studies sought to determine the main factor of amplitude reduction after a transient light pulse, in which either desynchrony alone (
      • Ukai H.
      • Kobayashi T.J.
      • Ueda H.R.
      • et al.
      Melanopsin-dependent photo-perturbation reveals desynchronization underlying the singularity of mammalian circadian clocks.
      ) or a combination of single-cell amplitude reduction and desynchrony (
      • Pulivarthy S.R.
      • Tanaka N.
      • Panda S.
      • et al.
      Reciprocity between phase shifts and amplitude changes in the mammalian circadian clock.
      ) 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 (
      • Ukai H.
      • Kobayashi T.J.
      • Ueda H.R.
      • et al.
      Melanopsin-dependent photo-perturbation reveals desynchronization underlying the singularity of mammalian circadian clocks.
      ). The experimental data is denoised using a discrete wavelet transform (
      • Leise T.L.
      • Harrington M.E.
      Wavelet-based time series analysis of circadian rhythms.
      ). 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 (
      • Hirota T.
      • Lee J.W.
      • Kay S.A.
      • et al.
      Identification of small molecule activators of cryptochrome.
      ). 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.
      Figure thumbnail gr5
      Figure 5Method allows direct comparison between model results and experimental bioluminescence profiles. (A) Response curves for the parameter perturbation used to model the first light pulse in panel B. The initial phase of the system is chosen to coincide with the strong dip in the population-level ARC (left). (B) Denoised bioluminescence data from Ukai et al. (
      • Ukai H.
      • Kobayashi T.J.
      • Ueda H.R.
      • et al.
      Melanopsin-dependent photo-perturbation reveals desynchronization underlying the singularity of mammalian circadian clocks.
      ) (top) demonstrates population-level amplitude change resulting from two light pulses (
      • Ukai H.
      • Kobayashi T.J.
      • Ueda H.R.
      • et al.
      Melanopsin-dependent photo-perturbation reveals desynchronization underlying the singularity of mammalian circadian clocks.
      ). Circadian oscillations are suppressed by the first light pulse and rescued by the second. These perturbations were reproduced in silico (bottom) using the model from Hirota et al. (
      • Hirota T.
      • Lee J.W.
      • Kay S.A.
      • et al.
      Identification of small molecule activators of cryptochrome.
      ). The experimental results are qualitatively captured by the model, demonstrating the suitability of the method in capturing bioluminescence experiments. The relative amplitude of the model equations has been scaled to match that of the normalized bioluminescence profiles. (C) Example probability density functions for the model trajectory shown in panel A. The initial population (red) is effectively desynchronized (orange) after the first light pulse. After the second light pulse, stronger peaks are seen (teal), which gradually damp as time evolves. To see this figure in color, go online.
      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 (
      • Tischkau S.A.
      • Mitchell J.W.
      • Gillette M.U.
      • et al.
      Ca2+/cAMP response element-binding protein (CREB)-dependent activation of Per1 is required for light-induced signaling in the suprachiasmatic nucleus circadian clock.
      ). 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 (
      • Rohling J.H.T.
      • vanderLeest H.T.
      • Meijer J.H.
      • et al.
      Phase resetting of the mammalian circadian clock relies on a rapid shift of a small population of pacemaker neurons.
      ). 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 (
      • Chen Z.
      • Yoo S.-H.
      • Takahashi J.S.
      Small molecule modifiers of circadian clocks.
      ). 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 (
      • Abraham U.
      • Granada A.E.
      • Herzel H.
      • et al.
      Coupling governs entrainment range of circadian clocks.
      ), 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

      References

        • Lamia K.A.
        • Storch K.-F.
        • Weitz C.J.
        Physiological significance of a peripheral tissue circadian clock.
        Proc. Natl. Acad. Sci. USA. 2008; 105: 15172-15177
        • Kornmann B.
        • Schaad O.
        • Schibler U.
        • et al.
        System-driven and oscillator-dependent circadian transcription in mice with a conditionally active liver clock.
        PLoS Biol. 2007; 5: e34
        • Pittendrigh C.S.
        • Kyner W.T.
        • Takamura T.
        The amplitude of circadian oscillations: temperature dependence, latitudinal clines, and the photoperiodic time measurement.
        J. Biol. Rhythms. 1991; 6: 299-313
        • Abraham U.
        • Granada A.E.
        • Herzel H.
        • et al.
        Coupling governs entrainment range of circadian clocks.
        Mol. Syst. Biol. 2010; 6: 438
        • Marcheva B.
        • Ramsey K.M.
        • Bass J.
        • et al.
        Disruption of the clock components CLOCK and BMAL1 leads to hypoinsulinemia and diabetes.
        Nature. 2010; 466: 627-631
        • Hatori M.
        • Vollmers C.
        • Panda S.
        • et al.
        Time-restricted feeding without reducing caloric intake prevents metabolic diseases in mice fed a high-fat diet.
        Cell Metab. 2012; 15: 848-860
        • Chang H.-C.
        • Guarente L.
        SIRT1 mediates central circadian control in the SCN by a mechanism that decays with aging.
        Cell. 2013; 153: 1448-1460
        • St John P.C.
        • Hirota T.
        • Doyle 3rd, F.J.
        • et al.
        Spatiotemporal separation of PER and CRY posttranslational regulation in the mammalian circadian clock.
        Proc. Natl. Acad. Sci. USA. 2014; 111: 2040-2045
        • Welsh D.K.
        • Yoo S.-H.
        • Kay S.A.
        • et al.
        Bioluminescence imaging of individual fibroblasts reveals persistent, independently phased circadian rhythms of clock gene expression.
        Curr. Biol. 2004; 14: 2289-2295
        • Hirota T.
        • Lee J.W.
        • Kay S.A.
        • et al.
        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
        • Ramanathan C.
        • Xu H.
        • Liu A.C.
        • et al.
        Cell type-specific functions of period genes revealed by novel adipocyte and hepatocyte circadian clock models.
        PLoS Genet. 2014; 10: e1004244
        • Baggs J.E.
        • Price T.S.
        • Hogenesch J.B.
        • et al.
        Network features of the mammalian circadian clock.
        PLoS Biol. 2009; 7: e52
        • Pulivarthy S.R.
        • Tanaka N.
        • Panda S.
        • et al.
        Reciprocity between phase shifts and amplitude changes in the mammalian circadian clock.
        Proc. Natl. Acad. Sci. USA. 2007; 104: 20356-20361
        • Ukai H.
        • Kobayashi T.J.
        • Ueda H.R.
        • et al.
        Melanopsin-dependent photo-perturbation reveals desynchronization underlying the singularity of mammalian circadian clocks.
        Nat. Cell Biol. 2007; 9: 1327-1334
        • Leloup J.-C.
        • Goldbeter A.
        Toward a detailed computational model for the mammalian circadian clock.
        Proc. Natl. Acad. Sci. USA. 2003; 100: 7051-7056
        • Becker-Weimann S.
        • Wolf J.
        • Kramer A.
        • et al.
        Modeling feedback loops of the mammalian circadian oscillator.
        Biophys. J. 2004; 87: 3023-3034
        • Rand D.A.
        • Shulgin B.V.
        • Millar A.J.
        • et al.
        Design principles underlying circadian clocks.
        J. R. Soc. Interface. 2004; 1: 119-130
        • Daan S.
        • Pittendrigh C.S.
        A Functional analysis of circadian pacemakers in nocturnal rodents.
        J. Comp. Physiol. A Neuroethol. Sens. Neural Behav. Physiol. 1976; 106: 253-266
        • Taylor S.R.
        • Doyle 3rd, F.J.
        • Petzold L.R.
        Oscillator model reduction preserving the phase response: application to the circadian clock.
        Biophys. J. 2008; 95: 1658-1673
        • Pfeuty B.
        • Thommen Q.
        • Lefranc M.
        Robust entrainment of circadian oscillators requires specific phase response curves.
        Biophys. J. 2011; 100: 2557-2565
        • Kim J.K.
        • Kilpatrick Z.P.
        • Josić K.
        • et al.
        Molecular mechanisms that regulate the coupled period of the mammalian circadian clock.
        Biophys. J. 2014; 106: 2071-2081
        • Kramer M.A.
        • Rabitz H.
        • Calo J.M.
        Sensitivity analysis of oscillatory systems.
        Appl. Math. Model. 1984; 8: 328-340
        • Gunawan R.
        • Doyle 3rd, F.J.
        Isochron-based phase response analysis of circadian rhythms.
        Biophys. J. 2006; 91: 2131-2141
        • Taylor S.R.
        • Gunawan R.
        • Doyle III, F.J.
        • et al.
        Sensitivity measures for oscillating systems: application to mammalian circadian gene network.
        IEEE Trans. Automat. Contr. 2008; 53: 177-188
        • Jewett M.E.
        • Kronauer R.E.
        Refinement of a limit cycle oscillator model of the effects of light on the human circadian pacemaker.
        J. Theor. Biol. 1998; 192: 455-465
        • Achermann P.
        • Kunz H.
        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
        • van der Veen D.R.
        • Shao J.
        • Duffield G.E.
        • et al.
        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
        • Castejón O.
        • Guillamon A.
        • Huguet G.
        Phase-amplitude response functions for transient-state stimuli.
        J. Math. Neurosci. 2013; 3: 13
        • van Oosterhout F.
        • Lucassen E.A.
        • Meijer J.H.
        • et al.
        Amplitude of the SCN clock enhanced by the behavioral activity rhythm.
        PLoS ONE. 2012; 7: e39693
        • An S.
        • Harang R.
        • Herzog E.D.
        • et al.
        A neuropeptide speeds circadian entrainment by reducing intercellular synchrony.
        Proc. Natl. Acad. Sci. USA. 2013; 110: E4355-E4361
        • Rabitz H.
        • Kramer M.
        • Dacol D.
        Sensitivity analysis in chemical kinetics.
        Annu. Rev. Phys. Chem. 1983; 34: 419-461
        • Rougemont J.
        • Naef F.
        Collective synchronization in populations of globally coupled phase oscillators with drifting frequencies.
        Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 2006; 73: 011104
        • Kuramoto Y.
        Chemical Oscillations, Waves, and Turbulence.
        Springer, Berlin, Germany1984
        • Teramae J.N.
        • Tanaka D.
        Robustness of the noise-induced phase synchronization in a general class of limit cycle oscillators.
        Phys. Rev. Lett. 2004; 93: 204103
        • Stein R.B.
        A theoretical analysis of neuronal variability.
        Biophys. J. 1965; 5: 173-194
        • Rougemont J.
        • Naef F.
        Dynamical signatures of cellular fluctuations and oscillator stability in peripheral circadian clocks.
        Mol. Syst. Biol. 2007; 3: 93
        • Chirikjian G.S.
        Stochastic Models, Information Theory, and Lie Groups. Vol. 1. Birkhäuser, Boston, MA2009
        • Mardia K.V.
        • Jupp P.E.
        Directional Statistics.
        John Wiley, New York2009
      1. Andersson, J. 2013. A General-Purpose Software Framework for Dynamic Optimization. Ph.D. thesis. Faculty of Engineering, KU Leuven, Leuven, Belgium.

        • Sanft K.R.
        • Wu S.
        • Petzold L.R.
        • et al.
        STOCHKIT2: software for discrete stochastic simulation of biochemical systems with events.
        Bioinformatics. 2011; 27: 2457-2458
        • Novák B.
        • Tyson J.J.
        Design principles of biochemical oscillators.
        Nat. Rev. Mol. Cell Biol. 2008; 9: 981-991
        • Gillespie D.T.
        Exact stochastic simulation of coupled chemical reactions.
        J. Phys. Chem. 1977; 84: 2340-2361
        • Field R.J.
        • Noyes R.M.
        Oscillations in chemical systems. IV. Limit cycle behavior in a model of a real chemical reaction.
        J. Chem. Phys. 1974; 60: 1877
        • Hirota T.
        • Lee J.W.
        • Kay S.A.
        • et al.
        Identification of small molecule activators of cryptochrome.
        Science. 2012; 337: 1094-1097
        • Leise T.L.
        • Harrington M.E.
        Wavelet-based time series analysis of circadian rhythms.
        J. Biol. Rhythms. 2011; 26: 454-463
        • Tischkau S.A.
        • Mitchell J.W.
        • Gillette M.U.
        • et al.
        Ca2+/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
        • Rohling J.H.T.
        • vanderLeest H.T.
        • Meijer J.H.
        • et al.
        Phase resetting of the mammalian circadian clock relies on a rapid shift of a small population of pacemaker neurons.
        PLoS ONE. 2011; 6: e25437
        • Chen Z.
        • Yoo S.-H.
        • Takahashi J.S.
        Small molecule modifiers of circadian clocks.
        Cell. Mol. Life Sci. 2013; 70: 2985-2998
      Advertisement