Integrating oscillator-based circadian clocks with crop growth simulations

Abstract

Circadian rhythms play critical roles in plant physiology, growth, development and survival, and their inclusion in crop growth models is essential for high-fidelity results, especially when considering climate change. Commonly used circadian clock models are often inflexible or result in complex outputs, limiting their use in general simulations. Here we present a new circadian clock model based on mathematical oscillators that easily adapts to different environmental conditions and produces intuitive outputs. We then demonstrate its utility as an input to Glycine max development models. This oscillator clock model has the power to simplify the inclusion of circadian cycles and photoperiodic effects in crop growth models and to unify experimental data from field and controlled environment observations.

1. INTRODUCTION

Circadian rhythms were first discovered in a plant (probably Mimosa pudica) nearly three centuries ago (de Mairan 1729). Since that time, circadian clocks have been shown to play key roles in plant physiology, development and survival (Dodd et al. 2005; McClung 2006). Seemingly disparate behaviours such as hypocotyl elongation (Dowson-Day and Millar 1999; McClung 2001; Webb 2003), pulvinar leaf movements (McClung 2001; Webb 2003), jasmonic acid production (Hsu and Harmer 2014), stomatal opening (Webb 1998, 2003), root permeability (Clarkson et al. 2000; Webb 2003), starch consumption (Hsu and Harmer 2014; Izumi 2019), shade response (Salter et al. 2003) and flowering time (Song et al. 2015) have all been linked to clock function in various plants. In the particular case of Glycine max (soybean), transgenic plants with a light-sensing gene from Arabidopsis thaliana exhibit altered transcript abundances of several native clock genes, extended pod development and increased seed yields, demonstrating a clear link between circadian rhythms and food production (Preuss et al. 2012).

Given the importance of circadian clocks to plants, significant effort has been made to incorporate them into growth and production models, usually following one of two strategies. In the indirect approach, clocks are included by determining plant behaviour according to daily sunrise times and photoperiod lengths (seeSection S1 of the Supporting Information for a discussion of this terminology), often calculated from latitude and longitude using celestial mechanics (Goudriaan and Van Laar 1978; Penning de Vries and Laar 1982). Such models assume that plants are able to sense patterns in daily solar radiation. It is rare for these models to mention circadian clocks explicitly, yet it is well established that clocks are required for this sensitivity (Greenham and McClung 2015; Song et al. 2015). In the direct approach, clocks are included via systems of differential equations representing gene networks (Chew et al. 2014, 2017; Zardilis et al. 2019). The indirect method has an advantage over gene networks because its clock outputs—sunrise time and photoperiod length—can be readily understood and straightforwardly used as inputs to other model components. On the other hand, gene networks can be more adaptable because they represent the intrinsic response of a plant to illumination and are therefore independent of the light source. Furthermore, they provide the fundamental link between genetics and plant behaviour.

Here we present a method that combines the simple and intuitive clock outputs of the indirect approach with the adaptability of a gene network. We first describe the form of this clock model, which is based on mathematical oscillators. After testing its reliability and accuracy, we use the oscillator-based clock to estimate G. max flowering times from experimental weather data, demonstrating its utility for predicting an important determinant of soybean yield. More generally, since a clock can be used to control many behaviours within a comprehensive plant model, this simple clock provides a means to simultaneously reduce parameters and predict responses more flexibly by providing a straightforward and robust way to respond to the lighting environment.

2. MATERIALS AND METHODS

Surface radiation and temperature data were obtained from two publicly available data sources: WARM’s Champaign, IL weather monitoring station (40°05′02″N, 88°22′23″W) and SURFRAD’s Bondville, IL station (40°03′07″N, 88°22′23″W) (Illinois State Water Survey 2015; Augustine et al. 2000). WARM provides data at daily time steps and full yearly data sets are available from 1990 to 2019. SURFRAD provides data at 1- or 3-min time steps and full yearly data sets are available from 1995 to 2019. Data were converted to hourly time steps using the procedures described in Section S2 of the Supporting Information.

Numerical integration of systems of differential equations was performed using version 1.71 of the Boost odeint library (‘Boost Version 1.71.0’ 2019). Unless noted otherwise, fifth-order Runge-Kutta integration with fourth-order Cash-Karp error estimation was used with absolute and relative error tolerances set to 1 × 10⁻⁶ for all simulations, and all calculations were performed using the default parameters in Tables 1 and 2. While integrating, weather data were treated as a continuous piecewise function defined by linear interpolation between the hourly values.

3. USING OSCILLATORS TO MODEL A CIRCADIAN CLOCK

Real-world clocks use physical oscillators such as pendulums or piezoelectric tuning forks to keep time (Hebra 2010), and the same approach can be used to model circadian clocks with mathematical oscillators (Gonze 2011). Examples include phase-only oscillators (Roenneberg et al. 2010) or phase-amplitude oscillators such as the Poincaré oscillator (Glass and Mackey 1988; Winfree 2001). In general, these oscillators possess a natural period $T0$ but can be coupled to an external signal with a different period T. Under certain circumstances, the resulting oscillator trajectories will essentially exhibit the external signal’s periodicity. Thus, the oscillators exhibit the same important features that characterize a classical circadian rhythm: their cycles are endogenous in the sense that they oscillate with a natural period even without an external driver, and they are entrainable in the sense that they can adopt the period of a driving signal and become locked to its phase.

However, in addition to adopting the Earth’s 24-h day length, plant circadian clocks also respond to the daily photoperiod. At a molecular level, this sensitivity is accomplished via groups of genes that are primarily expressed in either the morning or evening (Greenham and McClung 2015). The core idea of our model is to mimic this functionality by including two oscillators that separately respond to the beginning and end of the daily light exposure. As shown schematically in Fig. 1A, the dawn oscillator receives a signal at the beginning of the photoperiod which alters its phase angle $φdawn$ (represented as an arrow rotating counterclockwise). When fully entrained, $φdawn$ resets to zero at dawn. Likewise, the dusk oscillator resets its phase angle (⁠$φdusk$⁠) at the end of the photoperiod.

Figure 1B shows the results from 1 day of a simulation where both oscillators have run long enough to become fully entrained to the daily cycle of solar radiation. Since the phase angles increase at a constant rate, $φdawn$ and $φdusk$ are proportional to the time interval since sunrise and sunset, while the difference $Δ=φdawn−φdusk$ is proportional to the clock’s estimate of the photoperiod length (⁠$ρclock$⁠). Figure 1C shows $ρclock$⁠, sunrise time and sunset time determined by the oscillator-based clock throughout an entire year, where phase angles have been converted to time intervals. (A phase change of $2π$ rad corresponds to 24 h, e.g. $ρclock=24⋅Δ/2π$⁠.)

4. MATHEMATICAL DESCRIPTION OF THE OSCILLATOR CLOCK MODEL

4.1 Poincaré oscillators

In principle, any entrainable oscillator model can be used as the clock base. Due to its simple form and ubiquity in the field of chronobiology (Glass and Mackey 1988; Winfree 2001), we have chosen to use the Poincaré oscillator, which is defined by the following set of differential equations:

where $R$ and $φ$ represent the radius and phase angle of the oscillator in phase space, respectively. The general behaviour of this oscillator can be understood directly from these equations. If $R$ exceeds $R0$⁠, its rate of change is negative. Likewise, $R$ increases when $R<R0$⁠. In other words, $R$ tends towards $R0$ regardless of its initial state, and its rate of approach is determined by the friction coefficient $γ$ (seeSection S4 of the Supporting Information for a justification of this terminology). $φ$ advances at a constant rate set by $ω0$⁠, the oscillator’s natural frequency. Consequently, changes in phase angle are proportional to changes in time. As an angle, the phase effectively resets after a change of 2π rad, which requires a time interval of $T0=2π/ω0$⁠, referred to as the oscillator’s natural period.

A Poincaré oscillator can be coupled to an external driving signal $Z(t)$ by projecting $R$ and $φ$ onto rectangular coordinates (⁠$a=Rcos⁡φ$ and $b=Rsin⁡φ$⁠) and adding $Z(t)$ to one of the resulting derivatives, which results in a new set of equations:

If $Z(t)$ is periodic with period $T$ and sufficiently strong, the resulting oscillator trajectory will be periodic with the same period rather than its natural period $T0$⁠.

As an illustrative example, we integrate Equations (3) and (4) to determine the phase space trajectory (Fig. 2A) and phase evolution (Fig. 2B) of a Poincaré oscillator before, during and after the application of a sinusoidal $Z(t)$⁠, where the oscillator begins with the initial state $φ=0$⁠, $R=R0/2$⁠. Since $R<R0$ to start, the radius gradually increases during the beginning of the simulation. Eventually, this transient behaviour fully decays and the oscillator reaches a stable orbit, moving counterclockwise around a circle of radius $R0=1.0$ at a rate of $ω0=π/15$ rad h⁻¹. When the driving signal is applied, the phase space radius increases and the trajectory becomes elliptical, but the phase advances at a constant rate of π/12 rad h⁻¹ (set by the driving force). When the signal is removed, the oscillator reverts back to its original path. The friction coefficient (⁠$γ=0.1$⁠) has been chosen so that the trajectory stabilizes within one cycle following these changes in the driving signal.

4.2 Coupling oscillators to light

In the clock model, two driving signals (⁠$Zdawn$ and $Zdusk$⁠) are derived from photosynthetic photon flux density (PPFD) values $Φ(t)$ and separately coupled to the dawn and dusk oscillators. To determine these signals, $Φ$ is first converted into a binary indicator of the presence of light (⁠$L$⁠) which smoothly but sharply transitions between 0 and 1 as $Φ$ increases from 0 to $ΦT$ according to a logistic formula

where $ΦT$ is a PPFD threshold. When $Φ=0$⁠, $L=1/(1+e10)≈e−10$ and when $Φ=ΦT$⁠, $L=1/(1+e−10)≈1−e−10$⁠.

Next, $L$ is used to determine the values of a day tracker $D$ and night tracker $N$⁠, quantities that evolve according to

and

During the day, $L=1$⁠, so $D$ and $N$ reach equilibrium values of 1 and 0, respectively. Likewise, during the night, they equilibrate to 0 and 1. The parameter $αtracker$ determines the length of the transition between equilibrium values, which occurs as an exponential approach; for example, if $αtracker=2.0$ h⁻¹, the transition is ~98 % complete after 2 h (⁠$1−e−2⋅2.0≈0.98$⁠). $D$ can be thought of as a protein produced only in the presence of light that undergoes exponential decay at all times, with a similar analogy for $N$⁠.

Finally, $Zdawn$ and $Zdusk$ are calculated according to

and

where $αZ$ determines the overall magnitude of the signals. $Zdawn$ is appreciably non-zero only at the start of the day, as that is the only time when both N and L are non-zero (Fig. 3A). In this sense, $Zdawn$ can be thought of as a molecular clock component activated by $N$ protein in the presence of light. Likewise, $Zdusk$ is only non-zero at the start of night (Fig. 3B), analogous to a molecular clock component activated by $D$ protein in darkness.

Apart from a small spike around dusk, $Zdawn$ is only non-zero around dawn and consists of a roughly triangular peak lasting around 2 h (Fig. 3A). Therefore, for a realistic $Φ$ input, $Zdawn$ is nearly periodic with an approximately 24-h period. Likewise, $Zdusk$ (Fig. 3B) has a similar shape and is non-zero for about 2 h near sunset. A critical feature of this scheme is that because both driving signals are only non-zero during transitions between light and dark, they both become zero under constant light (⁠$Φ≡Φ0>ΦT$⁠) or constant dark (⁠$Φ≡0$⁠) conditions, allowing the clock to run at its natural frequency in the absence of an environmental light cycle.

Table 1 lists defaults for the clock model parameters and initial states. They have been chosen for reliability and accuracy, although good results can be obtained for other values. In order to start the clock in phase with the environment and ensure its accuracy during the first days of the year, the clock model was run for the entirety of 2019 using WARM weather data and the final states of the oscillators and trackers at midnight on the last day of the year were chosen as the default initial states. Additional insight into the influence of the remaining parameters on model output can be gained through sensitivity analysis [seeSection S8 of the Supporting Information].

5. OSCILLATOR CLOCK ENTRAINMENT, RELIABILITY AND ACCURACY

For use in crop growth simulations, an oscillator-based clock must have a stable output independent of its starting conditions, respond quickly to changes in the starting time and length of the daily photoperiod and accurately estimate sunrise time and photoperiod length. There is a trade-off between these traits, so we examine the behaviour of this model with our parameterization.

In general, a driven Poincaré oscillator exhibits a transient response at the beginning of its motion that depends on its initial conditions rather than the driving signal (Fig. 2A). For an oscillator-based clock, this means that its initial outputs will not be reliable estimates for sunrise time or photoperiod length. To investigate the impact of this transient response, we use an artificial light signal where daily PPFD values are given in units of μmol s⁻¹ m⁻² at hour $h$ by the Gaussian function $Φ(h)=1000⋅e−(h−12)2/(2σ)2$⁠. By keeping the light profile width constant at $σ=2.75$ h and varying the initial phase of the dawn oscillator (⁠$φdawn,i$⁠) between 0 and $2π$⁠, it is possible to determine an approximate value for $tind$⁠, the time at which the transient response has fully decayed (Fig. 4A). Although the $ρclock$ values calculated with different $φdawn,i$ vary greatly during the first few days, the differences are small after 20 days and all values are virtually identical by Day 30, i.e. $tind$ = 30 days. For all subsequent times $ρclock$ is determined only by the light profile and therefore represents the clock’s true estimate of the photoperiod length.

The clock also exhibits an adjustment period following a change in the starting time or length of the daily photoperiod. The duration of this adjustment period can be thought of as the time to achieve entrainment (⁠$tentrain$⁠). By suddenly changing the width of the Gaussian light profile from $σshort=2.75$ h to $σlong≥σshort$ after $tind$⁠, it is possible to estimate $tentrain$ as the time when $ρclock$ first lies within 1 min of its final value (Fig. 4B). For large jumps the entrainment process can take up to 20 days, but for the smallest jump (~30 min), full entrainment occurs within 10 days.

For PPFD profiles derived from sunlight measurements, the length of the daily photoperiod changes slowly but continuously (by less than 3 min per day at 40° latitude), precluding a direct observation of entrainment in response to a single change as in Fig. 4B. Instead, $ρclock$ is expected to lag behind the true length with a delay determined by the clock’s sensitivity to changes in the light profile. For small changes, the clock is able to respond within a few days (see Fig. 4B), so errors in $ρclock$ that occur because of this delay are expected to be small.

To test the accuracy of the oscillator clock model, the model was run for all years in the WARM and SURFRAD weather data sets. Hourly values derived from the daily WARM data follow smooth profiles, while hourly values derived from the sub-hourly SURFRAD data include more realistic noise. For each day of each year, the difference between the clock’s photoperiod length output and the photoperiod length calculated by celestial mechanics (⁠$Δρ=ρclock−ρcelestial$⁠) was determined. The average daily difference across all years in each data set is shown in Fig. 5A. For both the WARM and SURFRAD data sets, $ρclock$ typically deviates from $ρcelestial$ by less than 0.5 h, with a bias towards shorter photoperiod lengths, especially earlier in the year. Although it is small, this error is not random and results from a systematic underestimation of photoperiod length (due to non-zero $ΦT$⁠) and a slight delay in the clock response, as discussed previously. Since this error is systematic, it is not expected to preclude using the oscillator clock output to reliably control other aspects of plant growth.

Although they are a small source of systematic error, $ΦT$ and the delayed clock response help to produce a reliable and smooth output. This is apparent in the median, interquartile range and range determined from daily photoperiod lengths across all years in the WARM data set (Fig. 5B). The interquartile range is always smaller than 0.8 h and is frequently smaller than 0.5 h, demonstrating that the clock is relatively insensitive to random variations in photosynthetically active radiation (PAR) data between years. An analogous plot based on the SURFRAD data exhibits similar interquartile range values [seeSupporting Information—Fig. S6], showing that the clock model is robust enough to function well with noisier input data, a conclusion also supported by the similarity of the error curves in Fig. 5A [see Section S5 of the Supporting Information].

6. APPLICATION EXAMPLE: G. MAX FLOWERING

The year-to-year repeatability of the clock model evident in Fig. 5B suggests that its outputs would be suitable as inputs to other established plant models, with the caveat that recalculation of some parameters may be necessary. One such model is the G. max development model described by Grimm et al. (1993) (referred to as the CROPGRO model by Piper et al. (1996)), which determines flowering dates (⁠$DAYR1$⁠) from temperature and night length throughout the growth period. This model was parameterized using data from soybean plots across North America and has been shown to typically reproduce experimental R1 dates in that region to within 3 days (Grimm et al. 1993; Piper et al. 1996). The original CROPGRO model uses celestial mechanics calculations to determine night length (⁠$η=24−ρ$⁠) and includes five parameters that are allowed to vary between cultivars in order to fit the experimental data: minimum and optimal night length (⁠$ηmin$ and $ηopt$⁠), minimum and optimal temperature (⁠$Tmin$ and $Topt$⁠) and threshold physiological age for flowering (⁠$THf$⁠).

To demonstrate the utility of the oscillator clock as an input to the CROPGRO model, flowering dates were predicted using either celestial mechanics (⁠$DAYcelestialR1$⁠; seeSection S3 of the Supporting Information) or an oscillator clock (⁠$DAYclockR1$⁠) as the source for the night length values. For each year, five different sowing dates were tested (Days 120, 130, 140, 150 and 160; corresponding to late April through early June) using the parameter values in Table 2. The clock was run beginning on Day 1 of each year to ensure that its initial transient response would not have an effect on the results. To achieve the best overall agreement between flowering dates, the values of $ηmin$ and $ηopt$ used to calculate $DAYclockR1$ were varied using a brute-force approach on a finite grid to minimize the sum of $(DAYcelestialR1−DAYclockR1)2$ across all sowing dates and years in the WARM data set [seeSection S6 of the Supporting Information].

Flowering time was also predicted using the SURFRAD data set without reparameterizing the models, giving an out-of-training data set prediction as a check for robustness. The results of this comparison are shown in Fig. 6, where $DAYclockR1$ is plotted against $DAYcelestialR1$ for each combination of year and sowing date in the WARM (Fig. 6A) and SURFRAD (Fig. 6B) weather data sets. The dashed lines are guides to the eye representing perfect agreement. For the WARM data set, most $DAYclockR1$ values lie within 2 days of $DAYcelestialR1$ (the mean and standard deviation of the differences $DAYclockR1−DAYcelestialR1$ are 0.0 days and 0.9 days, respectively), and even for the out-of-training SURFRAD data set, differences are generally within 3 days [seeSupporting Information—Fig. S9]. These discrepancies are not larger than the CROPGRO model’s inherent accuracy, and therefore are not expected to introduce significant errors into its predictions. As indicated in Table 2, changes in $ηmin$ and $ηopt$ were required for this good agreement. These changes can be understood by noting that the increase in $ηopt$ results from the small but systematic underestimation of photoperiod length in the oscillator clock model, while the larger range (⁠$ηopt−ηmin$⁠) reduces the impact of random fluctuations in $ρclock$⁠.

7. DISCUSSION AND CONCLUSIONS

Many physiological processes in plants depend on the starting time or length of the daily photoperiod, so neglecting photoperiodic effects in plant growth models can have deleterious effects on model accuracy. For example, growing degree days (GDDs) depend solely on temperature and are sometimes used to estimate soybean development, where developmental stages are assumed to occur at particular GDD thresholds. To illustrate the importance of photoperiodic sensitivity in soybeans, we compare flowering dates calculated by either a GDD threshold (⁠$DAYGDDR1$⁠) or the CROPGRO model (⁠$DAYcelestialR1$⁠) using the WARM data set (Fig. 7A). Growing degree day values were calculated throughout the growing season using the model and parameters in Ruiz-Vera et al. (2018) (although hourly air temperatures were used in place of daily average canopy temperatures) and the threshold value for flowering (520 °C d) was chosen to yield the best overall agreement with the CROPGRO model [seeSection S7 of the Supporting Information]. Here it is apparent that the GDD method often deviates from the more accurate CROPGRO model by more than 3 days. However, there is a more serious issue: the difference between $DAYGDDR1$ and $DAYcelestialR1$ is strongly correlated with the average summer temperature, with the GDD model predicting earlier flowering dates during hotter years (Fig. 7B).

In the face of global climate change, it is therefore essential to include photoperiodic effects in plant growth and development models in order to avoid these types of systematic errors. Consequently, it is critical to develop a circadian clock model that is both intuitive to use and applicable to a broad range of plant environments. The commonly used clock models—the indirect method based on celestial mechanics calculations and the direct method using gene networks—do not meet these two important criteria in isolation.

For example, experiments have shown that hypocotyl elongation in A. thaliana is controlled by its clock and occurs primarily around dawn (Niwa et al. 2009). With an indirect clock model, it is simple to represent this behaviour as a time-dependent hypocotyl growth rate with a peak centred at sunrise. If instead the clock output consists of the time-dependent concentrations of the various proteins involved in the clock gene network, then this task becomes more difficult, especially without detailed knowledge of the regulatory pathway between the clock and the relevant genes involved in hypocotyl elongation.

On the other hand, the gene network method can easily adapt to different conditions in the environment or in the plant. For example, when comparing plants grown in a field to others in a growth chamber, only the light profile (an input to the model) needs to be changed. In contrast, the indirect method requires a change to the model itself, since sunrise time and photoperiod length can no longer be calculated from properties of the Earth’s orbit. The indirect method also cannot account for any possible dependence of clock functionality on nutrient availability (Haydon et al. 2015).

For solar illumination, small systematic differences exist between the oscillator clock’s output and the results of celestial mechanics calculations. There are multiple approaches for dealing with this disagreement. For an independently validated model that requires photoperiod length or sunrise time as an input, one route is to simply reparameterize it to produce the same results when using an oscillator clock or celestial mechanics calculations, as demonstrated in Fig. 6. This is analogous to synchronizing two clocks. Alternatively, it would be possible to directly fit the model to experimental data when they are available, bypassing any comparison to celestial mechanics. As a concrete example of this route, one could simply vary $ηmin$ and $ηopt$ in the CROPGRO model to achieve the best agreement with measured $DOYR1$ values while using an oscillator clock to determine photoperiod length. Finally, we note that the oscillator clock parameters could be additionally constrained by performing an experiment similar to the one in Fig. 4B, i.e. by measuring the time-dependent response of real plants to sudden changes in photoperiod under artificial illumination. Ultimately, accuracy in predicting observable clock-controlled plant behaviours such as $DOYR1$ is more important than reproducing celestial mechanics calculations, especially when considering that plants determine day length by sensing light, not calculating celestial mechanics.

The oscillator clock model presented here combines the best aspects of the two approaches commonly used in the literature. As with the indirect clock, the oscillator clock’s outputs are sunrise time and photoperiod length. Similar to the gene network clock, the oscillator clock functions independently of the light source and its parameters can be modified in response to nutrient availability as in a real plant. Small systematic differences exist between the oscillator clock and celestial mechanics calculations. However, these discrepancies can be addressed in several straightforward ways. A major advantage of the oscillator clock model is that it has the power to simplify the integration of experimental data from fields illuminated by the sun and growth chambers illuminated by artificial light sources.

SUPPORTING INFORMATION

The following additional information is available in the online version of this article—

Table S1. Numbers of weather data points determined from PAR, solar irradiance, and interpolation in the SURFRAD data set.

Table S2. Numbers of weather data points determined from insolation and interpolation in the WARM data set.

Table S3. Error values determined for different combinations of ηmin and ηopt, where the entry with the lowest error is highlighted with a gray background.

Figure S1. Percentage of weather data points determined from PAR, solar irradiance, and interpolation in the SURFRAD data set.

Figure S2. Percentage of weather data points determined from insolation and interpolation in the WARM data set.

Figure S3. Hourly PAR data from SURFRAD and WARM data sets throughout 2019.

Figure S4. Hourly temperature data from SURFRAD and WARM data sets throughout 2019.

Figure S5. Comparison between magic and Poincaré oscillators. Top: phase space trajectory for magic (left) and Poincaré (right) oscillators. Bottom: time evolution of the phase angle for magic (left) and Poincaré (right) oscillators. Integration was performed numerically using the Runge-Kutta Cash-Karp method from the Boost ODEINT library (“Boost Version 1.71.0,” 2019) using the following parameters: ω0 = π/12 rad/hour, fexternal and Z(t) = A cos(ωt), A = 1.0, ω = π/15 rad/hour, γ = 0.2, R0 = 1.5, and (a,b) = (1,0) at t = 0. A was set to 0 for t > toff. For clarity, early times of the simulation (where transients are significant) have been omitted from the plots.

Figure S6. Median, interquartile range, and range of photoperiod length values calculated at each day of the year based on the full (a) SURFRAD (n = 24 years) and (b) WARM (n = 29 years) data sets. Dashed line shows reference photoperiod from celestial mechanics calculations. Data in (b) is identical to Figure 5(b) in the main manuscript.

Figure S7. Interquartile range values calculated at each day of the year based on the full SURFRAD (n = 24 years) and WARM (n = 29 years) data sets.

Figure S8. Median, interquartile range, and range of photoperiod length values calculated at each day of the year based on the SURFRAD data set, excluding the year 2015. Dashed line shows reference photoperiod from celestial mechanics calculations.

Figure S9. Differences in flowering day predicted by the CROPGRO model using either an oscillator clock or celestial mechanics to determine the night length. Multiple sowing dates (120, 130, 140, 150, and 160) were tested for each year in the WARM (left) and SURFRAD (right) weather data sets. Average summer temperature was determined by the mean of hourly temperature values between days 120 and 210.

Figure S10. Total errors determined for different values of the GDD threshold GDDR1 (solid circles), along with a parabolic fit (dashed line) determined by least-squares regression.

ACKNOWLEDGEMENTS

The authors wish to thank Stephen P. Long for helpful discussion and feedback regarding this manuscript. Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the authors and do not necessarily reflect the views of the US Department of Agriculture. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the US Department of Agriculture. USDA is an equal opportunity provider and employer.

SOURCES OF FUNDING

This work was supported by the project Realizing Increased Photosynthetic Efficiency (RIPE) that is funded by the Bill & Melinda Gates Foundation, Foundation for Food and Agriculture Research (FFAR) and the UK Foreign, Commonwealth and Development Office (FCDO) under grant number OPP1172157.

DATA AVAILABILITY

The weather data and all R scripts used to produce the figures in this manuscript are available as a public GitHub repository: https://github.com/eloch216/oscillator-based-circadian-clock-analysis/tree/v1.0.0.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

CONTRIBUTIONS BY THE AUTHORS

E.B.L. contributed to the conceptualization, modeling, data analysis and writing (original draft). J.M.M contributed to conceptualization and revising the manuscript.

LITERATURE CITED