Integrating Oscillator-Based Circadian Clocks with Crop Growth Simulations

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.

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 behavior according to daily sunrise times and photoperiod lengths (see Section S1 of the Supplemental 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., 2017(Chew et al., , 2014Zardilis 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 behavior.
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 behaviors 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.

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) (Augustine et al., 2000;Illinois State Water Survey, 1995). WARM provides data at daily time steps and full yearly data sets are available from 1990-2019. SURFRAD provides data at 1 or 3 minute time steps and full yearly data sets are available from 1995-2019. Data were converted to hourly time steps using the procedures described in Section S2 of the Supplemental 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, 5 th -order Runge-Kutta integration with 4 th -order Cash-Karp error estimation was used with absolute and relative error tolerances set to 1x10 -6 for all simulations, and all calculations were performed using the default parameters in Tables 1 and 2. While integrating, weather data was treated as a continuous piecewise function defined by linear interpolation between the hourly values.

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 but can be coupled to an external signal with a different period . 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 twenty-four hour 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 Figure 1a, the dawn oscillator receives a signal at the beginning of the photoperiod which alters its phase angle (represented as an arrow rotating counterclockwise). When fully entrained, resets to zero at dawn. Likewise, the dusk oscillator resets its phase angle ( ) at the end of the photoperiod. Figure 1b shows the results from one 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, and are proportional to the time interval since sunrise and sunset, while the difference is proportional to the clock's estimate of the photoperiod length ( ). Figure 1c shows , 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 rad corresponds to 24 hr, e.g. .)

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 and represent the radius and phase angle of the oscillator in phase space, respectively. The general behavior of this oscillator can be understood directly from these equations. If exceeds , its rate of change is negative. Likewise, increases when . In other words, tends towards regardless of its initial state, and its rate of approach is determined by the friction coefficient (see section S4 of the Supplementary Information for a justification of this terminology). advances at a constant rate set by , 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 ⁄ , referred to as the oscillator's natural period.
A Poincaré oscillator can be coupled to an external driving signal ( ) by projecting and onto rectangular coordinates ( and ) and adding ( ) to one of the resulting derivatives, which results in a new set of equations: If ( ) is periodic with period and sufficiently strong, the resulting oscillator trajectory will be periodic with the same period rather than its natural period .
As an illustrative example, we integrate Equations 3 and 4 to determine the phase space trajectory ( Figure 2a) and phase evolution ( Figure 2b) of a Poincaré oscillator before, during, and after the application of a sinusoidal ( ), where the oscillator begins with the initial state , . Since to start, the radius gradually increases during the beginning of the simulation. Eventually this transient behavior fully decays and the oscillator reaches a stable orbit, moving counterclockwise around a circle of radius at a rate of rad/hr. 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/hr (set by the driving force). When the signal is removed, the oscillator reverts back to its original path. The friction coefficient ( ) has been chosen so that the trajectory stabilizes within one cycle following these changes in the driving signal.

Coupling Oscillators to Light
In the clock model, two driving signals ( and ) are derived from photosynthetic photon flux density (PPFD) values ( ) 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 ( ) which smoothly but sharply transitions between 0 and 1 as increases from 0 to according to a logistic formula where is a PPFD threshold. When , ( ) and when , ( ) .
Next, is used to determine the values of a day tracker and night tracker , quantities that evolve according to During the day, , so and reach equilibrium values of 1 and 0, respectively. Likewise, during the night, they equilibrate to 0 and 1. The parameter determines the length of the transition between equilibrium values, which occurs as an exponential approach; for example, if hr -1 , the transition is ~98% complete after two hours ( ). 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 .
Finally, and are calculated according to where determines the overall magnitude of the signals. is appreciably nonzero only at the start of the day, as that is the only time when both N and L are nonzero ( Figure 3a). In this sense, can be thought of as a molecular clock component activated by protein in the presence of light. Likewise, is only nonzero at the start of night (Figure 3b), analogous to a molecular clock component activated by protein in darkness.
Apart from a small spike around dusk, is only nonzero around dawn and consists of a roughly triangular peak lasting around two hours ( Figure 3a). Therefore, for a realistic input, is nearly periodic with an approximately twenty-four hour period. Likewise, ( Figure  3b) has a similar shape and is nonzero for about two hours near sunset. A critical feature of this scheme is that because both driving signals are only nonzero during transitions between light and dark, they both become zero under constant light ( ) or constant dark ( ) 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 (Section S8 of the Supplemental Information).

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 behavior 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 ( Figure 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 -1 •m -2 at hour by the Gaussian function ( ) ( ) ( ) ⁄ . By keeping the light profile width constant at hr and varying the initial phase of the dawn oscillator ( ) between 0 and , it is possible to determine an approximate value for , the time at which the transient response has fully decayed ( Figure 4a). Although the values calculated with different 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., = 30 days. For all subsequent times 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 ( ). By suddenly changing the width of the Gaussian light profile from hr to after , it is possible to estimate as the time when first lies within one minute of its final value (Figure 4b). For large jumps the entrainment process can take up to 20 days, but for the smallest jump (~30 minutes), full entrainment occurs within ten days.
For PPFD profiles derived from sunlight measurements, the length of the daily photoperiod changes slowly but continuously (by less than three minutes per day at 40° latitude), precluding a direct observation of entrainment in response to a single change as in Figure 4b. Instead, 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 Figure 4b), so errors in 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 ( ) was determined. The average daily difference across all years in each data set is shown in Figure 5a. For both the WARM and SURFRAD data sets, typically deviates from by less than 0.5 hours, 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 nonzero ) 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, 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 (Figure 5b). The interquartile range is always smaller than 0.8 hr and is frequently smaller than 0.5 hr, demonstrating that the clock is relatively insensitive to random variations in PAR data between years. An analogous plot based on the SURFRAD data exhibits similar interquartile range values (Supplemental Figure 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 Figure 5a.

APPLICATION EXAMPLE: GLYCINE MAX FLOWERING
The year-to-year repeatability of the clock model evident in Figure 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 ( ) 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 three days (Grimm et al., 1993;Piper et al., 1996). The original CROPGRO model uses celestial mechanics calculations to determine night length ( ) and includes five parameters that are allowed to vary between cultivars in order to fit the experimental data: minimum and optimal night length ( and ), minimum and optimal temperature ( and ), and threshold physiological age for flowering ( ).
To demonstrate the utility of the oscillator clock as an input to the CROPGRO model, flowering dates were predicted using either celestial mechanics ( ; see Section S3 of the Supplemental Information) or an oscillator clock ( ) 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 and used to calculate were varied using a brute-force approach on a finite grid to minimize the sum of ( ) across all sowing dates and years in the WARM data set (see Section S6 of the Supplemental Information).
Flowering time was also predicted using the SURFRAD data set without reparametrizing the models, giving an out-of-training data set prediction as a check for robustness. The results of this comparison are shown in Figure 6, where is plotted against for each combination of year and sowing date in the WARM (Figure 6a) and SURFRAD (Figure 6b) weather data sets. The dashed lines are guides to the eye representing perfect agreement. For the WARM data set, most values lie within two days of (the mean and standard deviation of the differences are 0.0 days and 0.9 days, respectively), and even for the out-of-training SURFRAD data set, differences are generally within three days (Supplemental Figure 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 and were required for this good agreement. These changes can be understood by noting that the increase in results from the small but systematic underestimation of photoperiod length in the oscillator clock model, while the larger range ( ) reduces the impact of random fluctuations in .

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 ( ) or the CROPGRO model ( ) using the WARM data set (Figure 7a). GDD 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 (see section S7 of the Supplemental Information). Here it is apparent that the GDD method often deviates from the more accurate CROPGRO model by more than three days. However, there is a more serious issue: the difference between and is strongly correlated with the average summer temperature, with the GDD model predicting earlier flowering dates during hotter years ( Figure  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 behavior as a time-dependent hypocotyl growth rate with a peak centered 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 reparametrize it to produce the same results when using an oscillator clock or celestial mechanics calculations, as demonstrated in Figure 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 and in the CROPGRO model to achieve the best agreement with measured 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 Figure 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 behaviors such as 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.   (Grimm et al., 1993). Parameters used with celestial mechanics calculations as the night length source are from the entry for the Williams cultivar (maturity group III, as is typical in central Illinois) in Table 3 of Grimm et al. (1993) with the exception of the value for , which is implied in Keisling (1982). Parameters used solely with the oscillator clock model were determined as described in the text.