Appropriate and timely changes in gene expression are essential for cell life. The transcriptome is finely regulated by a large number of molecular mechanisms able to adjust the balance between mRNA production and degradation. Every aspect of transcript life is subject to elaborate control but, traditionally, the focus of the research has been on transcriptional regulation 
. However, whereas mRNA abundance results from the dynamic interplay between transcription and degradation, the speed by which cells can adjust their mRNA levels is critically dependent on the rate of mRNA turnover 
. As a result, small changes in mRNA stability may dramatically drive rapid variations of transcript abundance. Efforts to understand the underlying principles of mRNA decay and transcription co-ordination are very important since the balance between transcription and decay influences most, if not all, the cell responses to endogenous and exogenous signals 
The current widespread interest in this topic has been fostered by the finding of specific regulatory mechanisms of mRNA stability such as, for example, RNA binding proteins 
and small RNAs 
. Regulation of transcript stability cannot be considered a simple “disposal system” but a sophisticated tool for the proper orchestration of the global cell response to internal and external stimuli 
. Remarkably, a key role of mRNA stability has been reported in cancer, inflammatory diseases and Alzheimer's 
. In recent years there has been a surge in empirical studies that measured, on a genome-wide scale in a variety of environmental conditions, messenger half-lives of many organisms, including plants 
and fungi 
. The discovery of such new regulatory layer has clarified that, in order to obtain a clear picture of the underlying regulatory machinery, it is necessary to complement the traditional time-course experiment measuring the cell transcriptional response under certain conditions far from steady state with decay rates data under the same
Experimental procedures for the evaluation of mRNA decay rates are based on measuring gene expression upon inhibition of transcription 
or on pulse-chase RNA labeling protocols 
. Such protocols are very critical (see for a comparison among different studies summarized in ), since, for instance, transcriptional shut-off blocks growth and has a profound effect on cellular physiology, as well as on mRNA metabolism 
. In fact, Wang et al.
and Grigull et al.
datasets show a low value of the Pearson correlation (
), and no correlation at all can be found (
) between Munchel et al.
and Wang et al.
datasets (see and respectively). Despite the same experimental conditions (asynchronous growth), the two half-life independent measurements obtained by Wang et al.
and Munchel et al.
are uncorrelated, probably due to differences in the shut off protocol (pulse chase for 
and thermal inactivation for 
), whereas Grigull et al.
and Wang et al.
appears significantly correlated, probably due to the same shut off protocol used (thermal inactivation).
Basic comparison statistics among the yeast S. cerevisiae mRNA half-lives during asynchronous growth measured by three independent laboratories.
Half-lives experimental measurements.
It has been shown that genes having the same biological function 
are likely to share similar half-life values. Consistently, by averaging using functional groups, we found an increase in correlation between Wang et al.
and Grigull et al.
, see Supplementary Figure S1A
), and still no correlation between Munchel et al.
and Wang et al.
, see Supplementary Figure S1B
Here, we developed a stochastic computational model of the expression kinetics to identify condition-specific mRNA stabilities which makes use only of experimental mRNA time profiles. We also assumed that degradation rates are gene-specific but approximately constant over the experiment time course. Predictions of our algorithm, termed DRAGON (Decay RAtes from Gene expressiON), were validated on experimental mRNA abundance 
and turnover 
data, both collected during the Intraerythrocytic Developmental Cycle (IDC) of Plasmodium falciparum
. The estimations were in line with the experimental measurements. Remarkably, the DRAGON estimated half-lives were consistent with the finding of a peculiar pattern of mean half-life values along the wave of sequentially induced genes in subsequent stages of P. falciparum
development. We also applied our methodology to public time-series datasets for which half-lives data, under the same experimental conditions, have not been experimentally measured. In particular we focused on budding yeast reproductive 
and metabolic cycle data 
. In fact, for the yeast Saccharomyces cerevisiae
, only half-life data under asynchronous growth are publicly available 
. Our study showed the presence of the same periodic pattern of mean half-life values in all datasets, thus suggesting that such behavior may be a general feature, not limited to the Plasmodium falciparum
mRNA kinetics and half-life
Experimental evidence suggest that the majority of mRNAs are degraded with a first-order decay rate 
. This allows to characterize mRNA disappearance time profiles by a first-order rate equation
is the decay rate (or half-life
is the mRNA concentration and
is the promoter activity (the rate of production of new mRNAs). It is worth noting that, the degradation rate
cannot be estimated from the concentration time profile
for a single gene, since the term
is not usually available in the typical time-course microarray experiment. The measurement of the promoter activity time profile would require additional experiments (such as those described in 
) but, in this paper, we will assume that only mRNA abundance time-series data are available. At steady-state
are constant so that
and, consequently, we obtain
From the above equation, it is clear that at steady-state an increase (decrease) in mRNA concentration can be produced either by an increase (decrease) of transcription or by a decreased (increased) value of the decay rate: the two regulatory strategies have therefore an equivalent outcome. As a result, from steady-states measurements, it is hopeless to reveal the relative contribution of transcription and degradation and, most importantly, their co-ordinated activity as well. By contrast, the whole kinetics of induction and relaxation, as measured by time-courses experiments, depends on the degradation and production rate in different ways: increasing (decreasing) the production rate results in a proportionally increased (decreased) mRNA abundance, whereas the rise time (i.e.
the time required for the response to rise from
of its final value) is not affected. Increasing the decay rate results in a faster rise time both in the induction and relaxation phases, whereas a decrease results in slower rise time 
. This key point is illustrated in and in Supplementary Figures S2A and S2B
Another important consequence of half-life specificity is the regulation of the timing of gene induction, as pointed out by Elkon et al.
. In fact, an expression wave, i.e.
the sequential activation of genes, is usually interpreted as resulting from the corresponding activation of a multi-step transcription factors cascade (as illustrated in ). Whereas such mechanism is certainly very important, there is also an alternative way to obtain an expression wave by means of a “stability gradient”. As illustrated in , a single transcription factor may initiate transcription of a set of target genes and their peak of induction can be modulated by a stability “gradient”, i.e.
by specifically adjusted decay rates. More precisely, early induced genes would have short half-lives and late responding genes would have long half-lives. Clearly, both mechanism may well act in cells, thus generating a wide spectrum of responses.
Time-courses are a very common design for microarray analysis, which allows researchers to follow the dynamics of the cellular response to perturbations 
. Such data are available for a very large number of experimental conditions and organisms: only the Stanford Microarray Database includes to date 1545 time course data sets. Among the examples later illustrated in the paper, it is worth mentioning the genome-wide gene expression time-series obtained during the reproductive cell cycle 
, the metabolic cycle 
and the P. falciparum
. The time-series datasets used in this paper are summarized in .
Gene expression time-series experimental measurements.
DRAGON–an algorithm for half-life estimation
The goal of the DRAGON methodology is to derive a robust estimate of each mRNA species half-life starting from all available gene expression pairs. The rationale for the algorithm mainly draws on properties of pairs exhibiting a certain degree of common promoter activity (as in 
). Besides, DRAGON infers common promoter activity using a statistical model that simulates both gene-specific and common effects.
The rate of change of mRNA concentration for a generic pair of genes, say gene
where the symbols
represent the mRNA time profiles of the gene pair
are the promoters activity, and
are the degradation rate of mRNA of gene
, respectively. The terms
are not known since we considered the case in which only mRNA abundance is measured.
We modeled promoter activities as the sum of two terms, the first one common to the pair and the other one specific for each gene:
is the common part, scaled by constants
are gene specific independent stochastic processes with zero mean, that is
. Equations (3)
encompass the case of:
- , fully correlated (correlation ) for which and
- , partially correlated (correlation ) for which and
- , un-correlated (correlation ) for which either or .
can therefore be written for all available gene pairs; thus, for a set of
genes, we have
pairs to analyze. For each gene pair DRAGON provides an estimate of the time profile of
, of all the parameters,
, and the covariance matrix of the stochastic processes. For each gene we therefore have
estimates of the decay rate
, one for each pair containing that gene.
Notice that equations (3)
yield a couple of linear stochastic differential equations. Since measurements of mRNA concentrations are available only at given time points, it is necessary to transform (3)–(4) in a couple of discrete stochastic equations. The exact discretization of (3)–(4) is possible since they are linear 
. The Kalman filter
is used on the resulting discrete equations and a maximum likelihood algorithm is exploited to generate the best possible estimate of the parameters.
A complete description of the mathematical model and of the discretization and parameters estimation procedure is given in the paragraph Stochastic modeling of expression kinetics and Kalman filtering
of the Materials and Methods