A spatial model of the plant circadian clock reveals design principles for coordinated timing

We first implemented a reduced network model of the A. thaliana circadian clock, hereafter referred to as the De Caluwé model (De Caluwé et al, 2016). To decrease the complexity of the model, the authors grouped functionally similar genes into single entities (Fig 1A). The compact network model incorporates known light inputs to the network including “acute” activation of CIRCADIAN CLOCK ASSOCIATED 1 (CCA1)/LATE ELONGATED HYPOCOTYL (LHY) and PSEUDO‐RESPONSE REGULATOR 9 (PRR9)/PRR7 at dawn, a constant increase in the synthesis of CCA1/LHY and ELF4/LUX ARRHYTHMO (LUX) under light, and altered degradation of several components in the light or dark (Fig 1A). The model qualitatively recapitulates clock dynamics under both LD cycles and LL, and at only 9 equations and 34 parameters is also computationally tractable for spatial simulations. We modified the De Caluwé model to include a repression rather than activation interaction between CCA1/LHY and PRR9/PRR7 and included a term for CCA1/LHY repressing its own transcription, as these interactions have recently been shown experimentally (Adams et al, 2015). We also adjusted the degradation rate of PRR5/TOC1 to be higher in the dark, as several studies have established that the degradation of both proteins is increased in the dark (Más et al, 2003; Kiba et al, 2007; Kim et al, 2007; Fujiwara et al, 2008) (Appendix Table S1, Materials and Methods). The identified parameter values were robust to at least 5% variation (Appendix Fig S1, Materials and Methods).

We implemented the modified De Caluwé model on a simplified template of a seedling. The template consisted of approximately 800 cells, classified into cotyledon, hypocotyl, root, and root tip regions (Fig 1B, Materials and Methods). Although a number of studies have demonstrated local cell‐to‐cell coupling between clocks in A. thaliana (Fukuda et al, 2007, 2012; Wenden et al, 2012; Endo et al, 2014; Takahashi et al, 2015; Gould et al, 2018; Greenwood et al, 2019), the identity of the coupling signal, or signals, is unclear. Initially, to model the signal, the level of CCA1/LHY mRNA in one cell (Fig 1C, black squares) was assumed to be coupled to the average level of the cell’s neighbors (Fig 1C, white squares). The coupling strength, J_local, determines the extent that molecules are shared (Materials and Methods). Although we did not model diffusion directly, this coupling function was designed to simulate the passive movement of molecules which commonly occurs through plasmodesmata in plants (Faulkner, 2018). To simulate period variation between cells, at each time step, we multiplied the level of the mRNA and protein by a “scaling” parameter. For each cell, this parameter was randomly selected from a normal distribution to give a unique value for each cell through the simulation (Materials and Methods). This approach generates between‐cell but not within‐cell period differences, allowing us to focus on the one source of variation. Further, we could set the range of the distribution differently for each organ, as informed by experimental data (Appendix Fig S2).

We next attempted to recapitulate in our model the differences in clock period and phase in different organs that have been observed in experiments using a single‐cell CCA1‐YFP reporter (Gould et al, 2018) and a GIGANTEA luciferase reporter (Greenwood et al, 2019). In these experiments, faster rhythms were observed in the cotyledon, hypocotyl, and root tip, with slower rhythms in the rest of the root. We first reanalyzed existing luciferase data (Greenwood et al, 2019) and confirmed that these relationships held for several of the core clock genes in our model, PSEUDO‐RESPONSE REGULATOR 9 (PRR9), TIMING OF CAB EXPRESSION 1 (TOC1), and ELF4 (Fig 2A). Whereas in our previous “phase‐only” model, we fixed the periods to be different in each part of the plant, with our spatial network model, we could now investigate what causes the differences in periods. Previously it has been hypothesized that varied sensitivities to light alter periods across the plant (Bordage et al, 2016; Greenwood et al, 2019; Nimmo et al, 2020). This was based on previous experiments demonstrating: (i) period differences between organs when exposed to equal amounts of light (Bordage et al, 2016; Greenwood et al, 2019); (ii) the loss of some period differences in the light‐sensing mutant phyb‐9 (Greenwood et al, 2019; Nimmo et al, 2020); and (iii) the spatial expression pattern of PHYB and other light‐sensing genes (Somers & Quail, 1995; Bognár et al, 1999; Tóth et al, 2001). Experiments, however, are confounded by changes in metabolism and development caused by light (Nozue et al, 2011). We therefore tested the light‐sensitivity hypothesis in our spatial model, by setting the sensitivity to light to differ depending on the region, and examining whether this could generate the period differences observed across the plant.

We entrained the cells in our simulations to LD cycles for 4 days before releasing them into LL for a further 6 days and measured the periods, as carried out in previous experiments (Greenwood et al, 2019). We initially set the local cell‐to‐cell coupling parameter, J_local, to 2. When assuming high sensitivity to light in the cotyledon (L_sens = 1.6) and hypocotyl (L_sens = 1.0), but lower in the root (L_sens = 0.65) and root tip (L_sens = 0.95), all regions entrained to the LD cycles (Appendix Fig S3). Upon transfer to LL, we were able to generate different periods (Fig 2B and C) and phases (Fig 2D–F and Appendix Fig S4) in each organ, matching those observed experimentally. This is due to higher light sensitivity causing the clock to run faster in our simulations (Appendix Fig S5), as expected for a diurnal organism (Aschoff & Pohl, 1978). Finally, by setting the sensitivity to light in our locally coupled model to be zero in all regions of the seedling, we were also able to simulate the loss of period differences observed in the light‐sensing mutant phyb‐9 (Greenwood et al, 2019) (Fig EV1A and B). Thus, our results revealed that different sensitivities to environmental inputs are sufficient to generate the experimentally observed spatial differences in period and phase across the plant.

Previously we observed two waves of clock gene expression, one traveling up, and one down the root, in CCA1, PRR9, and GI reporters (Gould et al, 2018; Greenwood et al, 2019). These waves could be explained in a “phase‐only” model by local cell‐to‐cell coupling (Fukuda et al, 2007, 2012; Gould et al, 2018; Greenwood et al, 2019). We next tested whether a plausible mechanism for cell‐to‐cell coupling, sharing of mRNA between cells (Maizel et al, 2020), can recapitulate the experimental observations. To provide a benchmark, we first analyzed the seedlings carrying transcriptional reporters for PRR9, TOC1, and ELF4 (Greenwood et al, 2019) at the sub‐tissue level (Materials and Methods). As in previous studies (Fukuda et al, 2007, 2012; Wenden et al, 2012; Gould et al, 2018; Greenwood et al, 2019), space‐time plots revealed spatial waves of gene expression within and between organs (Fig 3A, Appendix Fig S6A and C and Movie EV1). The direction of the waves could be clearly observed in plots of the final peaks of expression (Fig 3B and Appendix Fig S6E and G). For each gene, the wave directions appeared similar, traveling from the faster oscillating regions to the slower oscillating regions. Particularly distinct were the waves converging from the hypocotyl and root tip to the slower oscillating root region.

We next analyzed at the sub‐tissue level the simulations performed with varied light sensitivity and local coupling (Appendix Fig S7B, Materials and Methods), to see if they capture the spatial dynamics. For each gene, the wave patterns appeared similar to experiments, traveling from the fast oscillating regions that are more sensitive to light, into the slower regions with lower sensitivity to light (Fig 3C and D, Appendix Fig S6B, D, F and H, and Movie EV2). These waves required cell‐to‐cell coupling through the local sharing of clock components, as we only observed waves with coupling strengths, J_local, above approximately 1 (Fig 3E and Appendix Fig S8). Similar simulation results were found when keeping regional differences of light sensitivity and local cell‐to‐cell coupling, but assuming no cell‐to‐cell variability within regions (Appendix Fig S9). We also found that our simulation results were qualitatively similar when using different clock gene mRNA (Fig EV2A) or protein (Fig EV2B) as the coupling component, suggesting that any cell‐to‐cell sharing of clock components could explain the experimentally observed spatial dynamics. Additionally, we ran simulations assuming local coupling between 8 rather than 4 neighbor cells, and simulations assuming global (all‐to‐all) coupling. Increasing the local coupling to be between 8 neighbor cells gave similar results (Fig EV3A and B). However, with global coupling all cells adopted the same phase at higher coupling strengths, regardless of position in the plant (Fig EV3C). Finally, we further analyzed our simulations in which we set the sensitivity to light input to be equal in all regions of the seedling. At the sub‐tissue level, we observed the loss of spatial waves (Fig EV1D and F) seen in the light‐sensing mutant phyb‐9 (Greenwood et al, 2019) (Fig EV1C and E). Taken together, these results show that the assumptions of local cell‐to‐cell coupling and differential light sensitivity between regions are the key aspects of our model that allow a match to experimental data.

In addition to local coupling, the long‐distance shoot‐to‐root movement of ELF4 has been proposed to couple clocks in different organs (Chen et al, 2020). We simulated seedlings assuming that ELF4 protein expressed from cells in the shoot is shared with cells in the root tip region. We did this to approximate a known destination for phloem signals in plants (Oparka et al, 1994), and the location that ELF4 is observed in experiments (Chen et al, 2020). In these simulations, the peak times of cells in the root tip became closer to those in the shoot, but we did not see the spatial waves of gene expression observed in experiments (Fig 4A). We then simulated seedlings with combinations of local and long‐distance coupling, to see if together they improved the match to experimental data. In these simulations, we saw a clear spatial structure with waves of clock gene expression (Fig 4B), matching previous simulations with only local coupling (Fig 3). These simulations confirm that it is the local, not long‐distance, coupling that drives the spatial waves of gene expression.

We hypothesized that long‐distance coupling instead contributes by matching clock periods between different organs of the plant, as has been proposed for ELF4 (Chen et al, 2020) and light piping from the shoot to the root (Nimmo, 2018). To test this in our model, we removed the assumption of high light sensitivity at the root tip (Materials and Methods), and re‐simulated the model with combinations of local and long‐distance coupling. Without long‐distance coupling (J_long = 0), as expected, cells were slower at the root tip, and the spatial wave up the root was lost (Fig 4C, black). However, with J_long above approximately 0.1, we observed an increase in the speed of cells at the root tip, and the recovery of spatial waves up the root (Fig 4C, yellow). This approximated the typical spatial structure observed with higher light sensitivity at the root tip (Fig 4B). Together, this suggests that although long‐distance coupling does not generate spatial waves, it can create period differences between regions of the plant. Further experiments will be required to understand the extent of long‐distance clock coupling during plant development, so henceforth we focus on our model assuming different sensitivities to light and only local coupling.

We simulated our locally coupled model, with the sharing of CCA1/LHY mRNA, under LD cycles to investigate whether local flexibility persisted with rhythmic input. We first simulated idealized LD cycles, where light is fully on during the daytime and off at night (Fig 5A, left; Materials and Methods). The sensitivity of cells to light during the daytime varied between regions, as with LL simulations, and we included local cell‐to‐cell coupling. Under LD cycles, cells had more similar timings than under LL; however, inspection of the peaks of PRR9/PRR7 expression revealed differences between regions (Fig 5B, black dots). These phase differences were qualitatively the same as observed under LL (Fig 3) and persisted even at high coupling strengths (Fig 5B, black dots). Also as observed under LL, with J_local above approximately 1, spatial waves could be observed (Fig 5B, black dots). These simulations qualitatively match the experimentally observed spatial waves under LD cycles (Greenwood et al, 2019). A similar structure was observed in the other clock genes (Appendix Fig S10, black dots), although spatial waves of PRR5/TOC1 and ELF4/LUX were altered because their peaks coincided with the LD transition which had a synchronizing effect (Appendix Fig S10C and D, black dots). Simulations also captured whole‐plant experimental data under short or long days (Fig EV4) and predicted that a spatial structure remains under these conditions (Appendix Fig S11).

Our modeling suggested that different sensitivities to environmental inputs across the plant allow the clock to be locally flexible, by allowing it to adopt different phases across the plant even under LD cycles. However, it was not clear whether a locally flexible clock would still entrain to environmental cycles and maintain a stable spatial structure of phase differences when the environment is noisy, as is the case outside of the laboratory. To test this, we simulated more realistic LD cycles containing fluctuations in the light intensity (Fig 5A, right). We first considered the scenario where the fluctuations in light intensity are uncorrelated for neighboring cells, to simulate microenvironmental differences (Materials and Methods). Without cell‐to‐cell coupling (J_local = 0) the peak times of PRR9/PRR7 expression were highly variable within a region and any spatial structure was difficult to resolve (Fig 5B, red dots). However, with increasing strengths of coupling we observed a decrease in variability, revealing a spatial structure (Fig 5B, red dots) similar to that observed under idealized LD cycles (Fig 5B, black dots). A similar effect was observed for the remaining clock genes (Appendix Fig S10, red dots). Overall, local cell‐to‐cell coupling appears to reduce the within‐region variability of cellular rhythms. Yet, despite the stabilizing effect of the coupling, between‐region differences persist, due to the between‐region differences in sensitivity to the environment.

This emergence of spatial structure suggests that local coupling improves the timing of individual cells under noisy environments. To quantify this effect, we measured the timing of the peaks (Fig 5C and Appendix Fig S12) and troughs (Fig 5D and Appendix Fig S12) of expression relative to simulations without noise in the LD cycle (Materials and Methods). The timing of each gene was mis‐timed under noisy LD cycles compared to idealized LD cycles. For example, with J_local = 0, there was a timing error of 0.48 ± 0.02 h (mean ± SD) for the peaks, and 0.46 ± 0.02 h for the troughs of PRR9/PRR7 expression. Each gene showed a differing degree of error in their peak (Fig 5C and Appendix Fig S12) and trough (Fig 5D and Appendix Fig S12) times, which reflects differences in their regulation under light and dark conditions (Appendix Fig S12 and Appendix Table S1). For example, the troughs of ELF4/LUX expression were synchronized (Fig 5D, blue line) due to rapid degradation in the dark (Appendix Fig S12D). However, for each gene, we observed a decrease in the timing error with increasing strengths of local cell‐to‐cell coupling (Fig 5C and D). This is consistent with local coupling having a stabilizing effect on cellular oscillations, increasing their robustness to perturbation from the noisy environment. Local cell‐to‐cell coupling similarly reduced the timing error when we assumed weak correlations in the LD cycles that cells are exposed to, to simulate shared microenvironments (Appendix Fig S13A–C, Materials and Methods). However, local coupling no longer reduced the timing error when we assumed stronger correlations in the LD cycles, as neighboring cells experienced more similar LD cycles (Appendix Fig S13D and E). Taken together, our results show that by averaging fluctuations between neighboring cells, local cell‐to‐cell coupling can improve timing accuracy of individual cells under noisy environments. Because coupling acts locally, between‐region differences in phase, due to between‐region differences in light sensitivities, are able to persist despite the averaging effect of cell‐to‐cell coupling.

As a model for the plant circadian clock, we utilized the compact model introduced by De Caluwé et al (2016). The original compact model consists of 9 ordinary differential equations. Among them, 8 equations describe the temporal evolution of the mRNA and protein levels of the core clock genes. The clock genes are grouped into four sets of pairs labeled as: CL (CCA1 and LHY), P97 (PRR9 and PRR7), P51 (PRR5 and TOC1), and EL (ELF4 and LUX). The 9‐th equation is for protein P, a yet‐identified protein which mediates the acute activation of transcription by light observed in experiments (Ito et al, 2005; Locke et al, 2005b).

In our modeling, we modified the original De Caluwé model. According to the experiments reported in Para et al (2007), the compact model assumed that the P97 variable is activated by CL. However, more recent work has shown that LHY acts as a repressor of PRR9 and PRR7 (Fogelmark & Troein, 2014; Adams et al, 2015). We therefore replaced the activation term that represents the connection from CL to P97 with a repression. In addition, we added a term for CL repressing its own transcription, as this has also recently been observed for CCA1 and LHY (Adams et al, 2015).

We also made modifications to the light input mechanism. In the original De Caluwé model, light signals feed into the network at multiple levels: (i) The transcription of CL and P97 is acutely activated at dawn by the dark‐accumulating protein P; (ii) The transcription of EL is increased continuously by light; (iii) The translation of CL is increased by light; (iv) The degradation rate of CL mRNA is increased by light; (v) The degradation rate of P51 protein is increased in the light; and (vi) The degradation rates of P97 and EL proteins are increased in the dark. In our model, we adjusted the degradation rate of P51, as several studies have established that the degradation of both PRR5 and TOC1 are also increased in the dark (Más et al, 2003; Kiba et al, 2007; Kim et al, 2007; Fujiwara et al, 2008) (Appendix Table S1).

Additionally, we introduced L_sens, to represent the cellular sensitivity to light. In the presence of light L = 1, a cell senses light according to its sensitivity, L_sens L > 0. In the absence of light L = 0, the cell does not sense light, L_sens L = 0. We include L_sens in the terms describing protein P degradation, EL and CL translation, and CL mRNA degradation, as in these terms light increases the rates of reaction. We additionally include L_sens in the term describing the degradation of the EC, as it has previously been suggested that ELF3 degradation (an EC component not explicitly modeled in De Caluwé) is increased in the light (Pokhilko et al, 2012). L_sens was set to 1 in the single cell implementation of the model, but varied between cells in the multicellular implementation (see “Spatial molecular model”). The revised single‐cell model is now described by the following differential equations:

cαm and cαp represent the concentrations of the α‐th mRNA and protein (or protein complex) respectively, for α = CL, P97, P51, and EL. P represents the activity of the dark accumulated protein P. The parameters v_j, k_j, p_j, and d_j represent the reaction rates for transcription, mRNA degradation, translation, and protein degradation respectively, for the various reactions j = 1… (Appendix Table S1). K_j represent a scaling factor for the protein concentrations. D represents a dark signal, which is equal to 1 when L = 0 and equal to 0 when L > 0. This explicit representation of a dark signal is commonly used in circadian clock models to capture the different reaction rates observed under dark conditions in experiments (Pokhilko et al, 2012; Fogelmark & Troein, 2014; De Caluwé et al, 2016).

The original De Caluwé model contains 34 parameters, the values of which have been obtained previously through automated optimization (De Caluwé et al, 2016). Firstly, we adjusted the value of d_3L to be less than d_3D. We did this to account for the fact that several studies have demonstrated that the degradation of both PRR5 and TOC1 are increased in the dark (Más et al, 2003; Kiba et al, 2007; Kim et al, 2007; Fujiwara et al, 2008). We then used automated optimization to find the values related to the modified interaction between CL and P97 (K0,K4,K5, and K5b). We optimized by minimizing a cost function, combining Sobol search and simulated annealing to optimize within the range Ki∈[0.1,10] for i=0,4,5,5b.

We computed the cost function as follows. First, we simulated the model under 12‐h light–12‐h dark conditions for a total of 60 days and then released to LL for 60 days, followed by constant darkness (DD) for 60 days. In each light condition, we discarded the first 55 days as a transient dynamic. To ensure detectable rhythmicity under LL and DD conditions, we required that all variables must have a minimum value of 0.1, as well as a minimum difference of 10% between their minimum and maximum values. We penalized any solution that did not meet these criteria with an arbitrarily large score. Then, we calculated the free‐running period from CL gene expression using the chi‐square periodogram (Sokolove & Bushell, 1978) at a significance level of 1%. We gave a score of 0 to a solution having a free‐running period between 24 and 25 h under LL and between 25 and 28 h under DD. For solutions with free‐running periods outside these ranges, we allocated the scores of 2 (τ_LL − 24.5)²/(0.1⋅24.5)² and 2 (τ_DD − 26.5)²/(0.1⋅26.5)², where τ_LL and τ_DD represent free‐running periods under LL and DD respectively. Under 12‐h light–12‐h dark cycles, we penalized solutions that did not entrain to the LD cycles with an arbitrarily large score. We assigned entrained solutions a score of 0 for each gene that attained peak expression within ± 1 h of the expected ZT, which were as follows: CL = 1.5, P97 = 6, P51 = 12, and EL = 9. We scored expression peaks lying outside these intervals as (ZT_CL − 1.5)²/(0.1⋅24)², (ZT_P97 − 6)²/(0.1⋅24)², (ZT_P51 − 12)²/(0.1⋅24)², (ZT_EL − 9)²/(0.1⋅24)², where ZT represents the peak time measured from the corresponding genes expression. Finally, we penalized solutions that showed a large phase shift, Δ, when transferred from LD cycles to LL. The phase shift was scored as Δ²/(0.1⋅24)². The optimized parameter values together with the original values determined by De Caluwé et al (2016) are listed in Appendix Table S1. The free‐running period of the revised model was 25.5 h under LL (L = 1; D = 0) and 27.4 h under constant darkness (L = 0; D = 1). Under LD cycles, the peak times (ZT) of the gene expressions were CL = −0.3 h, P97 = 6.0 h, P51 = 11.8 h, and EL = 9.7 h.

To test the robustness of the single cell model we carried out a parameter sensitivity analysis. We distributed all 34 parameters uniformly in the ± 5% range of their optimized value. For 50,000 sets of randomized parameter values, the free running period under LL (25.5 ± 0.4 h), free running period under DD (27.4 ± 0.5 h), amplitude under LL (0.9 ± 0.1), and amplitude under DD (1.4 ± 0.2) of CL expression were calculated. The values were distributed in a narrow range centered at their target values, indicating that they are not sensitive to the selected values of the 34 parameters (Appendix Fig). Our sensitivity analysis could in future be extended (e.g., by using Approximate Bayesian Computation) to calculate confidence intervals for the parameter estimates. S1

We analyzed experimental luciferase and confocal data to characterize the periods within a seedling. We previously completed an organ‐level analysis of period and phase for PRR9::LUC, ELF4::LUC, and TOC1::LUC expression (Greenwood et al, 2019). In this analysis, 315 µm diameter regions of interest (ROI) were defined to represent the cotyledon, hypocotyl, root, and root tip regions of the seedling (Greenwood et al, 2022). These ROI were used to generate time series, from which period estimates were made across a number of individual seedlings. Here, we pooled estimates made from individual seedlings carrying PRR9::LUC, ELF4::LUC, or TOC1::LUC reporter genes, and calculated the mean for each region (Fig 2A). Only periods from time series classified as rhythmic (as defined previously (Greenwood et al, 2019)) were used in the calculation.

To characterize the variability of the periods, we measured the between‐cell variation of clock periods in each region of the plant (Appendix Fig S2). To do this we analyzed a dataset of single cell time‐lapse microscopy, CCA1::CCA1‐YFP fusion protein reporter data (Gould et al, 2018). We defined cells as being within the root tip region if they are less than 315 µm from the most distal cell of the root tip. Where there are multiple imaging sections in the same organ, cells are pooled to give one region for the cotyledon, hypocotyl, and root. Variability was estimated using the coefficient of variation (CV; SD/mean) of all cells within these regions. Only periods from time series classified as rhythmic (as defined previously (Gould et al, 2018)) were used in the calculation.

We implemented our optimized molecular model on a simplified template of a seedling. The template consisted of approximately 800 cells, classified into cotyledon, hypocotyl, root, and root tip regions (Fig 1B). Each cell contained an implementation of the revised compact model. To simulate growth of the seedlings, we added a row of cells to the root tip every 24 h. During this growth, we kept the root tip region of the template fixed in size. To do this, after the addition of new cells, the previously uppermost root tip cells became root cells instead.

We set the light sensitivity of the cell, L_sens, to vary depending on the position within the template. In doing so, we found that modest differences in L_sens were sufficient to generate period differences between regions. To approximate the period differences that we observed between regions in experiments (Fig 2A) we set L_sens = 1.6 for cotyledon cells, L_sens = 1.0 for hypocotyl cells, L_sens = 0.65 for root cells, and L_sens = 0.95 for root tip cells. As the template grows, the uppermost root tip cells become root cells, so the light sensitivity of these cells decreases so that they have a sensitivity characteristic of root cells, L_sens = 0.65. This caused the clock period of these cells to slow.

As the coupling agents have yet to be clearly identified by experimental studies, we initially assumed individual cells to be coupled through CL. Our model for coupled plant cells is described as follows:

As in the revised single‐cell model, cαm and cαp represent the concentration of the α‐th mRNA and protein (or protein complex) respectively for α = CL, P97, P51, and EL, in the i‐th cell (i=1,2,…,N). P again represents activity of the dark accumulated protein P. To generate cell‐to‐cell variability in periods, we multiply the derivative with respect to the time of each molecular species by the time scaling parameter τ. τ is a real number randomly selected from the normal distribution N(1, 0.059), N(1, 0.028), N(1, 0.073), N(1, 0.089) for cells from the cotyledon, hypocotyl, root, or root tip respectively, as informed by the analysis of single‐cell experimental data (see “Characterizing clock periods from experiments”).

In the first equation, the CL gene is locally coupled to its neighboring cells, where J_local is the coupling strength and c¯CL,im=14∑j∈NicCL,jm represents averaged expression level over 4 neighboring cells (left, right, above, and below). In the case of local coupling between 8 neighboring cells (left, right, above, below, left above, left below, right above, and right below) or global coupling, the averaged expression becomes c¯CL,im=18∑j∈NicCL,jm or c¯CL,im=1N∑j=1NcCL,j(m) respectively. In Fig EV2 we varied the identity of the local coupling component. Here, the averaged expression becomes c¯α,im=14∑j∈Nicα,jm or c¯α,ip=14∑j∈Nicα,jp for coupling through mRNA or protein respectively, for α = CL, P97, P51, EL. The coupling term is appended to the equation of the corresponding molecular species. We set the local coupling strength as J_local = 0, 0.01, 0.1, 1, 2, or 4.

In Fig 4, we simulated long‐distance coupling between the hypocotyl and root tip cells through the movement of ELF4. To implement this, a long‐distance coupling term is appended to the equations describing EL protein in root tip cells, where J_long is the long‐distance coupling strength and c¯EL,longp=1NR∑j∈RcEL,jp represents the averaged expression level over cells in the hypocotyl, R. We set the long‐distance coupling strength as J_long = 0, 0.01, 0.1, 1, 2, or 4.

To assess our model, we quantified the times of the peaks and peak‐to‐peak times of simulated gene expression. We firstly did this from ROI averages to allow for comparison to experimental data (e.g., Fig 2) (Greenwood et al, 2019). We defined 5‐by‐3 (width‐by‐height) pixel ROI on the seedling template, approximating the positions used in the experimental analysis (Appendix Fig S7A). We took the mean of these ROI at each time point to give the simulated time series. We identified the peaks of the time series within 24–144 h after transfer to LL using the “findpeaks” MATLAB (MathWorks, UK) function. We imposed that peaks must be more than 19 h apart and greater than the mean of the time series. Next, we measured the period as the mean of the times between consecutive peaks. In Figs 2 and 3, we plotted the analysis of the simulated data with previously collected experimental data (Greenwood et al, 2019) for comparison. For visual clarity, only peaks in which all organs completed the full cycle within the time window were plotted. Additionally, only peaks from the experimental time series classified as rhythmic (as defined previously (Greenwood et al, 2019)) were plotted.

To visualize the spatial clock dynamics, we created space‐time intensity plots of experimental and simulated gene expression (e.g., Fig 3A and C). We made plots of PRR9::LUC, ELF4::LUC, and TOC1::LUC expression as described previously (Greenwood et al, 2019). In brief, we average expression across sections that are perpendicular to the primary axis of the seedling. For simulations, we take the mean expression from 1‐by‐5 cell cross‐sections through the primary axis of the template (Appendix Fig S7B, yellow lines). We additionally plotted the times of the final peaks of expression of the space‐time plots (e.g., Fig 3B and D). We detected peaks as before (see “Analysis of peaks and periods of expression”), but restricted the detection to within a 24 h window of the expected time of the final peak of expression. If more than one peak was detected we plotted the peak with the greatest height. However, under some conditions, we consistently observed two peaks of TOC1 expression within 24 h. We therefore plotted the earliest peak of TOC1 expression.

We simulated light‐dark cycles with a day length of 8, 12, or 16 h. We termed these cycles without noise “idealized LD” (Fig 5A, left). To simulate fluctuations in the light cycle we utilized an algorithm developed previously (Pittayakanchit et al, 2018), which approximates meteorological data (Gu et al, 2001). Fluctuations lasted an interval of time drawn randomly from an exponential distribution of mean 2.4 h. The intensity of light during each fluctuation was drawn randomly from a second function of uniform distribution in the range 0.5 to 1.5. We termed this condition “noisy LD” (Fig 5A, right). LD cycles were input to the molecular model through the parameter L.

Individual plant cells may not experience a unique noisy LD cycle but instead have some correlation in the light cycle between one another. To examine the dependence of our results on the correlation, we extended our algorithm to consider two types of fluctuations, global and individual. The global fluctuation, ξ_global, is applied commonly to all cells, while the individual fluctuation, ξ_individual, is applied only to individual cells. The global and individual fluctuations are drawn randomly from the same distributions as before, but are statistically independent from each other. By combining the two fluctuations, each cell receives a light input as Lt=1‐{λξglobal(t)+(1‐λ)ξindividual(t)} with a mixture ratio of 0≤λ≤1. If λ=1, the same fluctuation is applied to all the cells. If on the other hand λ=0, fluctuations applied to individual cells are all independent from each other. We simulate noisy LD cycles with λ = 0, 0.1, 0.25, 0.5, and 1.

We introduced the cell timing error as a relative measure of clock precision under noisy LD cycles. For all cells of a simulation (h = 1, 2 ,..., N), we computed the times of the final peaks of expression T. We computed the peaks under the noisy LD condition T_noisy, and the idealized LD condition, T_idealized, with a matching strength of local coupling. The latter condition acted as a benchmark for the cell’s ideal time of expression, in order to estimate the error caused by the noise. We then calculate the timing error E at each coupling strength as

In the experimental analysis, we excluded seedlings from the analysis if they were not rhythmic, as described in the previous experimental studies (Mockler et al, 2007; Gould et al, 2018; Greenwood et al, 2019). We plotted the mean and SD when the results of the analysis were normally distributed, and otherwise plotted the median, 25‐th, and 75‐th percentile. We chose representative seedlings for the space‐time intensity plots. For simulations, we performed 9 simulations under each condition. We performed each simulation with a different distribution of the time scaling parameter, τ, (see “Spatial molecular model”). We matched these distributions between different conditions to allow for comparison. Summary statistics and representative space‐time intensity plots were presented, as for experiments. We did not use statistical tests to select the sample size or perform blinding during the analysis.