BACKGROUND
pattern recognition is central to many scientific disciplines and is often a first step in building a model that explains the data.
in particular, the study of periodic phenomena and frequency detection has received much attention, leading to the well-established field of spectral analysis.
biology is rich with  periodic behaviour, with sustained oscillations in the form of limit cycles playing important roles in many diverse phenomena such as glycolytic metabolism, circadian rhythms, mitotic cycles, cardiac rhythms, hormonal cycles, population dynamics, epidemiological cycles, etc.
<cit> .
a conventional method for frequency detection is fourier analysis.
it is based on the fact that it is possible to represent any integrable function as an infinite sum of sines and cosines.
the fourier transform  uses this property to reveal the underlying components that are present in a signal  <cit> .
fourier theory has given rise to a wide range of diverse developments and far-reaching applications, demonstrating the theory's undisputed importance and impact.
for frequency detection, however, it is known that the ft works optimally only for uniformly sampled, long, stationary time series.
furthermore, common procedures of pre-processing the data can cause problems.
time series can contain low frequency background fluctuations or drift that are unrelated to the signal of interest.
for the ft, it is then necessary to remove the trends using detrending techniques.
as has been shown previously, this detrending leads to convolution of the signal that can both remove evidence for periodicity and add false patterns  <cit> .
another known problem is aliasing.
if a signal containing high frequencies is recorded with a low sampling rate, peaks of high frequencies can fold back into the frequency spectrum, appearing as low frequencies  <cit> .
the gibbs phenomenon provides another example where spurious peaks appear in a ft. it occurs at points of discontinuity in a periodic function, and results in so-called ringing artefacts around the "true" frequency peak  <cit> .
as for the accuracy of the frequency estimate, no direct information of this is given by the output from a ft, since the sharpness of the peaks depends on a combination of factors such as noise levels and the length of the time series.
for further details, see the extensive ft literature .
wavelet transforms  <cit>  offer an attractive alternative to fourier transforms.
the main difference is that they are localised in both the time and frequency domain.
this property makes wavelets better adapted to problems with truncated data.
wavelets have found wide-ranging applications and have proven to be particularly powerful for image processing and data compression  <cit> .
bayesian inference provides another approach for analysing data .
it addresses additional aspects of the problem, such as the inherent uncertainty of the analysis and the effects of external noise.
using this framework  <cit> , the method of bayesian spectrum analysis  was developed by bretthorst  <cit>  and applied to nuclear magnetic resonance  signals and parameter estimation with great success  <cit> .
there are several advantages of the bayesian approach, including an inherent mechanism for estimating the accuracy of the result and all parameters, as well as the ability to compare different hypotheses directly.
focus is shifted to the question of interest by integrating out other parameters and noise levels.
initial knowledge of the system can be incorporated in the analysis and expressed in the prior probability distributions.
there has been a recent flood of bayesian papers with some convincing applications and promising developments in systems biology .
the bayesian approach to time series analysis has proven its value in fields such as nmr and ion cyclotron resonance analysis .
in this paper, we describe the development, implementation and testing of bayesian model development coupled with bsa and nested sampling, in a biological context.
we present a comparison of this approach with the ft, applied to a number of simulated test cases and two types of biological time series that present challenges to accurate frequency detection.
we first present some necessary background, upon which we build to develop to our approach.
bayesian inference
data is rarely available in sufficient quantity and quality to allow for exact scientific deduction.
instead we are forced to infer models from incomplete knowledge.
bayesian inference is based on bayes' rule, which evaluates a hypothesis, h, in light of some data, d, and some prior information, i. it is a method of assigning probabilities based on the current state of knowledge, allowing for subsequent re-evaluation as new data becomes available.
the goal is to determine p , the posterior probability distribution of the hypothesis, given the data and the prior information.
with bayes' rule, the posterior can be expressed in terms of other probabilities as
where p is the probability of observing the data given the hypothesis and the prior, p is the prior probability of the hypothesis, and p is the probability of the data given the prior.
when the hypothesis is the variable and the data is held constant, p is called the likelihood function, and when the hypothesis is constant it is called the probability of obtaining a specified outcome .
when evaluating only one hypothesis, p is a normalising constant, but when investigating more than one hypothesis this term plays a key role and is called the evidence  <cit> .
bayesian spectrum analysis
our presentation in this section follows closely that of bretthorst  <cit> .
the aim is to infer the most probable frequency , ω, from the given data.
we start by building a model  for the observed data, parameterised by angular frequency, ω, and amplitudes, c, and then use bayes' rule to compute the posterior probability of the parameters, p. by assigning priors to the model parameters c and integrating over these, we arrive at the posterior probability for the parameter of interest, ω, .
this is referred to as the marginal posterior probability of ω. we note that ω is an r-tuple, {ω <dig>  ω <dig>  ..., ωr}, with as many elements as there are distinct frequencies in the data.
a general model for observed data sampled at n discrete time points, d = {d, ..., d}, includes the signal of interest, s, a possible background function, g, and the noise present in the system, e,
the signal function will usually be unknown and may be complicated, but can be approximated by a linear combination of ms model functions, ψi, that we parameterise by the quantity of interest, ω:
in which  are the expansion coefficients.
similarly, the background function, g, can be approximated by a set of mg functions, ζi, that are independent of ω,
where  are the background model function expansion coefficients.
since a and b are not the main focus of the analysis, we will aim to integrate them out of the equations by marginalisation.
parameters that are treated in this manner as often referred to as nuisance parameters, which we denote here by .
although the signal function depends on ω, whereas the background function does not, for notational purposes we introduce the set of model functions, ϕi, which consists of both ψi and ζi.
this allows us to condense the model equation into
such that now d = f + e and m is now the total number of model functions, ms + mg. the model functions will typically not be orthogonal functions over the time series domain.
this, however, can be achieved by cholesky decomposition.
in all subsequent calculations an orthogonal basis is used.
from bayes' rule, the joint probability distribution of the model for the parameters ω and c is
the likelihood function, p, is calculated by comparing data produced by the model signal, equation , to real experimental data.
if the model perfectly captures the signal, the difference between the model data and the real data is simply the noise in the system.
the model of the data in equation  includes noise, e, which we assume to be time independent in the further developments.
the true noise level is unknown, but for a given noise power, σ <dig>  the principle of maximum entropy leads to the use of a normal distribution,
a noise model of this form ensures that the accuracy of the results is maximally conservative for a given noise power.
we will later integrate over all possible noise levels to remove the dependence on σ. with the described signal and this noise model, the likelihood was calculated by bretthorst  <cit>  to be
where n is the number of data points, and
where  is the mean-square of the data, , and Φjk is the matrix of the model functions, .
the goal of the analysis is to compute the posterior probability for frequencies in the data, i.e.
to go from the joint probability distribution to a posterior probability of ω, independent of the other parameters.
by integrating over all possible values of the parameters σ and c, the remainder is the marginal posterior of the parameters of interest, ω={ω <dig>  ω <dig>  ..., ωr}.
this is an essential advantage of the bayesian framework, allowing the analysis to focus on estimating the parameters of interest, regardless of the values of the others.
if necessary, the other parameters can be estimated at a later point.
to integrate over the σ and c values, priors must first be assigned to them.
we chose uniform priors for c and ω, representing complete lack of knowledge.
we know that σ is continuous and must be positive, and in such cases a jeffreys prior is appropriate, p = 1/σ.
both the uniform distribution over continuous variables and the jeffreys prior are known as improper priors if bounds are not specified as they cannot be normalised.
for more information on prior assignment see  <cit> .
using the general model, equation , assigning the priors, calculating the likelihood function, equation , and integrating out the amplitudes and noise parameters, the posterior probability distribution of ω is proportional to
where h is the projection of the data onto the orthonormal model functions, , and  is the mean-square of the hj, ,  <cit> .
this expression of the posterior allows us to identify the strongest frequencies present in the data.
for a good model, there will be a high probability peak in the posterior distribution at that ω = {ω <dig>  ω <dig>  ..., ωr}.
RESULTS
we employed the framework developed by jaynes  <cit>  and bretthorst  <cit>  to investigate the frequency components in a number of biological time series.
model comparison
after evaluating the probability of parameters in light of a certain hypothesis, it is important to question the validity of that hypothesis.
thus, the next step in bayesian inference is to compare the probability of different hypotheses.
the hypothesis is now a particular model of the signal, hi, out of a set of m models {h <dig>  ..., hm}, and using bayes' rule, the posterior probability of this model is
then two different models, hi and hj, can be compared by taking their ratios,
the probability of the data given our prior information, p, which was a normalisation constant in equation , will now vary between models, and is called the evidence.
it evaluates the fit of the data to the model, whilst penalising models that include more parameters.
each additional model parameter should be followed by a significant increase in probability, otherwise the simpler model is preferred.
thus, bayesian model comparison naturally follows the principle of occam's razor  <cit> .
model development
it will often not be obvious which function to choose to model trends in the data, so an approach using basis functions and expanding these to different orders will be of advantage, as in equation .
each expansion represents a different model, hi, and these can be compared using inference techniques.
likewise, different functions for capturing the signals in the data and modelling a different number of signals correspond to different models for data.
following  <cit> , we use the posterior ratio to evaluate different models.
this model ratio can be used to determine the number of background model functions for each time series.
the posterior probability ratio is calculated between model hn and hn+ <dig>  where hn is a model including n background functions.
to obtain the model ratio, priors are assigned to each of the models and their likelihood functions are calculated.
assigning equal prior probability to all models reduces this to the ratio of evidences.
to compute the evidences we need to integrate the likelihood, p from equation , over ω, σ, and c for each model hn.
by assigning proper normalised priors to all model parameters it is possible to integrate over them around the maximum likelihood estimate.
following bretthorst's derivation for location parameters  <cit> , we assign gaussian priors to the amplitudes with hyperparameters for the variances.
since the variances are scale parameters, they are subsequently assigned jeffreys' priors with an upper and lower bound.
this allows us to normalise them and integrate, leaving the defined bounds as parameters in the final equation.
for models with the same bounds these terms cancel out in the model ratio.
the evidence for a given model, hn, was calculated by bretthorst  <cit>  to be
where δ, γ and σ are the prior variances for amplitudes, frequencies and noise, respectively, rδ, rγ and rσ are the ratios of the integral bounds for these variances,  is the mean-square projection of the data onto the orthonormal model functions at the maximum likelihood point for model hn,  is the mean-square of the ω value that maximises the likelihood, , and r is the number of ω parameters, ω = {ω <dig>  ..., ωr}.
the jacobian  is obtained by orthogonalising the taylor-expansion of  around the maximum-likelihood point, .
see bretthorst for further details  <cit> .
for cases in which the number of frequencies in the data exceeds the dimension of omega, for instance multiple frequency data with a single frequency model, the above approximation for the evidence is ill-suited as the posterior will cease to be unimodal.
for such scenarios, either multiple expansions or mcmc offer attractive solutions to marginalisation.
for comparison we have included results from nested sampling  <cit>  as a means to perform the integration and compute the evidence.
nested sampling is a variant of mcmc that employs a likelihood based sorting of sample points to efficiently guide the search strategy of the posterior distribution  <cit> .
when the model ratio, hn/hn+ <dig>  becomes greater than  <dig>  the simpler model, hn, is favoured over hn+ <dig>  <cit> .
adding more background functions than are justified by the data  may lead to a lower probability for the frequency and in some cases possibly a location shift.
this model development approach used for the background functions above can also be used to decide on the number of underlying frequencies in the data.
the model ratios of a time series containing one frequency  and a time series containing two  are presented in table  <dig>  analysed with both a one- and two-frequency model.
the results show, as expected, a preference for the one-frequency model in case a, and for the two-frequency model in case b.
case a has only one frequency, nω =  <dig>  and case b has two, nω =  <dig>  the model ratio of a model including one frequency, h1ω, and a model with two, h2ω, is well above  <dig> for case a, showing that the simpler model should be chosen.
for case b, the first ratio is close to  <dig>  showing that the more complex model with two frequencies is far more likely.
values below 10- <dig> are listed as  <dig>  the second ratio is above  <dig>  preferring the model with two frequencies to the one with three.
we point out that the proposed method stops once the current best model has been found but is not guaranteed to find the global maximum from a predefined set of models.
the procedure is thus part of model development rather than model selection.
if the set of hypotheses are known in advance then the posterior ratios over the full set should be used to find the best model.
testing
we first show the power of the bsa approach on test cases using simulated data.
in these tests, we sought to recover known input parameters from the simulated data, to validate the bsa approach.
we employed sines and cosines as model functions ).
for comparison, discrete fourier transforms were computed using fast fourier transforms   <cit> .
in the test cases, we varied key parameters such as noise levels, trace length, sampling intervals, amplitudes, frequencies, background trends and shape of oscillations.
representative cases of noise levels and background trends are shown in table  <dig>  including fft results on the same data set.
a key observation is that the bayesian approach extracts the correct answer from the data with high precision.
bsa also computes the signal-to-noise ratios which is a useful indication of how much of the data cannot be accounted for in the model.
furthermore, the amplitudes do not impact the bsa results since they are integrated out.
each time series was generated with a sine function of angular frequency, ω, of  <dig>  rad/s with a level of noise in amplitude, ea, and phase, ep.
in some time series a background trend  was included, and in case number  <dig> an additional sine function of  <dig>  rad/s is present.
the resulting function was sampled  <dig> times at an interval of  <dig> s. results from fft are presented in the form of the angular frequency with the highest power, , and the estimated standard deviation, σfft.
the bayesian frequency estimate at the maximum posterior point is denoted by  and its standard deviation by σbsa.
for comparison, the expectation value of ω and its standard deviation computed using bsa and nested sampling  are denoted by  and σbsa-ns.
values of σ below 10- <dig> are listed as  <dig>  the estimated signal-to-noise ratio  from the bayesian analysis is given in the last column.
the bsa and bsa-ns approaches deliver the same results, apart from the case of multiple frequencies in a 1d search of ω , which for bsa-ns leads an intermediate estimate between the frequencies with a higher standard deviation.
bsa has a clear advantage over fft when the data is non-uniformly sampled.
fft requires uniform sampling, whilst bsa is less stringent and delivers the correct result with higher precision.
bretthorst also noted that non-uniformly sampled data removes aliases from the frequency domain, another significant advantage  <cit> .
five further distinct cases emerged from the tests in which bsa delivers superior results to fft: time series which have background trends, few data points, high noise levels, multiple frequencies, and non-harmonic oscillations.
background trends
additional file  <dig>  figure s <dig>  is an example of a time series with a strong background trend.
in table  <dig>  the model ratios for different background functions are shown.
the ratio is initially well below  <dig>  but the ratio of models with expansion orders of two and three legendre polynomials is above  <dig>  thus, background functions of legendre polynomials to expansion order two is more likely, and should be used in the estimation of ω.
posterior probability ratios of models including a different a number of background functions, ζ, in the analysis of the time series in additional file  <dig>  figure s <dig>  the first two ratios favour the more complex model .
the ratio between models with two and three background functions, h2ζ /h3ζ, is above  <dig> thus favouring the simpler model.
the analysis would then normally automatically stop at two background functions, but for demonstration we include more functions here.
the ratio stays above  <dig> thereby building a chain of decisions that always prefer the simpler model.
examples 8- <dig> in table  <dig> also include trends, and without pre-processing fft cannot pick out the correct frequency.
in contrast, bsa includes background functions in the model signal and delivers the desired result.
including background functions, however, results in over-estimation of the signal-to-noise ratio.
short time series
additional file  <dig>  figure s <dig>  shows the results from analysing a short time series.
the fft power spectrum is very broad , which comes as no surprise given the fft dependence on the number of data points.
bsa estimates the correct frequency sharply, but the maximum probability drops compared to longer time series .
this demonstrates the higher uncertainty associated with fewer time points.
high noise levels
bsa is also successful at handling high levels of noise, as highlighted in examples 1- <dig> in table  <dig>  the frequency estimates are correctly reproduced by the fft.
in these simple test scenarios, the bsa posterior probability distribution estimates the frequency with a significantly higher precision than the fft.
whereas the estimated uncertainty of parameter expectation values is a built-in aspect of any probabilistic treatment such as bsa, fft has no inherent mechanism for assessing the accuracy of the results.
the fft output is summarised by the average, , and the standard deviation, σfft, over the transformed data set.
we show that different noise levels influence the σfft more than the σbsa .
multiple frequencies
example  <dig> in table  <dig> has two frequencies present in the data.
both bsa and fft show these two peaks in the resulting plots.
although bsa can be used in this manner with a one-dimensional ω to scan through frequency space and estimate the number of frequencies in the data and their location, if more than one frequency is present, the model should be extended to reflect this.
without this extension the integration procedure around a single point is not well suited, so we employed nested sampling to compute the marginalisation in these cases.
for the extension approach, when the posterior probability over ω = {ω1} reveals two strong frequencies, then a better model would be ω = {ω <dig>  ω2}.
for example, additional file  <dig>  figure s <dig>  shows bsa and fft results for a test case that includes higher harmonics which give rise to multiple peaks in the log p plot.
if more than one peak in the resulting posterior probability emerges, then the model can be extended further.
one peak in the posterior probability over the number of modelled frequencies signifies that the correct number of frequencies has been captured.
as another example, additional file  <dig>  figure s <dig>  shows the result of a two-frequency search.
the bsa posterior probability distribution is now two-dimensional with a peak at the two correct frequencies .
the fft results are also shown .
additional file 5a, figure s5a, shows a time series with a high noise level and two very close frequencies of  <dig>  and  <dig>  rad/s.
fft cannot distinguish them and shows only one peak .
bsa breaks the resolution and precision limitations inherent to fft by introducing a continuous probability distribution instead of the fixed number of points and can therefore sample the posterior more finely in areas of high probability.
this approach gives rise to a high-resolution probability plot in which two distinct frequencies emerge .
the peaks have a larger variance at this local level, but the qualitative information of two underlying frequencies is revealed.
to develop bsa further, we used windowing of the time series to compute the posterior probability distribution of ω at each time point.
we call this bsa local .
the robustness and negligible peak broadening of bsa with fewer time points allows for this windowing to proceed without the introduction of artefacts due to truncation.
this local bsa captures changes in frequency, as shown in figure 2b and figure 3b.
the bsal was compared to short-time fourier transform , which is a windowed fourier transform, and to wavelets .
for the wavelet power spectrum a morlet mother wavelet was used  <cit> .
the advantages of bsal are that it remains within the same bsa framework, has high accuracy, and does not require pre-processing of the data.
non-harmonic oscillations
bsa results for oscillations with a non-harmonic shape are superior to the fft.
it highlights an essential difference in the two methods since biological data is often repetitive, but with a wide range of oscillatory patterns.
to demonstrate this further, figure 4a shows a time series simulated from an ordinary differential equation  model of cellular calcium  signals  <cit> .
such time series presents two potential problems: the time series is chaotic and thus not perfectly periodic, and the signal shape is non-harmonic.
the calculation of interspike intervals  of the time series show that multiple intervals are present .
the highest peak of the fft plot  suggests that the entire time series is one period, while bsa suggests a strong angular frequency around  <dig>  rad/s .
the bsa suggestion is similar to the second fft peak.
this highlights the differences between frequencies in the data and spike intervals.
isi are a common way of characterizing spike data, however, multiple isi need not correspond to multiple frequencies in the data.
of the four strong isi shown here, both bsa and fft identify only one of these as a regular period.
summary
after extensive test cases we find that bsa delivers superior results in cases where the fft assumptions are too constraining, most notably in the five cases above.
bsa is a flexible method allowing the underlying hypothesis to be changed depending on the focus of the analysis, and to directly compare the validity of different hypotheses.
it can handle non-uniformly sampled data and has no need for pre-processing procedures.
the price of these superior results comes at a computational cost that ranged from tens to hundreds of seconds for the examples shown here.
calcium spiking data
the first biological data set comes from intracellular signalling in plant-microbe interactions.
symbiotic bacteria induce calcium oscillations, called ca2+ spiking, in legume root cells .
these are non-stationary and often noisy time series, causing problems in identifying periodicity.
one hypothesis for signal transduction in this system is via frequency encoding  <cit> , so concluding whether there is underlying periodicity, and at what frequency, is of great interest.
the ca2+ spiking has background trends present due to fluorescence bleaching and cell movements, which are assumed to be unrelated to the underlying signal in the cell.
therefore, accounting for the background functions plays a key role in the analysis.
example time series are shown in figure 5a.
nine spiking cells from the model legume medicago truncatula were analysed for an underlying period.
the data is obtained by microinjecting a root hair cell with the calcium indicator dyes oregon green  and texas red , and exposing the plant to the bacterial signal molecule that induces the oscillations.
the data is a ratio of the fluorescence from the two dyes, showing changes in ca2+ concentrations.
the data has been published in  <cit> .
the fft of the ca2+ data results in a very broad periodogram, due to multiple frequencies and high noise levels .
also, the spiking produces a non-harmonic signal, which might be another problem for the fft.
for comparison, we also present results from the multitaper method .
the mtm is a non-parametric method of spectral analysis that uses tapers to minimize the variance in the power estimate, and is targeted at short and noisy time series  <cit> .
the mtm results were very similar to the fft .
these periodograms do not address the question of interest: is there a key period in the ca2+ signal?
in the bsa analysis , the ca2+ spiking data used the legendre background functions to an expansion order of 1- <dig>  depending on the individual trace.
nested sampling was used to compute the evidences.
frequencies with high probabilities were picked out, but varied in the interval of approximately 50- <dig> s .
however, the strongest periods were in the interval of 75- <dig> s. if periodicity plays a role in the signal transduction of this system, then the key period should be in this interval.
the signal-to-noise ratios were relatively high, between 100- <dig>  possibly as a consequence of including several background functions.
analysed calcium oscillations in m. truncatula root hair cells in response to bacterial symbionts.
the bsa shows strong underlying signals in the data, but the cell-to-cell variability is high.
bsa and the bsa-ns disagree when there are multiple peaks present, as can be seen by the large standard deviation, σ, of the bsa-ns average period estimate.
circadian data
the second biological data set shows gene expression of so-called clock genes.
many processes in plants follow a circadian rhythm .
a number of genes in arabidopsis thaliana have been shown to regulate circadian rhythms, and time series of rna levels show how these clock genes are expressed in cycles  <cit> .
time series with only a couple of cycles are common in biology and provide another suitable test case.
for these circadian rhythms, we chose to analyse rt-pcr data from four clock genes in two genotypes of a. thaliana.
the plants are either wild type, fri;flc, or mutants, fri;flc, of the genes fri and flc.
the rna was extracted from seedlings, and each time series is an average of two biological replicates.
an example of the rna levels of a clock gene is shown in figure 6a.
the data has been published in  <cit> .
fft on the rna levels of these clock genes did not give any clear periods, either having only a vague peak or none at all .
this is caused by the fft's dependence on the length of the time series, which in this case was only 1- <dig> cycles.
the mtm method had more of a peak in the 20- <dig> h period, but still lacking in precision .
bsa on the other hand provides a clear peak close to  <dig> h , consistently for all eight time series .
nested sampling was used to compute the evidences.
the assigned probabilities are relatively low, but the signal-to-noise ratios were between 2- <dig>  and similar probabilities were obtained using simulated data with few data points and high noise levels.
the period peaks are very stable over all the time series, and suggests a probable period which is unaffected by the mutations .
this is in agreement with the original conclusions of the experiment  <cit> .
the bsa results of rna levels of four so-called clock-genes in a. thaliana, measured in two different genotypes.
the plants are either wild type, fri;flc, or mutants, fri;flc, of the genes fri and flc.
all eight genotypes displayed oscillations in rna levels with a period of approximately  <dig> hours.
both bsa methods identify similar periods.
CONCLUSIONS
bayesian inference offers a powerful way of analysing biological time series.
despite the undisputed value of fourier theory, there are cases when the necessary requirements for its optimality for time series analysis are not met.
this is a consequence of the underlying assumptions of a fourier transform, causing it to work optimally only for uniformly sampled, long, stationary, harmonic signals that have either no or white noise.
in biology these requirements are rarely fulfilled, requiring pre-processing of the data, such as noise reduction and detrending techniques, with the risk of convoluting the signal and losing valuable information.
by placing the problem of frequency extraction in the framework of bayesian inference, the known and well-documented problems of fourier analysis can be overcome.
this approach also breaks the resolution and precision limitations inherent to the fft by introducing a continuous probability distribution instead of the fixed number of points maintained by the discrete fourier transform.
as we demonstrated here, bsa coupled with automated model development can give superior results to the fft when faced with short, noisy time series, non-stationarity and non-harmonic signals.
the suggested automated model development worked well in our hands but must be used with caution in practice as the approach is not guaranteed to find a global optimum in model space.
alternate models should be explored and compared using posterior probability ratios or approximations thereof.
we found nested sampling  <cit>  to provide a powerful means of estimating evidences for cases in which a single peak could not be identified.
other mcmc techniques such as simulated annealing running in parameter exploration mode or standard metropolis-hastings algorithms offer attractive alternatives  <cit> .
bsa calculates signal-to-noise ratios, provides parameter precision estimates, and can handle high noise levels as well as background trends and therefore has no need for pre-processing.
more importantly, the bayesian framework offers flexibility in the underlying model and enables direct comparison of hypotheses.
the work presented here is a merely a first step in this direction.
we have employed conservative priors  that make an analytical treatment tractable but in some cases more information could warrant a different choice of prior that might require substantial alternations to our approach to handle the numerics of marginalisation.
there are many known examples in biology in which oscillations play a key role and methods for their detection will be of value, especially in cases where subtle differences are of importance and for short, noisy time series.
in the presented examples, we demonstrated the improvements that can be gained from employing this approach.
although in these cases, the biological conclusions would not have changed, one can envision scenarios in which a higher accuracy in frequency detection may allow subtle changes to be detected, which may otherwise have been swamped by noise and less powerful techniques.
we believe that the presented methodology offers an attractive alternative to other approaches and will be a useful addition to the toolbox of systems biologists.

 pattern recognition is central to many scientific disciplines and is often a first step in building a model that explains the data. in particular, the study of periodic phenomena and frequency detection has received much attention, leading to the well-established field of spectral analysis. biology is rich with periodic behaviour, with sustained oscillations in the form of limit cycles playing important roles in many diverse phenomena such as glycolytic metabolism, circadian rhythms, mitotic cycles, cardiac rhythms, hormonal cycles, population dynamics, epidemiological cycles, etc.. a conventional method for frequency detection is fourier analysis. it is based on the fact that it is possible to represent any integrable function as an infinite sum of sines and cosines. the fourier transform uses this property to reveal the underlying components that are present in a signal. fourier theory has given rise to a wide range of diverse developments and far-reaching applications, demonstrating the theory's undisputed importance and impact. for frequency detection, however, it is known that the ft works optimally only for uniformly sampled, long, stationary time series. furthermore, common procedures of pre-processing the data can cause problems. time series can contain low frequency background fluctuations or drift that are unrelated to the signal of interest. for the ft, it is then necessary to remove the trends using detrending techniques. as has been shown previously, this detrending leads to convolution of the signal that can both remove evidence for periodicity and add false patterns. another known problem is aliasing. if a signal containing high frequencies is recorded with a low sampling rate, peaks of high frequencies can fold back into the frequency spectrum, appearing as low frequencies. the gibbs phenomenon provides another example where spurious peaks appear in a ft. it occurs at points of discontinuity in a periodic function, and results in so-called ringing artefacts around the "true" frequency peak. as for the accuracy of the frequency estimate, no direct information of this is given by the output from a ft, since the sharpness of the peaks depends on a combination of factors such as noise levels and the length of the time series. for further details, see the extensive ft literature. wavelet transforms offer an attractive alternative to fourier transforms. the main difference is that they are localised in both the time and frequency domain. this property makes wavelets better adapted to problems with truncated data. wavelets have found wide-ranging applications and have proven to be particularly powerful for image processing and data compression. bayesian inference provides another approach for analysing data. it addresses additional aspects of the problem, such as the inherent uncertainty of the analysis and the effects of external noise. using this framework, the method of bayesian spectrum analysis was developed by bretthorst and applied to nuclear magnetic resonance signals and parameter estimation with great success. there are several advantages of the bayesian approach, including an inherent mechanism for estimating the accuracy of the result and all parameters, as well as the ability to compare different hypotheses directly. focus is shifted to the question of interest by integrating out other parameters and noise levels. initial knowledge of the system can be incorporated in the analysis and expressed in the prior probability distributions. there has been a recent flood of bayesian papers with some convincing applications and promising developments in systems biology. the bayesian approach to time series analysis has proven its value in fields such as nmr and ion cyclotron resonance analysis. in this paper, we describe the development, implementation and testing of bayesian model development coupled with bsa and nested sampling, in a biological context. we present a comparison of this approach with the ft, applied to a number of simulated test cases and two types of biological time series that present challenges to accurate frequency detection. we first present some necessary background, upon which we build to develop to our approach. bayesian inference data is rarely available in sufficient quantity and quality to allow for exact scientific deduction. instead we are forced to infer models from incomplete knowledge. bayesian inference is based on bayes' rule, which evaluates a hypothesis, h, in light of some data, d, and some prior information, i. it is a method of assigning probabilities based on the current state of knowledge, allowing for subsequent re-evaluation as new data becomes available. the goal is to determine p, the posterior probability distribution of the hypothesis, given the data and the prior information. with bayes' rule, the posterior can be expressed in terms of other probabilities as where p is the probability of observing the data given the hypothesis and the prior, p is the prior probability of the hypothesis, and p is the probability of the data given the prior. when the hypothesis is the variable and the data is held constant, p is called the likelihood function, and when the hypothesis is constant it is called the probability of obtaining a specified outcome. when evaluating only one hypothesis, p is a normalising constant, but when investigating more than one hypothesis this term plays a key role and is called the evidence. bayesian spectrum analysis our presentation in this section follows closely that of bretthorst. the aim is to infer the most probable frequency, ω, from the given data. we start by building a model for the observed data, parameterised by angular frequency, ω, and amplitudes, c, and then use bayes' rule to compute the posterior probability of the parameters, p. by assigning priors to the model parameters c and integrating over these, we arrive at the posterior probability for the parameter of interest, ω,. this is referred to as the marginal posterior probability of ω. we note that ω is an r-tuple, {ω ω..., ωr}, with as many elements as there are distinct frequencies in the data. a general model for observed data sampled at n discrete time points, d = {d,..., d}, includes the signal of interest, s, a possible background function, g, and the noise present in the system, e, the signal function will usually be unknown and may be complicated, but can be approximated by a linear combination of ms model functions, ψi, that we parameterise by the quantity of interest, ω: in which are the expansion coefficients. similarly, the background function, g, can be approximated by a set of mg functions, ζi, that are independent of ω, where are the background model function expansion coefficients. since a and b are not the main focus of the analysis, we will aim to integrate them out of the equations by marginalisation. parameters that are treated in this manner as often referred to as nuisance parameters, which we denote here by. although the signal function depends on ω, whereas the background function does not, for notational purposes we introduce the set of model functions, ϕi, which consists of both ψi and ζi. this allows us to condense the model equation into such that now d = f + e and m is now the total number of model functions, ms + mg. the model functions will typically not be orthogonal functions over the time series domain. this, however, can be achieved by cholesky decomposition. in all subsequent calculations an orthogonal basis is used. from bayes' rule, the joint probability distribution of the model for the parameters ω and c is the likelihood function, p, is calculated by comparing data produced by the model signal, equation, to real experimental data. if the model perfectly captures the signal, the difference between the model data and the real data is simply the noise in the system. the model of the data in equation includes noise, e, which we assume to be time independent in the further developments. the true noise level is unknown, but for a given noise power, σ the principle of maximum entropy leads to the use of a normal distribution, a noise model of this form ensures that the accuracy of the results is maximally conservative for a given noise power. we will later integrate over all possible noise levels to remove the dependence on σ. with the described signal and this noise model, the likelihood was calculated by bretthorst to be where n is the number of data points, and where is the mean-square of the data,, and Φjk is the matrix of the model functions,. the goal of the analysis is to compute the posterior probability for frequencies in the data, i.e. to go from the joint probability distribution to a posterior probability of ω, independent of the other parameters. by integrating over all possible values of the parameters σ and c, the remainder is the marginal posterior of the parameters of interest, ω={ω ω..., ωr}. this is an essential advantage of the bayesian framework, allowing the analysis to focus on estimating the parameters of interest, regardless of the values of the others. if necessary, the other parameters can be estimated at a later point. to integrate over the σ and c values, priors must first be assigned to them. we chose uniform priors for c and ω, representing complete lack of knowledge. we know that σ is continuous and must be positive, and in such cases a jeffreys prior is appropriate, p = 1/σ. both the uniform distribution over continuous variables and the jeffreys prior are known as improper priors if bounds are not specified as they cannot be normalised. for more information on prior assignment see. using the general model, equation, assigning the priors, calculating the likelihood function, equation, and integrating out the amplitudes and noise parameters, the posterior probability distribution of ω is proportional to where h is the projection of the data onto the orthonormal model functions,, and is the mean-square of the hj,,. this expression of the posterior allows us to identify the strongest frequencies present in the data. for a good model, there will be a high probability peak in the posterior distribution at that ω = {ω ω..., ωr}.  we employed the framework developed by jaynes and bretthorst to investigate the frequency components in a number of biological time series. model comparison after evaluating the probability of parameters in light of a certain hypothesis, it is important to question the validity of that hypothesis. thus, the next step in bayesian inference is to compare the probability of different hypotheses. the hypothesis is now a particular model of the signal, hi, out of a set of m models {h..., hm}, and using bayes' rule, the posterior probability of this model is then two different models, hi and hj, can be compared by taking their ratios, the probability of the data given our prior information, p, which was a normalisation constant in equation, will now vary between models, and is called the evidence. it evaluates the fit of the data to the model, whilst penalising models that include more parameters. each additional model parameter should be followed by a significant increase in probability, otherwise the simpler model is preferred. thus, bayesian model comparison naturally follows the principle of occam's razor. model development it will often not be obvious which function to choose to model trends in the data, so an approach using basis functions and expanding these to different orders will be of advantage, as in equation. each expansion represents a different model, hi, and these can be compared using inference techniques. likewise, different functions for capturing the signals in the data and modelling a different number of signals correspond to different models for data. following, we use the posterior ratio to evaluate different models. this model ratio can be used to determine the number of background model functions for each time series. the posterior probability ratio is calculated between model hn and hn+ where hn is a model including n background functions. to obtain the model ratio, priors are assigned to each of the models and their likelihood functions are calculated. assigning equal prior probability to all models reduces this to the ratio of evidences. to compute the evidences we need to integrate the likelihood, p from equation, over ω, σ, and c for each model hn. by assigning proper normalised priors to all model parameters it is possible to integrate over them around the maximum likelihood estimate. following bretthorst's derivation for location parameters, we assign gaussian priors to the amplitudes with hyperparameters for the variances. since the variances are scale parameters, they are subsequently assigned jeffreys' priors with an upper and lower bound. this allows us to normalise them and integrate, leaving the defined bounds as parameters in the final equation. for models with the same bounds these terms cancel out in the model ratio. the evidence for a given model, hn, was calculated by bretthorst to be where δ, γ and σ are the prior variances for amplitudes, frequencies and noise, respectively, rδ, rγ and rσ are the ratios of the integral bounds for these variances, is the mean-square projection of the data onto the orthonormal model functions at the maximum likelihood point for model hn, is the mean-square of the ω value that maximises the likelihood,, and r is the number of ω parameters, ω = {ω..., ωr}. the jacobian is obtained by orthogonalising the taylor-expansion of around the maximum-likelihood point,. see bretthorst for further details. for cases in which the number of frequencies in the data exceeds the dimension of omega, for instance multiple frequency data with a single frequency model, the above approximation for the evidence is ill-suited as the posterior will cease to be unimodal. for such scenarios, either multiple expansions or mcmc offer attractive solutions to marginalisation. for comparison we have included results from nested sampling as a means to perform the integration and compute the evidence. nested sampling is a variant of mcmc that employs a likelihood based sorting of sample points to efficiently guide the search strategy of the posterior distribution. when the model ratio, hn/hn+ becomes greater than the simpler model, hn, is favoured over hn+. adding more background functions than are justified by the data may lead to a lower probability for the frequency and in some cases possibly a location shift. this model development approach used for the background functions above can also be used to decide on the number of underlying frequencies in the data. the model ratios of a time series containing one frequency and a time series containing two are presented in table analysed with both a one- and two-frequency model. the results show, as expected, a preference for the one-frequency model in case a, and for the two-frequency model in case b. case a has only one frequency, nω = and case b has two, nω = the model ratio of a model including one frequency, h1ω, and a model with two, h2ω, is well above for case a, showing that the simpler model should be chosen. for case b, the first ratio is close to showing that the more complex model with two frequencies is far more likely. values below 10- are listed as the second ratio is above preferring the model with two frequencies to the one with three. we point out that the proposed method stops once the current best model has been found but is not guaranteed to find the global maximum from a predefined set of models. the procedure is thus part of model development rather than model selection. if the set of hypotheses are known in advance then the posterior ratios over the full set should be used to find the best model. testing we first show the power of the bsa approach on test cases using simulated data. in these tests, we sought to recover known input parameters from the simulated data, to validate the bsa approach. we employed sines and cosines as model functions ). for comparison, discrete fourier transforms were computed using fast fourier transforms. in the test cases, we varied key parameters such as noise levels, trace length, sampling intervals, amplitudes, frequencies, background trends and shape of oscillations. representative cases of noise levels and background trends are shown in table including fft results on the same data set. a key observation is that the bayesian approach extracts the correct answer from the data with high precision. bsa also computes the signal-to-noise ratios which is a useful indication of how much of the data cannot be accounted for in the model. furthermore, the amplitudes do not impact the bsa results since they are integrated out. each time series was generated with a sine function of angular frequency, ω, of rad/s with a level of noise in amplitude, ea, and phase, ep. in some time series a background trend was included, and in case number an additional sine function of rad/s is present. the resulting function was sampled times at an interval of s. results from fft are presented in the form of the angular frequency with the highest power,, and the estimated standard deviation, σfft. the bayesian frequency estimate at the maximum posterior point is denoted by and its standard deviation by σbsa. for comparison, the expectation value of ω and its standard deviation computed using bsa and nested sampling are denoted by and σbsa-ns. values of σ below 10- are listed as the estimated signal-to-noise ratio from the bayesian analysis is given in the last column. the bsa and bsa-ns approaches deliver the same results, apart from the case of multiple frequencies in a 1d search of ω, which for bsa-ns leads an intermediate estimate between the frequencies with a higher standard deviation. bsa has a clear advantage over fft when the data is non-uniformly sampled. fft requires uniform sampling, whilst bsa is less stringent and delivers the correct result with higher precision. bretthorst also noted that non-uniformly sampled data removes aliases from the frequency domain, another significant advantage. five further distinct cases emerged from the tests in which bsa delivers superior results to fft: time series which have background trends, few data points, high noise levels, multiple frequencies, and non-harmonic oscillations. background trends additional file figure s is an example of a time series with a strong background trend. in table the model ratios for different background functions are shown. the ratio is initially well below but the ratio of models with expansion orders of two and three legendre polynomials is above thus, background functions of legendre polynomials to expansion order two is more likely, and should be used in the estimation of ω. posterior probability ratios of models including a different a number of background functions, ζ, in the analysis of the time series in additional file figure s the first two ratios favour the more complex model. the ratio between models with two and three background functions, h2ζ /h3ζ, is above thus favouring the simpler model. the analysis would then normally automatically stop at two background functions, but for demonstration we include more functions here. the ratio stays above thereby building a chain of decisions that always prefer the simpler model. examples 8- in table also include trends, and without pre-processing fft cannot pick out the correct frequency. in contrast, bsa includes background functions in the model signal and delivers the desired result. including background functions, however, results in over-estimation of the signal-to-noise ratio. short time series additional file figure s shows the results from analysing a short time series. the fft power spectrum is very broad, which comes as no surprise given the fft dependence on the number of data points. bsa estimates the correct frequency sharply, but the maximum probability drops compared to longer time series. this demonstrates the higher uncertainty associated with fewer time points. high noise levels bsa is also successful at handling high levels of noise, as highlighted in examples 1- in table the frequency estimates are correctly reproduced by the fft. in these simple test scenarios, the bsa posterior probability distribution estimates the frequency with a significantly higher precision than the fft. whereas the estimated uncertainty of parameter expectation values is a built-in aspect of any probabilistic treatment such as bsa, fft has no inherent mechanism for assessing the accuracy of the results. the fft output is summarised by the average,, and the standard deviation, σfft, over the transformed data set. we show that different noise levels influence the σfft more than the σbsa. multiple frequencies example in table has two frequencies present in the data. both bsa and fft show these two peaks in the resulting plots. although bsa can be used in this manner with a one-dimensional ω to scan through frequency space and estimate the number of frequencies in the data and their location, if more than one frequency is present, the model should be extended to reflect this. without this extension the integration procedure around a single point is not well suited, so we employed nested sampling to compute the marginalisation in these cases. for the extension approach, when the posterior probability over ω = {ω1} reveals two strong frequencies, then a better model would be ω = {ω ω2}. for example, additional file figure s shows bsa and fft results for a test case that includes higher harmonics which give rise to multiple peaks in the log p plot. if more than one peak in the resulting posterior probability emerges, then the model can be extended further. one peak in the posterior probability over the number of modelled frequencies signifies that the correct number of frequencies has been captured. as another example, additional file figure s shows the result of a two-frequency search. the bsa posterior probability distribution is now two-dimensional with a peak at the two correct frequencies. the fft results are also shown. additional file 5a, figure s5a, shows a time series with a high noise level and two very close frequencies of and rad/s. fft cannot distinguish them and shows only one peak. bsa breaks the resolution and precision limitations inherent to fft by introducing a continuous probability distribution instead of the fixed number of points and can therefore sample the posterior more finely in areas of high probability. this approach gives rise to a high-resolution probability plot in which two distinct frequencies emerge. the peaks have a larger variance at this local level, but the qualitative information of two underlying frequencies is revealed. to develop bsa further, we used windowing of the time series to compute the posterior probability distribution of ω at each time point. we call this bsa local. the robustness and negligible peak broadening of bsa with fewer time points allows for this windowing to proceed without the introduction of artefacts due to truncation. this local bsa captures changes in frequency, as shown in figure 2b and figure 3b. the bsal was compared to short-time fourier transform, which is a windowed fourier transform, and to wavelets. for the wavelet power spectrum a morlet mother wavelet was used. the advantages of bsal are that it remains within the same bsa framework, has high accuracy, and does not require pre-processing of the data. non-harmonic oscillations bsa results for oscillations with a non-harmonic shape are superior to the fft. it highlights an essential difference in the two methods since biological data is often repetitive, but with a wide range of oscillatory patterns. to demonstrate this further, figure 4a shows a time series simulated from an ordinary differential equation model of cellular calcium signals. such time series presents two potential problems: the time series is chaotic and thus not perfectly periodic, and the signal shape is non-harmonic. the calculation of interspike intervals of the time series show that multiple intervals are present. the highest peak of the fft plot suggests that the entire time series is one period, while bsa suggests a strong angular frequency around rad/s. the bsa suggestion is similar to the second fft peak. this highlights the differences between frequencies in the data and spike intervals. isi are a common way of characterizing spike data, however, multiple isi need not correspond to multiple frequencies in the data. of the four strong isi shown here, both bsa and fft identify only one of these as a regular period. summary after extensive test cases we find that bsa delivers superior results in cases where the fft assumptions are too constraining, most notably in the five cases above. bsa is a flexible method allowing the underlying hypothesis to be changed depending on the focus of the analysis, and to directly compare the validity of different hypotheses. it can handle non-uniformly sampled data and has no need for pre-processing procedures. the price of these superior results comes at a computational cost that ranged from tens to hundreds of seconds for the examples shown here. calcium spiking data the first biological data set comes from intracellular signalling in plant-microbe interactions. symbiotic bacteria induce calcium oscillations, called ca2+ spiking, in legume root cells. these are non-stationary and often noisy time series, causing problems in identifying periodicity. one hypothesis for signal transduction in this system is via frequency encoding, so concluding whether there is underlying periodicity, and at what frequency, is of great interest. the ca2+ spiking has background trends present due to fluorescence bleaching and cell movements, which are assumed to be unrelated to the underlying signal in the cell. therefore, accounting for the background functions plays a key role in the analysis. example time series are shown in figure 5a. nine spiking cells from the model legume medicago truncatula were analysed for an underlying period. the data is obtained by microinjecting a root hair cell with the calcium indicator dyes oregon green and texas red, and exposing the plant to the bacterial signal molecule that induces the oscillations. the data is a ratio of the fluorescence from the two dyes, showing changes in ca2+ concentrations. the data has been published in. the fft of the ca2+ data results in a very broad periodogram, due to multiple frequencies and high noise levels. also, the spiking produces a non-harmonic signal, which might be another problem for the fft. for comparison, we also present results from the multitaper method. the mtm is a non-parametric method of spectral analysis that uses tapers to minimize the variance in the power estimate, and is targeted at short and noisy time series. the mtm results were very similar to the fft. these periodograms do not address the question of interest: is there a key period in the ca2+ signal? in the bsa analysis, the ca2+ spiking data used the legendre background functions to an expansion order of 1- depending on the individual trace. nested sampling was used to compute the evidences. frequencies with high probabilities were picked out, but varied in the interval of approximately 50- s. however, the strongest periods were in the interval of 75- s. if periodicity plays a role in the signal transduction of this system, then the key period should be in this interval. the signal-to-noise ratios were relatively high, between 100- possibly as a consequence of including several background functions. analysed calcium oscillations in m. truncatula root hair cells in response to bacterial symbionts. the bsa shows strong underlying signals in the data, but the cell-to-cell variability is high. bsa and the bsa-ns disagree when there are multiple peaks present, as can be seen by the large standard deviation, σ, of the bsa-ns average period estimate. circadian data the second biological data set shows gene expression of so-called clock genes. many processes in plants follow a circadian rhythm. a number of genes in arabidopsis thaliana have been shown to regulate circadian rhythms, and time series of rna levels show how these clock genes are expressed in cycles. time series with only a couple of cycles are common in biology and provide another suitable test case. for these circadian rhythms, we chose to analyse rt-pcr data from four clock genes in two genotypes of a. thaliana. the plants are either wild type, fri;flc, or mutants, fri;flc, of the genes fri and flc. the rna was extracted from seedlings, and each time series is an average of two biological replicates. an example of the rna levels of a clock gene is shown in figure 6a. the data has been published in. fft on the rna levels of these clock genes did not give any clear periods, either having only a vague peak or none at all. this is caused by the fft's dependence on the length of the time series, which in this case was only 1- cycles. the mtm method had more of a peak in the 20- h period, but still lacking in precision. bsa on the other hand provides a clear peak close to h, consistently for all eight time series. nested sampling was used to compute the evidences. the assigned probabilities are relatively low, but the signal-to-noise ratios were between 2- and similar probabilities were obtained using simulated data with few data points and high noise levels. the period peaks are very stable over all the time series, and suggests a probable period which is unaffected by the mutations. this is in agreement with the original conclusions of the experiment. the bsa results of rna levels of four so-called clock-genes in a. thaliana, measured in two different genotypes. the plants are either wild type, fri;flc, or mutants, fri;flc, of the genes fri and flc. all eight genotypes displayed oscillations in rna levels with a period of approximately hours. both bsa methods identify similar periods.  bayesian inference offers a powerful way of analysing biological time series. despite the undisputed value of fourier theory, there are cases when the necessary requirements for its optimality for time series analysis are not met. this is a consequence of the underlying assumptions of a fourier transform, causing it to work optimally only for uniformly sampled, long, stationary, harmonic signals that have either no or white noise. in biology these requirements are rarely fulfilled, requiring pre-processing of the data, such as noise reduction and detrending techniques, with the risk of convoluting the signal and losing valuable information. by placing the problem of frequency extraction in the framework of bayesian inference, the known and well-documented problems of fourier analysis can be overcome. this approach also breaks the resolution and precision limitations inherent to the fft by introducing a continuous probability distribution instead of the fixed number of points maintained by the discrete fourier transform. as we demonstrated here, bsa coupled with automated model development can give superior results to the fft when faced with short, noisy time series, non-stationarity and non-harmonic signals. the suggested automated model development worked well in our hands but must be used with caution in practice as the approach is not guaranteed to find a global optimum in model space. alternate models should be explored and compared using posterior probability ratios or approximations thereof. we found nested sampling to provide a powerful means of estimating evidences for cases in which a single peak could not be identified. other mcmc techniques such as simulated annealing running in parameter exploration mode or standard metropolis-hastings algorithms offer attractive alternatives. bsa calculates signal-to-noise ratios, provides parameter precision estimates, and can handle high noise levels as well as background trends and therefore has no need for pre-processing. more importantly, the bayesian framework offers flexibility in the underlying model and enables direct comparison of hypotheses. the work presented here is a merely a first step in this direction. we have employed conservative priors that make an analytical treatment tractable but in some cases more information could warrant a different choice of prior that might require substantial alternations to our approach to handle the numerics of marginalisation. there are many known examples in biology in which oscillations play a key role and methods for their detection will be of value, especially in cases where subtle differences are of importance and for short, noisy time series. in the presented examples, we demonstrated the improvements that can be gained from employing this approach. although in these cases, the biological conclusions would not have changed, one can envision scenarios in which a higher accuracy in frequency detection may allow subtle changes to be detected, which may otherwise have been swamped by noise and less powerful techniques. we believe that the presented methodology offers an attractive alternative to other approaches and will be a useful addition to the toolbox of systems biologists. 
mathematical relationships between the time and frequency domain, such as fourier transforms and wavelets, are commonly used to extract information about the underlying signal from a given time series. fourier transforms are a convenient tool for analysing the frequency domain of time series. the results show that the bayesian frequency detection approach can provide useful results in specific areas where fourier analysis can be uninformative or misleading. we show the value of placing spectral analysis in the framework of bayesian inference and demonstrate how model comparison can automate this procedure. we demonstrate further benefits of the bayesian approach for time series analysis, such as direct comparison of different hypotheses, inherent estimation of noise levels and parameter precision, and a flexible framework for modelling the data without pre-processing. in this paper we build on the method of bayesian spectrum analysis and demonstrate its advantages over conventional methods by applying it to a number of test cases, including two types of biological time series. however, for truncated, noisy time series with background trends this unique mapping breaks down and the question reduces to an inference problem of identifying the most probable frequencies. in this paper we present an alternative approach based on bayesian inference. modelling in systems biology often builds on the study of time-dependent phenomena. secondly, circadian rhythms in gene expression measured over only two cycles highlights the problem of time series with limited length. this one-to-one mapping from time points to frequencies inherently assumes that both domains contain the complete knowledge of the system. however, there are well-known limitations of this method, such as the introduction of spurious frequencies when handling short and noisy time series, and the requirement for uniformly sampled data. biological time series often deviate significantly from the requirements of optimality for fourier transformation.
 a first step in building a mathematical model of a biological system is often the analysis of the temporal behaviour of key quantities. mathematical relationships between the time and frequency domain, such as fourier transforms and wavelets, are commonly used to extract information about the underlying signal from a given time series. this one-to-one mapping from time points to frequencies inherently assumes that both domains contain the complete knowledge of the system. however, for truncated, noisy time series with background trends this unique mapping breaks down and the question reduces to an inference problem of identifying the most probable frequencies.  in this paper we build on the method of bayesian spectrum analysis and demonstrate its advantages over conventional methods by applying it to a number of test cases, including two types of biological time series. firstly, oscillations of calcium in plant root cells in response to microbial symbionts are non-stationary and noisy, posing challenges to data analysis. secondly, circadian rhythms in gene expression measured over only two cycles highlights the problem of time series with limited length. the results show that the bayesian frequency detection approach can provide useful results in specific areas where fourier analysis can be uninformative or misleading. we demonstrate further benefits of the bayesian approach for time series analysis, such as direct comparison of different hypotheses, inherent estimation of noise levels and parameter precision, and a flexible framework for modelling the data without pre-processing.  modelling in systems biology often builds on the study of time-dependent phenomena. fourier transforms are a convenient tool for analysing the frequency domain of time series. however, there are well-known limitations of this method, such as the introduction of spurious frequencies when handling short and noisy time series, and the requirement for uniformly sampled data. biological time series often deviate significantly from the requirements of optimality for fourier transformation. in this paper we present an alternative approach based on bayesian inference. we show the value of placing spectral analysis in the framework of bayesian inference and demonstrate how model comparison can automate this procedure. 