ISRN Computational BiologyVolume 2013 (2013), Article ID 719138, 11 pageshttp://dx.doi.org/10.1155/2013/719138

Research Article

Effective Single-Step Posttranscriptional Dynamics Allowing for a Direct Maximum Likelihood Estimation of Transcriptional Activity and the Quantification of Sources of Gene Expression Variability with an Illustration for the Hypoxia and TNFα Regulated Inflammatory Pathway

T. D. Frank,¹ A. J. F. Collins,² and A. Cheong²

¹Center for the Ecological Study of Perception and Action, Department of Psychology, University of Connecticut, 406 Babbidge Road, Storrs, CT 06269, USA²Systems Biology Ireland (SBI), University College Dublin, Belfield, Dublin 4, Ireland

Received 31 May 2013; Accepted 31 July 2013

Academic Editors: S. Kalyana-Sundaram and B. Oliva

Copyright © 2013 T. D. Frank et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Data analysis methods for estimating promoter activity from gene reporter data frequently involve the reconstruction of the dynamics of unobserved species and numerical search algorithms for determining optimal model parameters. In contrast, we argue that posttranscriptional dynamics effectively behave like a singlestep stochastic process when gene expression variability is relatively low and, half-lives of the unobserved species are relatively small compared to characteristic observation time scales. In this case, by means of maximum likelihood estimators, for which analytical expressions exist, transcriptional activity of gene promoters can be estimated directly from observed gene reporter data without the need for numerical search algorithms and the reconstruction of unobserved variables. In addition, the model-based data analysis approach yields a single variable that measures the effective strength of the sources that give rise to gene expression variability. The approach is applied to conduct a model-based analysis of the inflammatory pathway under hypoxia condition and stimulation with tumor necrosis factor alpha in HEK293 cells.

1. Introduction

A problem in the field of computational biology is how to model and determine quantitatively promoter activity from observed reporter data. Deterministic approaches suggest that when the activity of a promoter is constant over a period of time, then reporter data should be linearly increasing (see Section 2.1). The steepness of the increase is proportional to the activity of the promoter. Linear regression analysis may be applied to determine the rate of increase [1]. This deterministic perspective is limited in its scope, which becomes clear when considering stochastic approaches as alternatives. For example, deterministic approaches regard gene expression fluctuations as errors. In contrast, according to stochastic approaches gene expression fluctuations indicate that the transcriptional machinery functions properly because the machinery is based on biochemical reactions that are stochastic in the very nature. Stochastic accounts, for example, based on chemical Langevin equations, provide a mathematical framework to address both the deterministic and stochastic components of promoter activity [2]. In the literature, modeling of gene transcription often starts at the promoter level with the transcription event [3–6]. Accordingly, mRNA is produced at a certain rate . Subsequently, mRNA is translated into proteins at a rate . The proteins are finally exported out of the cell with an export rate . Measurement devices of gene reporter systems typically measure the protein abundance in this extracellular space [3, 7–10]. Taken together, the post-transcriptional dynamics is a multistep process that involves different types of entities such as the mRNA, the intracellular proteins, and the extracellular proteins. More detailed models for post-transcriptional dynamics featuring even more entities and species [11] and the impact of repressors [12] have been discussed in the literature as well. Several studies have exploited stochastic modeling of post-transcriptional dynamics by means of chemical Langevin equations in order to estimate promoter activity from gene reporter data (e.g., [3, 5]). In these studies, post-transcriptional dynamics was modeled in terms of the aforementioned multistep process. Consequently, it was necessary to reconstruct the dynamics of hidden, unobserved entities, for example, the mRNA and the intracellular proteins. This in turn makes it difficult to obtain analytical expressions for estimators of the relevant model parameters. Parameter spaces need to be searched for the optimal parameters by means of numerical methods [3]. Although nowadays numerical methods have reached high standards, the issue is to develop complementary approaches that are based on analytical expressions for parameter estimators.

The impact of hypoxia and tumor necrosis factor-alpha (TNF) on cellular responses is of particular interest for cell biology. While oxygen is of vital importance for all aerobic organisms, oxygen levels have to be maintained between certain thresholds. That is, oxygen homeostasis is required on the cellular level for proper cell functioning [13–16]. Diseases such as various types of cancers have been associated with the failure to regulate oxygen homeostasis in the cell [14, 17]. Moreover, the dysfunction of the oxygen homeostasis has adverse consequences for prenatal development [13]. Likewise, high oxygen levels (i.e., hyperoxia rather than hypoxia) immediately after birth may have negative consequences for long term development of preterm infants [18, 19]. Finally, oxygen homeostasis and inflammation are closely linked. Hypoxic tissue conditions increase the chance of inflammation and, vice versa, inflammatory diseases are likely to lead to hypoxic tissue conditions [20]. In this context, the hypoxia induced factor (HIF) has been identified as master regulator of oxygen homeostasis [14, 17, 21]. However, TNF levels seem to play a crucial role as well because it has been shown that hypoxic conditions induce increased TNF levels [9, 22–25]. Therefore, the transcriptional response of the human cyclooxygenase-2 (COX2) promoter to changes in HIF and TNF levels has been studied and a synergistic regulation was found [1]. COX2, in turn, is known to be involved in the cellular response mechanism to inflammation (e.g., [26]). Apart from its central role for inflammatory diseases, the TNF-NFκB pathway seems to be involved in other signaling networks such as the regulatory network for cardiac myocyte growth [27].

The goal of the present study is to demonstrate that post-transcriptional dynamics effectively behaves like a single-step stochastic process when gene expression variability is relatively low and the half lives of the unobserved entities, that is, the mRNA and the intercellular proteins, are short relative to the measurement period. The latter issue means that mRNA and proteins are degraded relatively quickly in the cell. In fact, this is frequently the case: half lives are typically below one hour, whereas reporter data may be collected over periods of several hours.

In Section 2, we will derive an effective stochastic single-step model for post-transcriptional dynamics that comes with analytical expressions for maximum likelihood estimators. We will demonstrate how these estimators can be used to reconstruct the deterministic and stochastic components of promoter activity from gene reporter data. The results of Section 2 will be applied in Section 3. In Section 3, we will analyze the activity of the COX2 promoter when stimulated via two pathways: the TNF pathway and the hypoxia pathway. We will show that a particular model prediction, namely, the proportionality of the mRNA production rate and the strength of the promoter activity fluctuations, holds when stimulated via either of the two pathways. However, the prediction is violated when the COX2 promoter is stimulated by both pathways simultaneously.

2. Effective Single-Step Model for the Reconstruction of Promoter Activity and Gene Expression Volatility from Reporter Data

2.1. Single Cell Level

2.1.1. The Stochastic Single Cell Multistep Model

Let denote the number of mRNA molecules at time (in what follows molecule numbers are assumed to be large such that they can be treated as real numbers). Let and denote the protein numbers in the cell and in the extracellular space. In the deterministic case, the post-transcriptional model sketched in the introduction reads as follows:

Accordingly, mRNA is produced at a rate and degraded with a degradation constant . Intracellular proteins are produced by the activity of mRNA (in a nonconsumable way) at a rate and degraded with a degradation constant . Proteins are exported outside the cell at a rate such that in the evolution equation for the term occurs. In contrast, the evolution equation for the extracellular protein, features the same term with the plus sign. Note that the degradation of proteins in the extracellular space is considered as negligible, which is the reason why there is no decay term in (2).

When and reach their stationary values and , respectively, we have

with

Equation (5) is an important result because it shows that the parameter is proportional to the transcription rate . That is, if can be estimated from experimental data under different conditions and if these conditions do not affect the parameters , , , and of the post-transcriptional system, then variations of reflect variations of the promoter activity . Furthermore, note that the derivation of (4) proves the claim made in the introduction that under certain conditions deterministic approaches predict that reporter data should linearly increase as a function of time and that the steepness of the increase is proportional to the promoter activity.

The chemical Langevin equations accounting for the stochastic aspects of the biochemical reactions (1) and (2) read [2, 3] as

Here , , and denote statistically independent Wiener processes that (without loss of generality) have variance 2. The volatilities , , and are given by [2, 3]

2.1.2. Identifying the Relevant Fluctuating Forces in the Low Variability Limit

We consider the case in which the system is almost deterministic. That is, we assume that the fluctuating forces described by the terms involving the Wiener processes have a small impact on the evolution of , , and . In this case, we substitute the stationary values defined by (3) into the volatilities defined by (7). Thus, the volatilities in (7) become independent of the mRNA and protein numbers:

Such constant volatilities are known to describe so-called additive noise systems, which is the reason why we have added the subindex “”. Exploiting the additive noise volatilities, the stochastic model defined by (6) can be written like

The error terms , , and account for effects due to moderate and high volatility.

2.1.3. Identifying the Dominant Part in the Case of Fast Relaxing Intercellular mRNA and Protein Dynamics

The multivariate Langevin equation defined by (6) describes a three-variable Markov diffusion process. It is known that if we consider only one variable of a multivariate Markov process, then the process in that variable may become non-Markovian in nature [28]. This is indeed the case for the process described by (6). If we solve the evolution equations for and formally and substitute the result into the evolution equation of , then is described by a non-Markovian process (for a similar workedout example see [28]). However, if the decay constants and are relatively large, such that and decay quickly to their stationary values, then is approximately given by a Markov process, which can be shown using the method of adiabatic elimination [29–33]. More precisely, let us consider the stochastic process

where , , and are constants and is a Wiener process. It can be shown that in the limiting case variations of from the fixed point are determined by the Wiener process like [29, 30]

For sake of clarity, we note that (13) reads in integral form

Equation (13) describes an approximation to the exact solutions of (12). The approximation (13) is called the adiabatic elimination of because (13) can be used to replace the expression in other formula. According to the adiabatic elimination method, from the stochastic evolution equations (9) and (10) we obtain

The error terms in (15) and (16) account for non-Markovian effects that occur when the limiting case , does not hold and account for effects due to moderate and high volatility. Substituting (15) into (16) and substituting the result into (11), we obtain

with defined by (5). The error term (volat., nonMarkov) comprises effects arising from the adiabatic elimination of and , on the one hand, and the impacts of moderate and high volatility, on the other. In particular, it includes the error terms occurring in (15) and (16). The net impact of the three independent Wiener processes , , and can be described by means of a single Wiener process when the volatility is adjusted appropriately. We have

with

Note that the equivalence in (18) holds in probability (as indicated). Equation (19) can be expressed in terms of the model parameters by substituting (8) into (19):

Having these results at our disposal, (17) can equivalently be expressed as

The stochastic model defined by (21) provides us with two key parameters. The parameter is a measure for promoter activity (as argued previously). The term describes sources that give rise to fluctuations in gene expression levels. Moreover, the squared expression []² formally corresponds to the diffusion coefficient of a Markov diffusion process when we put the error term equal to zero. Consequently, the parameter and the expression []² are two alternative measures that quantify the strength of sources that lead to gene expression variability.

2.1.4. On the Effective Stochastic Single-Step Process

Note that (20) can be written like with . Furthermore, note that for a classical chemical reaction with a constant production rate the chemical Langevin equation holds [2]. Equation (21) deviates from this form because we have and (21) exhibits the error term . Irrespectively, (21) describes a single-step stochastic process. That is, we have reduced the multistep process model to a model of a single-step process.

We compare next the effective single-step process defined by (21) with a Wiener process subjected to a constant drift :

In (22) the parameter denotes the volatility of the process. Equation (22) exhibits an exact analytical solution for the transition probability density (see, e.g., [34, 35]). For data given at discrete time points with the maximum likelihood estimators (MLEs) of and of are defined by

Here, the total observation time is . Equation (23) can be obtained from the aforementioned exact solution of the transition probability density.

As argued above, for low volatility and when the decay rates and are sufficiently large, then the error term in (21) should have only a small impact on the evolution of . In this case, the Wiener process with drift defined by (22) is a good approximation of the single-step post-transcriptional dynamics model (21). Likewise, and are useful estimates of the parameters and defined by (5) and (20).

In order to illustrate this line of argument, we simulated the original three-variable Markov diffusion process defined by (6) and (7) and subsequently estimated and from the simulated protein data using (23). We assumed that mRNA and the protein have half lives of 30 min, whereas observations are made every 2 hours (i.e., ) over a sufficiently long period. In this case, mRNA and intracellular proteins decay relatively quickly on the time scale of the observations. In this sense, the decay rates and (and consequently ) are relatively large. Explicitly, for a half-life of 30 min, the decay constant is about 1.4/h (note: the decay constant can be computed from the half-life using the formula: decay constant = log(2)/half life). Moreover, we used in the simulations relative large production rates of proteins/h. We found that for these rates the variability is small relatively to the mean values as required for the validity of our approach.

Figure 1 shows the sample mean values of and for a sample of 10 repetitions obtained for several values of . It also shows, for comparison purposes, the analytical results for and defined by (5) and (20). We see that and correspond fairly well to the mean values of the estimates and . Figure 1 illustrates that the effective single-step model for post-transcriptional dynamics given by the Wiener process with drift defined by (22) is a useful approximations of the original stochastic three-variable model in (6) and (7).

(a)

(b)

Figure 1: Drift and diffusion parameters as used in the stochastic stimulation of the multistep model defined by (6) and (7) and as estimated from the approximate single-step model defined by (22) as functions of the transcription rate . (a) Drift coefficients and computed from (5) and (23), respectively. (b) Diffusion coefficients and as obtained from (20) and (23). Sample mean values are presented. Error bars represent SDs. Number of repetitions: 10. Parameters: , , /h, /h, and  h. Stochastic simulations were performed by means of an Euler forward method with computational time step of 0.01 h, and trajectories were computed in the interval from 0 to 100 h.

2.2. Population Level

Reporter data are frequently recorded from cell populations rather than from single cells. Let us consider a population of cells. The cell population is exposed to a stimulus such that the cells produce stimulus-specific proteins and export them into the extracellular space. The number of proteins deposited by cell into the extracellular space will be denoted by . We assume that the measurement devise records a physical quantity that is proportional to the total number of stimulus-specific proteins in the extracellular space at time :

In (24) the parameter is a proportionality factor. The sum in (24) may be interpreted as the total number of proteins that have been deposited by the cell population in the extracellular space up to time .

Let us evaluate the right-hand side of (24). As argued previously, each cell satisfies an evolution equation of the form (21) such that

where denote statistically independent Wiener processes with variance 2. Substituting (25) into (24), we obtain

with

The sum of the statistically independent Wiener processes can be replaced by a single Wiener process like

where is a Wiener process with variance 2. Exploiting (28) and (29), (26) becomes

with . Alternatively, with the help of (20), we obtain

Just as in our discussion of (21), we point out that the stochastic model defined by (30) features two key parameters that are essential for the characterization of promoter activity. The parameters and reflect deterministic and stochastic aspects of promoter activity. The parameter is a measure proportional to the promoter activity (i.e., the rate of transcription initiation). If is estimated from experimental data under various conditions that do not affect the parameters of the post-transcriptional system (, , , , , and ), then variations observed in inform us about variations in that are induced by these different conditions. The parameter or alternatively the squared parameter are quantitative measures for the strength of the sources that give rise to gene expression fluctuations. Again, by estimating from experimental data under varying conditions, we can identify how these conditions affect the strength of the sources that lead to gene expression variability.

In this context, it is important to point out that and are not independent parameters. By comparing (27) and (31) for and , we recognize that they both depend on the same set of parameters: , , , , , , and . This property can be exploited to obtain additional insights into how deterministic and stochastic aspects of cellular signaling are regulated by extracellular stimuli. We will return to this issue in Section 3.

Equation (30) may be compared with a Wiener process with constant drift

whose coefficients can be estimated from data using the aforementioned MLEs:

In order to demonstrate that on the population level the multistep stochastic model defined by (6) and (7) with (24) behaves under appropriate conditions (low noise and fast decays of mRNA and intracellular proteins) like the Wiener process with drift defined by (32), we conducted a series of simulations. We simulated a small population of cells () using (6) and (7) for each variable and (24) to obtain the measurement variable . We downsampled the trajectory to obtain a sequence at discrete time points . We fitted the computer generated reporter data to the Wiener process with drift (32) by means of the MLEs of and . Finally, we compared and with the analytical expressions and given by (27) and (31). The results of this procedure are reported in Figure 2.

(a)

(b)

Figure 2: Drift and diffusion parameters as used in the stochastic stimulation of the multistep model defined by (6), (7), and (24) and as estimated from the approximate single-step model defined by (32) as functions of the transcription rate . (a) Drift coefficients and computed from (27) and (33), respectively. (b) Diffusion coefficients and as obtained from (31) and (33). Sample mean values are presented. SDs for were smaller than 1, which is the size of the symbol. SDs for are shown as error bars. Number of repetitions: 10. Parameters as in Figure 1. Additional parameters: and . Stochastic simulations were performed by means of an Euler forward method with computational time step of 0.01 h, and trajectories were computed in the interval from 0 to 100 h.

Inspection of Figure 2 reveals that the estimates and of the Wiener process model with drift correspond well with the analytical expressions for and . This indicates that for our choice of parameters the error term in (30) affects the evolution of only to a small degree. Consequently, for appropriate parameters the multistep post-transcriptional process for the cell population as defined by (6), (7) and (24) behaves effectively like a single-step stochastic process of the form (32).

3. Analysis of Synergistic Transcriptional Activity of the COX2 Promoter Using a Single-Step Chemical Langevin Model

We conducted a model-based analysis of COX2 transcriptional activity data reported by Bruning et al., 2012 [1].

3.1. Sketch of the Experimental Procedure

In Bruning et al. [1] the activity of the COX2 promoter in HEK293 cell populations was measured when exposing cells to either hypoxia or tumor necrosis factor-alpha (TNF) or both. To this end, HEK293 cells were transfected with a Gaussia luciferase reporter under the control of a fragment of the COX2 promoter. The transfected cells were exposed to either 21% O₂ (normoxia) or 1% O₂ (hypoxia) with or without stimulation by 1 ng/mL TNF such that promoter activity was measured under four different conditions in total. The Gaussia luciferase was measured every 3 h in a 12 h period. Data from 8 replicates were recorded. For more details, see Bruning et al. [1].

3.2. Model-Based Analysis

Comparing (27) and (31) we see that both and depend linearly on the production rate . Consequently, variations in due to experimental manipulations should affect the deterministic and stochastic components of the post-transcriptional system as quantified by and to the same degree provided that all other parameters (, , , and ) of the post-transcriptional machinery are invariant across the experimental manipulations. In order to address this issue from a slightly different angle, we substitute (27) into (31) such that

with depending on , , , and as defined in Section 2.1.4. This implies that the ratio is independent of :

For each replicate we calculated and from the five data points related to the time points , 3, 6, 9, 12 h using (33). For each replicate we also calculated the ratio . If the error term in the single-step stochastic model (30) for the post-transcriptional dynamics is relatively small, then the model (30) corresponds approximately to the Wiener process with drift defined by (32). In this case, the ratios (as estimated by ) obtained for the four different conditions should correspond to the ratios of the Wiener process model (32) and consequently should be on the same level (see (35)).

3.3. Results

Figure 3 reports sample mean values of , and the ratio . From Figure 3(a) it follows that the promoter activity (as estimated ) was minimal under baseline condition (normoxia with 0 ng/mL TNF) with (measured in arbitrary units, i.e., relative light units per hour; see [1]) and maximal when both pathways were stimulated (hypoxia with 1 ng/mL TNF) with . The single stimulation conditions produced transcriptional activity levels in the medium range of (normoxia with 1 ng/mL TNF) and again (hypoxia with 0 ng/mL TNF). In particular, under dual stimulation the activity was higher () than the sum of the activity levels observed under single stimulation (), which indicates that the COX2 promoter activity was synergistically regulated [36] and is consistent with the results obtained by Bruning et al., 2012 [1] using linear regression analysis. Two-way ANOVA analysis was conducted to test the statistical significance of the synergy effect [37]. Both main effects were statistically significant. That is, decreasing the oxygen level from 21% (normoxia) to 1% (hypoxia) increased significantly the transcription rate as measured by the population-level drift coefficient (). Likewise, stimulation with 1 ng/mL TNF induced a significant increase of (, ). Most importantly, the interaction effect was statistically significant indicating that the increase of the transcriptional activity due to TNF stimulation was stronger under hypoxia than under normoxia (, ).

(a)  

(b)

(c)

Figure 3: Transcription rate as measured by the drift coefficient (a), strength of the fluctuating force inducing gene expression variability as measured by the diffusion coefficient (b), and ratio (c) observed under four different stimulation conditions. Error bares indicated standard errors of the means.

Figure 3(b) reveals that the measure for the strength of the sources of gene expression variability (as estimated by ) followed a similar pattern. The sources had minimal strength under the baseline condition and maximal strength under the dual stimulation condition. A two-way ANOVA revealed no significant effects. Since the sample size was relatively small, several non-parametric sign tests were conducted as an alternative to the ANOVA. We compared the baseline and the two single stimulation conditions using three sign tests. No significant differences were found. We tested a contrast given by the average of the three aforementioned conditions versus the dual stimulation condition. The sign test showed that the diffusion coefficient was significantly higher for the dual stimulation condition than the average of the three other stimulation conditions ().

Figure 3(c) reports the sample mean values of the ratios as estimated by . The baseline condition and the two conditions in which only one pathway was activated (normoxia with 1 ng/mL TNF or hypoxia with 0 ng/mL TNF) exhibited approximately the same ratios . That is, although the promoter activity was higher in the conditions in which one pathway is activated than in the baseline condition (Figure 3(a)), the ratios of the diffusion coefficients versus the drift parameters were found to be approximately on the same level. This is because the diffusion coefficients were larger in the conditions in which one pathway was stimulated than in the baseline condition (Figure 3(b)), just as the drift parameters were larger in the conditions in which one pathway was stimulated. The dual stimulation condition violated the prediction defined by (35). Under dual stimulation, the ratio was much higher than the ratios obtained in the other three conditions. Although both and were maximal in the dual stimulation condition, if we compare, for example, the baseline condition with the dual stimulation condition, then and increased differently and, in particular, did not increase by the same percentage value as one would expect on the basis of (35). The two-way ANOVA showed no statistically significant effects. Paired sign tests between the baseline condition and the two single stimulation conditions confirmed that these three conditions showed no significant differences. However, a paired signed test showed that the ratio for the dual stimulation was significantly higher than the average ratio for the three other conditions (). We will return to this issue in Section 4.

3.4. Model Checking

Model checking was conducted by testing to what extent the approximation of a Wiener process (32) applied to the experimental data. Equation (32) can be written like

The increments are normally distributed with variance and reflect white noise. Consequently, they are statistically independent (see, e.g., [35]). From (36) it follows that the residuals defined by

are proportional to the increments (as indicated) and, in doing so, are statistically independent. From the time points , 3, 6, 9, 12 h, we calculated four residuals , 2, 3, 4. Model checking typically is conducted by calculating lag- autocorrelation coefficients and comparing them with the threshold , where is the number of data points [38]. In our experiment, only 4 residuals were given. Therefore, only the lag-1 autocorrelation coefficient could be tested. Rather than testing individual samples we tested the mean lag-1 autocorrelation obtained from samples. The threshold for the sample mean lag-1 autocorrelation is . Table 1 presents the threshold value and the lag-1 autocorrelation coefficients observed in the experiment by Bruning et al. [1] under the four different stimulation conditions.

Table 1: Correlation coefficients of lag-1 of the residuals defined by (37) under the four stimulation conditions. In the first column, the threshold given indicates the statistically significant correlations (at a significance level of 0.05). Values below that threshold are considered not to be significant.

Threshold with ,
Baseline (Normoxia and 0 ng/mL TNF)
TNFα stimulation (Normoxia and 1 ng/mL TNF)
Hypoxia stimulation (Hypoxia and 0 ng/mL TNF)
Dual (Hypoxia and 1 ng/mL TNF)

0.35
0.02
0.22
0.10
0.35

The correlation conditions obtained for all four stimulation conditions satisfied the white noise threshold indicating that the correlations were not statistically significant. However, the dual stimulation condition with the larger transcription rate exhibited a lag-1 correlation value that was right at the threshold.

4. Discussion

We argued that under certain conditions, multistep post-transcriptional dynamics effectively is captured by a single-step chemical Langevin equation that corresponds to a two-parametric Wiener process with drift. The drift parameter reflects the transcriptional activity, while the volatility parameter can be regarded as a bulk measure for the strength of fluctuating forces leading to gene expression variability. In order to derive the single-step model, we made use of adiabatic elimination procedure for stochastic processes suggested by Haken [29]. Numerical simulation supported the validity of our approach.

Interestingly, the model predicts that transcriptional activity and gene expression variability are not two independent parameters. Rather they are two different mappings of the same set of underlying biochemical parameters. Therefore, variations in one of the biochemical parameters by experimental manipulations should in general induce variations in both the transcriptional activity and the amount of gene expression variability.

As pointed out in the previous sections (Sections 2.1.4 and 2.2), the approach yields approximately good results under certain conditions. The error terms in (9)–(11) reflecting the impact of nonadditive noise should have only a negligible impact. This is the case when the mRNA production rate r is relatively large (Section 2.1.4). There should be a clear time scale separation (fast degradation of mRNA and proteins). If so, the approximate solution via adiabatic elimination is close to the exact solution and the error terms reflecting non-Markovian processes make only small contributions. Future efforts may be devoted to derive effective single step models under less severe assumptions. For example, it is known that, in principle, the adiabatic elimination method can also be applied to stochastic processes exhibiting nonadditive fluctuating forces [29]. While in Section 2.1.2 we introduced (9)–(11) as departure point for the adiabatic elimination method assuming that the error terms are relatively small such that they do not conflict with the elimination technique, a more general analysis may use (6) and (7) as starting point for the adiabatic elimination method.

The model presented in Section 2 captures several sources of variability related to the biochemical processes of the post-transcriptional dynamics that are stochastic in nature. However, there are other potential sources of variability that have not been included and that might be included in generalization of the single-step model. For example, in Section 2.2 we assumed that the population is composed of cells with identical parameters. In particular, the transcription rate parameter was assumed to be the same for all cells of a population subjected to a particular stimulation. However, due to cell-to-cell variability cell parameters might be distributed rather than fixed across the whole population. In particular, the degree of cell-to-cell variability may be affected by experimental manipulations such that the parameter distribution may change. Likewise, in Section 2.1 we assumed that the single cell parameters are constant over time. For example, the transcription rate parameter was conceptualized as a constant. However, single cell parameters may fluctuate over time. Cell-to-cell variability as well as fluctuating single cell parameters is likely to inflate the observed gene expressions variability relative to a transcriptional machinery that does not exhibit these effects. Finally, the model does not account for feedback loops (as addressed, e.g., in [21]) from the post-transcriptional side to the signaling network regulating the transcriptional activity. Such positive or negative auto-regulation mechanisms may inflate or suppress fluctuations in gene expression levels.

In our model, we only distinguished between intracellular and extracellular proteins. We did not take cellular compartments into account. As mentioned in the introduction, more detailed models for post-transcriptional dynamics have been studied in the literature. Such more detailed models as well as models involving cellular compartments require in general a more sophisticated mathematical model that involves a larger number of coupled differential equation as considered in Section 2.1.1. In this context, it is worth mentioning that the adiabatic elimination method can be applied to arbitrarily large sets of coupled evolution equations provided that the dynamical system as a whole features an appropriate time scale separation. For example, when the processes within compartments take place on relative short time scales, then in principle the adiabatic elimination technique applies and predicts that the transcriptional machinery including the comportment dynamics effectively behaves like a single-step process.

We applied the model-based analysis method to transcriptional activity data reported from the COX2 promoter under synergistic regulation by hypoxia and TNF [1]. In fact, the MLE parameter estimation method produced consistent results with the linear regression analysis conducted in Bruning et al. [1]. In particular, the model-based analysis supported the notion of a synergistic regulation of COX2 activity by hypoxia and TNF. Under the assumption that the stimulation type does not affect the biochemical parameters of the post-transcriptional machinery, we expected to see that the pattern in the volatilities mimics the pattern in the transcription rates. In fact, qualitatively, the same patterns were found (compare Figures 3(a) and 3(b)). However, quantitatively, only the single stimulation conditions exhibited the same ratio as the baseline condition. The dual stimulation condition exhibited a ratio that was about four times higher than for the three other stimulation conditions. Future research might be devoted to identify the reason for this observation. As discussed previously, there are various conditions that may lead to a violation of invariance of the “noise to signal ratio” . It might be the case that the cell-to-cell variability was increased under dual stimulation relative to the other stimulation conditions. Alternatively or in addition, the transcription rate might be subjected to fluctuations to a greater extent under dual stimulation. Furthermore, from (5) and (20) it follows that the transcription rate and the diffusion coefficient depend on the factor

However, depends linearly on this factor, whereas involves a component that depends in a quadratic way on the factor. An increase in this factor may result in an increase of the ratio . Consequently, another speculative explanation for the observed increased ratio under dual stimulation is that dual stimulation affected some of the biochemical constants entering into the factor defined by (38).

In the context of the study by Bruning et al. [1] as part of pilot work, cells were also stimulated by 0.1 ng/mL TNF (rather than 1 ng/mL TNF) under both normoxia and hypoxia. Only 3 samples were recorded (rather than 8). The data was not published in Bruning et al. [1]; however, it is available to the authors. In order to obtain insights into the emergence of the synergy effect on the transcription rate , we estimated the model parameter for these samples using the protocol described in Section 3. Table 2 presents the sample means values reported in Section 3 together with the sample means values for the additional conditions.

Table 2: Sample means of estimated drift parameters (transcription rates) under various experimental conditions. SDs in brackets are SDs from actual sample sizes (i.e., SDs calculated before mean substitution of missing data). Data presented in arbitrary units (relative light units per hour).

Normoxia21% oxygen
Hypoxia1% oxygen

0 ng/mL TNF
0.11 (0.05)
0.30 (0.05)

0.1 ng/mL TNF
0.27 (0.01)
0.45 (0.14)

1.0 ng/mL TNF
0.30 (0.13)
0.87 (0.25)

As can be seen from Table 2, when comparing the 0 ng/mL and the 0.1 ng/mL TNF stimulation conditions, the transcription rates increased due to stimulation with TNF and due to decreasing the oxygen level. However, there was no dramatic synergy effect. Both in the 0 ng/mL and the 0.1 ng/mL TNF stimulation conditions, lowering the oxygen level increased the activity by about 0.20 units. In contrast, when comparing the 0.1 ng/mL and the 1.0 ng/mL TNF stimulation conditions, the lowering of the oxygen level increased the activity by about 0.20 units for the 0.1 ng/mL TNF stimulated cells, whereas it increased the activity by almost 0.50 units for the 1.0 ng/mL TNF stimulated cells, suggesting that the transcriptional activity was synergistically regulated. That is, Table 2 suggests that the synergy effect discussed in Section 3 does not emerge at relative low TNF doses. It emerges in the range of higher TNF doses. We conducted an overall 2 (normoxia versus hypoxia) by 3 (three TNF levels) ANOVA on the drift parameter using mean substitution of the missing data. The interaction was significant (, ) indicating the presence of a synergy effect, as expected from the results presented in Section 3. In order to locate the interaction effect, we conducted two 2-by-2 ANOVA. The first 2-by-2 ANOVA considered only the 0 ng/mL and the 0.1 ng/mL TNF stimulation conditions. The interaction effect was not significant. The second 2-by-2 ANOVA considered only the 0.1 ng/mL and the 1.0 ng/mL TNF stimulation conditions. The interaction effect was significant (, ). These outcomes of the hypothesis testing procedure support our hypothesis. However, since we used mean substitution of missing data, future work is needed to verify this result.

Although in the application only a single transcription factor is explicitly mentioned, the modeling approach presented in Section 2 can be applied to transcriptional activity regulated by several transcription factors. To this end, the transcriptional activity parameter is considered as a function of the transcription factor concentrations. Explicit models, for example, using a thermostatistical approach have been proposed for that purpose (see [4, 36] and references therein). Unknown parameters of such thermostatistical models or functional dependencies between and transcription factor concentrations can be estimated when (in replacement for ) is estimated for dose responses of the set of transcription factors of interest.

Acknowledgment

A. J. F. Collins and A. Cheong were supported by Science Foundation Ireland (Grant no. 06/CE/B1129).

References

U. Bruning, S. F. Fitzpatrick, M. Birtwistle, T. Frank, C. T. Taylor, and A. Cheong, “NFκB and HIF display synergistic behaviour during hypoxic inflammation,” Cellular and Molecular Life Sciences, vol. 69, no. 8, pp. 1319–1329, 2012.

D. J. Wilkinson, Stochastic Modelling for Systems Biology, Chapman & Hall/CRC, London, UK, 2006.

B. Finkenstädt, E. A. Heron, M. Komorowski et al., “Reconstruction of transcriptional dynamics from gene reporter data using differential equations,” Bioinformatics, vol. 24, no. 24, pp. 2901–2907, 2008.

T. D. Frank, A. Cheong, M. Okada-Hatakeyama, and B. N. Kholodenko, “Catching transcriptional regulation by thermostatistical modeling,” Physical Biology, vol. 9, no. 4, Article ID 045007, 2012.

C. V. Harper, B. Finkenstädt, D. J. Woodcock et al., “Dynamic analysis of stochastic transcription cycles,” PLoS Biology, vol. 9, no. 4, Article ID e1000607, 2011.

P. J. Ingram, M. P. H. Stumpf, and J. Stark, “Network motifs: structure does not determine function,” BMC Genomics, vol. 7, article 108, 2006.

J. Alam and J. L. Cook, “Reporter genes: application to the study of mammalian gene transcription,” Analytical Biochemistry, vol. 188, no. 2, pp. 245–254, 1990.

A. J. Millar, S. R. Short, N.-H. Chua, and S. A. Kay, “A novel circadian phenotype based on firefly luciferase expression in transgenic plants,” Plant Cell, vol. 4, no. 9, pp. 1075–1087, 1992.

J. F. Schmedtje Jr., Y.-S. Ji, W.-L. Liu, R. N. DuBois, and M. S. Runge, “Hypoxia induces cyclooxygenase-2 via the NF-κB p65 transcription factor in human vascular endothelial cells,” Journal of Biological Chemistry, vol. 272, no. 1, pp. 601–608, 1997.

A. J. Millar, I. A. Carre, C. A. Strayer, N.-H. Chua, and S. A. Kay, “Circadian clock mutants in Arabidopsis identified by luciferase imaging,” Science, vol. 267, no. 5201, pp. 1161–1163, 1995.

T. Nakakuki, M. R. Birtwistle, Y. Saeki et al., “Ligand-specific c-fos expression emerges from the spatiotemporal control of ErbB network dynamics,” Cell, vol. 141, no. 5, pp. 884–896, 2010.

S. Kuttykrishnan, J. Sabina, L. L. Langton, M. Johnston, and M. R. Brent, “A quantitative model of glucose signaling in yeast reveals an incoherent feed forward loop leading to a specific, transient pulse of transcription,” Proceedings of the National Academy of Sciences of the United States of America, vol. 107, no. 38, pp. 16743–16748, 2010.

E. Maltepe and O. D. Saugstad, “Oxygen in health and disease: regulation of oxygen homeostasis-clinical implications,” Pediatric Research, vol. 65, no. 3, pp. 261–268, 2009.

G. L. Semenza, “Mechanisms of disease: oxygen sensing, homeostasis, and disease,” The New England Journal of Medicine, vol. 365, no. 6, pp. 537–547, 2011.

R. H. Wenger, “Mammalian oxygen sensing, signalling and gene regulation,” Journal of Experimental Biology, vol. 203, no. 8, pp. 1253–1263, 2000.

G. L. Semenza, “Hypoxia-inducible factors in physiology and medicine,” Cell, vol. 148, no. 3, pp. 399–408, 2012.

J. Pouysségur, F. Dayan, and N. M. Mazure, “Hypoxia signalling in cancer and approaches to enforce tumour regression,” Nature, vol. 441, no. 7092, pp. 437–443, 2006.

R. Deulofeut, A. Critz, I. Adams-Chapman, and A. Sola, “Avoiding hyperoxia in infants ≤1250g is associated with improved short- and long-term outcomes,” Journal of Perinatology, vol. 26, no. 11, pp. 700–705, 2006.

O. D. Saugstad, “Optimal oxygenation at birth and in the neonatal period,” Neonatology, vol. 91, no. 4, pp. 319–322, 2007.

H. K. Eltzschig and P. Carmeliet, “Hypoxia and inflammation,” The New England Journal of Medicine, vol. 364, no. 7, pp. 656–665, 2011.

L. K. Nguyen, M. A. S. Cavadas, C. C. Scholz et al., “A dynamic model of hypoxia-induced factor 1-alpha (HIF-1alpha) network,” Journal of Cell Sciences, vol. 126, pp. 1454–1463, 2013.

C. Culver, A. Sundqvist, S. Mudie, A. Melvin, D. Xirodimas, and S. Rocha, “Mechanism of hypoxia-induced NF-kappaB,” Molecular and Cellular Biology, vol. 30, no. 20, pp. 4901–4921, 2010.

E. P. Cummins, E. Berra, K. M. Comerford et al., “Prolyl hydroxylase-1 negatively regulates IκB kinase-β, giving insight into hypoxia-induced NFκB activity,” Proceedings of the National Academy of Sciences of the United States of America, vol. 103, no. 48, pp. 18154–18159, 2006.

S. F. Fitzpatrick, M. M. Tambuwala, U. Bruning et al., “An intact canonical NF-κB pathway is required for inflammatory gene expression in response to hypoxia,” Journal of Immunology, vol. 186, no. 2, pp. 1091–1096, 2011.

A. C. Koong, E. Y. Chen, and A. J. Giaccia, “Hypoxia causes the activation of nuclear factor κB through the phosphorylation of IκBα on tyrosine residues,” Cancer Research, vol. 54, no. 6, pp. 1425–1430, 1994.

H. Li, J. Alyce Bradbury, R. T. Dackor et al., “Cyclooxygenase-2 regulates Th17 cell differentiation during allergic lung inflammation,” American Journal of Respiratory and Critical Care Medicine, vol. 184, no. 1, pp. 37–49, 2011.

K. A. Ryall, D. O. Holland, K. A. Delaney, M. J. Kraeutler, A. J. Parker, and J. J. Saucerman, “Network reconstruction and systems analysis of cardiac myocyte hypertrophy signaling,” Journal of Biological Chemistry, vol. 287, pp. 42259–42268, 2012.

H. Risken, The Fokker-Planck Equation: Methods of Solutions and Applications, Springer, Berlin, Germany, 1989.

H. Haken, Advanced Synergetics, Springer, Berlin, Germany, 1987.

G. Schöner and H. Haken, “The slaving principle for stratonovich stochastic differential equations,” Zeitschrift für Physik B, vol. 63, no. 4, pp. 493–504, 1986.

B. C. Eu, Kinetic Theory and Irreversible Thermodynamics, John Wiley & Sons, New York, NY, USA, 1992.

T. D. Frank, A. Daffertshofer, P. J. Beek, and H. Haken, “Impacts of noise on a field theoretical model of the human brain,” Physica D, vol. 127, no. 3-4, pp. 233–249, 1999.

T. D. Frank, “A limit cycle oscillator model for cycling mood variations of bipolar disorder patients from cellular biochemical reaction equations,” Communications in Nonlinear Science and Numerical Simulation, vol. 18, pp. 2107–2119, 2013.

T. D. Frank, Nonlinear Fokker-Planck Equations: Fundamentals and Applications, Springer, Berlin, Germany, 2005.

C. W. Gardiner, Handbook of Stochastic Methods For Physics, Chemistry and the Natural Sciences, Springer, Berlin, Germany, 1985.

T. D. Frank, A. M. Carmody, and B. N. Kholodenko, “Versatility of cooperative transcriptional activation: a thermodynamical modeling analysis for greater-than-additive and less-than-additive effects,” PLoS One, vol. 7, no. 4, Article ID e34439, 2012.

B. K. Slinker, “The statistics of synergism,” Journal of Molecular and Cellular Cardiology, vol. 30, no. 4, pp. 723–731, 1998.

P. J. Diggle, Time Series: A Biostatistical Introduction, Clarendon Press, Oxford, UK, 1990.