An improved rhythmicity analysis method using Gaussian Processes detects cell-density dependent circadian oscillations in stem cells

Abstract Motivation Detecting oscillations in time series remains a challenging problem even after decades of research. In chronobiology, rhythms (for instance in gene expression, eclosion, egg-laying, and feeding) tend to be low amplitude, display large variations amongst replicates, and often exhibit varying peak-to-peak distances (non-stationarity). Most currently available rhythm detection methods are not specifically designed to handle such datasets, and are also limited by their use of P-values in detecting oscillations. Results We introduce a new method, ODeGP (Oscillation Detection using Gaussian Processes), which combines Gaussian Process regression and Bayesian inference to incorporate measurement errors, non-uniformly sampled data, and a recently developed non-stationary kernel to improve detection of oscillations. By using Bayes factors, ODeGP models both the null (non-rhythmic) and the alternative (rhythmic) hypotheses, thus providing an advantage over P-values. Using synthetic datasets, we first demonstrate that ODeGP almost always outperforms eight commonly used methods in detecting stationary as well as non-stationary symmetric oscillations. Next, by analyzing existing qPCR datasets, we demonstrate that our method is more sensitive compared to the existing methods at detecting weak and noisy oscillations. Finally, we generate new qPCR data on mouse embryonic stem cells. Surprisingly, we discover using ODeGP that increasing cell-density results in rapid generation of oscillations in the Bmal1 gene, thus highlighting our method’s ability to discover unexpected and new patterns. In its current implementation, ODeGP is meant only for analyzing single or a few time-trajectories, not genome-wide datasets. Availability and implementation ODeGP is available at https://github.com/Shaonlab/ODeGP.


qPCR data collection protocol
For the qPCR experiments to investigate the presence/absence of oscillations, we used the following protocol for Samples 1-4 (a schematic is presented in Figure 1).A frozen vial of mESCs was first thawed and grown on a 100 mm plastic dish coated with gelatin, under pluripotent conditions, till about 60-70% confluency.The cells were then split and grown on a six well plate for about a week (with passaging) under the conditions mentioned in Table 1.Finally, for Dex synchronization and Trizol collection, the cells were transferred to 24 well-plates under the same growth conditions as mentioned in Table 1.To ensure equal cell density at time of Trizol collection, during transfer to the 24-well plates two groups were created -"Day 1" cells were plated into 5 wells of a 24-well plate corresponding to eventual collection time-points 12, 15, 18, 21 and 24 hours, while "Day 2" cells were plated on 4 wells of a separate 24-well plate representing eventual collection time-points 0, 3, 6 and 9 hours (Figure 1).Thus one full 24-hour cycle with 3-hour sampling was represented.For the Day 1 cells, 25000 cells were seeded per well while for the Day 2 cells, 12500 cells were seeded per well.After allowing a few hours for attachment in the 24 well plates, the Day 1 cells were Dex synchronized beginning from 9 am in the following manner: time-point 24 synchronized at 9am, time-point 21 at 12pm, all the way till time-point 12 which was synchronized at 9 pm.Finally, these Day 1 cells were Trizol collected next day at 10 am.Two hours later, beginning from 12 pm, the Day 2 cells were synchronized every three hours, and finally collected using Trizol at 10 pm (Figure 1).This somewhat elaborate strategy was followed to ensure that at time of collection, all time points had similar cell densities.Assuming an average division time of 12 hours for mESCs, this protocol allowed for 2 doubling times for the Day 1 cells and 3 doubling times for the Day 2 cells.The final expected cell numbers at time of Trizol collection are shown in Table 1.For Samples 5-6, an identical protocol was followed, but with an extra day and an additional collection time, to allow sample RNA collection over 36 hours.Samples 1-4 are considered "low density" since at collection they were about 40-60% confluent while Samples 5-6 are considered "high density" since they were almost completely confluent at time of collection.Note that this high density was achieved only within the 2-3 days that the cells were plated in the 24-well plates, since prior to that the cells were always maintained at less than 60-70% confluency in the 100 mm dish or 6-well plate.

Growth conditions Seeding density
Expected cell number at time of Trizol collection  2 Fold-change and error analysis for qPCR datasets

Fold change computation
We collected Ct values for a gene of interest from three technical replicates of each sample.
For every sample, we also collected three Ct values for the reference gene (i.e.GAPDH).The calculation for fold change is as follows: • We take the average for three samples (⟨C s ⟩) and three reference samples(⟨C r ⟩): • We calculate the difference of these two: ∆ct = ⟨C s ⟩ − ⟨C r ⟩.
• ∆ct at time point zero is used as reference ∆ct.We named it ∆ct Ref .
• We calculate ∆∆ct as the difference between ∆ct and ∆ct Ref .
• The fold change is defined as 2 −∆∆ct

Error propagation
To calculate the error in fold change, we did basic error propagation computation as follows: • Uncertainty in fold change depends on the error in ∆∆ct as follows: where, δ f c/∆∆ct is uncertainty in fold change due to uncertainty in ∆∆ct and δ ∆∆ct is uncertainty in ∆∆ct • δ ∆∆ct depends on uncertainty in ∆∆ct due to uncertainty in ∆ct (δ ∆∆ct/∆ct ) and uncertainty in ∆∆ct due to uncertainty in ∆ct Ref (δ ∆∆ct/∆ct Ref ) as follows: • δ ∆∆ct/∆ct can be calculated as follows: where ∆∆ct = ∆ct − ∆ct Ref and δ ∆ct is the uncertainty in ∆ct • Similarly, δ ∆∆ct/∆ct Ref can be calculated as follows: where δ ∆ct Ref is the uncertainty in ∆ct Ref • δ ∆ct can be calculated as follows: • δ ∆ct Ref is just the special case of equation 7, calculated at time point zero.
• δ ∆ct/⟨Cs⟩ can be calculated as follows: where δ ⟨Cs⟩ standard error of mean for C s (i) • δ ∆ct/⟨Cr⟩ can be calculated as follows: where δ ⟨Cr⟩ standard error of mean for C r (i)

An overview of Gaussian Process Regression
Gaussian processes (GP) provide a framework for performing non-parametric regression to estimate a target function via Bayesian inference [1].A function-space rather than weight-space approach is taken in GP regression, where instead of fixing a function and placing a prior on its parameters (to learn the optimal parameters via Bayesian inference), the functional form itself is learned by placing a GP prior over functions.A GP is simply a collection of random variables, any finite combination of which has a multivariate normal distribution.If each random variable represents the value assumed by the target function f at a particular domain point x, then the joint distribution of these random variables corresponds to the prior distribution over possible target functions.The multivariate normal is described by a covariance matrix Σ, whose entries are calculated from a kernel function capturing the smoothness and length-scales associated with the data.Note that the "parameters" of the kernel are actually hyperparameters of the inference problem -they do not directly specify any one functional form to be learnt; rather they help prioritize a class of functions with certain smoothness properties and decay lengthscales.
Let X * be the vector of domain points x ∈ R at which the value of f is to be estimated.Then the prior distribution of f (X * ) (before any data is observed) is given by Here, µ : R → R is a 'mean' function such that µ(X * ) = E[f (X * )].The covariance matrix Σ X * X * is determined by a positive semi-definite kernel function K : R 2 → R as follows: K thus describes the similarity between values of f at different domain points.The parameters of K are considered to be hyperparameters in the context of the regression.With µ and K being well-defined on R and R 2 respectively, X * can contain any number of x ∈ R, allowing us to model f at a theoretically infinite number of domain points.
Given a vector of N observed data values Y at corresponding domain points X, Bayesian inference can be carried out to obtain a posterior distribution for f : where An optimal posterior that best models the given data while avoiding overfitting can be found by maximization of the data likelihood P (Y |f ) = N (Y |µ, K).The log marginalized form of this likelihood can be expressed as With the hyperparameters of K being the only variables here, this marginal log likelihood can be maximized by finding optimal values for the hyperparameters using standard techniques.

Implementation details of ODeGP
Pre-processing the time-series data Let X be a list of N points in the domain at which the data has been observed.In our case, X is a list of times.If d distinct replicates of the data have been obtained for the N time points in X, let these be Y 1 , Y 2 , . . ., Y d .The average and standard deviation of these d lists at each of the N domain points is calculated to obtain Y and S respectively.This step is not carried out if Y and S are directly supplied as input to the method.
Next, Y is detrended.This is done by performing simple linear regression on the points in Y , and setting the new Y as the difference between the initial Y and the values of the resulting best-fit line at the corresponding time points in X.
Finally, Y is 0-centered by setting Y = Y − mean(Y ).This allows us to set the GP prior distribution's mean as zero, and simplifies further computations required for GP regression.

Kernel definitions
The diagonal kernel K D and the non-stationary kernel K N S are positive semi-definite functions (R × R) → R defined as follows: where We use K N S with Q = 1, resulting in the following updated expression: where The diagonal kernel has a single hyperparameter ϵ, whereas the non-stationary kernel has ten hyperparameters (ϵ, along with three each for the given functional forms of w(x), l(x) and µ(x)).ϵ acts as a global noise term in these kernels.The local noise encoded by S is incorporated as follows: for s i ∈ S, s 2 i is added to the ith diagonal entry of the covariance matrix (irrespective of which kernel it is constructed with).

The optimization routine
Given the data, optimal hyperparameters for each of K D and K N S are found through maximization of the marginal log likelihood (MLL).This is the log likelihood of observing the data, given the kernel and the domain points at which the observations are made, i.e. log P (Y |K, X).
where K XX is K(X, X), i.e. the N × N covariance matrix generated by applying the kernel function K on all pairs of points (x i , x j ) such that x i , x j ∈ X, and then adding s 2 i to K ii for s i ∈ S.
The DEoptim [2] package, which uses a differential evolution algorithm, is used to carry out this optimization in ODeGP.The DEoptim method is designed to exclusively carry out minimizations, so the negative MLL is passed to this method as the function to be optimized.
The (negative) MLL, as seen from 17, is a function of X, Y and K. Since both X and Y are fixed by a given dataset, the optimization is done by searching over different possible K, or more specifically different sets of hyperparameters for K.
To carry out its optimization, DEoptim requires specification of the upper and lower bounds of the search space for each hyperparameter.For the purpose of defining these bounds, some measures are derived from X and Y : • Y step: the maximum difference between two consecutive values in the list obtained by arranging the values in Y in ascending order K D has a single hyperparameter, ϵ.Since this encodes global noise in ODeGP, we pass the lower and upper bounds for this hyperparameter to DEoptim as 0 and Y range respectively.The same bounds are used for the ϵ hyperparameter in K N S .
The bounds for the remaining nine hyperparameters in K N S are defined in Table 2: Table 2: List of the hyperparameters and their lower and upper bounds used in optimization of the marginal log likelihood.
These hyperparameters control the variation of w(x), l(x) and µ(x) through the functional form described in Equation 16, which is introduced in [3].Their bounds are set as above with the following reasoning: • w min and w max are respectively the minimum and maximum possible value taken by w(x), which controls the amplitude of the correlation between the GP posterior's values at different time points.Since Y range is the amplitude of the data as a whole, we take w max ≤ Y range.w min cannot exceed Y step, else w(x) will not be able to generate correlations with amplitudes comparable to the difference between closely spaced Y values.Y step thus also acts as a lower bound for w max .
• l min and l max are respectively the minimum and maximum possible value taken by l(x), which encodes the lengthscale of correlation between the GP posterior values at different time points.It is clear that 0 ≤ l(x) ≤ X range.l max cannot be smaller than timestep, else l(x) will not be able to generate correlations between consecutive time points.timestep thus also acts as an upper bound for l min .
• µ min and µ max are respectively the minimum and maximum possible value taken by µ(x), which controls the GP posterior's oscillation period at different time points.µ min ≥ the sampling interval of the data (i.e.timestep).µ max ≤ the total time range of the data (X range), since for time periods > X range a complete oscillation will not be present in the data.µ max cannot be smaller than X range/2, else larger time periods possibly present in the data will not be captured.X range/2 thus also acts as an upper bound for µ min .
• c w , c l and c µ control the 'curvature' of the functions w(x), l(x) and µ(x); i.e. the speed with which they respectively approach their respective maximum values [3].All of these are bounded between 0 and X range, with the assumption that this curvature should not exceed the total domain range of the data.
The minimization of negative MLL is done with a maximum of 50 iterations for K D , and 200 iterations for K N S .Further parameter values used that are specific to DEoptim's implementation can be found in the code for ODeGP.
When complete, the optimization returns the set of hyperparameter values h for the kernel in consideration that results in the minimum attained negative MLL value.Let this be h D and h N S for K D and K N S respectively.

Model selection with the Bayes factor
Given X and Y , h D and h N S , let the corresponding optimal MLL values for K D and K N S be MLL D and MLL N S respectively.
The Bayes factor is calculated as a ratio of the marginal likelihoods of two competing hypotheses.
In our case, it compares the suitability of the two kernels (models) for the given data.Since the MLL values previously calculated are log marginal likelihoods, the Bayes factor is evaluated by applying the exponential function to the difference of MLL N S and MLL D .
The calculated BayesF is returned by the method as output.The larger the resulting value of BayesF, the stronger the support for K N S being a better model for the data.Hence, if this metric exceeds a defined threshold then the data is judged to be oscillatory, else it is not.

The posterior distribution and further processing
The GP posterior distribution for kernel K at any given list of M time points W has a mean µ and covariance matrix Σ which can now be calculated as follows: where K XW is the N × M covariance matrix K(X, W ), K W W is the M × M covariance matrix K(W, W ) and K XX is the same as defined in 17.Here, while evaluating the kernel function, the hyperparameters of K are fixed as h, obtained through optimization of the MLL.
The confidence interval of µ at time point w i ∈ W is evaluated as 2 √ Σ ii .The strongest detected oscillation period in µ is obtained using R's spectrum function from the stats library.

Correction to the Bayes Factor for multiple-hypothesis testing
While ODeGP is primarily intended for detecting oscillations in single trajectories, a multiple hypothesis testing correction in the Bayesian setting is also provided in case a dataset has multiple trajectories that need to be simultaneously analyzed for the presence of rhythms.Here we provide a brief explanation of the correction factor, and how it compares with the frequentist Bonferroni approach of controlling the Family-Wise Error Rate.
As discussed above, ODeGP computes the Bayes factor as the ratio of the marginal likelihoods of the alternate (oscillatory; H A ) versus null (not-oscillatory; H 0 ) hypotheses.Using Bayes theorem, it is easy to show that: where the left hand side represents the posterior odds, and the right hand side is the product of the prior odds and the Bayes Factor.When using just the Bayes factor to classify a waveform as oscillatory or not, ODeGP implicitly assumes that the prior odds is 1, i.e. π(H A ) = π(H 0 ).This results in the posterior odds being identical to the Bayes factor.In this setting, multiple hypothesis corrections can be carried out by suitably modifying the prior odds to generate a posterior odds that is less than the Bayes factor [4], as described below.
Assume there are k trajectories that are to be tested for oscillations, with corresponding null and alternate hypotheses H 0i , H Ai respectively, where i = 1, ..., k.Suppose that initially π(H 0i ) is set to 0.5 for all i, implying that the assessor is truly objective about the prior probability of the hypotheses being true or false.Given independent null hypotheses, the joint probability of all the null hypotheses being true becomes 0.5 k , which is possibly much smaller than desired.
In this case, one could begin with assigning a value for the joint probability Π 0 = k i=1 π(H 0i ), for instance 0.5, and then modifying the π(H 0i ) accordingly.If there is no preference for any particular H 0i , then the revised π(H 0i ) becomes equal to Π 1/k 0 .In this case, the prior odds for each trajectory is no longer 1 and the posterior odds in Equation.21 gets modified to: Note that the Bayes factor for each hypothesis, P (Y |H Ai ) P (Y |H 0i ) , is calculated independently as explained previously.ODeGP reports this prior odds as a correction term to the Bayes factor, with a default value Π 0 = 0.5.It can be easily demonstrated (see [4] for the derivation) that under certain conditions (k → ∞ and 1/Bayes factor → 0), this method of multiple hypothesis correction is similar to the frequentist Bonferroni correction, since the modified posterior probability π(H Ai |Y ) is scaled by a factor approximately equal to k.

Generating simulated datasets
Simulated waves were generated as sets of 3 replicates, with a duration of 48 hours and an interval of 3 hours between consecutive observations (when no data was missing).Each such set of 3 replicates was produced as follows: a base waveform was generated at the designated timepoints -let these times be contained in a vector T , the corresponding waveform values be contained in a vector D, and the length of T = length of D = N .Taking the standard deviation of D to be σ D , each replicate was generated by adding to D, a vector of N independent samples drawn from U (−σ D , σ D ).
The base waveform D is generated by applying a function f (t) to the times in T (see Tables 3, 4, 5 and 6 for details).Following this, a fraction γ of the times in T are randomly chosen, and these times along with the corresponding values in D are deleted from T and D respectively.
Construction of an ROC curve requires a dataset consisting of both oscillatory and non-oscillatory waves.Each such dataset was generated as a collection of 500 oscillatory waves and 500 nonoscillatory waves, where each wave is a set of 3 replicates produced from a defined base waveform as described above.
For oscillatory waves, the definitions of f depend on the type of waves being simulated, and are provided in Tables 3, 4, 5 and 6 along with the values of γ used in different cases.For non-oscillatory waveforms, f was defined as f (t) = N (0, σ), and γ = the γ being used to generate the 500 oscillatory waves comprising the other half of that particular dataset of 1000.

Specification of parameters for running existing methods
Let X be the list of time points at which the data has been observed.We have earlier defined X range = max(X) − min(X).Considering that the individual elements of a list are accessed by 1-indexing in R, we additionally define MetaCycle's [5] meta2d function requires the minimum and maximum period length of interest to be passed as arguments.When running MetaCycle on our simulated datasets, these were respectively specified as • maxper = X range RAIN's [6] rain function requires the period to test for and/or the range of periods to test for to be passed as arguments.These were respectively specified during testing as: • period = X range/2 • period.delta= ( X range/2)− delta t For all methods tested that accept the raw data (without merging of replicates) as input, the replicates were passed to the method in an appropriate format.

AUC values of tested methods on simulated datasets
The numbers assigned to simulated datasets in the previous section describing their generation are used here.

Analysis of model complexity penalisation by MLL
The marginal likelihood 17 computed during GP regression automatically incorporates a measure of the complexity of the model being considered through the log|K| term.A larger log|K| for a given model generally results from higher covariance values between the data's values at different time points, which produce wider confidence bounds in the GP posterior.Wider bounds allow a greater flexibility of functions in the posterior distribution, indicating that the model is more complex.For this reason, we did not additionally penalise the number of hyperparameters of both models when comparing the two (something that is generally done in regression methods which use the BIC/AIC for model selection).Importantly, we found that the Bayes factor was distinctly different between datasets that were known to be oscillatory versus non-oscillatory, thus highlighting the usefulness of this metric in the classification problem.
The expression for K N S in 15 contains a variable Q, which denotes the number of components [7] of K N S .Every increase of 1 in Q adds 9 new hyperparameters to the defined K N S .We show below that increasing Q (and thus the number of hyperparameters) does not cause overfitting of the data in the posterior generated by K N S , and that the MLL does penalise the more 'complex' resulting model represented by K N S .
Non-oscillatory data was simulated by sampling from N (0, 1) at 3 hour intervals for a duration of 48 hours, and generating three replicates from this base waveform as described above.ODeGP was run on this data with multiple values of Q.The datafit and complexity term of the MLL for K N S were respectively calculated as −0.5Y  Figure 2 demonstrates that increasing the number of hyperparameters used in the non-stationary kernel does not enforce a closer fitting of the data by the posterior mean.This is also reflected in the datafit term of K N S in Table 7.The complexity term worsens as well (though not monotonically) and there is a clear trend of reduction in the Bayes factor.Overall, a larger Q leads to a stronger indication of K D representing a better model for the data.

Application of ODeGP to qPCR datasets
To understand how well ODeGP works on existing qPCR datasets, we analyzed gene expression profiles reported in a prior work [8].While the expression patterns and analysis of Rev-ERBβ and Per1 were shown in the main text Figure 4, here we present analysis of some of the other reported genes, Per2 and Npas2 along with the unsynchronized data for Per1.It is clear that between the Dex and vehicle treated datasets, the Bayes factor is always separated by an order of magnitude or more, demonstrating the power of ODeGP in classifying oscillatory datasets.
Similar to the analysis of the synchronized Per1 data shown in Table 4 of the main text, we also evaluated the Bayes factor bounds produced by each method for the unsynchronized Vehicletreated Per1 data in Figure 3C.The results of this are shown in Table 8.Barring MetaCycle and RAIN, all methods correctly classify the data as non-oscillatory with p-values greater than 0.05.The Bayes factor bounds obtained for the remaining methods are all sufficiently close to 1, strongly indicating that the data does not contain oscillations.ODeGP's non-stationary kernel posterior for this data in Figure 3C shows very wide confidence intervals, confirming that this kernel is not a good model for it.3C.p-values/Q-values were converted to corresponding Bayes factor bound [9] values where applicable to allow comparison with the Bayes factor returned by ODeGP.

Comparison of ODeGP with other GP based algorithms
Though not explored extensively, there are a few prior examples of the use of GPs in modeling biological data, for example in identifying differentially expressed genes [3,10], detecting oscillations [11,12] and the discovery of spatial patterns in gene expression [13].The method in [12] uses a principle similar to ours with the comparison of a non-oscillatory GP model to an oscillatory GP model.However, this method does not explicitly incorporate a non-stationary functional form to encode correlations between the data at different time points in its oscillatory model, and also compares its performance only with that of the LS Periodogram (which never falls among the best 3 performing methods in our analysis).GPrank [10] on the other hand uses GPs to model genome-wide time series data.Though it was not created for the specific purpose of detecting oscillations, the workflow followed by this method is similar to ours in terms of the combination of two GP models with Bayesian model selection.The significant difference between GPrank and ODeGP is in the choice of kernel used to define the alternate hypothesis (the RBF kernel is used in GPrank).Our results show that ODeGP almost always outperforms GPrank (Tables 2 and 3), demonstrating the importance of careful choice of the kernel for the specific problem at hand.This point also highlights the flexibility of GPs in modeling various kinds of datasets simply by appropriate choice of kernel functions.Similar to GPrank, the output metric used by ODeGP to make a decision on the presence of oscillations is the Bayes factor, which is a ratio of the marginal likelihoods of the competing models.
9 Modeling coupled oscillators to investigate the role of amplitude and coupling strength on the Bayes Factor We used a Poincare model of coupled phase-amplitude oscillators [14] to investigate how the Bayes Factor from ODeGP is affected by increasing the amplitude of single cell oscillators as well as coupling between them.Converting radial coordinates of the ith oscillator r i , θ i to Cartesian coordinates x i , y i with r i = x 2 i + y 2 i , we get the following equations for the Poincare oscillator: where Here N is the total number of single cell oscillators which are coupled via the mean field M , with the coupling strength parameterised by K.The free-running time period of the ith oscillator is τ i .All oscillators were assigned the same amplitude A i = A, and the amplitude relaxation rate of the ith oscillator is γ i .x i (t) and y i (t) were computed using the Euler-Maruyama method, solving the Equations 23 and 24 over time interval [0, T ] with timestep ∆t.

Figure 1 :
Figure 1: Schematic of qPCR experimental design to maintain the same cell density across all time points.Each rectangular block (brown and green representing Day 1 and Day 2 cells respectively) corresponds to a sample that was dexamethasone synchronized for 1 hour.Times of dex-synchronization are mentioned below each block.Trizol collection was performed at 10 am on Day 2 for Day 1 cells, and at 10 pm on Day 2 for Day 2 cells.The circadian time (time after synchronization) is mentioned inside each block.This protocol ensured that at the time of RNA collection, each sample across the 24 hours had approximately same cell density.

Figure 2 :
Figure 2: Posterior distributions of the non-stationary kernel with varying number of components Q, obtained by application of ODeGP on simulated non-oscillatory data.

Figure 3 :
Figure 3: Application of ODeGP on Per2, NPas2 and Per1 expression data measured using qPCR [8].(A) The posterior distributions (mean in red, standard deviation in blue) obtained by defining a GP prior using the non-stationary kernel, and performing GP regression on the relative gene expression of Per2 after treatment with DMSO (above) or Dexamethasone (below), measured at 4 hour intervals over a 48 hour duration and averaged over 3 biological replicates.(B) The GP posteriors of the non-stationary kernel applied to relative gene expression of Npas2 when treated with DMSO (above) or Dexamethasone (below), where the posterior mean and standard deviation are in red, blue respectively.(C) The GP posteriors of the non-stationary kernel applied to relative gene expression of Per1 when treated with DMSO (above) or Dexamethasone (below), where the posterior mean and standard deviation are in red, blue respectively.

Figure 4 :
Figure 4: Variation of the Bayes factor obtained by applying ODeGP on simulations of the Poincare model of coupled oscillators.(A)-(C) Bayes factor as a function of the amplitude A for fixed values of coupling strength K = 0.1, 0.15 and 0.2.(D)-(F) Bayes factor as a function of K for fixed values of A = 1, 2.5 and 5.In all panels, the different colours represent different noise strengths σ.

Table 1 :
Description of various samples, seeding densities and substrates used for the qPCR experiments.

Table 3 :
Functional forms and parameter values used to generate simulated stationary symmetric datasets.

Table 5 :
Functional forms and parameter values used to generate simulated non-stationary symmetric datasets.OR at the completion of each time period, a new time period τ is sampled from U (τ 1 , τ 2 ).If the total time elapsed upto the start of this new oscillation is ϕ, then

Table 6 :
Functional forms and parameter values used to generate simulated non-stationary asymmetric datasets.

Table 7 :
T K −1 XX Y and −0.5logK XX , as in 17.Comparison of the Bayes factor and individual terms of the MLL for K N S , obtained by testing ODeGP on the non-oscillatory data shown in 2 with different values of Q in K N S .A higher value of the datafit term indicates a closer fit of the data by the K N S posterior mean, and a higher value of the complexity term indicates a less complex (and more favourable) model generated by K N S .

Table 8 :
Comparison of output metrics of all methods tested on the Vehicle-treated Per1 data shown in Figure