Theory of time delayed genetic oscillations with external noisy regulation

We present a general theory of noisy genetic oscillators with externally regulated production rate and multiplicative noise. The observables that characterize the genetic oscillator are discussed, and it is shown how their statistics depend on the externally regulated production rate. We show that these observables have generic features that are observed in two different experimental systems: the expression of the circadian clock genes in fibroblasts, and in the transient and oscillatory dynamics of the segmentation clock genes observed in cells disassociated from zebrafish embryos. Our work shows that genetic oscillations with diverse biological contexts can be understood in a common framework based on a delayed negative feedback system, and regulator dynamics. Video Abstract: Theory of time delayed genetic oscillations with external noisy regulation


Introduction
Embryonic cells are highly dynamic systems with time dependent biochemical states [1]. These phenotypic states are characterized by dynamical gene expression profiles. A simple but important example of dynamical gene expression is when the concentration of a protein inside a cell evolves in a cyclic and stochastic manner, this is denoted as a genetic oscillation [2]. Such dynamic process acts in many instances as a clock that determines the timing between different physiological processes [3].
Genetic oscillations can emerge from the dynamic regulation between different components of a genetic network [4]. Several theoretical studies [5][6][7][8][9] suggest that the emergence of these oscillations is mainly due to a core element that exerts time-delayed negative feedback. Such delay appears naturally inside a cell since the time it takes to produce a protein through transcription and translation is non-negligible [5-7, 10, 11]. The feasibility of this mechanism has been shown in experiments with synthetic systems [10,11].
In this work we study a generic noisy time delayed negative feedback system W(t) with externally regulated production rate P(t) and multiplicative noise (figure 1). In this model W(t) correspond to protein concentration inside the cell. We develop a framework to understand the statistics that characterizes the time evolution of W(t). We are motivated to understand the similarities and differences in two different experimental cases, the single cell oscillations in the circadian clock and in the segmentation clock. Both systems switch stochastically between production and degradation phases with a characteristic time scale. Yet for circadian rhythms the process is sustained (figure 2(a)) and for the segmentation clock the amplitude increases and decreases transiently as the cells differentiate and cease oscillating (figure 2(e)). We show that the statistics characterizing W(t) have generic features found in both cases, and the differences reflects the type of regulation that P(t) exerts. Figure 1. Schematic of the noisy genetic oscillator considered in this work. W(t) is the protein level of a noisy oscillating gene that is regulated by P(t) which can also be noisy. The model is defined by equation (1): P(t) is the production rate of W(t), while a time-delayed negative feedback given by W(t − t d ) represses its own production.  (1) and (5) are simulated in a regime where W(t) exhibits oscillations with an amplitude that increases and decreases transiently, where P(t) evolves initially as a noisy ramp that decays back to its initial level once it reaches the threshold φ.

Model of genetic oscillations with external regulation
In our model the time evolution of W(t) is given by The parameter T 0 corresponds to the relaxation time, P(t) corresponds to the time dependent production rate and K 0 determines the mid-point to saturation of the Hill function H − h (X), which promotes repression via where h is the Hill coefficient. In a switch-like regime of production, h = ∞ and the Hill function becomes where θ(X) corresponds to the Heaviside function. The stochastic nature of gene expression in equation (1) is given by multiplicative noise, where we have chosen that the noise term Ξ(t) to be Gaussian and correlated in time as and with Ξ(t) = 0. Our choice of multiplicative noise ensures that the values of W(t) are always positive, and that the oscillation peaks are noisier than the troughs. This corresponds to a basic feature observed in the experimental data. When the production rate is constant (P(t) = P 0 ) and in the absence of fluctuations (ζ 2 = 0), W(t) becomes oscillatory when it undergoes a Hopf bifurcation as P 0 and t d cross specific thresholds (for details see [12]). The introduction of fluctuations (ζ 2 > 0) makes the component W(t) switch stochastically between production and degradation phases. This switching occurs with a characteristic time scale t d [13] (figures 2(a) and (d)), and can be interpreted as the switching on and off of transcription in a given gene [5]. When the production rate varies with time, it influences the dynamics of W(t) in a characteristic manner that we investigate below.
We find that two experimental cases are captured by our model with the parameters shown in table 1. Case I is the circadian clock, for which we analyzed the spatial average in luminescence intensity from single fibroblasts expressing PER2::LUC [14] (figure 2(a)). PER2 is a component of the mammalian circadian clock which indirectly represses its own production [16]. In this case the genetic oscillator W(t) corresponds to the concentration of PER2 in a single cell (figure 2(b)), which shows self-sustained oscillations with an amplitude that fluctuates around a constant mean value. In this case the production rate is constant (P(t) = P 0 ), and the amplitude fluctuations are given by the multiplicative noise (figures 2(b) and (c), parameters in table 1).
Case II is the segmentation clock, which determines the sequential segmentation of the pre-somitic mesoderm into body segments in vertebrate embryos [17]. We analyzed the spatial average in fluorescence intensity from single cells expressing Her1-YFP, that were previously dissociated from the presomitic mesoderm of zebrafish embryos (figure 2(d), see appendix G for experimental procedures). The protein Her1 is one of the components of the segmentation clock and also represses indirectly its own production. For the segmentation clock W(t) corresponds to the mean concentration of Her1 in a single cell, showing oscillations with an amplitude that increases and decreases transiently once per cell as the cells differentiate in these culture conditions. In this case the oscillatory regime has specific initiation and termination times (figure 2(e)). The simplest way to reproduce these features is by having a time-dependent production that initially increases as a noisy ramp, and then decays after reaching a threshold φ. For case II the time dependent production rate is given by where t o is the time where the production rate initiates and t φ is the time where P(t) reaches the threshold value φ (figure 2(f)). Fluctuations in the production rate are given by the noise term B P (t) which is Brownian with time correlation B P (t) B P (t ) = ζ 2 P min(t, t ) where min(t i , t j ) = t j when t i > t j > 0. In case II, the number of observed Her1 cycles, the amplitude envelope, and the duration of the oscillatory regime varies between cells. Our model reproduces such variability through the fluctuations in the trajectory of P(t), where as a result, the time to reach the threshold φ varies between different realisations (see appendix A).

Statistical analysis of the model
To understand the key statistics of this model, it is convenient to analyze equation (1) in the switch-like case and rewrite it in terms of normalized concentration w = W/K 0 and time τ = t/T 0 where χ(τ ) = P(τ )T 0 /K 0 , τ d = t d /T 0 is the normalized time delay and B(τ ) is the time integral of T 0 Ξ(τ ), which is Brownian noise with B(τ )B(τ ) = σ 2 0 min(τ , τ ) and σ 2 0 = ζ 2 0 T 2 0 . In our analysis we use the Stratonovich convention for stochastic integrals, which is used when the noise term T 0 Ξ(τ ) is interpreted as colored noise with a correlation time in the limit τ c → 0 [18].
First we study the statistics of equation (7) when χ is constant and σ 2 0 is low. Oscillatory dynamics in w(τ ) is observed when χ 1; in this regime the nonlinear term χH − ∞ (w(τ − τ d )) switches between two states: 0 and χ (figure 3(a)). We calculate the statistics of the values of w(τ ) at the transition points between these two states. The value at the transition 0 → χ is defined asŵ − i , and the value at the transition χ → 0 is defined asŵ + i (figure 3(a)), where the index i denominates the cycle number. In this case the values ofŵ − i are independent realizations of the stochastic variablê which conditional probability distribution is given by a log-normal distribution with where its mean and variance are given by (for derivation see appendix B). At the same timeŵ + i are independent realizations of the stochastic variablê the mean and variance forŵ + are given by (see black lines in figure A.4 and derivation details in appendix B). It can be shown that for τ d = ∞, the probability distribution function forŵ + is given by an inverse gamma distribution [15] ρ(ŵ + |χ) where Γ(X) corresponds to the gamma function.
. The valuesτ ± i are independent realizations ofτ ± determined by whereŵ − 1 andŵ + 1. Their mean values are approximated by  (7) the term χH − ∞ (w(τ − τ d )) alternates between 0 and χ (brown line) and we determine the values of w(τ ) at these transition points. The setŵ − corresponds to values at the transition 0 → χ, andŵ + to values at the transition In time series observed in experiment and simulation, we define proxies forŵ ± andτ ± with similar statistics. We use the topographic prominence to define the set of extrema w ± that are similar toŵ ± . The time interval between w + i−1 and w − i is the fall time τ − i , and between w − i and w + i is the rise time τ + i .
(black lines in figure A.4). Note that the mean ofτ − depends on the mean ofŵ + , while the mean ofτ + depends on the mean ofŵ − (details in appendix D). Also note that they do not depend equally on χ, on averageτ − becomes slower as χ increases, while the average ofτ + becomes faster χ increases. Numerical simulations show that their variance σ 2 τ ± behaves in a similar manner, such that forτ − the variance increases with χ, while the variance forτ + decreases (see in figure A.4(d)). In appendix F we analyze the case when equation (7) is forced with additive noise instead of multiplicative noise. For this case the expression for the mean values ofŵ ± and the varianceŵ − behave similarly to those shown in equations (10) and (11) and the means ofτ ± behave similarly to the ones in equations (18) and (19), but the variance ofŵ + is independent of χ in the additive noise case (see equation (A18)).
We can now consider the case when χ evolves stochastically over time, as in case II. If χ(τ ) evolves with a time scale much slower than τ d , then we can think of probability distributions forŵ ± andτ ± as approximated by where ρ(χ) is the probability distribution function for χ. Using these we proceed to analyze both experimental cases.

Analysis of experimental cases
To analyze the experimental cases we require observables that are proxies forŵ ± with similar statistics and that can be extracted directly from the observed time traces. We defined the set of minima w − and maxima w + that corresponds to the extrema at each observed oscillatory cycle ( figure 3(b)). These extrema are determined using their topographic prominence, an algorithmic measure that discerns between all local extrema and determines the most prominent peaks (see appendix C). We also define the rise τ + i and fall τ − i times that correspond to the time difference between points w − i and w + i , and between points w + i−1 and w − i respectively. The sets τ + and τ − serve as proxies forτ + andτ − . We have checked that the analytical expressions forŵ ± andτ ± are good approximations for w ± and τ ± in the low noise regime (see appendix E). We have been able to fit our model to both experimental cases I and II guided by equations (9)-(21). We numerically solved equation (1) (and (5) for case II), and did a parameter sweep to chose the parameter  (8) and (12) in case I, and equations (8), (12) and (20) for case II. (c) and (d) Distributions of the rise and fall times τ ± , lines and symbols as in the previous panels, the dashed lines were obtained from numerically solving equations (16) and (17) in case I, and equations (16), (17) and (21) for case II. set (table 1) that maximized the coefficient of determination R 2 . In figure 4 we plot the probability distribution functions ρ(w ± ) and ρ(τ ± ) from experimental data (symbols) and simulations (thick lines). In each case the distribution ρ(w − ) is single peaked (figures 4(a) and (b)). The distribution ρ(w + ) is always broader than ρ(w − ). By fitting the distribution ρ(w ± ) of the model to the experimentally measured distributions, we have obtained all the fitting parameter values in our normalized model (equation (7)). Now we turn to the resulting distributions ρ(τ ± ) for the rise and fall times and observe that the resulting distributions for τ ± are different between case I and II. In case I the distribution for τ + is slightly biased towards lower values compared to the distribution for τ − . Both distributions have a similar variance ( figure 4(c)). On the other hand, in case II the bias of the distribution for τ + towards lower values, and of τ − towards higher values is clear ( figure 4(d)). Also the variance τ + is lower than that for τ − .
To interpret the distributions of ρ(w ± ) and ρ(τ ± ) obtained numerically, we look at their corresponding distributions given by our analytical calculations (dashed lines, figure 4). To plot the distributions from our analytical expressions we solved numerically equations (8), (12), (16) and (17) for both cases I and II, and equations (20) and (21) for case II. The same parameters shown in table 1 were used to solve these equations. The distributions obtained from the analytical expressions do not fit perfectly those obtained from numerical simulations. Our analysis suggest that there are three main sources for this discrepancy: 1) our analytical expressions were obtained for h = ∞, while for the simulations h is finite; 2) as noise increases there is a higher discrepancy between in theŵ + and w + ; and 3) for the distributions of τ ± obtained from equations (16) and (17), there is a sharp edge to the left, failing to capture values of τ ± < τ d (dashed lines, figures 4(c) and (d)). This is because according to equations (16) and (17), any value w ∓ ≈ 1 will result in τ ± ≈ τ d . Nevertheless, there are important features between the distributions from analytical expressions and those from simulations that are conserved.
In both cases I and II the distribution of ρ(w − ) is always sharper than the distribution of ρ(w + )(dashed lines, figures 4(a) and (b)). This is because the mean and varianceŵ + depend on the values of χ (equation (12)), whileŵ − is independent of χ (equation (8)). In case I there is an overlap between the distributions for τ − and τ + (dashed and solid lines, figure 4(c)). Inspection of figure A.4 shows that at the value χ = 2 used for the fit, the mean values for μ τ ± are similar and hence give the overlap. In case II it is clear that for the distributions obtained for equations (16), (17), (20) and (21), the distribution for τ + is sharper and spans lower values compared to that for τ − (dashed lines, figure 4(d)). This feature is also conserved in the distributions obtained from numerical simulations, although for τ + it is less sharp and for τ − it is less broad. From figure A.4, this is a consequence of the fact that at larger values of χ the mean values and variance for τ ± show a clear deviation, where the values of τ + are faster on average compared to τ − , suggesting this is a generic feature of high amplitude genetic oscillations.
Finally, we fitted the same data substituting in equation (7) the multiplicative noise term with an additive noise term (appendix F). With additive noise the model cannot recapitulate the distributions for case I (figures (a) and (c)A.6), and for case II the fit is not as satisfactory as with multiplicative noise (figures A.6(b) and (d)). Taken together we conclude that a model with a time-delayed negative feedback influenced with multiplicative noise (equation (1)) is a good description for genetic oscillations.

Conclusions
In this work we have shown that complex genetic oscillations determined by many factors [19], can be represented by a simple stochastic time-delayed negative feedback oscillator with external regulation. Our model, which is motivated by reference [13], is based on a feedback loop that switches the production of a protein alternatingly on and off. This process captures details of the oscillations of two very different genetic oscillators, the circadian clock and the segmentation clock. These systems share similarities in the distributions of extrema of the oscillations w ± , as well as differences in the distributions of rise and fall times τ ± . The oscillations differ in the time dependence of their amplitudes, which we account for by differences in external regulation.
A key feature of genetic oscillations is that they are very noisy, with noise in both the oscillation period and amplitude. Our model describes the statistics of these fluctuations, and can help to understand the origin of the noises in frequency and amplitude, as well their correlation. There have been several attempts to characterise noisy genetic oscillations using their interpeak intervals [20][21][22][23], and to determine the bounds of extrema in deterministic genetic oscillations [13,24]. Our work has connected for the first time the statistics of the extrema with the statistics of timing for noisy genetic oscillations, as well their dependence on the model parameters. In addition, our work shows that the same statistics the can be used to characterize both stationary and non-stationary cases, i.e. where the external regulation P(t) is constant and when it evolves in time.
Our model highlights the importance of distinguishing between internal noise, which is noise of the oscillator process itself, and external noise, which is noise in the external regulatory process. Also, our model shows the importance of distinguishing between additive and multiplicative noise. With multiplicative noise, the fluctuations ofŵ + are stronger than those forŵ − (compare figures 4 and A.4). A biochemical origin for this fluctuations might come from the random bursting in protein production [25].
Our analysis can provide insight into the regulation of genetic oscillators. The circadian clock, which maintains a regular daily rhythm in individual cells over many years, was used as a reference case for steady state genetic oscillations. It is known that in mammals the circadian clock involves the production of CLOCK and BMAL1 which induce the production of CRY-PER dimers, which in turn exert negative feedback repression on the production on CLOCK and BMAL1 [16]. Our analysis reveals that amplitude fluctuations of the PER2::LUC signal that are slower than the oscillation cycles emerge directly from the multiplicative noise without the necessity for an external noisy regulation.
In contrast, the segmentation clock presents us with an opportunity to extend our insight into external regulation of genetic oscillators. The Her1 protein regulates its own production through negative feedback repression, but in contrast, in this case we find that external regulation has a strong effect on the Her1 amplitude profile. There is a history of modeling the transitions between oscillatory and non-oscillatory states in Hes/Her genetic circuits [26][27][28][29]. For the segmentation clock specifically, generic features of these transitions have been studied using canonical dynamical systems as the Stuart-Landau oscillator [27] and the Fitzhugh-Nagumo excitable system [29]. These studies sought to model the persistent noisy oscillations of the undifferentiated progenitor cells found in the tailbud of the embryo. In contrast, our work investigates the behavior of differentiating cells in the segmentation clock that have left the tailbud and show a characteristic, one-time developmental trajectory involving non-stationary period and amplitude. Our work extends these previous studies, first by using a model that has direct biological interpretation, and second by showing that a non-generic feature, the asymmetry between τ ± is dependent on amplitude. This asymmetry makes the oscillations slower as amplitude increases. The amplitude increase and the slowing has been observed in vivo across the antero-posterior axis of zebrafish embryos [30]. Our work suggests that it is the increasing amplitude that drives the oscillations to become slower, where this slowing is a necessary feature to observe traveling waves at a tissue level in the segmentation clock [31].
The slowing of oscillations in the embryo is poorly understood, with various mechanisms having been suggested [32,33]. Control of the amplitude and slowing is accounted for in our model by a time dependent change in the protein production rate. The model is agnostic about exactly which biochemical process are involved in the regulation of protein production, thus the model makes a general prediction about factors involved in production that can now be tested by experiment. Regulation of the production rate by the regulatory factor P(t) can be interpreted as the transition from an inactive to an active transcriptional state in a Her gene. This interpretation may be favored, because the transcriptional level is also where the Her transcriptional repressors exert negative feedback in the model. Furthermore, the transcription factor Tbx6 is known to be an activator of Her1 production [34], and is required for normal segmentation [35,36]. The levels of Tbx6 increase from posterior to anterior across the PSM [35]. Therefore in our model the activity of the external production rate P(t) could be analogous to concentrations of Tbx6 in the embryo.
Non-stationary dynamics in genetic oscillators are not limited to the segmentation clock. Examples in other cell types have been observed with roles in pluripotency and differentiation, as well as in response to pathological perturbations, and their regulation remains an open active question [37][38][39][40]. Application of this theory should help to disentangle the internal noise of genetic oscillators and the external noise of regulatory processes in naturally occurring systems, and help to optimise the design of synthetic systems.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.  (1) and (5). The variability emerges from the stochasticity in the production rate P(t). To determine the role of stochasticity in the production rate P(t), we analyzed the distribution of times taken to reach the last peak in a given time trace (t lastpeak ). In figure A.2 we show the distributions obtained from experiment, along with the distributions obtained for simulations in two different cases. In the first case, we set the production rate noise to zero, and increase the intrinsic noise of w(t) (red curve, figure A.2(a)). We note that very high intrinsic noise covers only a very limited range of the t lastpeak distribution observed in experiment, and the noise is so strong that the time traces no longer resemble those from experiment ( figure A.2(b)). On the other hand, reducing the intrinsic noise of w(t) and setting the noise of P(t) to a finite value, we note that we cover satisfactorily the range of (t lastpeak values observed in experiment (blue curve, figure A.2(a)), and the resulting distribution is similar to the experimental one. The resulting time traces are also similar to the ones observed in experiment figure A.2(c). The observed timing variability is related to the different times it takes P(t) to reach the threshold φ (see equation (5)).

Appendix B. Statistics at turning points
In the main text the dynamics of the genetic oscillator w(τ ) in the switch-like case are given by For the case when χ is constant we write an ansatz solution for a single realization of the form where G(X 1 , X 2 ) = θ(X − X 1 ) − θ(X − X 2 ), τ + thr(i) is the time when the threshold w(τ ) = 1 is crossed upwards, and τ − thr(i) is the time when the threshold w(τ ) = 1 is crossed downwards.  (7) with stochastic production P(τ ) defined in equation (5). To derive the statistics of w(τ ) at the turning points, we assume that once the threshold w = 1 is crossed, the signal does not return back. Each time the threshold is crossed, a turning point happens after a time τ d . With this in mind w(τ ) in the degradation phase with initial condition w o is Using the value of the threshold w o = 1 as an initial condition, the dynamics ofŵ − i are given bŷ Using the known property < e B(τ d ) >= e <B 2 (τ d )>/2 for geometric Brownian motion we arrive at equations (11) and (12). Since the dynamics ofŵ − behave as a geometric Brownian motion, its probability density is given by a log-normal distribution. The dynamics at the production phase are given by (see reference [15] for the solution and moments in the Ito convention). Using again the threshold value as the initial condition, the statistics ofŵ + are given bŷ the mean value shown in equation (11) is easily computed by taking the mean and evaluating the integral.
To determine the variance, the second moment is calculated as To compute the third term, special care needs to be taken for the integration limits since the correlations for the Brownian noise are different for s 2 > s 1 and s 1 > s 2 . The coefficients shown in equation (12) are

Appendix C. Topographic prominence
The topographic prominence p j is a measure used in topography to quantify the prominence between mountain peaks. The algorithm developed by Maizlish searches the elevation of a summit relative to the highest point to which one must descend before reascending to a higher summit [41]. In this work it is used to determine which of the members of the extrema v ± j of w(τ ) is significant enough to be part of the set w ± i . The extrema v ± j becomes part of w ± i if it fulfills the condition p j > ( figure A.3). As a first step, we assume that the signal w(τ ) is down sampled such that we observe w(qΔt) where the index 1 q N and Δt is the sampling time. In practice we eliminate the fastest fluctuations by using a Savitzky-Golay filter, and we take Δt as our integration time step for numerical simulations and as the inverse sampling rate for experimental time traces. Then the algorithm to determine the topographic prominence p j for each extrema v + j is as follows: 1. We draw a horizontal line that crosses v + j , extending until it either touches the signal (red line in 2. Two time windows are defined for the range that spans the line, one to the left of v + j and another one to the right. 3. For each window that corresponds to each side of v + j , we find the global minima u + bj where b = 1 corresponds to the left side and b = 2 to the right side. 4. We compare the values u + 1j and u + 2j and select the highest one to calculate the topographic prominence p j = v + j − u + bj . 5. Finally, we define w + i with the following rule, if p j > then v + j → w + i where the index i orders the maxima w + i as they appear chronologically in the time series. In figure A.3 we have two examples showing how p j is defined in each case. In these examples p j and p j share one minimum which is labeled u + 2j and u + 1j correspondingly. To define the set of w − i we run the same algorithm on the negative time series −w(τ ).
Note that the variance in the additive case is independent of χ. We have fitted the experimental data using the same parameters for P(τ ) as before, and the rest parameters are shown in table 2. The theory with  additive noise cannot fit the data from case I (figure A.6(a)), since in experiment the distributions for w ± have different variances, and in this alternative model they have the same. As for case II, where the values of χ are distributed, some features are reproduced with additive noise (figure A.6(b)). Yet, there are two problems by using additive noise: first, note that there are several values of w − that are very close to zero (figure A.6(b)), and second, by introducing additive noise, the base line of w(τ ) fluctuates around positive and negative values when the production rate P(τ ) = 0. The introduction of multiplicative noise is the proper choice to avoid these problems, and we get better agreement between experiment and simulations ( figure 4). Now we turn to obtain the mean values forτ ± following a similar procedure as in appendix D. An approximation for the mean ofτ ± is given by which are very similar to the means for the multiplicative noise case given by equations (14) and (15). Again, we note that the timing distributions ρ(τ ± ) cannot be reproduced for case I (figure A.6(c)). For case II the fit for τ − is in good agreement (figure A.6(d)), but for τ + the fit is better in the case of multiplicative noise (compare figure 4(d) with figure A.6(d)).