## Abstract

The phase response curve (PRC) of the circadian clock provides one of the most significant indices for anticipating entrainment of outer cycles, despite the difficulty of making precise PRC determinations in experiments. We characterized the PRC of the *Arabidopsis* *thaliana* circadian clock on the basis of its phase-locking property to variable periodic pulse perturbations. Experiments revealed that the PRC changed remarkably from continuous to discontinuous fashion, depending on the oscillation amplitude. Our hypothesis of amplitude-dependent adaptability to outer cycles was successfully clarified by elucidation of this transition of PRC as a change in the collective response of the circadian oscillator network. These findings provide an essential criterion against which to evaluate the precision of PRC measurement and an advanced understanding of the adaptability of plant circadian systems to environmental conditions.

## INTRODUCTION

Almost all species on Earth have internal circadian clocks as an adaptation to diurnal environmental cycles (*1*–*3*). The circadian rhythm is entrained systematically to outer cycles through indigenous responses, which depend upon the phase in which the zeitgeber stimulus was received. These phase-dependent responses are described by the phase response curve (PRC), which plots phase shifts of the circadian system elicited by the zeitgeber stimuli as a function of the phases of the stimuli (*4*–*9*). The PRC is quite essential for physiological adaptations to the environment, for example, recovery from jet lag in mammals (*10*) and photosynthetic entrainment in plants (*11*). The mechanisms of the PRC have been well studied and applied to various fields of science and engineering, including physics, neuroscience, and chronobiology. However, precise characterization of the circadian system’s PRC requires tedious and time-consuming measurements. The obtained PRC is often extremely noisy or is sometimes oversimplified because of environmental fluctuations, individual variability, and other factors (*7*–*9*). Here, we exploit an undiscovered property of the circadian clock to precisely characterize the PRC.

In the standard method for determining the PRC, phase shifts are measured by applying pulse stimuli at various phases of the free-running period (*4*, *7*). However, the method requires many plant individuals as well as an extremely long measurement time because only a single perturbation can be applied to each individual. Thus, we exploited an alternative method in which input pulses are injected to each system in a periodic fashion. The PRC can be reconstructed by measuring the phase shifts induced by the periodic inputs (see also schematic illustration in fig. S1). Our method can be considered as a special case of the weighted spike-triggered average of fluctuating stimuli, which has been recently developed to measure PRCs in neuroscience (*12*). Compared to the standard method that applies only a single stimulus to each plant, the present approach applies many stimuli to each individual, thus requiring measurements on fewer individuals while generating abundant information on the phase response properties. Numerical study demonstrated that our method is capable of measuring PRCs without losing accuracy of the standard method (fig. S2; see Materials and Methods for the detailed verification of our proposed method).

## RESULTS

We used *Arabidopsis thaliana* as our model plant and characterized the PRC of the circadian system. We monitored the expression of the clock gene *CIRCADIAN CLOCK ASSOCIATED 1* (*CCA1*) using the transgenic *A. thaliana CCA1::LUC* strain, the firefly luciferase (*LUC*) gene fused to a promoter sequence containing a part upstream of the *CCA1* open reading frame (*13*). Under constant light condition, 2-hour dark pulses were applied to the plant circadian system with zeitgeber periods in the range, *T* = 16 to 32 hours (figs. S3 and S4). Figure 1A shows examples of individual-level *CCA1::LUC* oscillations perturbed by periodic stimuli. For *T* = 24 hours, the circadian system immediately entrained to the stimuli, but for *T* = 20 hours, it started to entrain after a long transient period (*t* > 192 hours) under decreasing the amplitude (Fig. 1, A to C). For each perturbation, the phase shift Δθ caused by the pulse stimulus was calculated as the difference in peak timing of *CCA1::LUC* expression before and after the dark pulse administration (as detailed in Materials and Methods). On the basis of Δθ measured for individual plants perturbed at various phases, θ, PRC shows a broad structure with unexpectedly large dispersion (Fig. 1D and fig. S5). As shown in Fig. 1E, dependence of Δθ on the oscillation amplitude, *A*, shows that the range of the phase shifts is significantly reduced as the amplitude increases, implying that the system becomes less phase-sensitive as the oscillation gets larger. To distinguish cases of large amplitudes from those of small amplitudes, the PRCs for large (*A* ≥ 0.5) and small (*A* ≤ 0.2) amplitudes are drawn separately in Fig. 1 (F and G, respectively). Variability in phase shift was drastically reduced by sorting the data by amplitude. The large amplitude plot shows a type 1 PRC characterized by a continuous transition from phase delay to phase advance without any breakpoint and a relatively small phase shift (maximally 0.1 to minimally −0.15 rad/2π), indicating that the 2-hour dark pulse has a relatively weak effect on the large amplitude state. The phase shift switches from positive to negative around θ = 0.7 rad/2π (approximately 6 hours after the beginning of the subjective dusk, that is, midnight). This corresponds to a stable point toward which the circadian phase tends when the dark pulse is administered. In the sense that this PRC reproduces essential features of the PRC measured by the standard technique (*7*), our method of measuring the PRC by periodic stimuli is considered to be reliable (fig. S6, A and B). For the small amplitude (*A* ≤ 0.2), the phase shifts showed a type 0 PRC (Fig. 1G), with a breakpoint (discontinuity) at θ = 0.2 rad/2π. In the vicinity of the breakpoint, significantly large phase shifts are induced, implying that the circadian system is very phase-sensitive in the small amplitude state. The linear profile of the PRC with a slope of −1 (Fig. 1G) indicates that any phase was moved to the stable point (θ = 0.7 rad/2π) within one circadian step. In other words, the system can be reset to the stable point by a single pulse input at any timing in the phase in the small amplitude state.

Note that the form of the PRC gradually changes from type 1 (continuous) to type 0 (discontinuous) as a function of amplitude, whereas basic features, that is, positions of the stable and unstable points, are preserved (Fig. 1, F and G). The range of the phase shifts, namely sensitivity of the phase response, also increases gradually as the amplitude decreases (Fig. 1E). A rapid transition from type 1 to type 0 takes place at around *A* = 0.2. This transition was observed under different conditions of the dark pulse periods *T* (fig. S7), implying that it is independent of the details of our measurement. Moreover, consistent results were also obtained by the standard method (fig. S6).

A theoretical framework was developed by considering a network of *N* cellular oscillators with phase, ϕ_{j}. In accordance with previous experimental studies (*14*–*16*), the intercellular coupling in plant circadian system was assumed to be very weak and negligible. As the main variables determine the network amplitude, phase and amplitude of individual oscillators play a major role. Although individual amplitudes have a direct influence on the network amplitude, the level of the cellular synchrony also determines the amplitude of the network dynamics. For instance, a gradual decrease in the oscillation amplitude of plant circadian systems under constant conditions is known to be partially due to desynchronized oscillations of the individual circadian cells (*14*–*20*). Therefore, our study focused on network synchrony as the key parameter controlling the network amplitude. By assuming that individual cells have unified amplitude, the network phase, θ, and amplitude, *R*, are defined as (*21*).

For the sake of simplifying the model, we defined the cell networks as having fully synchronized and fully desynchronized populations. In the fully synchronized population, each cell is assumed to be in the same phase. On the other hand, uniformly distributed phases are assumed for the desynchronized population. Individual cells have the same phase response property, the function of which was derived by the least-squares fitting of the experimental PRC data for large amplitude state (Eq. 16; see Materials and Methods for detailed function of the cellular PRC). Under this condition, when all circadian cells are completely synchronized (*R* = 1), the population follows the same PRC. Thus, the collective PRC coincides with the cellular PRC. In contrast, when the circadian cells are fully desynchronized and their phases are uniformly distributed, a large portion of the cells moves toward the stable point by following the cellular PRC [that is, the stable point θ = 0.7 rad/2π in Fig. 1 (F and G)]. As a result, the fully desynchronized population is somewhat clustered toward the stable point (fig. S8), moving the phase of the network to the stable point (note that only the averaged phase is moved to the stable point, whereas the individual cellular phases can be still widely distributed). Thus, in the fully desynchronized state (*R* = 0), the collective PRC gives rise to a line function with a slope of −1, which maps the arbitrary phase to the stable state within one iteration. Between the two extremes (*R* = 1 and *R* = 0), a gradual transition of PRC from type 1 to type 0 is observed. As described in detail in Materials and Methods, the PRC can be described by an analytic formula (Eq. 15). As demonstrated in the three-dimensional space (Fig. 2B and fig. S9), an experimentally observed transition of cellular PRC from type 1 to type 0 can be well reproduced by the theory. In both experiment and theory, a large amplitude corresponds to the PRC of type 1. As soon as the network enters into a small amplitude regime (*R* < 0.2), the network PRC changes to type 0. The stable and unstable points are located around θ = 0.7 rad/2π and θ = 0.2 rad/2π, respectively, for all amplitudes.

Although the analytical solution is derived under simplified conditions, that is, uniform phase distribution for desynchronized population and uniformly of the same phase for the synchronized population, essentially the same results were obtained for the more relaxed conditions as far as the phase distribution was unimodal (for example, Gaussian), and the cellular PRC had a single peak. Our theory is quite general in this sense.

To take into account the effect of cellular heterogeneity and the effect of intercellular coupling, the following phase oscillator model (*15*) was further simulated(1)

Here, ω_{j} stands for natural frequency of the *j*th oscillator, *K* is the coupling strength, and *L*(*t*) represents the light intensity [*L*(*t*) = 0 for light off and *L*(*t*) = 1 for light on]. The phase sensitivity function *Z*(θ) for the dark pulse is described in Eq. 19, which was based on the PRC determined by experimental data (details are described in Materials and Methods). Essentially, the same PRC as the analytical ones has been reproduced for 2-hour dark pulse stimuli (see fig. S10).

Finally, we simulated the entrainment experiment. In a similar manner as the experiment, 2-hour dark pulses were applied to the phase oscillator model of Eq. 1 with zeitgeber periods ranging from *T* = 16 to 32 hours. As an index of entrainment, the difference *T* − <τ> between dark pulse period *T* and oscillation period <τ> was calculated (Fig. 3, C and D, and fig. S11). On the basis of experimental data, the *CCA1::LUC* oscillations for a large amplitude state [*t*_{p}^{(2)} ≤ *t* ≤ *t*_{p}^{(7)}, where *t* and *t*_{p}^{(k)} represent measurement time and *k*th peak time of the bioluminescence signal] and a small amplitude state [*t*_{p}^{(9)} ≤ *t* ≤ *t*_{p}^{(14)}] were shown in Fig. 3 (A and B, respectively). For the large amplitude state, an entrainment with zero period difference, that is, *T* − <τ> ≅ 0, is observed from *T* = 21 to 25 hours. In the outer region of this entrainment (*T* < 21 hours and 25 hours < *T*), the period difference is nearly equal to *T* − τ_{0}, where τ_{0} = 23 hours is the endogenous oscillation period of *A. thaliana* under continuous light. On the other hand, in the small amplitude state, the period difference tends to be small overall, and the entrainment range is extended to 18 hours < *T* < 30 hours. This is because small amplitude weak oscillators are generally more easily entrained to external cycles than the large amplitude strong oscillators. The model results agree quite well with the experimental results.

## DISCUSSION

In summary, the present study showed that the oscillation amplitude, which can vary among individuals and can also change over time, has a profound effect on the form of the PRC. This suggests that the large variability, which often makes estimations of PRC quite noisy in circadian clocks, could be due to an uncontrolled level of amplitude. For precise characterization of the PRC, it should be of significant importance to simultaneously measure the phase and the amplitude.

Here, we used periodic 2-hour dark pulses to measure the PRCs. Compared to the standard method based on a single pulse perturbation to individual plant, our approach has the advantage of requiring fewer individuals for the measurement because multiple stimuli applied to each individual provide abundant information on the phase response properties. Under the condition that the effect of the dark pulse is not too strong and the pulse period is not too short, our method was shown to be capable of precisely measuring the PRCs (see also Materials and Methods). Owing to the vast amount of phase shift data collected from the periodic pulse experiments, the transition between type 1 and type 0 PRCs has been successfully observed.

In our experiment, the main cause of increasing the variability of the PRC was the decrease in the oscillation amplitude of the plant circadian system. As the experiment proceeds, the amplitude gradually decreased, making the system more phase-sensitive and consequently enlarging the entrainment range. There are two scenarios of the decreased amplitude in plant circadian clocks: (i) desynchronization of the cellular oscillations and (ii) damped oscillations in the individual cells. The desynchronization scenario has been supported by many experimental results from imaging studies of clock gene expression in plant organs (*14*–*20*, *22*, *23*). Phase waves in leaf and root (*14*, *16*, *22*, *23*), differences of inherent periods among tissues (*24*, *25*), and spatiotemporal analytical data (*14*–*16*) are reported to result in desynchronized cellular oscillations. Our preceding study has also shown that the level of synchrony decreases as time grows, implying that desynchronization is one of the major factors in reducing the oscillation amplitude (*15*). On the other hand, damped oscillations in individual cells have also been observed by single-cell imaging (*26*). A combination of the two mechanisms should describe changes in the oscillation amplitude. Our working model is based on the scenario of desynchronized cellular oscillations, and additional effects of damping in individual cells, as found in experimental findings, may further refine this model.

In our theory, we assume that individual cells have exactly the same cellular PRC independent of network dynamics. This implies that a change in PRC at the network level does not result from the internal changes of the individual cells but from the desynchronized state of the cellular oscillators. Similar phenomena may occur not only in the plant circadian clocks but also in other systems of oscillator networks. In the circadian clock of mouse, it has been reported that intercellular coupling governs rigidity and entrainment range, and the individual-level PRC changes from type 1 to type 0 by a mutation that reduces the circadian oscillation amplitude (*27*–*29*). When the cell population is completely synchronized with a large oscillation amplitude, the network PRC becomes identical to the cellular PRC. This implies that cellular PRC can be assessed by measuring the phase responses of an associated cellular network that exhibits large oscillation amplitudes. For this reason, amplitude should be carefully examined simultaneously with quantification of PRC. These approaches should be of great importance for plant and animal circadian clocks for which it is difficult to observe phase responses at a cellular level.

Here, the plant circadian clock showed a weak type 1 response to the dark pulse when the system maintained a large oscillation amplitude. As the time proceeds, the plants transitioned to a small amplitude oscillation, making the system phase-sensitive and capable of entraining to dark pulses with an extended range of the zeitgeber period. This suggests that the plant clock may achieve both robustness to external perturbations and high responsiveness to weak stimulus by changing its PRC by tuning the level of synchrony as well as the amplitude of oscillations.

## MATERIALS AND METHODS

### Growth and bioluminescence monitoring conditions

Plants were grown on 4 ml of gellan gum–solidified Murashige and Skoog plant salt mixture with 2% (w/v) sucrose in a 40-mm dish under 12-hour light (100 μmol m^{−2} s^{−1} of fluorescent white light): 12-hour dark cycles at 22° ± 0.5°C for 14 days. The plants were treated with 500 ml of 1 mM luciferin (made up in water) 24 hours before the start of bioluminescence monitoring.

Bioluminescence measurements were carried out at 22° ± 0.5°C under a light-emitting diode (LED) illumination with red LED (80 μmol m^{−2} s^{−1}; λ_{p} = 660 nm) and blue LED (20 μmol m^{−2} s^{−1}; λ_{p} = 470 nm). Dark pulses with zeitgeber period *T*, dark period (2-hour darkness), and light period (*T* − 2-hour illumination) were repeatedly administered to 40 plants for each treatment. The zeitgeber period was set to *T* = 16, 18, 20, 21, 22, 23, 24, 25, 26, 27, 30, and 32 hours. Monitoring was carried out for 7 days for *T* = 18 and 27 hours and 14 days for all other treatments using an automated monitoring system, Kondotron, developed by Kondo *et al*. (*30*), with a measurement interval of 20 min.

### Data processing

Peak detection. For bioluminescence signal {*c*_{l} : *l* = 1, 2, … }, the slope *a*_{l} at the *l*th data point was calculated as(2)

Given a sampling time interval of 20 min and the shortest zeitgeber period in this study of 16 hours, the moving average was set to *dw* = 24. The peak point of the bioluminescence oscillation was defined as the point at which *a*_{l} ≥ 0 and *a*_{l + 1} < 0, whereas the bottom point was defined as the point at which *a*_{l} ≤ 0 and *a*_{l + 1} > 0.

Definition of phase. Setting the *k*th peak time of the bioluminescence signal as *t*_{p}^{(k)} and the natural period of *CCA1* to τ_{0}, the phase, θ, at time, *t*, was defined as(3)

Here, τ_{0} = 23.0 hours is the free-running period averaged over *t*_{p}^{(2)} ≤ *t* ≤ *t*_{p}^{(6)} under constant light. Because the average time from peak to bottom was 12.0 hours, the phase of the bottom points was θ = 12 of 23 × 2π rad.

Definition of phase shift. Let the time at which a dark pulse is administered be *t*_{dp}. For two consecutive peaks *t*_{p}^{(k)} and *t*_{p}^{(k + 1)} that satisfy *t*_{p}^{(k)} ≤ *t*_{dp} ≤ *t*_{p}^{(k + 1)}, the phase shift Δθ induced by the dark pulse was calculated as(4)

The peak points appearing from 2 hours before the dark pulse to 4 hours after the dark pulse were not used in the calculation to exclude the transient response. The phase shift from the bottoms was calculated in the same manner.

Definition of amplitude. Normalized amplitude *A*(*t*) was calculated as(5)where *p*(*t*) and *b*(*t*), respectively, represent the peak and bottom of the bioluminescence signal at time *t*, which were obtained by linear interpolation of the observed peaks and bottoms.

Normalized bioluminescence was then calculated as(6)where *t*_{l} represents the time at *l*th data point.

### Theoretical analysis of collective phase response change

We assumed that the number of plant circadian cells, *N*, is sufficiently large (*N*→*∞*). In this case, distribution of the cellular phases is given by a probability density function. For simplicity, the circadian cells were divided into two populations: fully synchronized population and fully desynchronized population. In the fully synchronized population, all cells are in the same phase, θ_{1}. In the fully desynchronized population, the phases are assumed to be uniformly distributed. Taking the two populations together, the phase density function is given by(7)

Here, *x* represents cellular phase, *R*_{1} determines the phase coherence of the cellular network, and δ(*x*) stands for the Dirac delta function. The density function, , satisfies the following normalization condition(8)

According to Kuramoto (*21*), for a network of *N* cellular oscillators, the network amplitude *R* and phase θ are defined as , where φ_{j} represents the phase of the *j*th cell and *i* is the imaginary unit (*i* = ). Applying the density function, , to this formula and taking the limit of *N*→*∞*, we obtained(9)

The above equation indicates that the coherence parameter, *R*_{1}, and the phase, θ_{1}, of the synchronized population correspond to the network amplitude and phase, respectively.

Let the phase response function for the 2-hour dark pulse be *g*(*x*). Then, the individual cellular phase *x* was shifted to *x* + *g*(*x*) by the dark pulse stimulus (no intercellular coupling is assumed). The dark pulse transforms the network amplitude and phase to *R*_{2} and θ_{2} as(10)

Note that the integral represents a constant value, which is determined by *g*(*x*). Here, we denote it as(11)

By using Eq. 10, we obtained the phase shift Δθ = θ_{2} − θ_{1} as(12)

Hence,(13)

By dividing the first equation by the second one, we obtained(14)

Finally, the phase shift was described by the function of the network amplitude and phase as(15)*g*(*x*) was estimated by the least-squares fitting of second-order Fourier expansion to the data points with *A* ≥ 0.6 as follows(16)

In this case, the *g*(*x*)-dependent constants were calculated as θ′ = 4.39 rad and *R*′ = 0.134.

### Numerical simulation of the phase oscillator network

Simulation of the phase oscillator network in Eq. 1 was performed, as shown in Fig. 3 (C and D) and figs. S10 and S11. Cell number and coupling strength were set to *N* = 1000 and *K* = 0.01, respectively (*15*). Natural frequencies of the cells were normally distributed as *N*[2π/τ, (0.05 × 2π/τ)^{2}], where τ is the natural period and τ = 23 hours. Under the initial conditions, all cells were in the same phase (*R* = 1), and the initial phases were uniformly distributed among all phases (θ = 0 to 2π rad). The zeitgeber period was set to *T* = 14, 15, …, 38 hours as was used in the experiments. In Fig. 3 (C and D) and fig. S11, the average of the initial phase was −2 × 2π/τ. In fig. S10, θ and *R* are the phase and the network amplitude at which a dark pulse was started, and Δθ was first calculated as the phase difference between the phases at the start and end times of the dark pulse and then the difference in those phases minus 2 × 2π/τ.

For the simulation of Eq. 1, the phase sensitivity function *Z*(θ) to the dark pulse was derived from the phase response function *g*(*x*), which is defined as having a 2-hour dark pulse and shows the following relationship(17)

Replacing *t* with *t* − 1 yields(18)

Because the 2-hour duration is sufficiently short compared to the plant natural period, τ, under the constant light condition, approximations were made as θ(*t* − 1) ≅ θ(*t*) − 2π/τ and . Finally, the phase sensitivity function *Z*(θ) was approximated as(19)

### Change of phase distribution by dark pulse perturbation

The phase distribution function before and after the application of a dark pulse is given by *f*_{1}(*x*) and *f*_{2}(*x*), respectively, for phase *x* ∈ [0, 2π). When the phase response function of a dark pulse is given by *g*(*x*), the phase after the application of a dark pulse is *x* + *g*(*x*). We assumed that the dark pulse is a weak perturbation, such that *dg*(*x*)/*dx* > −1 for all *x*. The integrated values of *f*_{1}(*x*) and *f*_{2}(*x*) are equal in the intervals before the dark pulse [*a*, *x*] and after the dark pulse [*a* + *g*(*a*), *x* + *g*(*x*)] for an arbitrary constant *a*, 0 ≤ *a* < *x*. Hence(20)

Differentiating both sides of the equation with respect to *x*(21)gives(22)

In the case that *f*_{1}(*x*) = 1 and *g*(*x*) are given by Eq. 16, the phase distribution is as shown as fig. S8.

### PRC measurement from extremely small amplitude data by the standard method

Our method based on periodic dark pulses has been applied to measure PRCs with various oscillation amplitudes. The result obtained for large amplitude was compared with the PRC measured by the standard method in our previous study (*7*). To compare the results for extremely small amplitude, the corresponding PRC was measured by the standard method in fig. S6C. Conditions on the growth of plants and their bioluminescence monitoring were set the same as those of the periodic pulse experiments. The monitoring was carried out for 14 days, and the 2-hour dark pulse was applied 262 hours after the start of measurement. In the analysis of bioluminescence signals, the peak points were detected in the same manner as in the case of periodic pulse experiments. The phases before and after the dark pulse were obtained by linear interpolation of the three (or two) peak points before and after the pulse by the least-squares method.

### Verification of our proposed method

Here, we applied periodic 2-hour dark pulses to measure PRCs. Basic requirements for the method of periodic pulse perturbation to precisely measure the PRC are summarized as follows. (i) Effect of the dark pulse perturbation should be weak enough so that the perturbed dynamics does not deviate much from the unperturbed limit cycle trajectory (*21*). (ii) Relaxation time of the system against perturbation should be short compared to the dark pulse period *T* (*12*). (iii) The phase should not be so tightly locked to that of the dark pulse so that the phases at which the dark pulses were injected may not too densely concentrate on a particular value.

Here, the experimental data were collected in such a way that these requirements were fulfilled. Concerning the requirement (i), the 2-hour dark pulse was chosen because we had confirmed that the 4-hour dark pulse was a little too strong, and the perturbed plant oscillations deviated from the unperturbed ones. The obtained PRC was relatively noisy compared to the theoretical PRC (fig. S12). On the other hand, the 1-hour dark pulse was too weak to detect the phase shift. The moderate strength of the 2-hour dark pulse was found optimal to draw the PRCs. Concerning the requirement (ii), the minimum period of the dark pulses was chosen to be *T* = 16 hours because any period shorter than *T* = 16 hours lowered signal-to-noise ratio of the PRC. Concerning the requirement (iii), we had confirmed that the phase values obtained in our PRCs were randomly distributed without much concentration on a particular phase value (fig. S5).

To further examine the reliability of the present method for measuring the PRC, we compared it with the standard method by numerical simulations. To quantify the accuracy of the measured PRC, we had used the van der Pol equation (*31*) as a prototypical model, the PRC of which has been known very well. The advantage of using the mathematical model is that, when the equation of motion is explicitly given, the PRC can be accurately calculated using the so-called adjoint method (*32*), and the calculated curve can be used as the “true” PRC. Our numerical study was carried out as follows.

In the standard method of estimating the PRC, a single pulse was injected to the van der Pol equation as(23)where μ stands for a constant parameter (μ = 1) and *F*(*t*) represents the single pulse input (red bars in fig. S2, A and B). The phases before and after the pulse injection were calculated from the simulated time waveform {*x*(*t*)}. The phase at which a pulse stimulus was injected was calculated by linear interpolation of the three zero-crossing points before the pulse input. In a similar manner, the phase after the pulse input was computed by the linear interpolation of the three zero-crossing points after the pulse input, where the difference between the phases before and after the pulse injection provides the phase shift. To discard the transient, the time interval of 2 after the pulse input was not used for the phase calculation.

In our proposed method, periodic stimuli were applied to the van der Pol equation as (fig. S2B)(24)

The zero phase is defined as the point at which the variable *x* turns from negative to positive. In a similar manner as the dark pulse experiment that used both peaks and bottoms, the phase shift was calculated by using both points at which *x* turns from negative to positive and from positive to negative. The forcing period *T* was varied from 6.8 to 10.0 with an increment of 0.2. Duration of the dark pulse stimulus was set to 0.2. For each forcing condition, the simulation time was 65 with a fixed initial condition of *x*(0) = 2 and *v*(0) = 0.

In the adjoint method (*32*), the phase sensitivity function {*Zx*(θ), *Zy*(θ)} was calculated as the solution to the adjoint problem(25)where *f*(*x*, *y*) = {*dx*/*dt*, *dy*/*dt*} and *Df* is the Jacobian matrix of *f*. Here, the van der Pol equation is rewritten by using variables *x* and *y* as(26)

*Zx*(θ) and *Zy*(θ) were calculated by numerically integrating Eq. 25 backward in time, and the PRC subject to pulse stimulus (stimulus duration of 0.2) was calculated as(27)where τ ≅ 6.666 represents the natural period of the van der Pol oscillator. Figure S2C shows the results. The two curves obtained by the standard method and our method agree quite well with the one obtained by the adjoint method. This indicates a high reliability of the present method to estimate the PRC.

### PRC to the 4-hour dark pulse

To examine the dependence of the PRC measurement on the dark pulse duration, the PRC was measured by periodic application of 4-hour dark pulses. The growth condition was set to be the same as that of the 2-hour dark pulse experiments. The bioluminescence measurements were carried out at 22° ± 0.5°C under a LED illumination with red LED (80 μmol m^{−2} s^{−1}; λp = 660 nm) and without blue LED. The zeitgeber period was set to *T* = 21, 24, 28 (two experiments; #1 and #2), 32, and 36 hours. The monitoring was carried out for 10 days for *T* = 21, 24, 28 (#1), and 32 hours; 11 days for *T* = 28 (#2) hours; and 14 days for *T* = 36 hours. The first dark pulse was applied at *t* = 12 hours for *T* = 21 hours; *t* = 16 hours for *T* = 28 (#2) hours; *t* = 18 hours for *T* = 24, 28 (#1), and 32 hours; and *t* = 32 hours for *T* = 36 hours. The free-running period under constant light was τ_{0} = 22.5 hours, and the average time from peak to bottom was 11.4 hours.

### Circadian phase at the first dark pulse

Examples of bioluminescence signals grown under LD12:12 cycles and then stimulated by periodic application of dark pulses are shown in fig. S7 (A and C). At the starting time *t* = 0 hour of the measurement, the circadian phase is few hours before the peak expression of *CCA1::LUC*. The first dark pulse was then injected at *t* = *T* − 2 hours, which was determined with the value *T* of the cycle period. The phase at which the first dark pulse was applied was therefore dependent on *T*.

The circadian phase at the first dark pulse is important when measuring the phase shift for large amplitude rhythms. This is because the large amplitude state appears only in the initial stage of the measurement due to the damped oscillations inherent in the plant. The phase response property of the large amplitude state can be measured mainly with the first dark pulses. If the first dark pulse is injected at the same timing *t* for all conditions of dark pulse period *T*, the PRC should concentrate densely on the phase of the first pulse injection. To avoid this biased estimate, we had varied the phase of the first dark pulse by injecting the pulse at *t* = *T* − 2 hours so that the phase shift at large amplitude can be measured for various phases (fig. S5).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/3/10/e1700808/DC1

fig. S1. Schematic illustration of the PRC measurement.

fig. S2. Verification of PRC measurement method by the van der Pol oscillator.

fig. S3. Normalized bioluminescence for various dark pulse periods.

fig. S4. Normalized amplitude for various dark pulse periods.

fig. S5. Experimental data of phase response by dark pulse periods.

fig. S6. PRCs for the 2-hour dark pulse obtained by the standard method (one-time stimulus) and our method (periodic stimuli).

fig. S7. PRCs constructed from periodic stimuli with short and long periods of dark pulses.

fig. S8. Phase distribution after the dark pulse was applied to the uniform distribution.

fig. S9. Theoretical PRCs.

fig. S10. PRC simulation using a phase oscillator model (Eq. 1).

fig. S11. Range of entrainment in simulation using a phase oscillator model (Eq. 1).

fig. S12. PRC of the 4-hour dark pulse.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We are grateful to T. Kondo for providing the bioluminescence monitoring device and N. Nakamichi for supplying the

*CCA1::LUC*plants.

**Funding:**This study was supported, in part, by Grants-in-Aid for Scientific Research (nos. 25712029 and 16H05011) and Precursory Research for Embryonic Science and Technology (no. JPMJPR15O4) of the Japan Science and Technology Agency (to H.F.).

**Author contributions:**K.M., R.K., and H.F. designed the experiments. K.M., R.K., and K.U. simulated the model and analyzed the data. K.M., H.F., and I.T.T. wrote the manuscript. All authors discussed the results and commented on the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2017 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).