- Research article
- Open Access
mRNA stability and the unfolding of gene expression in the long-period yeast metabolic cycle
BMC Systems Biologyvolume 3, Article number: 18 (2009)
In yeast, genome-wide periodic patterns associated with energy-metabolic oscillations have been shown recently for both short (approx. 40 min) and long (approx. 300 min) periods.
The dynamical regulation due to mRNA stability is found to be an important aspect of the genome-wide coordination of the long-period yeast metabolic cycle. It is shown that for periodic genes, arranged in classes according either to expression profile or to function, the pulses of mRNA abundance have phase and width which are directly proportional to the corresponding turnover rates.
The cascade of events occurring during the yeast metabolic cycle (and their correlation with mRNA turnover) reflects to a large extent the gene expression program observable in other dynamical contexts such as the response to stresses/stimuli.
Ultradian self-sustaining energy-metabolic oscillations arising spontaneously in high density Saccharomyces cerevisiae continuous cultures exposed to glucose-limited growth have been known and studied for decades [1, 2], and have more recently been observed to induce genome-wide periodic patterns in different series of microarray experiments [3, 4], although with widely different periodicities, ~40 min for  and ~300 min for .
Many studies aim at understanding the mechanisms inducing these sustained oscillations and the rigorous temporal compartmentalization they induce, see [5, 6] for surveys. Suggested causes range from a single critical pathway (like the feedback effect of cysteine on the sulfur assimilation pathway ) to the alternation of aerobic and anaerobic respiratory modes (as deduced by the fluctuations in the concentration of dissolved O2 and of other observed metabolites ), from the interaction with cell cycle [8, 9] to the mutual incompatibility of different redox biochemical processes [10, 11].
The scope of this work is to emphasize a different aspect, intrinsically dynamical and post-transcriptional, which is likely to play an important role in the coordination of the "slower" yeast metabolic cycle (YMC) of , namely mRNA stability. We will show that there is a roughly linear relationship between the average half life (HL) of the transcripts, clustered according to expression or function, and the phase at which their concentration peaks in the cycle. More generally, there seems to be a strong correlation between HL and the shape of the pulses of gene expression: genes with short HL have short and sharp (almost impulsive in the time scale considered) pulses, while genes with long HL have pulses that are not only delayed but also broader and with more gentle slopes.
In recent years, post-transcriptional control is being recognized as an important aspect of gene regulation, especially in eukaryotic DNA, which lacks operonal structure [12–15]. It can occur in many guises, through mRNA turnover [16–20], or through "RNA regulons"  i.e., groups of genes coordinately guided in the RNA processing, localization and protein synthesis by RNA-binding proteins (RBPs) [22, 23], or even through the mediation of a metabolic substrate (typically a nutrient [24–26] or an enzyme ). Our result confirms the importance of post-transcriptional control, and points at mRNA turnover as a regulatory mechanism at a genome-wide level. Its peculiarity consists in putting the time axis into the picture in an intrinsically dynamical way. Consequently, in order to be observed, it requires times series sampled at a sufficiently high frequency and dynamics in the right time window, a combination seldom occurring in current expression profiling datasets. So for example the correlation between HL and phase/shape of the oscillations cannot be observed in the much faster YMC of , where HL and the period are of comparable duration, hence the system has no time to decay before the arrival of the next wavefront.
In order to emphasize the dynamical aspects, we shall treat the YMC as the time response of a genome-wide dynamical system to a sequence of impulsive "inputs" of transcription activation. We will show that grouping genes in terms of progressively delayed and broadened responses to a sequence of "input pulses" of transcriptional activation allows to see in a remarkably fine detail the causal chain of events constituting the transcriptional program of the cell. The few ambiguities resulting from this classification can be interpreted in terms of some other annotation, typically compartmental localization.
In the following we shall proceed in two complementary ways: first the YMC time series are clustered in a completely unsupervised manner, only according to gene expression. The linear relationship between pulse phase (also pulse width) and HL then emerges in a straightforward way. Next, we consider families of genes whose products share some common annotation, for example genes on the same pathway or genes that are subunits of the same protein complex, and look at the type of time series they produce and at their "position" along the YMC.
Both approaches confirm that the YMC represents an organized cascade of events, in response to precisely equispaced bursts of transcriptional activation, with the temporal order reflecting the transcript turnover rate. Extrapolating from the specific YMC context, this cascade of events is observable to a good extent also in other gene expression time series (such as the response to a pulse of nutrient of , or the stress responses of ), suggesting it might reflect a prototypical dynamical mode of action of transcriptional response.
Results and discussion
The ~2000 genes labeled as periodic by a periodogram test are subdivided into 16 clusters, see Fig. 1. In Fig. 1(a) the clusters are sorted in increasing order of HL (computed as the average of the HLs of the cluster elements). It is immediately evident that the typical profiles, both in terms of the phase of the peaks (for each gene the phase is computed maximizing the correlation with respect to a train of shifted sinusoids) and of their width (although in a less regular way) is modified in an almost continuous manner as we move along the clusters figures. Notice in particular how the peaks of the first clusters match the "valleys" of the last ones. For the average phase on each cluster, the phase/HL relationship is almost linear (Fig. 1(b)). The scatter plot in (d) confirms this linear proportionality, but also shows a growing variance along the HL axis (see Table 1 for details). The deviations from linearity of clusters 6 and 9 admit a reasonable explanation, mostly in terms of compartmental localization. Cluster 6 is essentially composed of retrotransposons (all Ty1 and Ty2) and long term repeat mRNAs (mostly of δ type) for a total of 73 out of 102 genes. For most of these genes (59) an HL measure is missing. Hence the average HL for this cluster (and this cluster alone) may be biased or unreliable. Cluster 9 instead is almost entirely composed of cytoplasmic ribosomal subunits (109 out of 151 genes). In between, Clusters 7 and 8 contain to a large extent genes with mitochondrial localization and/or function (mitochondria organization and biogenesis, protein import into mitochondrial matrix, oxidoreductase activity for Cluster 7, mitochondrial ribosomes, envelope and membranes for Cluster 8). As is explained in detail in the next paragraph, the large deviation from linearity seen in Cluster 9 can be due to an extremely fast and short lived response of the mRNAs deputed to the biosynthesis of the cytoplasmic ribosomal complexes, not deducible from the available HL data, neither from the current literature (in  it is affirmed that cytoplasmic ribosomal genes tend to be stabilized by nutrient uptake). Although less precise, also the relation between HL and pulse width on each cluster (Fig. 1(c)) is approximately linear. Unlike the phase/HL proportionality, this last result is expected from simple dynamical considerations, as longer HL means longer "kernel width", see also the dynamical model explanation below. The emergence of a linear relation between HL and phase once the genes are arranged in classes according to profile similarity suggests that a corresponding cascade of causally organized events may be taking place during the YMC. To some extent this is already visible through an ontological analysis of the clusters of Fig. 1 (see Table 1), but in order to investigate more in detail the biological meaning and significance of such a genomic "assembly line" we computed HLs, phases and pulse widths along the main yeast pathways and for some of the annotated yeast protein complexes. The data for the pathways (see Fig. S2 in Additional file 1) are then lumped together into the 15 functional macrocategories shown in Fig. 2. In terms of these macrocategories (sorted by phase), the result is that the mRNAs activation reflects tightly the gene expression program expected to take place in the cell, especially for the "fast" categories, i.e., transcription, nucleotide metabolism and translation starting essentially synchronously in the time scale of the YMC, followed by DNA replication and repair and amino acid metabolism. Progressing further toward the slow processes, one encounters the metabolism of energy, carbohydrates and lipids. Also for this classification, the progression in terms of phase along the cycle is substantially faithful to the increase in HL (in the top plot of Fig. 2(b) the most significant outlier is still the category "translation" already mentioned, see also Fig. 3), and the progression in phase is paralleled by an increase in pulse width (see bottom plot of Fig. 2(b)).
A detailed functional analysis
Using the ordering by phase of pathways and protein complexes (see Fig. S2 and S5 in Additional file 1), we can zoom on these categories in much more detail. The first phase of this cascade consists of the activation of the transcription machinery with the synchronous bursts of transcription of the three RNA polymerases (see Fig. S1 in Additional file 1) and of most of the RNA processing components, like the tRNA processing complexes (RNase P) and rRNA processing complexes (exosome, RNase MRP, SIK1, NOP1), with the nuclear splicing complexes following closely. While the mRNAs for the polymerases are highly coordinated, the same cannot be said for the basal transcription factors (TFs) required for their initiation. Overall only a few of these genes follow the bursting trend of the RNA polymerases, notably, among them, SPT15, which forms the TATA-binding protein and is also a component of the polymerase I core factor and of TFIIIB. Most other genes involved with these general TFs do not show any periodic pattern, and their mRNA concentrations never surpass very low levels.
From Fig. 2, the peak of mRNA concentrations associated with the category "translation" seems to be synchronous with the RNA processing burst. However, a more careful analysis reveals that this phase is an average of two "compartmentalized" activations of the translation machinery, having fairly different phases: while cytoplasmic translation follows almost simultaneously the RNA machinery, the mitochondrial translation activation has a phase lag of more than one sixth of the period. In terms of time delay, this amounts to approximately 50 min, see Fig. 3. More in detail, most of the mRNAs of ribosomal small and large subunits for both cytoplasmic and mitochondrial localizations are highly correlated within their complex (average Pearson correlation for both is around 0.8) and correlated with the translation complexes at the corresponding location. In particular, among the cytoplasmic translation complexes, the initiation factors eIF and the termination factors eRF are very coordinated and respond very fast, while of the three elongation factors only eEF2 and eEF3 are well-coordinated, whereas the larger complex eEF1 shows a less-defined response pattern, with only the subunit eEF1-β clearly expressed. Overall for the class of translation complexes the pattern of activation of the response reflects closely the corresponding HL distributions  (eIF and eRF have short HL, eEF has not). Notice that a simple comparison of the HLs of the cytoplasmic and mitochondrial ribosomal and translation machineries (both approximately 24 min) does not show the significant difference which can be seen on the time series profiles and which is instead revealed by the phase delay analysis. For cytoplasmic ribosomal biogenesis, a similar anomaly is encountered also in the stress/stimuli responses analyzed below. For mitochondria, the same type of pattern is verified also by other complexes, for example by both the translocases located in the outer and inner mitochondrial membranes (TOM and TIM) which are known to mediate the protein import into the mitochondria, see Fig. 3.
A neat organization can be seen also in the phase of the nucleotide and amino acid metabolism: while pyrimidine and purine synthesis, as well as e.g. the CTP synthase enzyme involved in pyrimidine biosynthesis, are synchronous with the burst of transcription, the peaks for most of the enzymes involved in amino acid pathways tend to be in phase with the activation of the translational machinery. Also the synthesis of aminoacyl-tRNAs, necessary for the delivery of the amino acids to the ribosomes during translation has a similar phase. As expected, the "synthesis" pathway of an amino acid always anticipates its "degradation" pathway (see Fig. S2 in Additional file 1). In order to start translation, the initiator tRNA carrying methionine is required, and in fact, among the amino acid metabolic pathways, methionine is one of the fastest. As a matter of fact, the pathways of sulfur metabolism and of the sulfur-related amino acids (methionine, cysteine, as well as the closely related selenoamino acid metabolic pathway) present very similar and very compact time series (see Fig. S3 in Additional file 1), with an early (synchronous with the main burst) but long lasting activation (duration of the pulse is more than 100 min). This tight coordination may hint at a special role played by the sulfur pathways in the yeast population synchronization [31, 32].
To conclude the protein synthesis, the nascent polypeptide chains must fold into 3D structures. The molecular chaperonin-containing T-complex and the Gim complex, which help in the folding, behave synchronously with the main burst. On the contrary, ubiquitin and proteasome, that proceed to the recognition and degradation of anomalous proteins, as well as the SCF and anaphase promoting complexes, that cause the proteolysis of the cyclin-CDK complexes, have patterns of activation which are more delayed and broadened. Actually, this class of proteolytic processes (macrocategory "folding, sorting and degradation" in Fig. 2) has the highest values of phase i.e., it has the slowest response to the transcription bursts.
The macrocategory "DNA replication and repair" (see Fig. 2) contains what remains of the "fast" responses to a large extent synchronous (protein complexes: DNA damage checkpoint, DNA repair, pre-replication, replication, replication fork, which includes all DNA polymerases, helicases and ligases, cyclin-CDK) or within a short time delay from the initial bursts of transcription. The peculiarity of this class is that the pulses are more long lived than in the "transcription" and (cytoplasmic) "translation" categories. Also the complexes regulating the cohesion and separation of sister chromatids during the S-phase (nuclear cohesion family of complexes) follow the same pattern (see Fig. S5 in Additional file 1).
Moving to the core of the cell's metabolic activity, the average phase increases further (see Fig. 2), but the main qualitative difference is on the shape of the pulses, which are now broader and often with an asymmetric rise/decay profile: still sufficiently fast activation but slower and less abrupt decay. This difference is likely to reflect the longer HL associated to these categories (all have average HL ≥ 30 min), and implies metabolic functions more overlapping than sequential. Along each metabolic pathway, the degree of correlation among enzymes catalyzing neighboring reactions is higher than it is expected (the "expected value" is inferred from a large collection of yeast microarray experiments, see Fig. S4 in Additional file 1) implying a coherent and coordinated temporal behavior along the metabolic routes. Especially for mitochondrially localized pathways such as citric acid cycle and oxidative phosphorylation the pulses are very broad, with a neat downregulation only in correspondence of the bursts of transcription and an overall profile often exhibiting a double peak on each period (occurring with a phase lag of ~100° one from the other, see Fig. S8 in Additional file 1). The four respiratory chain complexes for example follow this pattern in a fairly precise manner. As shown in Additional file 1, this double peak characteristic is often associable with pairs of genes whose products are isoenzymes oscillating in antiphase, especially for enzymes involved in oxidoreductive processes (e.g. along the pentose phosphate pathway).
Regulation via TFs versus RBPs
In terms of regulatory influence, while the importance of transcription initiation via TFs is widely studied and a large amount of data (computational and experimental) is available about the binding of TFs to target genes, similar post-transcriptional systematic data on the regulation by means of RBPs are still sporadic . Notable examples are mRNAs associable to the nuclear export proteins Mex67 and Yra1 , the Puf family of RBPs , and the 3' UTR motif collection of . Inspired by , we applied these RBP lists as well as the list of TF binding sites from [36, 37] to the YMC time series comparing the average correlation among genes being common targets of a TF or of a RBP. The two distributions are shown in Fig. 4. For both TFs and RBPs, only a few motifs emerge as having a significantly high correlation. The number of genes regulated by the same TF varies between 1 and 226 with a mean of 35.2, while the number of genes with a common target mRNA motif varies between 6 and 1138 with a mean of 81.7. If we draw from a null distribution representing random grouping, increasing the number of genes in a group the probability of finding a high mean correlation obviously decreases, so we expect the distribution for the second set to be tighter around 0. In our case, on the contrary, there are 6 groups out of 110 with a mean correlation ≥ 0.4 for the TF target genes (versus an expected value of 1 for random groups of genes with the same cardinalities of these groups) and 7 groups out of 83 for the genes with a target mRNA motif (versus the expected 0 for random groups with same cardinalities). This suggests that post-transcriptional regulation is more significant that transcriptional regulation in the coordination of the metabolic cycle, although the evidence is not conclusive. When checking the groups of periodic genes with high correlation we found the following significant annotations:
• 44 genes out of 56 having Fhl1p as TF and 10 genes out of 12 having Sfp1p as TF are constituents of cytoplasmic ribosomes; notice that instead other cytoplasmic ribosomal TFs such as Rap1p do not correspond to a sufficiently high correlation;
• 22 genes out of 26 having Hap4p as TF code for subunits of respiration chain complexes;
• 62 out of 220 genes whose mRNA is bound by Puf3p are annotated for mitochondrial transcription/translation (56 are part of mitochondrial ribosomes, of which 47 are periodic), see Fig. 3.
Dynamical features of the unfolding cycle
Possible origins of the sustained oscillations are discussed at length in the literature [3, 5–8, 10, 11, 38]. Also Tu et al. explain the cycle and its time compartmentalization in terms of metabolism and redox balance [4, 32, 39]. Rather than adding to the list of mechanisms for metabolic regulation, by viewing each cycle as the dynamical response to a burst of transcriptional activation, this work aims at providing a characterization of the dynamics of the unfolding of the cycle, i.e., of how these "impulse responses" are progressively delayed and broadened with respect to the input pulses, and of how this correlates with the stability of the corresponding transcripts. The compactness in terms of phase and width of the early categories over repeated oscillatory cycles is an argument in favor of the existence of a single triggering event for each cycle, corresponding to the transcriptional activation bursts mentioned above. In fact, sharp, equispaced pulses are maintained in spite of the broader and less coordinated profiles of the events immediately preceding them. This hypothesis is not in contradiction with the observations about the metabolic origin of the YMC, neither with the observed alterations of the period following a genetic disruption [8, 32, 39] (which could in principle preserve the sequence of events described). On the contrary, it merges the metabolic control level described in  with an extra regulatory element which is known to play a role in dynamical contexts. In fact, the mRNA stability reflects known properties of the corresponding gene products: while mRNAs encoding transcriptional machinery or regulatory components tend to be short-lived and to turn over more quickly, transcripts encoding core enzymatic proteins are typically more stable [15, 19, 20]. For what is known, protein synthesis tends to follow the concentration of the corresponding mRNA  and to be at least as stable, if not longer-lived [41, 42]. Hence, it is expected that the concentration of the gene products follows profiles that are similar to those of the mRNAs. The observation that the dynamics through a metabolic pathway can be considered as a timed and sequential process at the level of gene expression appears in several papers in the literature, see [43, 44]. The same principle seems to be reflected in the YMC, although it is not observable at the level of detail investigated e.g. in , but more macroscopically and at genome-wide level.
An input-output dynamical model
In terms of dynamical models, the progressive broadening and smoothing of the response to a sequence of (transcriptional) pulses can be described by means of simple linear input-output models (i.e., transfer functions in the Laplace domain) of increasing order having "low-pass" characteristics. As the time constant of this low-pass filter is essentially given by the HL of the mRNA, this type of model naturally predicts the correlation HL-pulse width. In order to describe correctly also the phase along the cycle, a time delay is added to the response, see Methods for a thorough description and Fig. 5(a, b) for an example. If the order of such a fitted minimal dynamical model is used to sort the annotated macrocategories of Fig. 2, we still recover both the same expected cascade of events and the same direct proportionality with HL, see Fig. 5, meaning that even in terms of the simplest possible dynamical model the kernels providing the best fitting become increasingly complex as we progress through the cycle. This is of course expected as the mRNAs gradually pass from fast turnover to high stability.
A common dynamical gene expression program
As the YMC is obtained only in particular conditions (long-term continuous cultures in chemostats), an intriguing question is whether this highly organized unfolding of the dynamical response to pulses of transcriptional activation is peculiar only of the YMC or can be observed also in other experimental conditions. For this purpose, we consider the gene expression response of steady-state yeast to pulses of glucose described in . In this case, the yeast shows a transient dynamical response but no oscillatory behavior. Furthermore, the transient peaks are more or less synchronous for all genes, i.e., there is no time-ordering in the dynamics, unlike in the YMC.
However, if for a gene we compare the maximal signed amplitude of each expression profile on these time series with the corresponding phase and pulse width in the YMC, a sizable anticorrelation emerges, see Fig. 6(a). If, on the contrary, we consider the stress responses time series of , the YMC phase/pulse width turn out to be positively correlated (rather than anticorrelated) with amplitude, i.e., categories appearing early in the YMC tend to be downregulated in most stress responses, while "late phase" categories tend to be upregulated, see Fig. 6(b). It is known that in the stress responses genes annotated for ribosomal proteins and/or RNA metabolism are in general downregulated, while e.g. respiratory genes (such as those of the citric acid cycle and of the oxidative phosphorylation) become upregulated . On Fig. 6, notice that also in all these responses cytoplasmic ribosomes (cluster 9 in Fig. 1) are aligned with the rest of the (cytoplasmic) transcriptional/translational machinery rather than with the assigned HL values.
The conclusion of this analysis is therefore that in intrinsically dynamical contexts some form of common response might indeed be taking place, although exerted by different means. Such genome-wide coordinated response shows a graded ordering which reflects the degree of stability of the genes involved.
In [4, 39] the time compartmentalization of the cycle is interpreted in terms of the need to accumulate sufficient products from the metabolic reactions in order to move on to the next phase of the cycle and to autoinduce further cycles of oscillations. This picture is not contradicted by our observations.
If, as we do in this paper, rather than looking at the YMC merely as cyclic oscillations, we study it as a highly organized dynamical response to pulses of transcriptional activation, then this response can be analyzed in much more detail at genome-wide level and we can observe how an important role in the coordination seems to be played by the mRNA turnover rate. The self-sustained character of what we consider the most upstream event of the cycle, the transcriptional activation burst, can still be conditioned to the accumulation of the required metabolites, while the unfolding of the cycle, which from the analysis of  is already known to be functional to the distribution of e.g. the redox load of the cells, is enriched of an extra, intrinsically dynamical feature. This feature is a fine-graded detail of our notion that genes with a fast turnover are typically regulatory, while slow genes are enzymatic and metabolic [15, 19]. It can be used to describe the sequence of events occurring in the YMC as a "natural" gene expression program.
Extrapolating from the specific YMC context, the ordered pattern of events described for the YMC is to a good extent similar to that found on other intrinsically dynamical contexts such as the stress/stimuli responses. Whether the mRNA stability is the cause of this coherent behavior or is simply another effect of a more profound regulatory mechanism is a question to which we cannot provide a definitive answer at the moment.
The YMC time series of , the compendium of 790 gene profile experiments (all performed with the Affymetrix GeneChip Yeast Genome S98 platform) and the data series from  were downloaded from Gene Expression Omnibus . The time series of  are performed with cDNA, hence values of the area under the profiles are intended as relative (to the basal mRNA abundance). For each gene, the values obtained for the two different glucose stimuli are averaged. Five stress responses from  (two heat shocks of different amplitude, hydrogen peroxide, diamide, and sorbitol responses) are considered. The amplitudes are averaged over the five data series (the signs of these responses are known to be highly similar, see ).
The metabolic pathways used are those of the Kyoto Encyclopedia of Genes and Genomes (KEGG) . Also the assembling into the 15 macrocategories follows the KEGG hierarchy.
The HLs are computed averaging the values of the three experimental datasets [17, 18, 20]. While the magnitudes of the HLs in the three collections show some differences, in "normalized" terms (looking e.g. at rank-ordered values), the agreement between the three sets is sufficiently good, see  for a comparison. No turnover data specific for long-term continuous cultures are currently available. However, it is not unlikely that even in this setting the relative differences of HL rates (and also their ordering) remains more or less unchanged. In any case, we expect the correlation phase/HL to improve in presence of more tailored mRNA turnover data.
Time series analysis
To each of the genes labeled as periodic, we associated a phase, computed maximizing the correlation with respect to a train of 360 shifted sinusoids (resolution of 1°). The 0 phase was chosen so as to anticipate of ~30° the "crucial" transcription bursts [see Additional file 1]. Given that the period is approximately 287.5 minutes (see Fig. S1 in Additional file 1), the phase delay ϕ can be transformed into time delay τ by means of the relation . Under the convention for the 0 phase, each period "begins" approximately 24 min before the transcription bursts. For each gene, the pulse width is computed estimating on each period the interval in which the expression level stays above the median value across consecutive samples.
Least squares regressions
A minimal dynamical model: low-pass transfer functions and their dynamical system realizations
The aim of this Section is to set up a minimal dynamical model describing the response to the periodic bursts of transcriptional activation represented as "impulsive inputs" to the system. Such a model has to be able to reproduce the following features observable in the dataset:
• impulse responses get delayed and broadened in a way which is roughly proportional to HL;
• profile changes get progressively less steep with HL;
• the system "discharges" completely (i.e. the mRNA concentrations return to a basal level) in absence of further pulses.
At the same time, to be internally consistent a dynamical model has to:
• respect causality (i.e., be non-anticipating);
• preserve positivity of the mRNA concentrations.
In the Engineering practice of Systems Theory, one of the most elementary formalism that can be used to build dynamical models is the input-output design based on Laplace transform and elementary transfer functions , see e.g.  for an application to a transcriptional time series.
The concentration of mRNA of a gene y can be described as the response to the pulse of transcriptional activation u by the linear integral
In the Laplace domain, a convolution integral such as (1) corresponds to
where s is the Laplace variable and G(s) is called a transfer function. If u(t) is a perfect impulse δ0 (Dirac delta) then U(s) = ℒ[δ0(t)] = 1. When the transfer function G(s) represents a linear differential equation (i.e. it derives from a linear convolution such as (1)), it can be expressed as a rational polynomial in the Laplace variable s. A simple such polynomial is
where s = -d1 is called the pole of G1 and s = -n1 its zero. Choosing d1 > 0 the transfer function is stable (the pole is in the left half of the complex plane), i.e., a bounded input will always result in a bounded output. When n1 > 0 the system is said to be minimum phase. In this context this is an important condition in order to guarantee positivity of the output signal for all times. The requirements above can be translated into easy-to-handle design specifications on the values of the poles and zeros of the transfer function. For example, the first requirement (at least for what concerns pulse broadening) is met by the class of so-called low-pass filters, the most basic of which has the form given in (3), provided we choose 0 <d1 <n1. The term "low-pass" literally means that low frequencies in the input signal pass unchanged through the transfer function G1(s), while high frequencies get damped, hence the impulsive input exits from G1(s) smoothed and with more gentle slopes. Such a transfer function is proper and therefore respects causality; it discharges completely as required (since it has no integrator, i.e., no factors of the form 1/s in G1(s)). Strictly speaking, it is not a positive filter , however as long as u(t) > 0 and 0 <d1 <n1 it is also y(t) > 0. In the Laplace domain, a time delay T1 has Laplace transform equal to . This operator does not add poles or zeros to (3) but yields the irrational transfer function
In the time domain, each convolution integral (1) can be expressed as a linear input-output systems (of ordinary differential equations). For the transfer function in (3) and the delay operator in (4) this corresponds to
i.e., the pole d1 plays the role of "degradation rate" while the activation amplitude is proportional to n1 - d1 (> 0). The typical impulse response of a low-pass filter transfer function such as (3) is shown in the top plot of Fig. 5(b). Given a pulse shape, the capabilities of a single low-pass filter in terms of broadening and smoothing of the responses are limited, hence, in order to obtain a progressive effect of delayed and broadened impulse responses, several delayed low-pass filters should be put in cascade. For example the order-2 transfer function obtained concatenating 2 filters is
or, in the time domain,
In this case both d1 and d2 contribute to forming the degradation profile of the mRNA concentration y2(t). Likewise both dynamical variables x1 and x2 contribute to shape the pulse of a gene. Typically this model induces a steeper upregulation and a slower degradation front, coherently with what we observe on the YMC time series. The intermediate variables x i are only meant to describe the complexity of the input-output relationship. Qualitatively, they might reflect intermediate steps in the gene expression program. For example, the transcription of the genes of the central metabolism is activated downstream of the genes for translation and amino acid synthesis, which in their turn follow the principal bursts of transcription machinery (polymerases and other RNA processing components). Downstream activation of the genes of a category translates in this modeling framework into delayed and broadened pulses. Typical output responses for 1, 2, 3, and 4 such concatenated blocks are shown in Fig. 5(b).
A simple parameter search can be set up to identify values of n i , d i and T i , i = 1,...,4, that guarantee for each gene a sufficiently well-reproduced time course. The best transfer function order for each gene is identified as that maximizing the correlation between true and model-based time series.
HL and the short-period YMC of 
The HL of a gene is defined as the time needed to halve the concentration of mRNA in absence of new transcription. Hence in order for a "full" degradation of mRNA to be observed, the interval between two consecutive waves of transcription has to be at least twice or three times the HL. For yeast, the mean HL extrapolated from [17, 18, 20] is ~26 ± 17 min. Hence for the long-period YMC the response to bursts of transcription has the time to exhaust completely before the arrival of the next wavefront. On the contrary, for the short-period YMC described in  the period is approximately 40 min, meaning that excitation and degradation fronts are substantially overlapping.
Kyoto Encyclopedia of Genes and Genomes
yeast metabolic cycle.
Parulekar SJ, Semones GB, Rolf MJ, Lievense JC, Lim HC: Induction and elimination of oscillations in continuous cultures of Saccharomyces cerevisiae. Biotechnol and Bioeng. 1986, 28 (5): 700-710. 10.1002/bit.260280509
Porro D, Martegani E, Ranzi BM, Alberghina L: Oscillations in continuous cultures of budding yeast: A segregated parameter analysis. Biotechnol Bioeng. 1988, 32 (4): 411-417. 10.1002/bit.260320402
Klevecz RR, Bolen J, Forrest G, Murray DB: A genomewide oscillation in transcription gates DNA replication and cell cycle. Proc Natl Acad Sci USA. 2004, 101 (5): 1200-1205. 10.1073/pnas.0306490101
Tu BP, Kudlicki A, Rowicka M, McKnight SL: Logic of the Yeast Metabolic Cycle: Temporal Compartmentalization of Cellular Processes. Science. 2005, 310 (5751): 1152-1158. 10.1126/science.1120499
Patnaik PR: Oscillatory metabolism of Saccharomyces cerevisiae: an overview of mechanisms and models. Biotechnol Adv. 2003, 21 (3): 183-192. 10.1016/S0734-9750(03)00022-3
Reinke H, Gatfield D: Genome-wide oscillation of transcription in yeast. Trends Biochem Sci. 2006, 31 (4): 189-191. 10.1016/j.tibs.2006.02.001
Wolf J, Sohn HY, Heinrich R, Kuriyama H: Mathematical analysis of a mechanism for autonomous metabolic oscillations in continuous culture of Saccharomyces cerevisiae. FEBS Lett. 2001, 499 (3): 230-234. 10.1016/S0014-5793(01)02562-5
Chen Z, Odstrcil EA, Tu BP, McKnight SL: Restriction of DNA Replication to the Reductive Phase of the Metabolic Cycle Protects Genome Integrity. Science. 2007, 316 (5833): 1916-1919. 10.1126/science.1140958
Futcher B: Metabolic cycle, cell cycle, and the finishing kick to Start. Genome Biol. 2006, 7 (4): 107- 10.1186/gb-2006-7-4-107
Lloyd D, Murray DB: Redox rhythmicity: clocks at the core of temporal coherence. Bioessays. 2007, 29 (5): 465-473. 10.1002/bies.20575
Xu Z, Tsurugi K: A potential mechanism of energy-metabolism oscillation in an aerobic chemostat culture of the yeast Saccharomyces cerevisiae. FEBS J. 2006, 273 (8): 1696-1709. 10.1111/j.1742-4658.2006.05201.x
Beyer A, Hollunder J, Nasheuer HP, Wilhelm T: Post-transcriptional Expression Regulation in the Yeast Saccharomyces cerevisiae on a Genomic Scale. Mol Cell Proteomics. 2004, 3 (11): 1083-1092. 10.1074/mcp.M400099-MCP200
Brockmann R, Beyer A, Heinisch JJ, Wilhelm T: Posttranscriptional Expression Regulation: What Determines Translation Rates?. PLoS Comput Biol. 2007, 3 (3): e57- 10.1371/journal.pcbi.0030057
García-Martínez J, González-Candelas F, Pérez-Ortín JE: Common gene expression strategies revealed by genome-wide analysis in yeast. Genome Biol. 2007, 8 (10): R222- 10.1186/gb-2007-8-10-r222
Mata J, Marguerat S, Bähler J: Post-transcriptional control of gene expression: a genome-wide perspective. Trends Biochem Sci. 2005, 30 (9): 506-514. 10.1016/j.tibs.2005.07.005
Cheadle C, Fan J, Cho-Chung YS, Werner T, Ray J, Do L, Gorospe M, Becker KG: Stability Regulation of mRNA and the Control of Gene Expression. Ann N Y Acad Sci. 2005, 1058: 196-204. 10.1196/annals.1359.026
Grigull J, Mnaimneh S, Pootoolal J, Robinson MD, Hughes TR: Genome-Wide Analysis of mRNA Stability Using Transcription Inhibitors and Microarrays Reveals Posttranscriptional Control of Ribosome Biogenesis Factors. Mol Cell Biol. 2004, 24 (12): 5534-5547. 10.1128/MCB.24.12.5534-5547.2004
Kuai L, Das B, Sherman F: A nuclear degradation pathway controls the abundance of normal mRNAs in Saccharomyces cerevisiae. Proc Natl Acad Sci USA. 2005, 102 (39): 13962-13967. 10.1073/pnas.0506518102
Wilusz CJ, Wilusz J: Bringing the role of mRNA decay in the control of gene expression into focus. Trends Genet. 2004, 20 (10): 491-497. 10.1016/j.tig.2004.07.011
Wang Y, Liu CL, Storey JD, Tibshirani RJ, Herschlag D, Brown PO: Precision and functional specificity in mRNA decay. Proc Natl Acad Sci USA. 2002, 99 (9): 5860-5865. 10.1073/pnas.092538799
Keene JD: RNA regulons: coordination of post-transcriptional events. Nat Rev Genet. 2007, 8 (7): 533-543. 10.1038/nrg2111
Shalgi R, Lapidot M, Shamir R, Pilpel Y: A catalog of stability-associated sequence elements in 3' UTRs of yeast mRNAs. Genome Biol. 2005, 6 (10): R86- 10.1186/gb-2005-6-10-r86
Gerber AP, Herschlag D, Brown PO: Extensive Association of Functionally and Cytotopically Related mRNAs with Puf Family RNA-Binding Proteins in Yeast. PLoS Biol. 2004, 2 (3): 342-354. 10.1371/journal.pbio.0020079
Vallari RC, Cook WJ, Audino DC, Morgan MJ, Jensen DE, Laudano AP, Denis CL: Glucose Repression of the Yeast ADH2 Gene Occurs through Multiple Mechanisms, Including Control of the Protein Synthesis of Its Transcriptional Activator, ADR1. Mol Cell Biol. 1992, 12 (4): 1663-1673.
Cereghino GP, Scheffler IE: Genetic analysis of glucose regulation in Saccharomyces cerevisiae: control of transcription versus mRNA turnover. EMBO J. 1996, 15 (2): 363-374.
Kim JH, Brachet V, Moriya H, Johnston M: Integration of Transcriptional and Posttranslational Regulation in a Glucose Signal Transduction Pathway in Saccharomyces cerevisiae. Eukaryotic Cell. 2006, 5: 167-173. 10.1128/EC.5.1.167-173.2006
Hall DA, Zhu H, Zhu X, Gerstein M, Snyder M: Regulation of Gene Expression by a Metabolic Enzyme. Science. 2004, 306 (5695): 482-484. 10.1126/science.1096773
Ronen M, Botstein D: Transcriptional response of steady-state yeast cultures to transient perturbations in carbon source. Proc Natl Acad Sci USA. 2006, 103 (2): 389-394. 10.1073/pnas.0509978103
Gasch AP, Spellman PT, Kao CM, Carmel-Harel O, Eisen MB, Storz G, Botstein D, Brown PO: Genomic Expression Programs in the Response of Yeast Cells to Environmental Changes. Mol Biol Cell. 2000, 11 (12): 4241-4257.
Yin Z, Wilson S, Hauser NC, Tournu H, Hoheisel JD, Brown AJP: Glucose triggers different global responses in yeast, depending on the strength of the signal, and transiently stabilizes ribosomal protein mRNAs. Mol Microbiol. 2003, 48 (3): 713-724. 10.1046/j.1365-2958.2003.03478.x
Sohn HY, Murray DB, Kuriyama H: Ultradian oscillation of Saccharomyces cerevisiae during aerobic continuous culture: hydrogen sulphide mediates population synchrony. Yeast. 2000, 16 (13): 1185-1190. 10.1002/1097-0061(20000930)16:13<1185::AID-YEA619>3.0.CO;2-W
Tu BP, Mohler RE, Liu JC, Dombek KM, Young ET, Synovec RE, McKnight SL: Cyclic changes in metabolic state during the life of a yeast cell. Proc Natl Acad Sci USA. 2007, 104 (43): 16886-16891. 10.1073/pnas.0708365104
Hieronymus H, Silver PA: A systems view of mRNP biology. Genes Dev. 2004, 18 (23): 2845-2860. 10.1101/gad.1256904
Hieronymus H, Silver PA: Genome-wide analysis of RNA-protein interactions illustrates specificity of the mRNA export machinery. Nature Genetics. 2003, 33: 155-161. 10.1038/ng1080
Palumbo MC, Farina L, De Santis A, Giuliani A, Colosimo A, Morelli G, Ruberti I: Collective behavior in gene regulation: Post-transcriptional regulation and the temporal compartmentalization of cellular cycles. FEBS J. 2007, 275 (10): 2364-2371. 10.1111/j.1742-4658.2008.06398.x
MacIsaac KD, Wang T, Gordon DB, Gifford DK, Stormo GD, Fraenkel E: An improved map of conserved regulatory sites for Saccharomyces cerevisiae. BMC Bioinformatics. 2006, 7: 113- 10.1186/1471-2105-7-113
Balaji S, Babu MM, Iyer LM, Luscombe NM, Aravind L: Comprehensive Analysis of Combinatorial Regulation using the Transcriptional Regulatory Network of Yeast. J Mol Biol. 2006, 360: 213-227. 10.1016/j.jmb.2006.04.029
Xu Z, Tsurugi K: Role of Gts1p in regulation of energy-metabolism oscillation in continuous cultures of the yeast Saccharomyces cerevisiae. Yeast. 2007, 24 (3): 161-170. 10.1002/yea.1468
Tu BP, McKnight SL: Metabolic cycles as an underlying basis of biological oscillations. Nat Rev Mol Cell Biol. 2006, 7 (9): 696-701. 10.1038/nrm1980
Preiss T, Baron-Benhamou J, Ansorge W, Hentze MW: Homodirectional changes in transcriptome composition and mRNA translation induced by rapamycin and heat shock. Nat Struct Biol. 2003, 10 (12): 1039-1047. 10.1038/nsb1015
Belle A, Tanay A, Bitincka L, Shamir R, O'Shea EK: Quantification of protein half-lives in the budding yeast proteome. Proc Natl Acad Sci USA. 2006, 103 (35): 13004-13009. 10.1073/pnas.0605420103
Hargrove Jl, Schmidt FH: The role of mRNA and protein stability in gene expression. FASEB J. 1989, 3 (12): 2360-2370.
Zaslaver A, Mayo AE, Rosenberg R, Bashkin P, Sberro H, Tsalyuk M, Surette MG, Alon U: Just-in-time transcription program in metabolic pathways. Nat Genet. 2004, 36 (5): 486-491. 10.1038/ng1348
Chechik G, Oh E, Rando O, Weissman J, Regev A, Koller D: Activity motifs reveal principles of timing in transcriptional control of the yeast metabolic network. Nat Biotechnol. 2008, 26 (11): 1251-1259. 10.1038/nbt.1499
Gene Expression Omnibus.http://www.ncbi.nlm.nih.gov/geo/
Kyoto Encyclopedia of Genes and Genomes. http://www.genome.jp/kegg/
Doyle JC, Francis BA, Tannenbaum AR: Feedback Control Theory. 1992, New York: Macmillan Publishing Co
Becskei A, Boselli MG, van Oudenaarden A: Amplitude control of cell-cycle waves by nuclear import. Nat Cell Biol. 2004, 6 (5): 451-457. 10.1038/ncb1124
Anderson BDO, Deistler M, Farina L, Benvenuti L: Nonnegative realization of a linear system with nonnegative impulse response. IEEE Trans Circuits Syst I, Fundam Theory Appl. 1996, 43 (2): 134-142. 10.1109/81.486435.
NS and CA conceived the study, analyzed the data and drafted the manuscript. MZ contributed to the data analysis and the interpretation of the results. LF contributed to the interpretation of the results. All authors read and approved the final manuscript.