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

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.


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][4][5][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][8][9][10]. Taken together, the post-transcriptional dynamics is a multistep process that involves different types of entities such 2 ISRN Computational Biology 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, posttranscriptional 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][14][15][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][23][24][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 posttranscriptional dynamics effectively behaves like a singlestep 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 singlestep 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. In the deterministic case, the post-transcriptional model sketched in the introduction reads as follows:

Effective Single-Step
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 = .
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]

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.

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][30][31][32][33]. More precisely, let us consider the stochastic process where > 0, , and > 0 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 ,tot (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 4 ISRN Computational Biology 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 eff when the volatility is adjusted appropriately. We have 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 eff ⋅ eff describes sources that give rise to fluctuations in gene expression levels. Moreover, the squared expression [ eff ] 2 formally corresponds to the diffusion coefficient of a Markov diffusion process when we put the error term equal to zero. Consequently, the parameter eff and the expression [ eff ] 2 are two alternative measures that quantify the strength of sources that lead to gene expression variability.
We compare next the effective single-step process defined by (21) with a Wiener process subjected to a constant drift 0 : In (22) the parameter > 0 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 = Δ ⋅ ( − 1) with = 1, . . . , the maximum likelihood estimators (MLEs) 0, of 0 and of are defined by Here, the total observation time is = − 1 = ( − 1) ⋅ Δ . 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 ,tot 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 posttranscriptional dynamics model (21). Likewise, 0, and are useful estimates of the parameters and eff 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 0, 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., Δ = 2) 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 halflife 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 > 10 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 0, and for a sample of 10 repetitions obtained for several values of ≥ 10. It also shows, for comparison purposes, the analytical results for and eff defined by (5) and (20). We see that and eff correspond fairly well to the mean values of the estimates 0, 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).

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 stimulusspecific 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) (6) and (7) and as estimated from the approximate single-step model defined by (22) (20) and (23) 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 ,eff denote statistically independent Wiener processes with variance 2. Substituting (25) into (24), we obtain ,eff + ,tot (volat., non-Markov) ,tot (volat., non-Markov) , ,tot (volat., non-Markov) .
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) 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 [ ] 2 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 [ ] 2 , 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 ( = 50) 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 0, and [ ] 2 . Finally, we compared 0, and with the analytical expressions and given by (27) and (31). The results of this procedure are reported in Figure 2.
Inspection of Figure 2 reveals that the estimates 0, and [ ] 2 of the Wiener process model with drift correspond well with the analytical expressions for and [ ] 2 . This indicates that for our choice of parameters the error term ,tot in (30) affects the evolution of only to a small degree. Consequently, for appropriate parameters the multistep posttranscriptional process for the cell population as defined by (6), (7) and (24) behaves effectively like a single-step stochastic process of the form (32).

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

Model-Based Analysis.
Comparing (27) and (31) we see that both and [ ] 2 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 [ ] 2 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 [ ] 2 / is independent of : For each replicate we calculated 0, and [ ] 2 from the five data points 1 , . . . , 5 related to the time points = 0, 3, 6, 9, 12 h using (33). For each replicate we also calculated the ratio [ ] 2 / 0, . If the error term in the single-step stochastic ISRN Computational Biology 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 [ ] 2 / 0 (as estimated by [ ] 2 / 0, ) obtained for the four different conditions should correspond to the ratios [ ] 2 / of the Wiener process model (32) and consequently should be on the same level (see (35)). Figure 3 reports sample mean values of 0, , [ ] 2 and the ratio [ ] 2 / 0, . From Figure 3(a) it follows that the promoter activity (as estimated 0, ) was minimal under baseline condition (normoxia with 0 ng/mL TNF ) with = 0.11 (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 = 0.87. The single stimulation conditions produced transcriptional activity levels in the medium range of = 0.30 (normoxia with 1 ng/mL TNF ) and = 0.30 again (hypoxia with 0 ng/mL TNF ). In particular, under dual stimulation the activity was higher ( = 0.87) than the sum of the activity levels observed under single stimulation ( = 0.60), 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 ( (1, 7) = 28.967, = 0.001). Likewise, stimulation with 1 ng/mL TNF induced a significant increase of ( (1, 7) = 96.446, < 0.001). 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 ( (1, 7) = 14.740, < 0.01). Figure 3(b) reveals that the measure [ ] 2 for the strength of the sources of gene expression variability (as estimated by [ ] 2 ) 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 8 ISRN Computational Biology 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. 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 ( < 0.01). Figure 3(c) reports the sample mean values of the ratios [ ] 2 / as estimated by [ ] 2 / 0, . 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 [ ] 2 / . 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 [ ] 2 versus the drift parameters were found to be approximately on the same level. This is because the diffusion coefficients [ ] 2 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 [ ] 2 / was much higher than the ratios obtained in the other three conditions. Although both and [ ] 2 were maximal in the dual stimulation condition, if we compare, for example, the baseline condition with the dual stimulation condition, then and [ ] 2 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 [ ] 2 / for the dual stimulation was significantly higher than the average ratio [ ] 2 / for the three other conditions ( < 0.01). We will return to this issue in Section 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 2Δ 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 = 0, 3, 6, 9, 12 h, we calculated four residuals = 1, 2, 3, 4. Model checking typically is conducted by calculating lagautocorrelation coefficients and comparing them with the threshold 2/ 1/2 , 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 2/( ⋅ ) 1/2 . 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.
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.

Discussion
We argued that under certain conditions, multistep posttranscriptional dynamics effectively is captured by a singlestep chemical Langevin equation that corresponds to a twoparametric 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 autoregulation 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 [ ] 2 / as the baseline condition. The dual stimulation condition exhibited a ratio [ ] 2 / 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" [ ] 2 / . 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 [ ] 2 depend on the factor However, depends linearly on this factor, whereas [ ] 2 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 [ ] 2 / . Consequently, another speculative explanation for the observed increased ratio [ ] 2 / 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.
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 ( (2, 14) = 11.925, = 0.001) 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 ( (1, 7) = 11.078, < 0.05). 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.