# 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

Authors: T. D. Frank; A. J. F. Collins; A. Cheong
Journal: ISRN Computational Biology
 (2013)

---

## 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.


---

## Body

## 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 Section2.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 r. 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 Section2, 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 r 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

Letm(t) denote the number of mRNA molecules at time t (in what follows molecule numbers are assumed to be large such that they can be treated as real numbers). Let p(t) and x(t) 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:
(1)ddtm=r-k·m,ddtp=ω·m-(α+ε)·p,(2)ddtx=α·p.Accordingly, mRNAm(t) is produced at a rate r and degraded with a degradation constant k. Intracellular proteins p(t) 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 p(t) the term -α·p occurs. In contrast, the evolution equation for the extracellular protein, x(t) 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).Whenm and p reach their stationary values ms and ps, respectively, we have
(3)ms=rk,ps=1α+εωkr,(4)ddtx=b.with(5)b=αα+εωkr.Equation (5) is an important result because it shows that the parameter b is proportional to the transcription rate r. That is, if b can be estimated from experimental data under different conditions and if these conditions do not affect the parameters α, ε, ω, and k of the post-transcriptional system, then variations of b reflect variations of the promoter activity r. 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
(6)dm=(r-k·m)dt+σm·dWm,dp=[ω·m-(α+ε)·p]dt+σp·dWp,dx=α·pdt+σx·dWx.HereWm, Wp, and Wx denote statistically independent Wiener processes that (without loss of generality) have variance 2. The volatilities σm, σp, and σx are given by [2, 3]
(7)σm=r+k·m2,σp=ω·m+(α+ε)·p2,σx=α·p2.


#### 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 ofm, p, and x. 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:
(8)σm,a=r,σp,a=ω·rk,σx,a=α·ω·r2(α+ε)k.Such constant volatilities are known to describe so-called additive noise systems, which is the reason why we have added the subindex “a”. Exploiting the additive noise volatilities, the stochastic model defined by (6) can be written like
(9)dm=(r-k·m)dt+σm,a·dWm+em(volat.),(10)dp=[ω·m-(α+ε)·p]dt+σp,a·dWp+ep(volat.),(11)dx=α·pdt+σx,a·dWx+ex(volat.).The error termsem, ep, and ex 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 m(t) and p(t) formally and substitute the result into the evolution equation of x, then x is described by a non-Markovian process (for a similar workedout example see [28]). However, if the decay constants k and k′=α+ε are relatively large, such that m and p decay quickly to their stationary values, then x 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
(12)dz=(-A·z+B)dt+σz·dWz,whereA>0, B, and σz>0 are constants and Wz is a Wiener process. It can be shown that in the limiting case A→∞ variations uz(t)=z(t)-B/A of z from the fixed point B/A are determined by the Wiener process Wz like [29, 30]
(13)uz(t)dt=(z(t)-BA)dt=σzA·dWz.For sake of clarity, we note that (13) reads in integral form
(14)∫t1t2uz(t)dt=σzA·[Wz(t2)-Wz(t1)].Equation (13) describes an approximation to the exact solutions z(t) of (12). The approximation (13) is called the adiabatic elimination of z because (13) can be used to replace the expression zdt in other formula. According to the adiabatic elimination method, from the stochastic evolution equations (9) and (10) we obtain
(15)[m-rk]dt=σm,akdWm+em(volat., non-Markov),(16)[p-ω·mα+ε]dt=σp,aα+εdWp+ep(volat., non-Markov).The error terms in (15) and (16) account for non-Markovian effects that occur when the limiting case k, k′→∞ 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
(17)dx=bdt+αα+εωkσm,a·dWm+αα+εσp,a·dWpdx=+σx,a·dWx+ex,tot(volat., non  Markov)withb defined by (5). The error term ex,tot(volat., nonMarkov) comprises effects arising from the adiabatic elimination of m and p, 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 Wm, Wp, and Wx can be described by means of a single Wiener process Weff when the volatility is adjusted appropriately. We have
(18)σeff·dWeff=Probαα+εωkσm,a·dWmσeff·dWeff=Prob+αα+εσp,a·dWp+σx,a·dWxwith(19)σeff2=[αα+εωk]2σm,a2+[αα+ε]2σp,a2+σx,a2.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):
(20)σeff2=12αα+εωk[1+2αα+ε(1+ωk)]·r.Having these results at our disposal, (17) can equivalently be expressed as
(21)dx=bdt+σeff·dWeff+ex,tot(volat., non-Markov).The stochastic model defined by (21) provides us with two key parameters. The parameter b is a measure for promoter activity (as argued previously). The term σeff·dWeff 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.


#### 2.1.4. On the Effective Stochastic Single-Step Process

Note that (20) can be written like σeff=[b(1+λ)/2]1/2 with λ=2α·[1+ω/k]/(α+ε). Furthermore, note that for a classical chemical reaction with a constant production rate b the chemical Langevin equation dx=b·dt+(b/2)1/2dW holds [2]. Equation (21) deviates from this form because we have λ>0 and (21) exhibits the error term ex,tot. 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 b0:
(22)dx=b0dt+g·dW.In (22) the parameter g>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 Xk given at discrete time points tk=Δt·(k-1) with k=1,…,n the maximum likelihood estimators (MLEs) b0,e of b0 and ge of g are defined by
(23)b0,e=Xn-X1T,ge2=Δt2(n-1)∑k=1n-1(Xk+1-XkΔt-b0,e)2.Here, the total observation time isT=tn-t1=(n-1)·Δt. 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 ratesk and k′ are sufficiently large, then the error term ex,tot in (21) should have only a small impact on the evolution of x. 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, b0,e and ge are useful estimates of the parameters b 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 b0,e and ge from the simulated protein data Xk using (23). We assumed that mRNA and the protein have half lives of 30 min, whereas observations are made every 2 hours (i.e., Δt=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 k and ε (and consequently k′) 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 r>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.Figure1 shows the sample mean values of b0,e and ge for a sample of 10 repetitions obtained for several values of r≥10. It also shows, for comparison purposes, the analytical results for b and σeff defined by (5) and (20). We see that b and σeff correspond fairly well to the mean values of the estimates b0,e and ge. 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).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 r. (a) Drift coefficients b and b0,e computed from (5) and (23), respectively. (b) Diffusion coefficients [σeff]2 and [ge]2 as obtained from (20) and (23). Sample mean values are presented. Error bars represent SDs. Number of repetitions: 10. Parameters: α=0.1, ω=0.1, k=1.4/h, ε=1.4/h, and Δt=2 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.
(a)
(b)


### 2.2. Population Level

Reporter data are frequently recorded from cell populations rather than from single cells. Let us consider a population ofN 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 k into the extracellular space will be denoted by xk. We assume that the measurement devise records a physical quantity y that is proportional to the total number of stimulus-specific proteins in the extracellular space at time t:
(24)y(t)=c∑k=1Nxk(t).In (24) the parameter c>0 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 t.Let us evaluate the right-hand side of (24). As argued previously, each cell k satisfies an evolution equation of the form (21) such that
(25)dxk=bdt+σeff·dWk,eff+ek,x,tot(volat., non-Markov),whereWk,eff denote statistically independent Wiener processes with variance 2. Substituting (25) into (24), we obtain
(26)dy=Bdt+c·σeff∑k=1NdWk,eff+ey,tot(volat., non-Markov)with(27)B=c·N·b=c·Nαα+εωkr,(28)ey,tot(volat., non-Markov)=c∑k=1Nek,x,tot(volat., non-Markov).The sum of the statistically independent Wiener processes can be replaced by a single Wiener process like(29)∑k=1NWk,eff=ProbN·Wy,whereWy is a Wiener process with variance 2. Exploiting (28) and (29), (26) becomes
(30)dy=Bdt+σy·dWy+ey,tot(volat., non-Markov)withσy=c·N1/2·σeff. Alternatively, with the help of (20), we obtain
(31)σy2=c2·N2αα+εωk[1+2αα+ε(1+ωk)]·r.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 B and σy reflect deterministic and stochastic aspects of promoter activity. The parameter B is a measure proportional to the promoter activity r (i.e., the rate of transcription initiation). If B is estimated from experimental data under various conditions that do not affect the parameters of the post-transcriptional system (c, N, α, ω, ε, and k), then variations observed in B inform us about variations in r that are induced by these different conditions. The parameter σy or alternatively the squared parameter [σy]2 are quantitative measures for the strength of the sources that give rise to gene expression fluctuations. Again, by estimating σy 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 thatB and σy are not independent parameters. By comparing (27) and (31) for B and [σy]2, we recognize that they both depend on the same set of parameters: r, c, N, α, ω, ε, and k. 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
(32)dy=B0dt+G·dWwhose coefficients can be estimated from dataYk using the aforementioned MLEs:
(33)B0,e=Yn-Y1T,Ge2=Δt2(n-1)∑k=1n-1(Yk+1-YkΔt-B0,e)2,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 (N=50) using (6) and (7) for each variable xk(t) and (24) to obtain the measurement variable y(t). We downsampled the trajectory y(t) to obtain a sequence Yk at discrete time points tk. We fitted the computer generated reporter data Yk to the Wiener process with drift (32) by means of the MLEs of B0,e and [Ge]2. Finally, we compared B0,e and Ge with the analytical expressions B and σy given by (27) and (31).The results of this procedure are reported in 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 r. (a) Drift coefficients B and B0,e computed from (27) and (33), respectively. (b) Diffusion coefficients [σy]2 and [Ge]2 as obtained from (31) and (33). Sample mean values are presented. SDs for B0,e were smaller than 1, which is the size of the symbol. SDs for [Ge]2 are shown as error bars. Number of repetitions: 10. Parameters as in Figure 1. Additional parameters: c=1 and N=50. 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.
(a)
(b)Inspection of Figure2 reveals that the estimates B0,e and [Ge]2 of the Wiener process model with drift correspond well with the analytical expressions for B and [σy]2. This indicates that for our choice of parameters the error term ey,tot in (30) affects the evolution of y 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).


## 2.1. Single Cell Level

### 2.1.1. The Stochastic Single Cell Multistep Model

Letm(t) denote the number of mRNA molecules at time t (in what follows molecule numbers are assumed to be large such that they can be treated as real numbers). Let p(t) and x(t) 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:
(1)ddtm=r-k·m,ddtp=ω·m-(α+ε)·p,(2)ddtx=α·p.Accordingly, mRNAm(t) is produced at a rate r and degraded with a degradation constant k. Intracellular proteins p(t) 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 p(t) the term -α·p occurs. In contrast, the evolution equation for the extracellular protein, x(t) 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).Whenm and p reach their stationary values ms and ps, respectively, we have
(3)ms=rk,ps=1α+εωkr,(4)ddtx=b.with(5)b=αα+εωkr.Equation (5) is an important result because it shows that the parameter b is proportional to the transcription rate r. That is, if b can be estimated from experimental data under different conditions and if these conditions do not affect the parameters α, ε, ω, and k of the post-transcriptional system, then variations of b reflect variations of the promoter activity r. 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
(6)dm=(r-k·m)dt+σm·dWm,dp=[ω·m-(α+ε)·p]dt+σp·dWp,dx=α·pdt+σx·dWx.HereWm, Wp, and Wx denote statistically independent Wiener processes that (without loss of generality) have variance 2. The volatilities σm, σp, and σx are given by [2, 3]
(7)σm=r+k·m2,σp=ω·m+(α+ε)·p2,σx=α·p2.


### 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 ofm, p, and x. 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:
(8)σm,a=r,σp,a=ω·rk,σx,a=α·ω·r2(α+ε)k.Such constant volatilities are known to describe so-called additive noise systems, which is the reason why we have added the subindex “a”. Exploiting the additive noise volatilities, the stochastic model defined by (6) can be written like
(9)dm=(r-k·m)dt+σm,a·dWm+em(volat.),(10)dp=[ω·m-(α+ε)·p]dt+σp,a·dWp+ep(volat.),(11)dx=α·pdt+σx,a·dWx+ex(volat.).The error termsem, ep, and ex 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 m(t) and p(t) formally and substitute the result into the evolution equation of x, then x is described by a non-Markovian process (for a similar workedout example see [28]). However, if the decay constants k and k′=α+ε are relatively large, such that m and p decay quickly to their stationary values, then x 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
(12)dz=(-A·z+B)dt+σz·dWz,whereA>0, B, and σz>0 are constants and Wz is a Wiener process. It can be shown that in the limiting case A→∞ variations uz(t)=z(t)-B/A of z from the fixed point B/A are determined by the Wiener process Wz like [29, 30]
(13)uz(t)dt=(z(t)-BA)dt=σzA·dWz.For sake of clarity, we note that (13) reads in integral form
(14)∫t1t2uz(t)dt=σzA·[Wz(t2)-Wz(t1)].Equation (13) describes an approximation to the exact solutions z(t) of (12). The approximation (13) is called the adiabatic elimination of z because (13) can be used to replace the expression zdt in other formula. According to the adiabatic elimination method, from the stochastic evolution equations (9) and (10) we obtain
(15)[m-rk]dt=σm,akdWm+em(volat., non-Markov),(16)[p-ω·mα+ε]dt=σp,aα+εdWp+ep(volat., non-Markov).The error terms in (15) and (16) account for non-Markovian effects that occur when the limiting case k, k′→∞ 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
(17)dx=bdt+αα+εωkσm,a·dWm+αα+εσp,a·dWpdx=+σx,a·dWx+ex,tot(volat., non  Markov)withb defined by (5). The error term ex,tot(volat., nonMarkov) comprises effects arising from the adiabatic elimination of m and p, 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 Wm, Wp, and Wx can be described by means of a single Wiener process Weff when the volatility is adjusted appropriately. We have
(18)σeff·dWeff=Probαα+εωkσm,a·dWmσeff·dWeff=Prob+αα+εσp,a·dWp+σx,a·dWxwith(19)σeff2=[αα+εωk]2σm,a2+[αα+ε]2σp,a2+σx,a2.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):
(20)σeff2=12αα+εωk[1+2αα+ε(1+ωk)]·r.Having these results at our disposal, (17) can equivalently be expressed as
(21)dx=bdt+σeff·dWeff+ex,tot(volat., non-Markov).The stochastic model defined by (21) provides us with two key parameters. The parameter b is a measure for promoter activity (as argued previously). The term σeff·dWeff 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.


### 2.1.4. On the Effective Stochastic Single-Step Process

Note that (20) can be written like σeff=[b(1+λ)/2]1/2 with λ=2α·[1+ω/k]/(α+ε). Furthermore, note that for a classical chemical reaction with a constant production rate b the chemical Langevin equation dx=b·dt+(b/2)1/2dW holds [2]. Equation (21) deviates from this form because we have λ>0 and (21) exhibits the error term ex,tot. 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 b0:
(22)dx=b0dt+g·dW.In (22) the parameter g>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 Xk given at discrete time points tk=Δt·(k-1) with k=1,…,n the maximum likelihood estimators (MLEs) b0,e of b0 and ge of g are defined by
(23)b0,e=Xn-X1T,ge2=Δt2(n-1)∑k=1n-1(Xk+1-XkΔt-b0,e)2.Here, the total observation time isT=tn-t1=(n-1)·Δt. 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 ratesk and k′ are sufficiently large, then the error term ex,tot in (21) should have only a small impact on the evolution of x. 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, b0,e and ge are useful estimates of the parameters b 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 b0,e and ge from the simulated protein data Xk using (23). We assumed that mRNA and the protein have half lives of 30 min, whereas observations are made every 2 hours (i.e., Δt=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 k and ε (and consequently k′) 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 r>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.Figure1 shows the sample mean values of b0,e and ge for a sample of 10 repetitions obtained for several values of r≥10. It also shows, for comparison purposes, the analytical results for b and σeff defined by (5) and (20). We see that b and σeff correspond fairly well to the mean values of the estimates b0,e and ge. 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).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 r. (a) Drift coefficients b and b0,e computed from (5) and (23), respectively. (b) Diffusion coefficients [σeff]2 and [ge]2 as obtained from (20) and (23). Sample mean values are presented. Error bars represent SDs. Number of repetitions: 10. Parameters: α=0.1, ω=0.1, k=1.4/h, ε=1.4/h, and Δt=2 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.
(a)
(b)


## 2.1.1. The Stochastic Single Cell Multistep Model

Letm(t) denote the number of mRNA molecules at time t (in what follows molecule numbers are assumed to be large such that they can be treated as real numbers). Let p(t) and x(t) 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:
(1)ddtm=r-k·m,ddtp=ω·m-(α+ε)·p,(2)ddtx=α·p.Accordingly, mRNAm(t) is produced at a rate r and degraded with a degradation constant k. Intracellular proteins p(t) 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 p(t) the term -α·p occurs. In contrast, the evolution equation for the extracellular protein, x(t) 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).Whenm and p reach their stationary values ms and ps, respectively, we have
(3)ms=rk,ps=1α+εωkr,(4)ddtx=b.with(5)b=αα+εωkr.Equation (5) is an important result because it shows that the parameter b is proportional to the transcription rate r. That is, if b can be estimated from experimental data under different conditions and if these conditions do not affect the parameters α, ε, ω, and k of the post-transcriptional system, then variations of b reflect variations of the promoter activity r. 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
(6)dm=(r-k·m)dt+σm·dWm,dp=[ω·m-(α+ε)·p]dt+σp·dWp,dx=α·pdt+σx·dWx.HereWm, Wp, and Wx denote statistically independent Wiener processes that (without loss of generality) have variance 2. The volatilities σm, σp, and σx are given by [2, 3]
(7)σm=r+k·m2,σp=ω·m+(α+ε)·p2,σx=α·p2.


## 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 ofm, p, and x. 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:
(8)σm,a=r,σp,a=ω·rk,σx,a=α·ω·r2(α+ε)k.Such constant volatilities are known to describe so-called additive noise systems, which is the reason why we have added the subindex “a”. Exploiting the additive noise volatilities, the stochastic model defined by (6) can be written like
(9)dm=(r-k·m)dt+σm,a·dWm+em(volat.),(10)dp=[ω·m-(α+ε)·p]dt+σp,a·dWp+ep(volat.),(11)dx=α·pdt+σx,a·dWx+ex(volat.).The error termsem, ep, and ex 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 m(t) and p(t) formally and substitute the result into the evolution equation of x, then x is described by a non-Markovian process (for a similar workedout example see [28]). However, if the decay constants k and k′=α+ε are relatively large, such that m and p decay quickly to their stationary values, then x 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
(12)dz=(-A·z+B)dt+σz·dWz,whereA>0, B, and σz>0 are constants and Wz is a Wiener process. It can be shown that in the limiting case A→∞ variations uz(t)=z(t)-B/A of z from the fixed point B/A are determined by the Wiener process Wz like [29, 30]
(13)uz(t)dt=(z(t)-BA)dt=σzA·dWz.For sake of clarity, we note that (13) reads in integral form
(14)∫t1t2uz(t)dt=σzA·[Wz(t2)-Wz(t1)].Equation (13) describes an approximation to the exact solutions z(t) of (12). The approximation (13) is called the adiabatic elimination of z because (13) can be used to replace the expression zdt in other formula. According to the adiabatic elimination method, from the stochastic evolution equations (9) and (10) we obtain
(15)[m-rk]dt=σm,akdWm+em(volat., non-Markov),(16)[p-ω·mα+ε]dt=σp,aα+εdWp+ep(volat., non-Markov).The error terms in (15) and (16) account for non-Markovian effects that occur when the limiting case k, k′→∞ 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
(17)dx=bdt+αα+εωkσm,a·dWm+αα+εσp,a·dWpdx=+σx,a·dWx+ex,tot(volat., non  Markov)withb defined by (5). The error term ex,tot(volat., nonMarkov) comprises effects arising from the adiabatic elimination of m and p, 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 Wm, Wp, and Wx can be described by means of a single Wiener process Weff when the volatility is adjusted appropriately. We have
(18)σeff·dWeff=Probαα+εωkσm,a·dWmσeff·dWeff=Prob+αα+εσp,a·dWp+σx,a·dWxwith(19)σeff2=[αα+εωk]2σm,a2+[αα+ε]2σp,a2+σx,a2.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):
(20)σeff2=12αα+εωk[1+2αα+ε(1+ωk)]·r.Having these results at our disposal, (17) can equivalently be expressed as
(21)dx=bdt+σeff·dWeff+ex,tot(volat., non-Markov).The stochastic model defined by (21) provides us with two key parameters. The parameter b is a measure for promoter activity (as argued previously). The term σeff·dWeff 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.


## 2.1.4. On the Effective Stochastic Single-Step Process

Note that (20) can be written like σeff=[b(1+λ)/2]1/2 with λ=2α·[1+ω/k]/(α+ε). Furthermore, note that for a classical chemical reaction with a constant production rate b the chemical Langevin equation dx=b·dt+(b/2)1/2dW holds [2]. Equation (21) deviates from this form because we have λ>0 and (21) exhibits the error term ex,tot. 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 b0:
(22)dx=b0dt+g·dW.In (22) the parameter g>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 Xk given at discrete time points tk=Δt·(k-1) with k=1,…,n the maximum likelihood estimators (MLEs) b0,e of b0 and ge of g are defined by
(23)b0,e=Xn-X1T,ge2=Δt2(n-1)∑k=1n-1(Xk+1-XkΔt-b0,e)2.Here, the total observation time isT=tn-t1=(n-1)·Δt. 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 ratesk and k′ are sufficiently large, then the error term ex,tot in (21) should have only a small impact on the evolution of x. 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, b0,e and ge are useful estimates of the parameters b 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 b0,e and ge from the simulated protein data Xk using (23). We assumed that mRNA and the protein have half lives of 30 min, whereas observations are made every 2 hours (i.e., Δt=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 k and ε (and consequently k′) 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 r>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.Figure1 shows the sample mean values of b0,e and ge for a sample of 10 repetitions obtained for several values of r≥10. It also shows, for comparison purposes, the analytical results for b and σeff defined by (5) and (20). We see that b and σeff correspond fairly well to the mean values of the estimates b0,e and ge. 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).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 r. (a) Drift coefficients b and b0,e computed from (5) and (23), respectively. (b) Diffusion coefficients [σeff]2 and [ge]2 as obtained from (20) and (23). Sample mean values are presented. Error bars represent SDs. Number of repetitions: 10. Parameters: α=0.1, ω=0.1, k=1.4/h, ε=1.4/h, and Δt=2 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.
(a)
(b)


## 2.2. Population Level

Reporter data are frequently recorded from cell populations rather than from single cells. Let us consider a population ofN 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 k into the extracellular space will be denoted by xk. We assume that the measurement devise records a physical quantity y that is proportional to the total number of stimulus-specific proteins in the extracellular space at time t:
(24)y(t)=c∑k=1Nxk(t).In (24) the parameter c>0 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 t.Let us evaluate the right-hand side of (24). As argued previously, each cell k satisfies an evolution equation of the form (21) such that
(25)dxk=bdt+σeff·dWk,eff+ek,x,tot(volat., non-Markov),whereWk,eff denote statistically independent Wiener processes with variance 2. Substituting (25) into (24), we obtain
(26)dy=Bdt+c·σeff∑k=1NdWk,eff+ey,tot(volat., non-Markov)with(27)B=c·N·b=c·Nαα+εωkr,(28)ey,tot(volat., non-Markov)=c∑k=1Nek,x,tot(volat., non-Markov).The sum of the statistically independent Wiener processes can be replaced by a single Wiener process like(29)∑k=1NWk,eff=ProbN·Wy,whereWy is a Wiener process with variance 2. Exploiting (28) and (29), (26) becomes
(30)dy=Bdt+σy·dWy+ey,tot(volat., non-Markov)withσy=c·N1/2·σeff. Alternatively, with the help of (20), we obtain
(31)σy2=c2·N2αα+εωk[1+2αα+ε(1+ωk)]·r.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 B and σy reflect deterministic and stochastic aspects of promoter activity. The parameter B is a measure proportional to the promoter activity r (i.e., the rate of transcription initiation). If B is estimated from experimental data under various conditions that do not affect the parameters of the post-transcriptional system (c, N, α, ω, ε, and k), then variations observed in B inform us about variations in r that are induced by these different conditions. The parameter σy or alternatively the squared parameter [σy]2 are quantitative measures for the strength of the sources that give rise to gene expression fluctuations. Again, by estimating σy 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 thatB and σy are not independent parameters. By comparing (27) and (31) for B and [σy]2, we recognize that they both depend on the same set of parameters: r, c, N, α, ω, ε, and k. 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
(32)dy=B0dt+G·dWwhose coefficients can be estimated from dataYk using the aforementioned MLEs:
(33)B0,e=Yn-Y1T,Ge2=Δt2(n-1)∑k=1n-1(Yk+1-YkΔt-B0,e)2,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 (N=50) using (6) and (7) for each variable xk(t) and (24) to obtain the measurement variable y(t). We downsampled the trajectory y(t) to obtain a sequence Yk at discrete time points tk. We fitted the computer generated reporter data Yk to the Wiener process with drift (32) by means of the MLEs of B0,e and [Ge]2. Finally, we compared B0,e and Ge with the analytical expressions B and σy given by (27) and (31).The results of this procedure are reported in 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 r. (a) Drift coefficients B and B0,e computed from (27) and (33), respectively. (b) Diffusion coefficients [σy]2 and [Ge]2 as obtained from (31) and (33). Sample mean values are presented. SDs for B0,e were smaller than 1, which is the size of the symbol. SDs for [Ge]2 are shown as error bars. Number of repetitions: 10. Parameters as in Figure 1. Additional parameters: c=1 and N=50. 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.
(a)
(b)Inspection of Figure2 reveals that the estimates B0,e and [Ge]2 of the Wiener process model with drift correspond well with the analytical expressions for B and [σy]2. This indicates that for our choice of parameters the error term ey,tot in (30) affects the evolution of y 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% O2 (normoxia) or 1% O2 (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 B and [σy]2 depend linearly on the production rate r. Consequently, variations in r due to experimental manipulations should affect the deterministic and stochastic components of the post-transcriptional system as quantified by B and [σy]2 to the same degree provided that all other parameters (α, ω, ε, and k) 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
(34)σy2=12c·B·(1-λ)withλ depending on α, ω, ε, and k as defined in Section 2.1.4. This implies that the ratio [σy]2/B is independent of r:
(35)σy2B=12c·(1-λ).For each replicate we calculatedB0,e and [Ge]2 from the five data points Y1,…,Y5 related to the time points tk=0, 3, 6, 9, 12 h using (33). For each replicate we also calculated the ratio [Ge]2/B0,e. 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 [G]2/B0 (as estimated by [Ge]2/B0,e) obtained for the four different conditions should correspond to the ratios [σy]2/B of the Wiener process model (32) and consequently should be on the same level (see (35)).


### 3.3. Results

Figure3 reports sample mean values of B0,e, [Ge]2 and the ratio [Ge]2/B0,e. From Figure 3(a) it follows that the promoter activity B (as estimated B0,e) was minimal under baseline condition (normoxia with 0 ng/mL TNFα) with B=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 B=0.87. The single stimulation conditions produced transcriptional activity levels in the medium range of B=0.30 (normoxia with 1 ng/mL TNFα) and B=0.30 again (hypoxia with 0 ng/mL TNFα). In particular, under dual stimulation the activity was higher (B=0.87) than the sum of the activity levels observed under single stimulation (B=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 B (F(1,7)=28.967,P=0.001). Likewise, stimulation with 1 ng/mL TNFα induced a significant increase of B (F(1,7)=96.446, P<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 (F(1,7)=14.740, P<0.01).Transcription rate as measured by the drift coefficientB (a), strength of the fluctuating force inducing gene expression variability as measured by the diffusion coefficient [σy]2 (b), and ratio [σy]2/B (c) observed under four different stimulation conditions. Error bares indicated standard errors of the means.
(a)
B
(b)
[
σ
y
]
2
(c)
[
σ
y
]
2
/
BFigure3(b) reveals that the measure [σy]2 for the strength of the sources of gene expression variability (as estimated by [Ge]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 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 (P<0.01).Figure3(c) reports the sample mean values of the ratios [σy]2/B as estimated by [Ge]2/B0,e. 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 [σy]2/B. That is, although the promoter activity B 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 [σy]2 versus the drift parameters B were found to be approximately on the same level. This is because the diffusion coefficients [σy]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 B 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 [σy]2/B was much higher than the ratios obtained in the other three conditions. Although both B and [σy]2 were maximal in the dual stimulation condition, if we compare, for example, the baseline condition with the dual stimulation condition, then B and [σy]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 [σy]2/B for the dual stimulation was significantly higher than the average ratio [σy]2/B for the three other conditions (P<0.01). 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
(36)y(tk+1)=y(tk)+B0Δt+G·ΔW.The incrementsΔW are normally distributed with variance 2Δt and reflect white noise. Consequently, they are statistically independent (see, e.g., [35]). From (36) it follows that the residuals defined by
(37)ξ(k)=y(tk+1)-y(tk)-B0ΔtGΔt=ΔWΔtare proportional to the incrementsΔW (as indicated) and, in doing so, are statistically independent. From the time points tk=0, 3, 6, 9, 12 h, we calculated four residuals k=1, 2, 3, 4. Model checking typically is conducted by calculating lag-p autocorrelation coefficients and comparing them with the threshold 2/N1/2, where N 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 R samples. The threshold for the sample mean lag-1 autocorrelation is 2/(R·N)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.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.
Threshold2/(R·N)1/2 with R=8, N=4
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.35The 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.


## 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% O2 (normoxia) or 1% O2 (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 B and [σy]2 depend linearly on the production rate r. Consequently, variations in r due to experimental manipulations should affect the deterministic and stochastic components of the post-transcriptional system as quantified by B and [σy]2 to the same degree provided that all other parameters (α, ω, ε, and k) 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
(34)σy2=12c·B·(1-λ)withλ depending on α, ω, ε, and k as defined in Section 2.1.4. This implies that the ratio [σy]2/B is independent of r:
(35)σy2B=12c·(1-λ).For each replicate we calculatedB0,e and [Ge]2 from the five data points Y1,…,Y5 related to the time points tk=0, 3, 6, 9, 12 h using (33). For each replicate we also calculated the ratio [Ge]2/B0,e. 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 [G]2/B0 (as estimated by [Ge]2/B0,e) obtained for the four different conditions should correspond to the ratios [σy]2/B of the Wiener process model (32) and consequently should be on the same level (see (35)).


## 3.3. Results

Figure3 reports sample mean values of B0,e, [Ge]2 and the ratio [Ge]2/B0,e. From Figure 3(a) it follows that the promoter activity B (as estimated B0,e) was minimal under baseline condition (normoxia with 0 ng/mL TNFα) with B=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 B=0.87. The single stimulation conditions produced transcriptional activity levels in the medium range of B=0.30 (normoxia with 1 ng/mL TNFα) and B=0.30 again (hypoxia with 0 ng/mL TNFα). In particular, under dual stimulation the activity was higher (B=0.87) than the sum of the activity levels observed under single stimulation (B=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 B (F(1,7)=28.967,P=0.001). Likewise, stimulation with 1 ng/mL TNFα induced a significant increase of B (F(1,7)=96.446, P<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 (F(1,7)=14.740, P<0.01).Transcription rate as measured by the drift coefficientB (a), strength of the fluctuating force inducing gene expression variability as measured by the diffusion coefficient [σy]2 (b), and ratio [σy]2/B (c) observed under four different stimulation conditions. Error bares indicated standard errors of the means.
(a)
B
(b)
[
σ
y
]
2
(c)
[
σ
y
]
2
/
BFigure3(b) reveals that the measure [σy]2 for the strength of the sources of gene expression variability (as estimated by [Ge]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 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 (P<0.01).Figure3(c) reports the sample mean values of the ratios [σy]2/B as estimated by [Ge]2/B0,e. 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 [σy]2/B. That is, although the promoter activity B 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 [σy]2 versus the drift parameters B were found to be approximately on the same level. This is because the diffusion coefficients [σy]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 B 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 [σy]2/B was much higher than the ratios obtained in the other three conditions. Although both B and [σy]2 were maximal in the dual stimulation condition, if we compare, for example, the baseline condition with the dual stimulation condition, then B and [σy]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 [σy]2/B for the dual stimulation was significantly higher than the average ratio [σy]2/B for the three other conditions (P<0.01). 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
(36)y(tk+1)=y(tk)+B0Δt+G·ΔW.The incrementsΔW are normally distributed with variance 2Δt and reflect white noise. Consequently, they are statistically independent (see, e.g., [35]). From (36) it follows that the residuals defined by
(37)ξ(k)=y(tk+1)-y(tk)-B0ΔtGΔt=ΔWΔtare proportional to the incrementsΔW (as indicated) and, in doing so, are statistically independent. From the time points tk=0, 3, 6, 9, 12 h, we calculated four residuals k=1, 2, 3, 4. Model checking typically is conducted by calculating lag-p autocorrelation coefficients and comparing them with the threshold 2/N1/2, where N 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 R samples. The threshold for the sample mean lag-1 autocorrelation is 2/(R·N)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.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.
Threshold2/(R·N)1/2 with R=8, N=4
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.35The 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 (Sections2.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 Section2 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 r 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 r 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 Section2.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 [G]2/B as the baseline condition. The dual stimulation condition exhibited a ratio [G]2/B 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” [G]2/B. 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 r might be subjected to fluctuations to a greater extent under dual stimulation. Furthermore, from (5) and (20) it follows that the transcription rate B and the diffusion coefficient [G]2 depend on the factor
(38)Θ=αα+εωk.However,B depends linearly on this factor, whereas [G]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 [G]2/B. Consequently, another speculative explanation for the observed increased ratio [G]2/B 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 B, we estimated the model parameter B 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)B 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 Table2, 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 B using mean substitution of the missing data. The interaction was significant (F(2,14)=11.925, P=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 (F(1,7)=11.078, P<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 Section2 can be applied to transcriptional activity regulated by several transcription factors. To this end, the transcriptional activity parameter r 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 r and transcription factor concentrations can be estimated when B (in replacement for r) is estimated for dose responses of the set of transcription factors of interest.



---

*Source: 719138-2013-09-09.xml*