BACKGROUND
various statistical approaches have been specifically developed to summarise the vast quantities of data that are produced in microarray studies  <cit> , employing analysis of variance , clustering and network modelling.
analysis of variance  has been used to identify those gene expression responses that are most affected by different treatments, often taking account of particular forms of treatment structure, such as the correlations between sample times in a time-course study  <cit> .
approaches for clustering genes with similar responses range from simple methods for observed data, the calculation of correlations between genes  <cit> , through to clustering based on linear  <cit>  or polynomial regression  <cit>  or spline models  <cit> .
network models are used to reconstruct transcription factor activity  <cit>  or infer regulatory networks  <cit> , assuming a particular mechanistic model for the behaviour of each regulation function based on observed microarray gene expression data.
this paper aims to use standard statistical non-linear regression models to enhance the biological interpretation of individual gene expression profiles.
such regression models provide accessible methods to describe the shape of each gene expression profile as a function of time, thus providing an insight into the underlying processes rather than simply identifying significant differences.
for example, non-linear models can be used to identify the time of a particular event in a gene expression profile, such as the time of rapid up- or down-regulation.
similarly, modelling transcript changes using parametric equations that allow biological interpretation can further allow the comparison or clustering of the shapes of the expression profiles based on biological interpretable parameters.
such non-linear regression techniques are commonly used in agronomic studies to describe responses to a range of quantitative input variables, but are not commonly used in the examination of gene expression data.
the initial model system used to investigate the potential of statistical parametric, non-linear regression approaches for gene profiling was fungal morphogenesis with data provided by quantitative reverse transcriptase pcr  which provides a more precise method than either northern analysis or microarrays  <cit> .
this system encompasses a range of growth forms from vegetative mycelium to multicellular organs which enable the fungus to respond to changes in nutrition and environment, and undergo pathogenesis or reproduction.
the fruiting bodies of the basidiomycete fungus agaricus bisporus, the cultivated mushroom, are ideal for studying fungal morphogenesis as they are macroscopic, the tissues are clearly de-lineated  and the initiation of fruiting body morphogenesis is controlled environmentally.
differential screening and targeted gene cloning procedures have identified genes up-regulated post-harvest in a. bisporus fruiting bodies, based on northern analysis  <cit> .
genes have also been identified which are expressed in developing fruiting bodies of several other fungi, including lentinula edodes  <cit> , pleurotus ostreatus  <cit> , flammulina velutipes  <cit>  and coprinus cinereus   <cit> .
this study investigated how expression profiles, generated from qrt-pcr data of differentially regulated genes in a. bisporus fruiting bodies, could be statistically modelled both to estimate the time of up-regulation and determine similar temporal expression patterns.
the five genes chosen for profiling are functionally distinct, and therefore unlikely to have obvious common regulatory mechanisms: they are cruciform dna binding protein, cytochrome p450ii, β  glucan synthase, glucuronyl hydrolase and riboflavin aldehyde-forming enzyme  <cit> .
whilst being functionally distinct, these genes were expected to show broadly similar patterns of expression following harvest, allowing the fitting of a single model form to all five responses.
this enables the comparison of the profiles via biologically interpretable parameters rather than simply clustering genes based on the observed data.
to determine the time when transcription first increased for each gene, transcript levels of each gene were examined at  <dig> h intervals for  <dig> h. transcript levels were also measured at  <dig> hour intervals over  <dig> days.
these profiles were modelled to provide directly interpretable parameter values, offering an insight into the regulation of the genes.
spatial control of gene expression was also assessed by comparing transcription in tissues of the harvested mushroom during the first  <dig> hours post-harvest storage.
furthermore, this study has applied these regression modelling approaches to publicly-available microarray data sets from published studies, to identify groups of genes showing similar regulatory patterns.
the potential application of this approach to fully exploit the large quantities of high-throughout data was discussed.
RESULTS
the qrt-pcr-derived data were initially assessed to ensure suitability for further study.
standard curves show that the pcr reactions were operating at close to 100% efficiency.
melting curve analyses showed that in all cases a single product was obtained showing that the specificity of the reaction was high.
the pcr control treatments showed no evidence of dna contamination or active reverse transcription after the 85°c treatment .
comparison of methods of measuring gene expression in a. bisporus
gene transcript levels were obtained by qrt-pcr and northern analysis  from the same rna extracts obtained from the mushroom tissues.
these data were then compared by calculating correlation coefficients and fitting simple linear regression relationships for northern analysis response as a function of the qrt-pcr response, for each gene separately within each of the three experiments.
northern analysis responses increased with increasing values of the qrt-pcr response, and in most cases linear regression lines fitted well .
however, the northern analysis response appears to reach an upper threshold at high values of the qrt-pcr response, e.g.
for cruciform dna-binding enzyme and glucuronyl hydrolase , and β glucan synthase .
exponential regression curves were then fitted to the responses and for ten of the  <dig> data sets these curves provided a better fit than the simple linear regression models , with the curvature of the response suggesting asymptotic regression, for example, glucuronyl hydrolase, 0–5d experiment .
transcription in the first  <dig> hours
the different expression patterns for each gene over the first  <dig> hours following harvest, as obtained using qrt-pcr, are shown in figure  <dig>  during this period, transcript levels increased significantly, by approximately  <dig> fold or more, for all genes except riboflavin aldehyde-forming enzyme .
transcript levels of cruciform dna-binding protein increased shortly after harvest to very high levels, while for cytochrome p450ii and glucuronyl hydrolase the profiles showed an initial period of little change followed by increases in transcript levels .
in contrast the transcript levels of β glucan synthase increased continuously during the first  <dig> hours post-harvest, with a notable sudden increase at about  <dig> h.
'split-line' or 'broken-stick' analysis was used to model the transcription profile data from the first  <dig> hours to estimate the time at which an initial increase in transcription occurred.
the model was applied successfully to the transcript data for cytochrome p450ii and glucuronyl hydrolase, with these data sets each described by two separate linear regression line segments, the first line having a slope of zero .
more precise estimates of when transcription first increased were determined in these analysis as  <dig>  hours ) for cytochrome p450ii  and  <dig>  hours ) for glucuronyl hydrolase .
this analytical approach was not successful in describing the transcription profiles of cruciform dna-binding protein and β glucan synthase, as there were insufficient early time points to establish an initial baseline response .
transcripts of riboflavin aldehyde-forming enzyme showed no clear upward trend in the first  <dig> hours, and so the 'broken-stick' analysis approach was again not successful.
transcription profiling over  <dig> days post-harvest
the 5-day transcript profiles, as obtained using qrt-pcr, showed an initial increase in transcript levels for all five genes for 1– <dig> days, followed by a plateau or a decline .
the extent of increase in transcript levels over the first two days ranged from  <dig> fold for cruciform dna-binding protein, to  <dig> fold for riboflavin aldehyde-forming enzyme and  <dig> fold for glucuronyl hydrolase.
gene expression over five days storage, as measured using qrt-pcr, was modelled as a function of storage time by non-linear regression analysis of log10-transformed data using a critical exponential curve  = a + rt + ε).
this curve was selected in preference over either the simpler single exponential curve or more complex double exponential curve, based both on the shape of the response, and a formal comparison of the goodness of fit of each model.
this was achieved by comparing each more complex model to the simpler alternative .
the improvement in fit was assessed using an f-test to identify the significance of the ratio of the change in residual variance between models to the residual variance for the more complex model.
for all genes, there was a significant improvement in the fit when choosing the critical exponential rather than the single exponential, but no significant improvement when choosing the double exponential over the critical exponential .
simultaneous fitting of the critical exponential curve to the data for all five genes allowed a comparison of the fitted parameters, and hence the detection of degrees of commonality between the patterns of expression for the genes.
this analysis identified that there was no significant improvement to the fit when allowing parameter r  to be different for the five genes, and so this parameter could be constrained to be the same.
however, constraining the other three parameters to be the same across all five genes resulted in significantly worse fits, though for some pairs of genes the fitted  parameters did suggest some similarities .
the different parameters  can be interpreted in terms of particular features of the shape of the response.
the size of parameter c indicates the magnitude of the decline from the maximum expression response, whilst the ratio of b over c is related to the time to the maximum response.
parameter a measures the asymptotic response level after lengthy storage, a + b is the response at time t =  <dig>  and the maximum response is dependent on all four parameters, obtained by inserting the time of maximum response into the critical exponential equation .
the analysis identified three distinct regulatory patterns .
cruciform dna-binding protein and cytochrome p450ii were identified as having similar shaped curves  and similar times of peak transcript levels , despite large differences in transcript magnitudes during this five-day period .
similarly, the transcript profiles of glucuronyl hydrolase and riboflavin aldehyde-forming enzyme shared a common curve shape  and similar transcript values for both the initial and asymptotic parts of the curve .
the time of maximum transcript level of each of these genes was later than for the other genes examined .
the  <dig> day β glucan synthase profile had a different curve shape with the earliest peak of transcript level  followed by a rapid decline .
maximum
response
response
fitted parameters for five genes  glucan synthase, and rafe = riboflavin aldehyde-forming enzyme) with parameter r constrained to be common across all five genes ).
additional descriptors of response shape derived from the fitted parameters.
transcript expression between tissues
transcript levels of each gene were measured in stipe, cap and gill tissues of harvested mushrooms over  <dig> days storage, using qrt-pcr .
differences in transcript levels, between tissues, between storage times and due to the interaction of these two factors, were assessed using anova.
average transcript levels of cruciform dna-binding protein, cytochrome p450ii, glucuronyl hydrolase and β  glucan synthase in stipe and cap tissues were similar and significantly higher than in the gills.
for riboflavin aldehyde-forming enzyme, transcript levels were significantly higher in the stipe tissue compared with the cap or gills, which had similar levels.
despite differences in transcript levels observed between the tissues, all genes showed an increased level of expression in all tissues from day  <dig> to  <dig>
application of regression approaches to microarray datasets
regression analyses were applied to time course profiles obtained from two published microarray datasets  <cit> .
gene responses of e. coli to treatment with paraquat  <cit>  were analysed by fitting an exponential function.
from over  <dig>  features on the microarray the profiles for 11% were described well by the exponential function, with significance levels of p <  <dig> ; 2% of the gene profiles fitted with significance levels of p <  <dig>  .
the critical exponential function was fitted to the gene time-course responses in the rat liver tissue to the corticosteroid, methylprednisolone,  <cit> .
the r. norvegicus microarray  consisted of  <dig> full length sequences and over  <dig> expressed sequence tagged clusters.
from over  <dig> features on the microarray, the profiles for over 25% were described well by the critical exponential , with over 9% with significant levels of p <  <dig>  .
for both studies the chosen functions allow the description of a number of distinct forms of response.
the fitted responses for the  <dig> most significant fits from each study demonstrate the variety of profiles that can be described by each of these models .
discussion
gene expression studies were conducted to investigate the benefit of applying statistical regression approaches to the fitting of mathematical models for the analysis of transcriptional data.
in this study these approaches have been developed using precisely measured transcript levels for a small number of genes and application of the approaches has been further demonstrated for microarray datasets.
as microarray technologies continue to be developed, the variability of gene expression data from such technologies will be reduced, leading to the widespread application of the regression modelling techniques developed in this study, thus allowing the comparative analysis of larger numbers of gene expression profiles.
a range of statistical regression techniques, both linear and non-linear, are readily available.
the selection of appropriate techniques is critically dependent on the specific question being addressed and the data that are collected.
for a. bisporus, the increases in transcript levels in the first  <dig> hours following harvest are largely due to transcription rather than losses due to mrna turnover.
'split-line or 'broken-stick' regression analysis was used to calculate the time when transcription was initiated, at  <dig>  h post-harvest for cytochrome p450ii and  <dig>  h post-harvest for glucuronyl hydrolase.
the novel application of this approach to transcript profiling demonstrates the potential value of this simple mathematical model which has been used previously in such diverse applications as estimating the thresholds for patch size in ecological studies  <cit> , humidity levels in plant pathogen germination experiments  <cit>  and the mineral density of bones  <cit> .
successful application of the 'split-line' or 'broken-stick' model for cruciform dna-binding protein and β  glucan synthase would require more sampling points to be made in the first  <dig> hours to provide sufficient time points to allow the fitting of the baseline response, zero-slope line, and hence allow estimation of the time of initial transcription.
for the riboflavin aldehyde-forming enzyme further data are needed beyond  <dig> hours  to allow estimation of the second linear regression segment, and to again allow calculation of the time of initial transcription.
the time of increased transcription for at least  <dig> of the  <dig> tested genes occurs at different times in the first  <dig> hours post-harvest.
this suggests that the response of the mushroom to harvest is not under the control of a single regulatory pathway, such as the signal transduction pathways described for fungal oxidative and osmotic stress responses  <cit> .
the controlling events in the mushroom are likely to be affected by a range of stimuli, such as stress, nutrient limitation, continued maturation and spore formation, which might illicit both primary or immediate responses, and secondary responses.
here we observed increased transcription first of cruciform dna-binding protein , followed by β  glucan synthase , cytochrome p450ii , glucuronyl hydrolase  and riboflavin aldehyde-forming enzyme .
cruciform dna-binding protein and β  glucan synthase may be part of a primary response, while cytochrome p450ii, glucuronyl hydrolase and riboflavin aldehyde-forming enzyme are from a secondary response or caused by a later stimulus.
statistical regression modelling showed similar expression patterns for two a. bisporus genes, glucuronyl hydrolase and riboflavin aldehyde-forming enzyme, the latter of which is known to be up-regulated during the development of the non-harvested mushroom  <cit> .
further study is required to determine whether the common pattern observed between the genes post-harvest is also observed in the morphogenesis of the non-harvested fruiting body.
the pattern of β  glucan synthase gene expression during  <dig> days storage was different from the other genes studied, with transcript levels falling greatly after  <dig> d following an initial increase in gene expression.
the gene has been hypothesised to be involved in cell wall synthesis  <cit> .
the initial increase in gene expression coincides with the period when hyphae in the cap and stipe elongate, i.e.
in the first  <dig> days post-harvest, a process for which cell wall synthesis is important.
similarly, the reduction in expression after day  <dig> coincides with the cessation of cell wall synthesis following the full extension of the cap  <cit> .
the fitted critical exponential models for the 5-day profiles of cruciform dna-binding protein and cytochrome p450ii had similar shapes, whilst the times of initial transcriptional increase for these two genes, as determined by split-line regressions, were markedly different .
this apparent paradox illustrates the importance of considering transcript responses over a number of different time-scales.
the precise function of cruciform dna-binding protein is not known, however, it is unlikely to be involved in recombination as transcript levels are low in the gill tissue where meiosis occurs.
in other organisms, proteins with similar cruciform dna-binding activity  cause the increased and decreased transcription of genes  <cit> .
it is possible that the early and abundant transcription of cruciform dna-binding protein in a. bisporus fruiting bodies acts to regulate the expression of other genes.
the spatial transcript levels were different between tissues and for all genes were significantly lower in the gill tissue, indicating a physiological difference from the stipe and cap.
the gills are actively growing, respiring and producing spores via meiosis.
they are a nutrient 'sink' and are therefore physiologically different from the stipe and cap tissues, which export nutrients and are subject to stress from cell damage.
however the increased respiration in gill tissue  <cit>  may result in increased ribosomal, and therefore 18s rrna, levels, to which transcript levels from qrt-pcr are normalised, so the quantitative differences between gills and stipe/cap tissues may be influenced by this.
while pairs of genes showed similar patterns of response when described using the critical exponential model, the spatial distribution of transcripts between different tissues varied in some cases.
for example, cytochrome p450ii and cruciform dna-binding protein showed similar overall patterns, but cytochrome p450ii showed a delayed initial increase in transcription in the gill tissue.
further and more detailed studies of the expression of these genes in the separate tissues are needed to understand the potentially different regulatory mechanisms acting in each tissue.
the use of a statistical non-linear regression approach to model the gene expression profiles over an extended period offers the opportunity to compare the shapes of the response curves between genes.
the critical exponential curve, selected by the observed response shapes and goodness of fit, can be explained in terms of the combination of two processes, in this case rna synthesis and degradation.
a more complex model for such a situation would be the double exponential curve, which is a natural function for a two-timescale process.
the critical exponential is the degenerate case of this model and occurs when two processes in a system have the same timescale.
in this case, however, there was no evidence for choosing the more complex model.
choice of an appropriate model can be important, but many common non-linear models are based on functional forms derived from observed biological processes.
increased gene transcription is responsible for the initial rapid increase in transcript levels between days  <dig> to  <dig>  transcript levels continue to rise until a maximum is reached, followed by a decline towards a steady level, possibly as a consequence of a balance between transcription and degradation .
the critical exponential model was successful in identifying genes with quantitatively similar patterns of response, which could not have been predicted from their putative protein functions.
this approach, therefore, offers a new method by which a large number of genes could be classified according to their initial transcription regulation and subsequent turnover.
our approach to first model and then cluster allows more genes to be considered and potentially a greater insight to understand the system.
the application of this approach to a microarray dataset allows the screening of genes to identify responses that can be described by a particular mathematical function.
thus each gene profile is reduced from noisy observations to a smaller set of biologically-interpretable parameters.
in the analysis of the microarray datasets, different shapes of profiles fitted by the exponential or critical exponential functions can be identified , allowing the grouping of genes based on the parameter values.
eliminating the inherent variability in the data through the regression modelling approach allows a more precise comparison of gene profiles and thus improved clustering.
for example, the 9% of r. norvegicus genes identified that followed the critical exponential curve at the p <  <dig>  level represents approx  <dig> genes compared with the  <dig> genes identified and modelled by jin et al   <cit> .
the groupings of genes then generated by the improved clustering propose hypotheses of regulatory association between genes.
for example, the aim of the e. coli microarray study was to identify those genes co-regulated with the main regulatory gene, soxs, which was demonstrated to have an exponential-type response following the application of paraquat  <cit> .
by using an initial regression analysis to identify the subset of genes that can be described by an exponential function, subsequent cluster analyses can focus on this subset of genes with similar, but not identical , shapes of expression profiles.
the fitted exponential function parameters for this gene subset could then be used to better identify those genes most closely co-regulated with soxs.
whilst our study demonstrates the application of a regression modelling approach to describe gene expression profiles, this approach can be expanded to fully exploit microarray datasets.
application of a wider range of functional forms  offers the potential to develop regulatory networks based on relationships between the shapes of expression profiles as captured by the fitted parameters.
further interpretation of the parameters, alongside knowledge of gene function, might allow the identification of the stimuli driving the observed gene expression responses.
compared with standard clustering for gene profiles of microarray data, the statistical regression modelling of mathematical functions to describe these profiles eliminates the inherent variability in the data and allows the direct comparison of profile shapes.
CONCLUSIONS
this study has illustrated how the use of standard statistical modelling approaches , linear regression modelling, non-linear regression modelling) commonly used in plant, microbial and ecological sciences, can be used to aid and extend the interpretation of gene expression profiles obtained from qrt-pcr.
these approaches have been applied to model profiles of larger numbers of genes obtained from expression microarray studies, and could be further applied to other high-throughput "-omic" technologies.
a wide range of statistical approaches have been specifically developed to analyse the vast quantities of data generated in microarray-based studies  <cit> , assessing both similarities and differences between genes and between treatments.
similarly a number of approaches have been developed to generate mathematical models for assumed networks of gene pathways   <cit> .
however, there appears to have been little statistical consideration for the detailed modelling of individual gene expression profiles.
the statistical regression modelling approaches applied in this study allow the estimation of parameters which succinctly describe the shapes of gene expression profiles.
these parameter estimates  can then be related directly to the processes stimulating and driving the expression of these genes.
comparison of the parameters and expression profiles for a set of genes could then indicate that a sub-set of these genes are co-regulated, with the potential to hypothesise a common regulatory mechanism.
hence, consideration of a wide range of non-linear regression models could provide building blocks for the development of more biologically realistic models of gene expression profiles.

 various statistical approaches have been specifically developed to summarise the vast quantities of data that are produced in microarray studies, employing analysis of variance, clustering and network modelling. analysis of variance has been used to identify those gene expression responses that are most affected by different treatments, often taking account of particular forms of treatment structure, such as the correlations between sample times in a time-course study. approaches for clustering genes with similar responses range from simple methods for observed data, the calculation of correlations between genes, through to clustering based on linear or polynomial regression or spline models. network models are used to reconstruct transcription factor activity or infer regulatory networks, assuming a particular mechanistic model for the behaviour of each regulation function based on observed microarray gene expression data. this paper aims to use standard statistical non-linear regression models to enhance the biological interpretation of individual gene expression profiles. such regression models provide accessible methods to describe the shape of each gene expression profile as a function of time, thus providing an insight into the underlying processes rather than simply identifying significant differences. for example, non-linear models can be used to identify the time of a particular event in a gene expression profile, such as the time of rapid up- or down-regulation. similarly, modelling transcript changes using parametric equations that allow biological interpretation can further allow the comparison or clustering of the shapes of the expression profiles based on biological interpretable parameters. such non-linear regression techniques are commonly used in agronomic studies to describe responses to a range of quantitative input variables, but are not commonly used in the examination of gene expression data. the initial model system used to investigate the potential of statistical parametric, non-linear regression approaches for gene profiling was fungal morphogenesis with data provided by quantitative reverse transcriptase pcr which provides a more precise method than either northern analysis or microarrays. this system encompasses a range of growth forms from vegetative mycelium to multicellular organs which enable the fungus to respond to changes in nutrition and environment, and undergo pathogenesis or reproduction. the fruiting bodies of the basidiomycete fungus agaricus bisporus, the cultivated mushroom, are ideal for studying fungal morphogenesis as they are macroscopic, the tissues are clearly de-lineated and the initiation of fruiting body morphogenesis is controlled environmentally. differential screening and targeted gene cloning procedures have identified genes up-regulated post-harvest in a. bisporus fruiting bodies, based on northern analysis. genes have also been identified which are expressed in developing fruiting bodies of several other fungi, including lentinula edodes, pleurotus ostreatus, flammulina velutipes and coprinus cinereus. this study investigated how expression profiles, generated from qrt-pcr data of differentially regulated genes in a. bisporus fruiting bodies, could be statistically modelled both to estimate the time of up-regulation and determine similar temporal expression patterns. the five genes chosen for profiling are functionally distinct, and therefore unlikely to have obvious common regulatory mechanisms: they are cruciform dna binding protein, cytochrome p450ii, β glucan synthase, glucuronyl hydrolase and riboflavin aldehyde-forming enzyme. whilst being functionally distinct, these genes were expected to show broadly similar patterns of expression following harvest, allowing the fitting of a single model form to all five responses. this enables the comparison of the profiles via biologically interpretable parameters rather than simply clustering genes based on the observed data. to determine the time when transcription first increased for each gene, transcript levels of each gene were examined at h intervals for h. transcript levels were also measured at hour intervals over days. these profiles were modelled to provide directly interpretable parameter values, offering an insight into the regulation of the genes. spatial control of gene expression was also assessed by comparing transcription in tissues of the harvested mushroom during the first hours post-harvest storage. furthermore, this study has applied these regression modelling approaches to publicly-available microarray data sets from published studies, to identify groups of genes showing similar regulatory patterns. the potential application of this approach to fully exploit the large quantities of high-throughout data was discussed.  the qrt-pcr-derived data were initially assessed to ensure suitability for further study. standard curves show that the pcr reactions were operating at close to 100% efficiency. melting curve analyses showed that in all cases a single product was obtained showing that the specificity of the reaction was high. the pcr control treatments showed no evidence of dna contamination or active reverse transcription after the 85°c treatment. comparison of methods of measuring gene expression in a. bisporus gene transcript levels were obtained by qrt-pcr and northern analysis from the same rna extracts obtained from the mushroom tissues. these data were then compared by calculating correlation coefficients and fitting simple linear regression relationships for northern analysis response as a function of the qrt-pcr response, for each gene separately within each of the three experiments. northern analysis responses increased with increasing values of the qrt-pcr response, and in most cases linear regression lines fitted well. however, the northern analysis response appears to reach an upper threshold at high values of the qrt-pcr response, e.g. for cruciform dna-binding enzyme and glucuronyl hydrolase, and β glucan synthase. exponential regression curves were then fitted to the responses and for ten of the data sets these curves provided a better fit than the simple linear regression models, with the curvature of the response suggesting asymptotic regression, for example, glucuronyl hydrolase, 0–5d experiment. transcription in the first hours the different expression patterns for each gene over the first hours following harvest, as obtained using qrt-pcr, are shown in figure during this period, transcript levels increased significantly, by approximately fold or more, for all genes except riboflavin aldehyde-forming enzyme. transcript levels of cruciform dna-binding protein increased shortly after harvest to very high levels, while for cytochrome p450ii and glucuronyl hydrolase the profiles showed an initial period of little change followed by increases in transcript levels. in contrast the transcript levels of β glucan synthase increased continuously during the first hours post-harvest, with a notable sudden increase at about h.'split-line' or'broken-stick' analysis was used to model the transcription profile data from the first hours to estimate the time at which an initial increase in transcription occurred. the model was applied successfully to the transcript data for cytochrome p450ii and glucuronyl hydrolase, with these data sets each described by two separate linear regression line segments, the first line having a slope of zero. more precise estimates of when transcription first increased were determined in these analysis as hours ) for cytochrome p450ii and hours ) for glucuronyl hydrolase. this analytical approach was not successful in describing the transcription profiles of cruciform dna-binding protein and β glucan synthase, as there were insufficient early time points to establish an initial baseline response. transcripts of riboflavin aldehyde-forming enzyme showed no clear upward trend in the first hours, and so the'broken-stick' analysis approach was again not successful. transcription profiling over days post-harvest the 5-day transcript profiles, as obtained using qrt-pcr, showed an initial increase in transcript levels for all five genes for 1– days, followed by a plateau or a decline. the extent of increase in transcript levels over the first two days ranged from fold for cruciform dna-binding protein, to fold for riboflavin aldehyde-forming enzyme and fold for glucuronyl hydrolase. gene expression over five days storage, as measured using qrt-pcr, was modelled as a function of storage time by non-linear regression analysis of log10-transformed data using a critical exponential curve = a + rt + ε). this curve was selected in preference over either the simpler single exponential curve or more complex double exponential curve, based both on the shape of the response, and a formal comparison of the goodness of fit of each model. this was achieved by comparing each more complex model to the simpler alternative. the improvement in fit was assessed using an f-test to identify the significance of the ratio of the change in residual variance between models to the residual variance for the more complex model. for all genes, there was a significant improvement in the fit when choosing the critical exponential rather than the single exponential, but no significant improvement when choosing the double exponential over the critical exponential. simultaneous fitting of the critical exponential curve to the data for all five genes allowed a comparison of the fitted parameters, and hence the detection of degrees of commonality between the patterns of expression for the genes. this analysis identified that there was no significant improvement to the fit when allowing parameter r to be different for the five genes, and so this parameter could be constrained to be the same. however, constraining the other three parameters to be the same across all five genes resulted in significantly worse fits, though for some pairs of genes the fitted parameters did suggest some similarities. the different parameters can be interpreted in terms of particular features of the shape of the response. the size of parameter c indicates the magnitude of the decline from the maximum expression response, whilst the ratio of b over c is related to the time to the maximum response. parameter a measures the asymptotic response level after lengthy storage, a + b is the response at time t = and the maximum response is dependent on all four parameters, obtained by inserting the time of maximum response into the critical exponential equation. the analysis identified three distinct regulatory patterns. cruciform dna-binding protein and cytochrome p450ii were identified as having similar shaped curves and similar times of peak transcript levels, despite large differences in transcript magnitudes during this five-day period. similarly, the transcript profiles of glucuronyl hydrolase and riboflavin aldehyde-forming enzyme shared a common curve shape and similar transcript values for both the initial and asymptotic parts of the curve. the time of maximum transcript level of each of these genes was later than for the other genes examined. the day β glucan synthase profile had a different curve shape with the earliest peak of transcript level followed by a rapid decline. maximum response response fitted parameters for five genes glucan synthase, and rafe = riboflavin aldehyde-forming enzyme) with parameter r constrained to be common across all five genes ). additional descriptors of response shape derived from the fitted parameters. transcript expression between tissues transcript levels of each gene were measured in stipe, cap and gill tissues of harvested mushrooms over days storage, using qrt-pcr. differences in transcript levels, between tissues, between storage times and due to the interaction of these two factors, were assessed using anova. average transcript levels of cruciform dna-binding protein, cytochrome p450ii, glucuronyl hydrolase and β glucan synthase in stipe and cap tissues were similar and significantly higher than in the gills. for riboflavin aldehyde-forming enzyme, transcript levels were significantly higher in the stipe tissue compared with the cap or gills, which had similar levels. despite differences in transcript levels observed between the tissues, all genes showed an increased level of expression in all tissues from day to application of regression approaches to microarray datasets regression analyses were applied to time course profiles obtained from two published microarray datasets. gene responses of e. coli to treatment with paraquat were analysed by fitting an exponential function. from over features on the microarray the profiles for 11% were described well by the exponential function, with significance levels of p <; 2% of the gene profiles fitted with significance levels of p <. the critical exponential function was fitted to the gene time-course responses in the rat liver tissue to the corticosteroid, methylprednisolone,. the r. norvegicus microarray consisted of full length sequences and over expressed sequence tagged clusters. from over features on the microarray, the profiles for over 25% were described well by the critical exponential, with over 9% with significant levels of p <. for both studies the chosen functions allow the description of a number of distinct forms of response. the fitted responses for the most significant fits from each study demonstrate the variety of profiles that can be described by each of these models. discussion gene expression studies were conducted to investigate the benefit of applying statistical regression approaches to the fitting of mathematical models for the analysis of transcriptional data. in this study these approaches have been developed using precisely measured transcript levels for a small number of genes and application of the approaches has been further demonstrated for microarray datasets. as microarray technologies continue to be developed, the variability of gene expression data from such technologies will be reduced, leading to the widespread application of the regression modelling techniques developed in this study, thus allowing the comparative analysis of larger numbers of gene expression profiles. a range of statistical regression techniques, both linear and non-linear, are readily available. the selection of appropriate techniques is critically dependent on the specific question being addressed and the data that are collected. for a. bisporus, the increases in transcript levels in the first hours following harvest are largely due to transcription rather than losses due to mrna turnover.'split-line or'broken-stick' regression analysis was used to calculate the time when transcription was initiated, at h post-harvest for cytochrome p450ii and h post-harvest for glucuronyl hydrolase. the novel application of this approach to transcript profiling demonstrates the potential value of this simple mathematical model which has been used previously in such diverse applications as estimating the thresholds for patch size in ecological studies, humidity levels in plant pathogen germination experiments and the mineral density of bones. successful application of the'split-line' or'broken-stick' model for cruciform dna-binding protein and β glucan synthase would require more sampling points to be made in the first hours to provide sufficient time points to allow the fitting of the baseline response, zero-slope line, and hence allow estimation of the time of initial transcription. for the riboflavin aldehyde-forming enzyme further data are needed beyond hours to allow estimation of the second linear regression segment, and to again allow calculation of the time of initial transcription. the time of increased transcription for at least of the tested genes occurs at different times in the first hours post-harvest. this suggests that the response of the mushroom to harvest is not under the control of a single regulatory pathway, such as the signal transduction pathways described for fungal oxidative and osmotic stress responses. the controlling events in the mushroom are likely to be affected by a range of stimuli, such as stress, nutrient limitation, continued maturation and spore formation, which might illicit both primary or immediate responses, and secondary responses. here we observed increased transcription first of cruciform dna-binding protein, followed by β glucan synthase, cytochrome p450ii, glucuronyl hydrolase and riboflavin aldehyde-forming enzyme. cruciform dna-binding protein and β glucan synthase may be part of a primary response, while cytochrome p450ii, glucuronyl hydrolase and riboflavin aldehyde-forming enzyme are from a secondary response or caused by a later stimulus. statistical regression modelling showed similar expression patterns for two a. bisporus genes, glucuronyl hydrolase and riboflavin aldehyde-forming enzyme, the latter of which is known to be up-regulated during the development of the non-harvested mushroom. further study is required to determine whether the common pattern observed between the genes post-harvest is also observed in the morphogenesis of the non-harvested fruiting body. the pattern of β glucan synthase gene expression during days storage was different from the other genes studied, with transcript levels falling greatly after d following an initial increase in gene expression. the gene has been hypothesised to be involved in cell wall synthesis. the initial increase in gene expression coincides with the period when hyphae in the cap and stipe elongate, i.e. in the first days post-harvest, a process for which cell wall synthesis is important. similarly, the reduction in expression after day coincides with the cessation of cell wall synthesis following the full extension of the cap. the fitted critical exponential models for the 5-day profiles of cruciform dna-binding protein and cytochrome p450ii had similar shapes, whilst the times of initial transcriptional increase for these two genes, as determined by split-line regressions, were markedly different. this apparent paradox illustrates the importance of considering transcript responses over a number of different time-scales. the precise function of cruciform dna-binding protein is not known, however, it is unlikely to be involved in recombination as transcript levels are low in the gill tissue where meiosis occurs. in other organisms, proteins with similar cruciform dna-binding activity cause the increased and decreased transcription of genes. it is possible that the early and abundant transcription of cruciform dna-binding protein in a. bisporus fruiting bodies acts to regulate the expression of other genes. the spatial transcript levels were different between tissues and for all genes were significantly lower in the gill tissue, indicating a physiological difference from the stipe and cap. the gills are actively growing, respiring and producing spores via meiosis. they are a nutrient'sink' and are therefore physiologically different from the stipe and cap tissues, which export nutrients and are subject to stress from cell damage. however the increased respiration in gill tissue may result in increased ribosomal, and therefore 18s rrna, levels, to which transcript levels from qrt-pcr are normalised, so the quantitative differences between gills and stipe/cap tissues may be influenced by this. while pairs of genes showed similar patterns of response when described using the critical exponential model, the spatial distribution of transcripts between different tissues varied in some cases. for example, cytochrome p450ii and cruciform dna-binding protein showed similar overall patterns, but cytochrome p450ii showed a delayed initial increase in transcription in the gill tissue. further and more detailed studies of the expression of these genes in the separate tissues are needed to understand the potentially different regulatory mechanisms acting in each tissue. the use of a statistical non-linear regression approach to model the gene expression profiles over an extended period offers the opportunity to compare the shapes of the response curves between genes. the critical exponential curve, selected by the observed response shapes and goodness of fit, can be explained in terms of the combination of two processes, in this case rna synthesis and degradation. a more complex model for such a situation would be the double exponential curve, which is a natural function for a two-timescale process. the critical exponential is the degenerate case of this model and occurs when two processes in a system have the same timescale. in this case, however, there was no evidence for choosing the more complex model. choice of an appropriate model can be important, but many common non-linear models are based on functional forms derived from observed biological processes. increased gene transcription is responsible for the initial rapid increase in transcript levels between days to transcript levels continue to rise until a maximum is reached, followed by a decline towards a steady level, possibly as a consequence of a balance between transcription and degradation. the critical exponential model was successful in identifying genes with quantitatively similar patterns of response, which could not have been predicted from their putative protein functions. this approach, therefore, offers a new method by which a large number of genes could be classified according to their initial transcription regulation and subsequent turnover. our approach to first model and then cluster allows more genes to be considered and potentially a greater insight to understand the system. the application of this approach to a microarray dataset allows the screening of genes to identify responses that can be described by a particular mathematical function. thus each gene profile is reduced from noisy observations to a smaller set of biologically-interpretable parameters. in the analysis of the microarray datasets, different shapes of profiles fitted by the exponential or critical exponential functions can be identified, allowing the grouping of genes based on the parameter values. eliminating the inherent variability in the data through the regression modelling approach allows a more precise comparison of gene profiles and thus improved clustering. for example, the 9% of r. norvegicus genes identified that followed the critical exponential curve at the p < level represents approx genes compared with the genes identified and modelled by jin et al. the groupings of genes then generated by the improved clustering propose hypotheses of regulatory association between genes. for example, the aim of the e. coli microarray study was to identify those genes co-regulated with the main regulatory gene, soxs, which was demonstrated to have an exponential-type response following the application of paraquat. by using an initial regression analysis to identify the subset of genes that can be described by an exponential function, subsequent cluster analyses can focus on this subset of genes with similar, but not identical, shapes of expression profiles. the fitted exponential function parameters for this gene subset could then be used to better identify those genes most closely co-regulated with soxs. whilst our study demonstrates the application of a regression modelling approach to describe gene expression profiles, this approach can be expanded to fully exploit microarray datasets. application of a wider range of functional forms offers the potential to develop regulatory networks based on relationships between the shapes of expression profiles as captured by the fitted parameters. further interpretation of the parameters, alongside knowledge of gene function, might allow the identification of the stimuli driving the observed gene expression responses. compared with standard clustering for gene profiles of microarray data, the statistical regression modelling of mathematical functions to describe these profiles eliminates the inherent variability in the data and allows the direct comparison of profile shapes.  this study has illustrated how the use of standard statistical modelling approaches, linear regression modelling, non-linear regression modelling) commonly used in plant, microbial and ecological sciences, can be used to aid and extend the interpretation of gene expression profiles obtained from qrt-pcr. these approaches have been applied to model profiles of larger numbers of genes obtained from expression microarray studies, and could be further applied to other high-throughput "-omic" technologies. a wide range of statistical approaches have been specifically developed to analyse the vast quantities of data generated in microarray-based studies, assessing both similarities and differences between genes and between treatments. similarly a number of approaches have been developed to generate mathematical models for assumed networks of gene pathways. however, there appears to have been little statistical consideration for the detailed modelling of individual gene expression profiles. the statistical regression modelling approaches applied in this study allow the estimation of parameters which succinctly describe the shapes of gene expression profiles. these parameter estimates can then be related directly to the processes stimulating and driving the expression of these genes. comparison of the parameters and expression profiles for a set of genes could then indicate that a sub-set of these genes are co-regulated, with the potential to hypothesise a common regulatory mechanism. hence, consideration of a wide range of non-linear regression models could provide building blocks for the development of more biologically realistic models of gene expression profiles. 
previous studies have summarised gene expression profiles using simple descriptive statistics, basic analysis of variance and the clustering of genes based on simple models fitted to their expression profiles over time. we report the novel application of statistical non-linear regression modelling techniques to describe the shapes of expression profiles for the fungus agaricus bisporus, quantified by pcr, and for e. coli and rattus norvegicus, using microarray technology. the use of parametric non-linear regression models provides a more precise description of expression profiles, reducing the "noise" of the raw data to produce a clear "signal" given by the fitted curve, and describing each profile with a small number of biologically interpretable parameters. the statistical non-linear regression approaches presented in this study provide detailed biologically oriented descriptions of individual gene expression profiles, using biologically variable data to generate a set of defining parameters. through careful choice of appropriate model forms, such statistical regression approaches allow an improved comparison of gene expression profiles, and may provide an approach for the greater understanding of common regulatory mechanisms between genes. this approach then allows the direct comparison and clustering of the shapes of response patterns between genes and potentially enables a greater exploration and interpretation of the biological processes driving gene expression. the vast quantities of gene expression profiling data produced in microarray studies, and the more precise quantitative pcr, are often not statistically analysed to their full potential. five-day profiles were modelled using the biologically-oriented, critical exponential curve, y = a + rt + ε. this non-linear regression approach allowed the expression patterns for different genes to be compared in terms of curve shape, time of maximal transcript level and the decline and asymptotic response levels. these approaches have application to the modelling and greater interpretation of profiles obtained across a wide range of platforms, such as microarrays. applying the regression modelling approach to microarray-derived time course data allowed 11% of the escherichia coli features to be fitted by an exponential function, and 25% of the rattus norvegicus features could be described by the critical exponential model, all with statistical significance of p <. "split-line" or "broken-stick" regression identified the initial time of gene up-regulation, enabling the classification of genes into those with primary and secondary responses. three distinct regulatory patterns were identified for the five genes studied.
 the vast quantities of gene expression profiling data produced in microarray studies, and the more precise quantitative pcr, are often not statistically analysed to their full potential. previous studies have summarised gene expression profiles using simple descriptive statistics, basic analysis of variance and the clustering of genes based on simple models fitted to their expression profiles over time. we report the novel application of statistical non-linear regression modelling techniques to describe the shapes of expression profiles for the fungus agaricus bisporus, quantified by pcr, and for e. coli and rattus norvegicus, using microarray technology. the use of parametric non-linear regression models provides a more precise description of expression profiles, reducing the "noise" of the raw data to produce a clear "signal" given by the fitted curve, and describing each profile with a small number of biologically interpretable parameters. this approach then allows the direct comparison and clustering of the shapes of response patterns between genes and potentially enables a greater exploration and interpretation of the biological processes driving gene expression.  quantitative reverse transcriptase pcr-derived time-course data of genes were modelled. "split-line" or "broken-stick" regression identified the initial time of gene up-regulation, enabling the classification of genes into those with primary and secondary responses. five-day profiles were modelled using the biologically-oriented, critical exponential curve, y = a + rt + ε. this non-linear regression approach allowed the expression patterns for different genes to be compared in terms of curve shape, time of maximal transcript level and the decline and asymptotic response levels. three distinct regulatory patterns were identified for the five genes studied. applying the regression modelling approach to microarray-derived time course data allowed 11% of the escherichia coli features to be fitted by an exponential function, and 25% of the rattus norvegicus features could be described by the critical exponential model, all with statistical significance of p <.  the statistical non-linear regression approaches presented in this study provide detailed biologically oriented descriptions of individual gene expression profiles, using biologically variable data to generate a set of defining parameters. these approaches have application to the modelling and greater interpretation of profiles obtained across a wide range of platforms, such as microarrays. through careful choice of appropriate model forms, such statistical regression approaches allow an improved comparison of gene expression profiles, and may provide an approach for the greater understanding of common regulatory mechanisms between genes. 