Vijayanarasimha H Pakka

Adam Prügel-Bennett

Srinandan Dasmahapatra

Correlated fluctuations carry signatures of gene regulatory network dynamics

Keywords: transcriptional regulation, gene regulatory networks, correlated fluctuations, fluctuation-dissipation relations, master equation

The dynamics of transcriptional control involve small numbers of molecules and result in significant fluctuations in protein and mRNA concentrations. The correlations between these intrinsic fluctuations then offer, via the fluctuation dissipation relation, the possibility of capturing the systems response to external perturbations, and hence the nature of the regulatory activity itself. We show that for simple regulatory networks of activators and repressors, the correlated fluctuations between molecular species show distinct characteristics for changes in regulatory mechanism and for changes to the topology of causal influence. Here, we do a stochastic analysis and derive time-dependent correlation functions between molecular species of regulatory networks and present analytical and numerical results on peaks and delays in correlations between proteins within networks. Upon using these values of peaks and delays as a 2-dimensional feature space, we find that different regulatory mechanisms separate into distinct clusters. This indicates that experimentally observable pairwise correlations can distinguish between gene regulatory networks.

Introduction

A gene regulatory network (GRN) summarizes the causal influences between genes that regulate the expression of proteins in cells. While the set of genes constitute the nodes of this network, the regulatory activity between genes is represented by its edges. The biochemical steps involved in this activity is inherently stochastic in nature, and fluctuating numbers of mRNA and protein levels have been observed in clonal populations. Many interesting biological phenomena have been attributed to these fluctuations (Arkin et al., 1998;[START_REF] Raser | Noise in gene expression: Origins, consequences, and control[END_REF]Raj and van Oudenaarden, 2008), which include phenotypic variation within an isogenic population of cells or their influence on switching between cell fates. Such examples include the lysis/lysogeny switch in bacteriophage-λ (Arkin et al., 1998) driven by fluctuations in Cro 2 and CI 2 ; competent-vegetative switching in B. subtilis cells [START_REF] Maamar | Noise in gene expression determines cell fate in bacillus subtilis[END_REF]; and in the regulatory circuit responsible for colour vision in Drosophila, where cells commit to either yellow(70%) or pale(30%) ommatidial subtypes [START_REF] Wernet | Stochastic spineless expression creates the retinal mosaic for colour vision[END_REF]. Stochastic effects in transcription factor proteins affect the regulation of the genes whose cis-regulatory upstream sequence it binds to; this is exemplified in the variability in the time of onset of early meiosis genes in diploid yeast cells subjected to nutritional deprivation linked to fluctuations in the level of Ime1 proteins [START_REF] Nachman | Dissecting timing variability in yeast meiosis[END_REF]. Cell-cell variability due to transcriptional noise has been shown to be advantageous in adverse environmental conditions in [START_REF] Blake | Phenotypic consequences of promoter-mediated transcriptional noise[END_REF]. While such examples establish the salience of stochasticity in cell biology, how fluctuations are organized and propagate within regulatory networks remains an important issue, and whether different networks leave identifiable traces on the measured fluctuations in protein numbers.

Simple models of gene regulation [START_REF] Kepler | Stochasticity in transcriptional regulation: origins, consequences, and mathematical representations[END_REF] identify experimental interactions that outline the contributions of the core processes in variability of gene expression [START_REF] Ozbudak | Regulation of noise in the expression of a single gene[END_REF]. Quantifying molecular fluctuations by such formulations allows one to study the effect of various network topologies, regulatory mechanisms and parameters on the stationary statistics of variables in GRNs [START_REF] Tao | Effect of feedback regulation on stochastic gene expression[END_REF][START_REF] Tomioka | Multivariate analysis of noise in genetic regulatory networks[END_REF][START_REF] Pedraza | Noise propagation in gene networks[END_REF]. Such observations make it possible for the quantification of noise to reveal the regulatory links between genes [START_REF] Austin | Gene network shaping of inherent noise spectra[END_REF][START_REF] Cox | Using noise to probe and characterize gene circuits[END_REF]. For example, [START_REF] Cox | Using noise to probe and characterize gene circuits[END_REF] consider single-gene models seeking to characterize regulatory networks by their specific features of their stationary auto-correlations, viz. the coefficient of variation and the time at which auto-correlations decay to half their maximum value. However, a crucial element missing in such studies is the temporal aspect. Since the entire process of gene expression is dynamic in nature, tracking the fluctuations in species numbers over time would yield significant information. For example, due to the inherent time-delay in the transcription and translation steps, a perturbation in the upstream process of gene activation say, would result in a delayed response in the protein fluctuations. Such a response could be evaluated by temporal correlations that are more likely to be characteristic of the regulatory mechanisms, structure and parameter values of the system. Such studies are possible by tracking of mRNA/protein numbers over time in individual cells, e.g., (Raj et al., 2008). With the aid of time-lapse fluorescence measurements, recent reports suggest that the time dependent correlations of fluctuations in protein levels indicate the presence of regulatory activity in simple GRNs [START_REF] Dunlop | Regulatory activity revealed by dynamic correlations in gene expression noise[END_REF][START_REF] Sigal | Variability and memory of protein levels in human cells[END_REF]. That is the object of our study.

A number of papers have studied the single-gene case in great detail [START_REF] Paulsson | Summing up the noise in gene networks[END_REF][START_REF] Lei | Stochasticity in single gene expression with both intrinsic noise and fluctuation in kinetic parameters[END_REF][START_REF] Shahrezaei | Analytical distributions for stochastic gene expression[END_REF] and have investigated the noise properties of multi-gene networks [START_REF] Tomioka | Multivariate analysis of noise in genetic regulatory networks[END_REF][START_REF] Pedraza | Noise propagation in gene networks[END_REF]. An important difference between these investigations and our current work lies in our reliance on time-delayed correlation functions, not equal time ones and Appendix C presents the simplest case of a single gene in some detail. The motivation behind looking at time delayed correlations lies in linear response theory [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry, 3rd Edition[END_REF].

Given the stochasticity of the very processes that constitute regulatory responses to stimuli, it is likely that correlations between fluctuations will illuminate the dynamics of regulatory interactions. Indeed the intuition behind Onsager's regression hypothesis [START_REF] Onsager | Reciprocal relations in irreversible processes[END_REF] is that the regression to an equilibrium state after a short-lived external perturbation or an intrinsic fluctuation are identical, which has been developed further in several fluctuation-dissipation theorems even away from equilibrium [START_REF] Keizer | Statistical Thermodynamics of Nonequilibrium Processes[END_REF][START_REF] Speck | Restoring a fluctuation-dissipation theorem in a nonequilibrium steady state[END_REF]. Regulation of expression determines responses to signals which are organized by the cell to form functional pathways. In this paper, we consider correlated fluctuations in elementary fragments of GRNs to illustrate what they can tell us about the nature of the regulatory function enacted by the network.

To achieve that, we study the time-dependent correlation functions between proteins (such as a transcription factor and a downstream target) in a network expressed in a sum-of-exponentials form. We expect distinct behavioural patterns in the correlated fluctuations for variations in the regulatory networks. Variations are introduced in elementary fragments of regulatory networks by employing three strategies: (i) adopting different regulatory mechanisms; activation/repression with or without co-operative mechanisms for the same topology of dependence between genes, (ii) introduction of additional genes into the network, and (iii) adopting different sets of values for the reaction rate constants. They amount to varying the regulatory mechanisms, the network structure and the parameter values respectively, and hence form the basis for deriving qualitative relationships between these attributes and features extracted from the fluctuation properties of the network, which we summarize in Figure 3 below. Therefore, by tracking the species numbers in single cells over time, one could predict the type of regulatory activity present in a GRN, as demonstrated here for the case of small networks.

After setting up the machinery of stochastic processes that we need to compute correlation functions in Section 2, we present specific details of the models and the values that we choose for the rate constants in Section 3. Section 4 establishes the dependences of correlations on parameter values in different networks to set the stage for the comparison presented in Section 5.

Dynamic Correlation Functions

The relation between the macroscopic dynamics and the fluctuation properties of the system is formalized by the Fluctuation-Dissipation Theorem (FDT). In statistical physics, the FDT states that the linear response of an equilibrium (or steady state) [START_REF] Keizer | Statistical Thermodynamics of Nonequilibrium Processes[END_REF] system to external perturbations is obtained from the correlation functions of the fluctuations in the unperturbed system. In other words, it proposes that there is an explicit relationship between the internal fluctuating force that is random in nature, and the observed macroscopic response of the system to an external force that governs the dynamics of the averages. In the present context of biochemical reacting systems, this relation emerges naturally out of a system-size expansion of the Chemical Master equation (CME) [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry, 3rd Edition[END_REF], where the CME is simplified to the Linear Noise Approximation (LNA). This LNA is equivalent to the Fokker-Planck (F-P) for these processes [START_REF] Gillespie | The chemical langevin equation[END_REF], whose solution are the distributions of the fluctuations around the average concentrations of the molecules, or the Langevin equation, which generates system trajectories. Here, we derive the dynamic correlations between these molecular fluctuations, via the Langevin reduction of the CME, and on the way observe the emergence of a relationship between the dynamics of the averages and the fluctuating part of the distributions.

Supposing that X(t) is a vector of the number of molecular species in the system and that the evolution of this vector is a jump-type Markov process, then the time-evolution of the joint probability distribution of all species of a spatially homogeneous chemical system is given by the CME

∂ ∂t P(X, t) = M j=1 P(X -ν j , t)a j (X -ν j ) -P(X, t)a j (X) , (1)
where the propensity function a j gives probability for one out of a possible M reactions r j (1 ≤ j ≤ M) to occur in time [t, t + dt) in volume Ω, thus changing the number of s i molecules by an amount ν i j called the stoichiometry. The initial condition for the system to be P(X, t 0 ) = δ(X -X 0 ). Upon multiplying eq. ( 1) by X and summing over all configurations X after making the substitution X = X -ν j + ν j , we get the following expressions for the evolution of the mean values X :

d X i dt = M j=1 ν i j a j (X) . (2)
Similarly, by multiplying eq. ( 1) by (X -X )(X T -X T ) and shifting variables in the sum over configurations X, we get the following expression for the covariances:

d Cov(X i , X j ) dt = M k=1 (X i -X i )ν jk a k (X) + (X j -X j )ν ik a k (X) + ν ik ν jk a k (X) (3) where Cov(X i , X j ) = (X i -X i )(X j -X j ) .
These equations do not close on themselves if the dependence of the propensities on particle numbers is non-linear. The rate equations for the averages of the concentrations x i = X i /Ω obey mass action kinetics of the biochemical system if we ignore fluctuations

d x i dt = M j=1 ν i j R j ( x ) ( 4 )
where R j (x) := 1 Ω a j ( X Ω ) are the deterministic rates of the M reactions, with mean molecular numbers substituted by deterministic concentrations. It must be noted that while the same functional forms are retained for the propensities a j and reaction rates R j , the probabilities in the stochastic description are replaced by rate constants appropriately scaled by powers of Ω to translate between molecule numbers and concentrations [START_REF] Gillespie | The chemical langevin equation[END_REF][START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry, 3rd Edition[END_REF].

Linearizing about the steady states x s which satisfy d dt x = 0 in (4), we find that the deviations δx(t) := x(t)x s are given by

d dt δx = Aδx, ( 5
)
where A is the Jacobian matrix with elements:

A ik = ∂ ∂ x k ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ M j=1 ν i j R j ( x ) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ x = x s . ( 6
)
At steady state, the means in the kinetic rates R j are time-independent, and so are the elements of A and

BB T = ν diag (R( x s )) ν T ,
yielding the closed form equation for the covariances C(t) := δx(t) δx(t) T (from eq. 3):

∂ ∂t C = AC + CA T + BB T .
The deviations from steady state are thus governed by a multi-variate Ornstein-Uhlenbeck process [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry, 3rd Edition[END_REF]

d δx(t) = A δx(t) dt + B dW(t),
where dW(t) is the Wiener measure. Introducing the Green's function Y(t) which satisfies the homogeneous equation (from eq. ( 5)) d dt Y(t) = AY(t) with initial condition Y(0) = 1, we obtain δx(t + τ) as a sum of a homogeneous solution and an inhomogeneous term

δx(t + τ) = Y(τ) δx(t) + t+τ t Y(t -t ) B dW(t ) dt , τ ≥ 0,
from which we obtain the two-time correlations C(t, t + τ):

C(t, t + τ) := δx(t + τ) δx(t) T = Y(τ)C(t, t) =: Y(τ)C(t). (7)
The variable t takes values in the regime where steady state has been reached and hence the stationary covariance C(t) is t-independent. In the rest of the document, all correlation functions computed and referred to are in the steady state regime, and two-point correlations C(t, t + τ) are independent of t, and depend only on the time interval τ. Using the independence and bi-orthogonality properties of the eigenvectors of the non-symmetric matrix A, it is possible [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry, 3rd Edition[END_REF][START_REF] Elf | Fast evaluation of fluctuations in biochemical networks with the linear noise approximation[END_REF] to write the time-covariance function between two molecular species x i and x j as:

δx i (t + τ) δx j (t) = k Y ik (τ)C k j = k l e λ l τ U T lk V il C k j . (8)
Here, the elements of V and U T , the matrices comprising of the right and left eigenvectors of A respectively and the eigenvalues λ i are explicitly written out.

For effective comparison of results from different GRNs, we use time-delayed correlations by normalizing the above covariances. Use of stationary auto-covariances for normalization helps in retaining the dynamic character along the τ-axis, though the magnitudes are rescaled between 0 and 1.

Corr [δx i (t + τ), δx j (t)] = δx i (t + τ) δx j (t) (δx i (t + τ)) 2 (δx j (t)) 2 (9)
The form of the covariance function brings out the existing relation between the covariances and the averaged responses of the system. The system Jacobian A and diffusion matrix BB T are derived through the stoichiometry ν and reaction rates R that are in turn obtained from the deterministic rate equations. Therefore these rate equations that are responsible for the time-evolution of averages also influence the internal fluctuations of the system. These have been separated into contributions characteristic of individual genes in the network (intrinsic noise) and those from various upstream or system-wide elements such as polymerase fluctuations (extrinsic noise) in the literature [START_REF] Raser | Noise in gene expression: Origins, consequences, and control[END_REF].

From the dynamic correlations we shall extract features for possible identification of the type of regulation present between pairs of genes for likely signatures of network structure from the time-series data from a GRN. To demonstrate the above, we consider regulatory mechanisms such as activation and repression via monomers/dimers and network mechanisms such as cascades, dual activators and feedforward loops. A simple and effective way to model these networks is by choosing a set of elementary reactions that capture the significant steps of the regulatory process. If X and Y denote the regulator and regulated elements, the constituent molecular species of such a model would then be the regulated gene G y , mRNAs M x , M y and the proteins P x , P y . By employing fewer molecular species we intend to reduce the complexity of the model, while retaining the salient features of the complete model.

Features of time-dependent protein correlations

Positive regulation between any two genes of a network not only induces positive covariances between their respective proteins but also gives a characteristic shape to them. This can be seen in the covariance function in the case of an elementary activator X → Y as shown in Figure 1(b), and detailed derivations of which are given in Appendix A. The normalized Corr[P x (t), P y (t + τ)] function has exactly the same shape as the covariance. The defining features of these functions are (i) the stationary covariance Cov(τ = 0), (ii) the maximum or peak covariance Cov * = Cov(τ = τ * ) or, the maximum value of the time-dependent correlation Corr * and (iii) τ * , which is the time taken for the covariance to reach this peak, as indicated in Figure 1(b). (Note that in the following, we assume a single maximum for the covariance as a function of time; it is true for the examples presented in this paper.) It shall be seen in the coming sections that these features are sensitive to the various attributes of the system such as the regulatory mechanism, network structure and the rate constants, and therefore act as tools for the identification of these attributes. For example, in Appendix C we show how the feature τ * is sensitive in the rate constants by deriving its expression in the case of a single gene.

In the case of X → Y, Corr[P x (t), P y (t + τ * )] = 0.8 and τ * = 49 minutes. Physically, this means that at equilibrium, the deviations δP y from the average protein concentration P y are influenced mostly by δP x occuring 49 minutes earlier to it. This time delay is the result of the causal link between the two genes. On the other hand, the stationary correlation (at τ = 0), which is of non-zero value (0.6) would not indicate of any such causality. This is supported by the fact that the Corr[P y (t), P x (t + τ)] function has the same non-zero value at time τ = 0, while for τ > 0, is a monotonically decreasing function. Therefore, stationary correlation alone cannot in effect predict the direction of the causal link between genes. Ideally it would not make much sense for the stationary correlation to have a nonzero value, since that would mean that δP x would have an instantaneous effect on δP y . However the reason behind the non-zero value of the stationary correlations is the slow decaying auto-correlation function of P x , due mainly to proteolysis. Therefore, if δP y at say time (t) is influenced by δP x at say time (tt 1 ), the stationary correlation between the two proteins at time (t) could still be non-zero if the auto-correlation function of δP x is non-zero for a time-period (> t 1 ).

Mathematically, the shape of the dynamic covariance curve is due to its functional form N l=1 h l e λ l τ , where the h's are some functions in the eigenvectors of the Jacobian. In the case of a single activator (Appendix A), the eigenvalues are the four decay rates, hence the covariance function is the sum of exponentials in time raised to the negative power of these decay rates. However, since the co-efficients of these exponentials are different functions of the eigenvectors, their values are different, leading to the characteristic shape of their resulting sum. These exponentials normalized between ±1 and 0 are shown in Figure 1(a).

Figure 1: (a) The normalized exponentials in time raised to the negative power of the decay rates, in the case of the elementary activator system. The normalization is performed by dividing the exponentials with their respective co-efficients which are 4.18 × 10 8 for k - P y , 1.0 × 10 7 for k - M x , -6.34 × 10 5 for k - M y , and -4.23 × 10 8 for k - P x . (b) The sum of the four exponentials along with their co-efficients results in the dynamic covariances.

Comparison with Results from Stochastic Simulations

The objective here is to compare the dynamic correlations obtained through the analytical formulation with those obtained by Monte Carlo simulations [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF]. We do this in the case of the elementary activator system X → Y, where correlations are drawn between the proteins P x and P y at steady-state conditions. At low copy numbers of molecules the question of the applicability of analytical results come to fore, which we address here by assuming the following mean concentration levels of: M x = 1.0 nM, M y = 2.0 nM, P x = 20.0 nM and P y = 100.0 nM. The values for the rate constants are chosen accordingly:

k + M x = 0.0407 nM min -1 , k + P x = 0.367 min -1 , k + M y = 0.308 min -1 , k + P y = 0.88 min -1 , k on = 1.0 nM -1 min -1 , k o f f = 20.0 min -1 , k - M x = 0.0407 min -1 , k - M y = 0.0770 min -1 , k - P x = 0.0184 min -1 , k - P y = 0.0176 min -1 .
The ensemble size of the system is equivalent to the number of runs of the simulation program. In Figure 2 we notice that as the ensemble size increases, the correlations obtained from simulations match closely with that of the analytical one. The evaluation of correlations from simulated runs, rests on the equivalence between the ensembleaveraging and time-averaging of the stochastic variables.

Clustering of regulatory networks by correlation features

We investigate the effect that various rate constants have on the features Corr * and τ * of the dynamic correlation functions between the regulated P y and regulator P x proteins. While stationary correlations, especially in the case of elementary networks, have been the object of study in recent years, we probe the variation in the temporal aspect of the correlations and show their effectiveness as network signatures. In section 5 we show the effect of rate constants on stationary correlation in the case of X → Y system. In section 4 we demonstrate the effect of rate constants on dynamic correlations in the case of different regulatory mechansims. The results of this section and section 4 provide sufficient evidence to show that the features of dynamic correlation functions could be effectively used as signatures in identifying regulatory mechanisms of regulatory networks.

Features of the dynamic correlations such as the maximum value of the 2point correlation, Corr * , and the time τ * taken to reach it are plotted in the case of different networks. Clustering of these features for variations in rate constants enables one to identify the regulatory mechanism and structure of the network. As a summary of the results obtained, the effectiveness of dynamic correlations in characterizing small gene regulatory networks is seen in the Figure 3. The following section goes into each cluster of points for a network in more detail. Here we point out that, for example, the top two panels indicate that the movement along the time-delay feature τ * is influenced not by the decay rate of the transcription factor but that of its target. In the case of co-operative activation X → Y ← Z variations were observed in the peak correlation Corr * whilst leaving unaltered the time delay τ * , for changes in the rate constants. The clusters corresponding to the cascading activator networks (cascade of activators and cascade of repressors) have a larger value of τ * , and differ from each other in the magnitude of peak correlations Corr * . By combining the above networks one could compute correlations that capture signatures of more complex regulatory mechanisms such as the coherent and incoherent feedforward networks that have been much studied in the literature (Alon, 2007). Once again these complex networks are clearly identifiable with the aid of the correlation features. Finally, the effect of dimerization is clearly noticed with increased Corr * values. Due to changes in the shape of these correlation plots, the use of integral over time τ could also act as an additional feature for clustering purposes. However, this is a topic for further investigation.

We have demonstrated that at steady-state conditions dynamic correlations of protein fluctuations could be used to gain insight into the microscopic mechanisms of gene regulation. This means that a lot could be learnt by tracking the evolution of mRNA or protein numbers in single cells over time. The most revealing information that can be expected out of such time-series data is the type of regulatory activity present between pairs of genes. As an example, given the information about the lifetimes of proteins and mRNAs of X and Y, the two-time correlations could infer the presence or absence of an intermediary gene. We also demonstrated the effect that various reaction rate constants have on the dynamic correlations between the proteins P x and P y , in the case of different regulatory systems. By adopting a learning methodology, the framework of sensitivity analysis could be used for estimating the parameters of a regulatory system given the dynamic correlations between its species. Experimentally this would be equivalent to altering the reaction rates, as done in [START_REF] Ozbudak | Regulation of noise in the expression of a single gene[END_REF] where transcription and translation rates were varied synthetically, and tracking the time-evolution of species in the altered system to yield data for such novel learning approaches.

Effect of parameters on dynamic correlations in the case of simple regulatory networks

Fragments of a GRN that represent different regulatory mechanisms such as repression, regulation via dimers, cascading activation, co-operative activation, and so on, not only display different dynamic characteristics in the determinitic domain but also show unique characteristics in their molecular fluctuations. These Figure 3: The inset in (a) depicts the networks used in this paper: (X → Y) the single activator; (2X → Y) a single activator network, but with X coding for a protein which performs the activation after forming a homodimer; (X → Z → Y) which describes a cascade of interactions, [40, 50, basevalue, 90, 120] minutes, where the basevalue corresponds to a half-life of 65 minutes. (b) Similarly, (t 1/2 ) P y is varied in steps of [40, 50, basevalue, 90, 120] For time-independent (equal time at steady-state) fluctuations, it has been shown that some rate constants exert influence on both means and variances of protein levels [START_REF] Thattai | Intrinsic noise in gene regulatory networks[END_REF]. Our intention here is to study the effect on fluctuations while the mean concentration levels are held constant at the values given in Table A.2. For example, supposing that the pair [k + M x , k - M x ] are at their base values, deviations in them does not affect

X upregulating Z which upregulates Y; (X-Z-Y) cascading repressors; (X → Y ← Z) the combinatorial action of genes X and Z on Y; (coherent FFL) which combines X → Y ← Z and X → Z; (incoherent FFL) which combines X → Y -Z and X → Z (a) (t 1/2 ) P x is varied in steps of
M x = k + M x /k - M x .
However there is a profound effect on the correlations for equal variations in [k + M x , k - M x ] as seen in Figure 5(a). This is suggestive of possible ways in which biological networks could alter their fluctuation properties without affecting the mean concentration levels of proteins and mRNAs.

From Figures 5 and6 we observe that the feature τ * is sensitive not only to the decay rates but also to the transcription and translation rates. However, the following analysis makes it clear that τ * is predominatly sensitive only to the four decay rates. The system has five eigenvalues, of which three are the decay rates [-k - M x , -k - M y , -k - P y ] whilst the remaining two are roots to the quadratic equation

x 2 + bx + c = 0, where b = k o f f + P x k on + G y k on + k - P x c = k - P x (k o f f + k on P x ) Figure 5: (a)
The dynamical correlation at steady-state between P x and P y is sensitive to changes in the production and decay rates of M x . The base value for the half-life of M x is 21 minutes. k - Mx and k + Mx are varied simultaneously so that their ratio

M x = k + Mx /k
= k o f f k o f f +k on P x . Therefore b = G y + K D G y k on +k - P x and c = k o f f k - P x /G y . Since K D := k o f f /k on is usually of the order of hundreds of nM whereas G y is less than 1 nM, b ≈ k o f f G y +k - P x .
The roots of the quadratic equation are now -k - P x and the much larger -k o f f /G y . The contribution of the large eigenvalues to the sum of exponentials form for the covariance function decays rapidly, which was why the exponential corresponding to -k o f f /G y , in the Figure 1(a) was not shown.

Figures 6(a) and 6(b) show that only those decay rates that correspond to species Y induce greater sensitivity in τ * as shown in . As the decay rates (halflife) of P y and M y decrease (increase), so does τ * increases. Rapid turnover of proteins P y makes it more sensitive to respond to fluctuations in P x which shortens the response times τ * . This explains why at decreased decay rates of k - M y , k - P y , the peak in the correlation function is delayed (τ * shifts to the right in the figures).

With regard to changes in K D , the protein correlations remain unchanged, except for K D < 200 nM when there is progressive reduction only in Corr * while τ * still remains at its base value of 49 minutes. The insensitivity of τ * to the binding/unbinding reactions is because they are much faster that the other elementary reactions. Further, the correlations are totally insensitive to cases where k on and k o f f take on different values such that their ratio K D = k o f f /k on remains invariant. In the elementary fragment X Y the sensitivity of the dynamic correlations between the repressor P x and the repressed protein P y to different values of protein half-life is shown in Figure 8. Note the close similarity between these correlations and those in Figure 6(b). Whilst a change in the regulatory mechanism has influenced the dynamic correlations such that they are now negative in value, the sensitivities of the correlations w.r.t the parameters, are very similar to those of the activator case. This similarity is due to the fact that the reaction structure of the systems are essentially the same, with the only change being the way M y is transcribed.

Elementary Repressor

¤ ¥ ¦ ¥ ¤ § ¦ § ¤ § ¦ §

Activation via Dimerization

In the case of activation via dimers 2X → Y, for short half-lives of P x (about 8 minutes), the effect of dimerization begins to wear off (Figure 10 P y , k - P y ] are varied simultaneously so that the mean concentration level of P y = M y k + P y /k - P y remains constant. While there is an increase in Corr * , τ * decreases from 102 to 36 minutes for step-wise increase in these parameters. molecules are now being produced and degraded faster (due to increase in values of [k + P x , k - P x ]) and thus not able to contribute to dimerization. On the other hand, in Figure 10(b) we do not see any such effect upon reduction in the half-life of P y due to its non-participation in the dimerization process. P x is varied in steps from its base value of 0.0184 min -1 that corresponds to a half-life of (t 1/2 ) P x = 65 minutes. k + P x is simultaneously varied such that P x remains unchanged. For increase in the values of [k + P x , k - P x ], there is an observable increase in Corr * , whilst τ * decreases from 53 to 46 minutes. (b) Likewise, for step-wise increase in the parameters [k + P y , k - P y ], τ * once again reduces from 81 to 20 minutes, whilst Corr * is nearly constant. The effect of different values of K Dim on the protein correlations is shown in Figure 11. The overall effect of dimerization is that there is a higher level of correlation between the proteins and more importantly, the shape of the dynamic correlations is different than that in the monomer case. In fact, for a low value of K Dim = 40 nM, corresponding to P x 2 = 3913 nM, the correlation function loses its characteristic localized peak and its magnitude remains nearly constant at around 0.98 for a long duration of time. Under such conditions, any perturbation in P x has an corresponding effect on P y over a long period of time, i.e., the effect of correlation is sustained over time. Further, (a) only the dissociation constants K D and K Dim affect the correlations and not the corresponding individual forward and reverse rates, (b) though Cov[P x 2 (t), P y (t + τ)] and Cov[P x (t), P y (t + τ)] exhibit distinct behaviours, their normalized dynamic correlations are exactly the same, due to normalization by their respective stationary auto-covariances. Therefore either P x or P x 2 could be considered as the output variable.

Cascading Activation

In the case of the three-gene cascaded activation X → Z → Y, the mRNA and protein decay rates of the intermediate gene Z have an effect on the τ * of the correlation plots (Figure 12(b)). This is due to the fact that while four of the eight eigenvalues of the system are

[-k - M x , -k - M z , -k - M y , -k - P y ]
, the other four are the roots to two quadratic equations and are approximated as [-k

- P x , -k - P z , - k r 1 G z , - k r 2
G y ]. Therefore the element Z further sensitizes the temporal characteristics of the correlations between the input and output elements of the GRN. In this regard it is interesting to note that, for the hypothetical case of mRNA and protein half-lives of Z of around 1 minute each, τ * new ≈ 49 as in the case of X → Y, suggesting that for rapidly decaying mRNA/proteins, the intermediate gene-node ceases to influence the dynamic character of the correlations. Another form of cascaded activation X Z Y throws up similar behaviour except that now the magnitude of correlations are higher.

Effect of parameters on stationary correlations in the case of X → Y

The stationary covariance between molecular species is a function of system parameters that are the reaction rate constants. We show this in Appendix C in the case of a simple model of transcription. [START_REF] Ozbudak | Regulation of noise in the expression of a single gene[END_REF] experimentally demonstrated the validity of the expression of protein variance of one such model by altering the transcription and translation rates. Variations in parameter values therefore provide us with some knowledge with regard to the system at hand. By varying each of the ten rate constants individually, over a vast range, we note the change in value of the stationary correlation Corr[P x (t), P y (t)].

Figure 12: Cascaded Activation X → Z → Y. (a) For higher values of protein half-life (t 1/2 ) P x , τ * increases as [63,87,97,109,118] minutes along with an increase in Corr * . (b) Similarly, for changes in k - P z or rather the half-life of P z , τ * increases as [71,84,97,118,133] minutes along with a decrease in both Corr * and the stationary correlation (at τ = 0).

From Figure 13(a), we see that the correlations are smaller for lower values of the decay rates. As the proteins/mRNAs decay slowly, the number of molecules present at any given time remain mostly constant and hence any fluctuations in P x would not be prominently correlated with the fluctuations in P y . On the other hand, for larger decay rates, due to the fewer number of molecules present, any fluctuations in them are highly correlated. The increase in the stationary correlation is greater in the range where the decay rates are smaller in value. For continuously increasing k - M x and k - P x , the stationary correlation decreases at a later stage. This could possibly point towards the network structure which in this case is X → Y.

Similar reasoning suffices for the variation in the correlation for changes in the other parameters. Firstly, for lower values of k on and for higher values of k o f f , stationary correlations are high. For higher values of K D = k o f f /k on , the transcription factor binds less to the DNA-complex thereby reducing the rate of production of M y and inturn of P y . However, the stationary correlation depicts a different picture, mainly due to the fact that the fluctuations in P x are now less rapid and therefore correlate more with other species of the system. On the other hand, for higher values of the transcription rate k + M y , as expected, protein correlations are high. On introducing the basal transcription rate k basal as an additional parameter to the system, we notice that it has influence only over the magnitude of the correlations, whether stationary or dynamic, and not on the time-delay of the peak τ * itself. Finally, in these feedforward networks translation rates have no effect on the stationary correlations since they do not directly contribute to the transcription process and the proportionate increases in magnitudes of the covariances are normalized away. Analyzing the effects of these rate constants on the fluctuation properties of the molecular species is important for two reasons. One is that the molecular fluctuations as observed in the correlations provide significant information regarding the rate constants of the reaction processes. The other complementary reason is that, such analyses act as analytical tools that aid in designing new synthetic regulatory networks whose fluctuation properties and performance in general, could be controlled.

Conclusions

In this paper, we set out to describe the temporal dependence of the propagation of fluctuations in small gene regulatory networks. In setting up the models, we have introduced the processes of promoter interactions, followed by transcription and translation and investigated the contribution of these elements to the fluctuation properties. There are further time dependent effects that we have disregarded, such as the time for the first complete transcription, following potential pauses of RNA polymerases, entry of the transcripts into the cytoplasm to be translated by the ribosomes and re-entry of transcription factors into the nucleus, and potential effects due to post-transcriptional (such as due to miRNA) or post-translational (such as phosphorylation of the proteins for activity priming). Some of these processes occur at smaller time scales, and just as we could eliminate effects of binding and unbinding of proteins to form dimers on the overall characteristics of the correlation profile, we expect some of them to have a small effect on the location of the peak for instance. It is clear, however, that for many regulatory networks, such as those underlying the circadian clock system [START_REF] Leloup | Modeling the mammalian circadian clock: sensitivity analysis and multiplicity of oscillatory mechanisms[END_REF], phosphorylation of proteins are critically linked to the time of response. It should be possible to incorporate such effects either via the introduction of further components or by explicitly incorporating time delays in the Jacobian. That is beyond the scope of this paper.

What we have shown in this paper is that for small motifs of transcriptional regulatory networks, the correlations between time-delayed fluctuations between pairs of proteins or mRNA at steady state can be projected onto a two-dimensional feature space. This space is spanned by the maximum value of the correlation between fluctuations in numbers of a pair of species Corr * , and the time interval τ * where this peak correlation is attained. Variations in parameters in one motif that preserve mean values of its constituent species in steady state dynamics traces out curves in this (τ * -Corr * ) feature space. We have plotted segments of these curves from different regulatory motifs and have observed that they exhibit a degree of clustering, thus indicating that there may be motif-identifying signatures in transcriptional noise. The extent to which the segments presented in the clustered feature space can be extended to these hypothesized curves, and how they are more globally separated, is linked to a sensitivity analysis of the features extracted which is part of an ongoing study.

A key assumption which motivates the explorations in this paper is the possibility of being able to tag and monitor in a continuous fashion, the quantities of proteins or mRNA in a single cell over time. Fluorescent reporters have increased our ability to observe individual cells, which is a major shift from the microarray realm where bulk averages were the measured quantities. Flow cytometry measurements of the relative fluorescence intensities of cell populations enable one to compute histograms of protein fluorescence levels and infer phenotypic variability in isogenic populations [START_REF] Ozbudak | Regulation of noise in the expression of a single gene[END_REF][START_REF] Elowitz | Stochastic gene expression in a single cell[END_REF][START_REF] Raser | Control of stochasticity in eukaryotic gene expression[END_REF] studied under varying parameter values. However, such measurements leave out the temporal aspects of the distributions obtained and potentially miss crucial elements of the dynamics of gene regulation, which technologies such as time-lapse fluorescence microscopy should be able to access. Variation in protein levels in human cells was observed by [START_REF] Sigal | Variability and memory of protein levels in human cells[END_REF], who tracked the fluorescently tagged proteins and concluded that the fluctuations varied slowly in time and that genes of the same pathway showed enhanced temporal correlations. Therefore such time-lapse measurements in single cells resulting in time-series of species could be of immense help in estimating the presence and nature of the regulatory activity between corresponding genes. Due to rapid technological advancements obtaining such fine measurements that are needed to track the causal dynamics as revealed by dynamic correlations is fast becoming a reality. Tagged transcripts as utilized in (Raj et al., 2008) or tagged proteins in [START_REF] Yu | Probing Gene Expression in Live Cells, One Protein Molecule at a Time[END_REF] and [START_REF] Cai | Stochastic protein expression in individual cells at the single molecule level[END_REF] are making single cell, single molecule observations possible, which makes the analytical techniques presented in this paper appropriate and useful. The suggestive clustering of motifs in the feature space defined by time dependent fluctuations provides scope for such optimism.

The rate equations are,

dG y dt = -k on P x G y + k o f f (α -G y ) ( A . 1 ) dM x dt = k + M x -k - M x M x (A.2) dM y dt = k + M y (α -G y ) -k - M y M y (A.3) dP x dt = k + M x M x -k - P x P x -k on P x G y + k o f f (α -G y ) ( A . 4 ) dP y dt = k + P y M y -k - P y P y (A.5) From the above, M y = k + My k - My α P x K D + P x , which follows Michaelis-Menten (hyperbolic) type kinetics. k +
M y is the transcription rate. Note that k basal is excluded in the model as it does not represent the regulatory link between the two genes. The quantities required to solve for the dynamic covariances are the deterministic rates R and the stoichiometry ν which are obtained from the above rate equations.

R = k on P x G y , k o f f (α -G y ), k + M y (α -G y ), k + M x , k + P x M x , k + P y M y , k - M x M x , k - M y M y , k - P x P x , k - P y P y T ν = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -1 +1 0 0 0 0 0 0 0 0 0 0 0 +1 0 0 -1 0 0 0 0 0 +1 0 0 0 0 -1 0 0 -1 +1 0 0 +1 0 0 0 -1 0 0 0 0 0 0 +1 0 0 0 -1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
The relation between A, and the above quantities is given in equation ( 6), where x is the vector of variables [G y , M x , M y , P x , P y ].

A = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -(k o f f + k on P x ) 0 0 -k on G y 0 0 -k - M x 0 0 0 -k + M y 0 -k - M y 0 0 -(k o f f + k on P x ) k + P x 0 -(k - P x + k on G y ) 0 0 0 k + P y 0 -k - P y ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
It now becomes essential to work within the range of the parameter values that is biologically plausible. For this purpose we adopt the {CHA4 → CHA1} regulatory link [START_REF] Bornaes | A regulatory element in the cha1 promoter which confers inducibility by serine and threonine on saccharomyces cerevisiae genes[END_REF]. Information on TFs regulating genes in yeast is available through projects such as YEASTRACT1 [START_REF] Monteiro | YEASTRACT-DISCOVERER: new tools to improve the analysis of transcriptional regulatory associations in Saccharomyces cerevisiae[END_REF]. Cha4p protein is a TF activating the regulation of CHA1. In our model, CHA4 is the regulator gene X and CHA1 is the regulated gene Y, whose mean steady state protein and mRNA levels, and their respective decay rates are obtained from the experimental datasets of (Arava et al., 2003;[START_REF] Ghaemmaghami | Global analysis of protein expression in yeast[END_REF][START_REF] Wang | Precision and functional specificity in mrna decay[END_REF][START_REF] Belle | Quantification of protein half-lives in the budding yeast proteome[END_REF] and are reproduced in Table A.2. The time taken for the yeast cells to divide, denoted as t double is assumed to be 90 minutes.

The mRNA and protein decay rates are therefore:

k - (M,P) = ln(2) (t 1/2 ) (M,P) + ln(2) t double (A.6)
The values of the other parameters, as given in Table A.1, are derived on solving for the deterministic rate equations at steady-state conditions. For example, at steady-state

k + M x = k - M x M x = 0.0407 × 0.2 = 8.14 × 10 -3 nM min -1 .
The OFF and ON rates of the DNA are chosen such that their ratio K D = 200 nM. by the action of repressor P x . M y is transcribed when the promoter region is free of this repressor.

φ k + Mx --→ M x P x + G y k on ------ k o f f C y M y k + Py --→ M y + P y M x k + Px --→ M x + P x G y k + My --→ G y + M y M x , M y k - Mx ,k - My -----→ φ P x , P y k - Px ,k - Py -----→ φ
The deterministic rate equations describing the evolution of average concentration levels of species are:

dG y dt = -k on P x G y + k o f f (α -G y ) dM x dt = k + M x -k - M x M x dM y dt = k + M y G y -k - M y M y dP x dt = k + M x M x -k - P x P x -k on P x G y + k o f f (α -G y ) dP y dt = k + P y M y -k - P y P y
The Jacobian matrix, the stoichiometric matrix and the vector of deterministic rates turn out to be same as in the case of an elementary activator system, with the main exception being the change in the sign of the transcription rates in the Jacobian matrix. This is denoted by an underbrace. This sign change in the offdiagonal element of A is what causes the correlations to go negative, due to the corresponding change in the signs of the eigenvectors.

A = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -(k o f f + k on P x ) 0 0 -k on G y 0 0 -k - M x 0 0 0 k + M y 0 -k - M y 0 0 -(k o f f + k on P x ) k + P x 0 -(k - P x + k on G y ) 0 0 0 k + P y 0 -k - P y ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ν = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -1 +1 0 0 0 0 0 0 0 0 0 0 0 +1 0 0 -1 0 0 0 0 0 +1 0 0 0 0 -1 0 0 -1 +1 0 0 +1 0 0 0 -1 0 0 0 0 0 0 +1 0 0 0 -1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ R = k on P x G y , k o f f (α -G y ), k + M y G y , k + M x , k + P x M x , k + P y M y , k - M x M x , k - M y M y , k - P x P x , k - P y P y T Appendix B.2

. Activation via Dimers

The reaction scheme is the same as that of elementary activation with the additional reactions of dimerization 2P

dG y dt = k o f f (α -G y ) -k on G y P x 2 dM x dt = k + M x -k - M x M x dM y dt = k + M y (α -G y ) -k - M x M x dP x dt = k + P x M x -k - P x P x + 2k b P x 2 -2k a P 2 x dP x 2 dt = k o f f (α -G y ) -k on G y P x 2 -k b P x 2 + k a P 2 x dP y dt = k + P y M y -k - P y P y A = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ (-k o f f -k on P x 2 ) 0 0 0 -k on G y 0 0 -k - M x 0 0 0 0 -k + M y 0 -k - M y 0 0 0 0 k + P x 0 (-4k a P x -k - P x ) 2k b 0 (-k o f f -k on P x 2 ) 0 0 2k a P x (-k b -k on G y ) 0 0 0 k + P y 0 0 -k - P y ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ν = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝
-1 +1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +1 0 0 -1 0 0 0 0 0 +1 0 0 0 0 0 0 -1 0 0 0 0 0 +2 -2 0 +1 0 0 0 -1 0 -1 +1 0 -1 +1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +1 0 0 0 -1

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ R = k on G y P x 2 , k o f f (α -G y ), k + M y (α -G y ), k b P x 2 , k a P 2 x , k + M x , k + P x M x , k + P y M y , k - M x M x , k - M y M y , k - P x P x , k - P y P y T Appendix B.3.

Cascading Activation

The set of reactions describing the cascaded activation system, that employs dual activators, is as follows:

φ k + Mx --→ M x M x k + Px --→ M x + P x P x + G z k on 1 ------ k o f f 1 C z C z k + Mz --→ C z + M z P z + G y k on 2 ------ k o f f 2 C y C y k + My --→ C y + M y M z k + Pz --→ M z + P z M y k + Py --→ M y + P y M x , M z , M y k - Mx ,k - Mz ,k - My --------→ φ P x , P z , P y k - Px ,k - Pz ,k - Py -------→ φ
In the following set of rate equations, the redundant variables C z and C y are eliminated and consequently the terms α and β denote the initial concentrations G z (t 0 )

and G y (t 0 ) respectively:

dG z dt = -k on 1 G z P x + k o f f 1 (α -G z ) dG y dt = -k on 2 G y P z + k o f f 2 (β -G y ) dM x dt = k + M x -k - M x M x dM z dt = k + M z (α -G z ) -k - M z M z dM y dt = k + M y (β -G y ) -k - M y M y dP x dt = k + P x M x -k - P x P x -k on 1 G z P x + k o f f 1 (α -G z ) dP z dt = k + P z M z -k - P z P z -k on 2 G y P z + k o f f 2 (β -G y ) dP y dt = k + P y M y -k - P y P y A = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -(k o f f 1 + k on 1 P x ) 0 0 0 0 -k on 1 G z 0 0 0 -(k o f f 2 + k on 2 P z ) 0 0 0 0 -k on 2 G y 0 0 0 -k - Mx 0 0 0 0 0 -k + Mz 0 0 -k - Mz 0 0 0 0 0 -k + My 0 0 -k - My 0 0 0 -(k o f f 1 + k on 1 P x ) 0 k + Px 0 0 -(k - Px + k on 1 G z ) 0 0 0 -(k o f f 2 + k on 2 P z ) 0 k + Pz 0 0 -(k - Pz + k on 2 G y ) 0 0 0 0 0 k + Py 0 0 -k - Py ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ν = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -1 +1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 +1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +1 -1 0 0 0 0 -1 +1 0 0 0 0 0 0 0 0 0 0 +1 -1 0 0 0 0 -1 +1 0 0 0 0 0 0 0 0 0 0 +1 -1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ R = [k on 1 G z P x , k o f f 1 (α -G z ), k + My (α -G z ), k on 2 G y P z , k o f f 2 (β -G y ), k + My (β -G y ), k + Mx , k + Px M x , k + Pz M z , k + Py M y , k - Mx M x , k - Mz M z , k - My M y , k - Px P x , k - Pz P z , k - Py P y ] T Appendix B.

Co-operative Activation

The following is the set of reactions describing the regulatory system, where two activators P x and P z act on the regulatory region upstream of G y , either individually or through combined action:

φ k + Mx ,k + Mz -----→ M x , M z M x k + Px --→ M x + P x M z k + Pz --→ M z + P z P x + G y k on 1 ------ k o f f 1 C y 1 P z + G y k on 2 ------ k o f f 2 C y 2 P z + C y 1 k on 3 ------ k o f f 3 C y 3 C y i k + My i ---→ C y i + M y (i = 1, 2, 3) M y k + Py --→ M y + P y M x , M z , M y k - Mx ,k - Mz ,k - My --------→ φ P x , P z , P y k - Px ,k - Pz ,k - Py -------→ φ
Once again in the rate equations, G y is chosen as the redundant variable, while the TF-DNA complexes C y i are retained in the model.

dC y 1 dt = k on 1 P x (α -C y 1 -C y 2 -C y 3 ) -k o f f 1 C y 1 + k o f f 3 C y 3 -k on 3 C y 1 P z dC y 2 dt = k on 2 P z (α -C y 1 -C y 2 -C y 3 ) -k o f f 2 C y 2 dC y 3 dt = k on 3 C y 1 P z -k o f f 3 C y 3 dM x dt = k + M x -k - M x M x dM z dt = k + M z -k - M z M z dM y dt = k + M y 1 C y 1 + k + M y 2 C y 2 + k + M y 3 C y 3 -k - M y M y dP x dt = k o f f 1 C y 1 -k on 1 P x (α -C y 1 -C y 2 -C y 3 ) + k + P x M x -k - P x P x dP z dt = k o f f 2 C y 2 -k on 2 P z (α -C y 1 -C y 2 -C y 3 ) + k o f f 3 C y 3 -k on 3 C y 1 P z + k + P z M z -k - P z P z dP y dt = k + P y M y -k - P y P y
As in previous cases, it is quite straightforward to derive the Jacobian matrix A, the vector of deterministic rates R and the stoichimetric matrix ν from the above set of deterministic dynamical equations.

Appendix B.5. Feedforward Loops

Another class of well-defined regulatory networks are the feedforward loops consisting of three genes where the downstream gene is regulated by two different TFs, out of which one TF acts as the regulator of the other TF. Thus a loop is formed. The type of regulation can either be that of activation or repression. If for example, both the TFs act as activators of the downstream protein, such a network is called as the coherent feedforward loop. On the other hand, if one is a repressor while the other an activator, the incoherent FFL is formed. It has been reported that the feedforward loop is a recurring network in organisms such as both Escherichia coli and the yeast [START_REF] Shen-Orr | Network motifs in the transcriptional regulation network of Escherichia coli[END_REF][START_REF] Kalir | A coherent feed-forward loop with a sum input function prolongs flagella expression in Escherichia coli[END_REF][START_REF] Mangan | The incoherent feedforward loop accelerates the response-time of the gal system of Escherichia coli[END_REF]. is:

k on G * + 2(-k on -k o f f )C 11 + k o f f (α -G * ) = 0 which on solving gives the auto-variance of G * , Var[G * (t), G * (t)] = (k on -k o f f ) G * + G(t 0 )k o f f 2(k on + k o f f ) = G * k o f f k on + k o f f
Continuing on similar lines, the covariances between all the variables are evaluated. Whiled doing so, we make approximations that the decay rates are much smaller when compared to the ON and OFF rates of the gene (

k on + k o f f ) k - (M,P)
and that k o f f > k ON . Both these approximations are valid for biologically plausible values of rate constants [START_REF] Bundschuh | The role of dimerization in noise reduction of simple genetic networks[END_REF].

Cov[G * (t), M(t)] = M k - M k o f f (k on + k o f f + k - M )(k on + k o f f ) ≈ M k - M k o f f (k on + k o f f ) 2 ≈ M k - M k o f f Var[M(t), M(t)] = M + M k + M k o f f (k on + k o f f + k - M )(k on + k o f f ) ≈ M + M k + M k o f f ≈ M Cov[M(t), P(t)] = M k + P k - M + k - P 1 + k + M k o f f (k on + k o f f )(k on + k o f f + k - M ) 1 + k - M (k on + k o f f + k - P ) ≈ M k + P k - M + k - P 1 + k + M k o f f (k on + k o f f ) 2 1 + k - M (k on + k o f f ) ≈ M k + P k - M + k - P 1 + k + M k o f f 1 + k - M k o f f ≈ M k + P k - M + k - P 1 + k + M k o f f ≈ M k + P k - M + k - P Var[P(t), P(t)] = P + k + P k - P Cov[M(t), P(t)] ≈ P + P k + P k - M + k - P
It has been observed in various real biological systems that the noise in the mRNA distribution which is nothing but its variance over mean squared has sometimes the characteristics of a Poissonian process while in other cases is far from it. A Poisson process being reflected from the fact that the distribution resulting out of such a process has its variance equal to its mean. Therefore in the case where the gene is constantly in the ON state, the transcription process results in mRNA molecules that are Poissonian distributed, while in the above described gene ON/OFF model the distribution of mRNAs and proteins are skewed and are non-Poissonian in character. However, with the approximation of (k on + k o f f ) k - (M,P) and k o f f > k on the expressions for the stationary covariance terms are exactly as would be obtained in the model where the gene is constantly in the ON state and the mRNAs are being produced by random birth-death events leading to a Poissonian distribution. Previous works [START_REF] Thattai | Intrinsic noise in gene regulatory networks[END_REF][START_REF] Paulsson | Summing up the noise in gene networks[END_REF] have involved investigating the stationary auto-covariance which is, as seen above, a function of the system parameters. However, as our aim is to include the additional factor of time into the analyses, the dynamic covariance between the mRNA and the protein variables are employed as representative of the system's internal fluctuations. In order to determine the time-covariance between the mRNA and protein Cov[M(t), P(t + τ)] we need the eigenvectors of the system. The right eigenvectors of the Jacobian are the column vectors of the following matrix:

V = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 (k on +k o f f -k - M )(k on +k o f f -k - P ) k + M k + P 0 -(k - M -k - P ) k + P -(k on +k o f f -k - P ) k + P 0 1 1 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
and by the bi-orthogonal property of the right and left eigenvectors U T V = I, we get the left eigenvectors as:

U T = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -k + M k + P (k on +k o f f -k - M )(k - M -k - P ) -k + P (k - M -k - P ) 0 k + M k + P (k on +k o f f -k - M )(k on +k o f f -k - P ) 0 0 k + M k + P (k on +k o f f -k - P )(k - M -k - P ) k + P (k - M -k - P ) 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
The corresponding eigenvalues are:

λ = [-k - M , -(k on + k o f f ), -k - P ]
Therefore the time-covariance between the mRNA and protein species is:

Cov[M(t), P(t + τ)] = Cov[X 2 (t), X 3 (t + τ)]

= 3 l=1 e λ l τ V 3l 3 k=1 U T lk C k2 = e -k - P τ k + M k + P (k on + k o f f -k - M )(k - M -k - P ) C 12 + k + P (k - M -k - P ) C 22 -e -k - M τ k + M k + P (k on + k o f f -k - M )(k - M -k - P ) C 12 + k + P (k - M -k - P ) C 22 + C 32 + e -(k on +k o f f )τ k + M k + P (k on + k o f f -k - M )(k on + k o f f -k - P ) C 12
Since (k on + k o f f ) is assumed to be much larger than the decay rates, we neglect the fast decaying exponential e -(k on +k o f f )τ and its allied term leading to, Cov[M(t), P(t + τ)] ≈ e -k -

P τ k + M k + P (k on + k o f f -k - M )(k - M -k - P ) C 12 + k + P (k - M -k - P ) C 22 -e -k - M τ k + M k + P (k on + k o f f -k - M )(k - M -k - P ) C 12 + k + P (k - M -k - P ) C 22 + C 32
On substituting for the C terms and doing approximations, Cov[M(t), P(t

+ τ)] ≈ M k + P (k - M -k - P ) e -k - P τ 2k - M k - M + k - P -e -k - M τ
This function is a sum or rather difference of two exponentials which in all possibility could be a non-monotonic function. To find if this is truly the case, differentiating the above covariance term partially w.r.t τ and equating to 0, yields the τ = τ * at which the covariance reaches peak magnitude. Since this time delay τ is shown to be a function of the decay rates, therefore the temporal character of the covariances can be controlled via these parameters. The interesting aspect of this time-delay is that it is related to the time delay in the response of the mean value of the regulated variable to a perturbation in the mean value of the regulator. For a perturbation of ΔM in the mRNAs, the response in the proteins is,

P(t) = P(t ss ) + k + P (ΔM) k - M -k - P e -k - P t -e -k - M t
This response in the protein level also has the characteristic peak after a time-delay of (t resp ), which is evaluated by equating the time-evolution of P to zero,

∂ P(t) ∂t = 0 k - M e -k

Figure 2 :

 2 Figure 2: Dynamic correlations between the proteins of the X → Y regulatory network, obtained through analytics and simulations.

Figure 4 :

 4 Figure3: The inset in (a) depicts the networks used in this paper: (X → Y) the single activator; (2X → Y) a single activator network, but with X coding for a protein which performs the activation after forming a homodimer; (X → Z → Y) which describes a cascade of interactions,X upregulating Z which upregulates Y; (X-Z-Y) cascading repressors; (X → Y ← Z) the combinatorial action of genes X and Z on Y; (coherent FFL) which combines X → Y ← Z and X → Z; (incoherent FFL) which combines X → Y -Z and X → Z (a) (t 1/2 ) P x is varied in steps of[40, 50, basevalue, 90, 120] minutes, where the basevalue corresponds to a half-life of 65 minutes. (b) Similarly, (t 1/2 ) P y is varied in steps of[40, 50, basevalue, 90, 120] minutes, where the basevalue corresponds to a half-life of 70 minutes. (c) The transcription rate or the basal transcription rate (in case of repression) k + M y is varied as [0.6, 0.8, 1.0, 1.2, 1.4] × basevalue. (d) The rate of protein translation k + P x is varied in steps as [0.6, 0.8, 1.0, 1.2, 1.4] × basevalue.

Figure 6 :

 6 Figure 6: (a) Transcription rate constant k + My , and k - My are varied simultaneously in steps, such that M y remains constant. While Corr * does change by an appreciable amount, there is a huge decrease in τ * , from 87 to 43 minutes. (b) Similar is the case for variations in k + Py and k - Py . There is again a large reduction in τ * , from 102 to 36 minutes accompanied by noticable change in the shape of the plots, as was in the case of [k + Px , k - Px ]. The observed variation in τ * is due to the decay rate rather than the translation rate.

Figure 7 :

 7 Figure 7: Schematic represention of the elementary repressor system. The basal transcription levels on the left is diminished by the action of the repressor.

  (a)), as P x

Figure 8 :

 8 Figure 8: Elementary repression X Y. [k +P y , k - P y ] are varied simultaneously so that the mean concentration level of P y = M y k + P y /k - P y remains constant. While there is an increase in Corr * , τ * decreases from 102 to 36 minutes for step-wise increase in these parameters.

Figure 9 :

 9 Figure 9: Schematic represention of the activator system where dimers are involved. As before, the figure on the left indicates basal transcription.

Figure 10 :

 10 Figure 10: Activation via Dimers 2X → Y. (a) k -P x is varied in steps from its base value of 0.0184 min -1 that corresponds to a half-life of (t 1/2 ) P x = 65 minutes. k + P x is simultaneously varied such that P x remains unchanged. For increase in the values of [k + P x , k - P x ], there is an observable increase in Corr * , whilst τ * decreases from 53 to 46 minutes. (b) Likewise, for step-wise increase in the parameters [k + P y , k - P y ], τ * once again reduces from 81 to 20 minutes, whilst Corr * is nearly constant.

Figure 11 :

 11 Figure 11: Activation via dimers 2X → Y. The base value of the dimer dissociation constant K Dim is fixed at 200 nM, and is further assigned values of 40, 80, 600 and 1400 nM. An equivalent variation in the value of k + M y is done such that M y remains unchanged. For increase in the values of [k + M y , K Dim ], τ * decreases marginally from 53 to 47 minutes.

Figure 13 :

 13 Figure 13: (a) Protein stationary correlation varies for changes in the decay rates whilst keeping the other parameters fixed at their base values. The decay rates are varied from {0.05 × base value} to {30 × base value}, where base values are those of Table A.1. (b) The transcription, translation and binding/unbinding rates too have a corresponding influence on the stationary correlation.

  Figure B.14: Coherent and Incoherent FeedForward Loops. Gene activation is denoted by an arrowhead and repression by a dashhead.

-

  Mx remains constant. While Corr * decreases significantly for increase in the two parameters, τ * registers miminal change from 56 to 46 minutes. (b) The translation and protein decay rates k + Px and k - Px are varied simultaneously so that P x = M x k + Px /k - Px remains constant. While there is not much difference in Corr * , τ * decreases from 55 to 45 minutes for increase in these parameters. The half-life of M x and the rate constant k - Mx are related as follows: k -

Mx = ln(2) × 1/(t 1/2 ) Mx + 1/t double , where t double = 90 minutes. We use this definition (as in eq. (A.6)) in the other figures as well.

At steady state, mean value of G y

Table A .

 A 2: The values are of the CHA4 → CHA1 activator link in yeast. Nuclear volume is taken to be around 1 μm 3 , resulting in nano-Molar concentrations for the mRNA and protein species. The mRNA and protein half-lives t 1/2 are in minutes.

	Std Name Sys Name mRNA protein mRNA t 1/2 protein t 1/2
		(nM)	(nM)	(min)	(min)
	CHA4	YLR098C 0.2	395.6	21	65
	CHA1	YCL064C 4.5	36602 10	70

  - M t respk - P e -k - P t resp = 0 Alon, U., Jun 2007. Network motifs: theory and experimental approaches. Nature Reviews Genetics 8, 450-461. Arava, Y., Wang, Y., Storey, J. D., Liu, C. L., Brown, P. O., Herschlag, D., 2003. Genome-wide analysis of mrna translation profiles in saccharomycescerevisiae. PNAS 100 (7), 3889-3894. Arkin, A., Ross, J., McAdams, H. H., Aug 1998. Stochastic kinetic analysis of developmental pathway bifurcation in phage λ-infected Escherichia coli cells. Genetics 149, 1633-1648.

			t resp =	ln k -M -k -k -M k -P P
	This is smallest for the case where k -M = k -P and is evaluated using l'Hopital's rule. Letting k -M /k -P = x, the expression for t resp is now,
	lim x→1	t resp = lim x→1	1 k -P	ln(x) (x -1)
		= lim x→1	1 k -P	1/x 1
		=	1 k -P

http://www.yeastract.com

Appendix A. Elementary Activator

The reaction set of an elementary activator X → Y is given in Table A.1. The regulatory sequence upstream of gene G y is bound by the transcription factor P x , which is an activator in the present case, resulting in complex C y which is in effect the active or ON state of the gene. Transcription of M y is initiated from this complex. The deterministic dynamics of this system are governed by a set of

Table A.1: Reaction set describing the process of activation between X and Y. Production of M x is assumed to be of Poissonian birth-death process. φ → M x therefore denotes spontaneous creation of M x from a constant source φ. nM stands for nano-molar concentration and min is minutes.

coupled ODEs that describe the time-evolution of the mean concentration levels of the species. The total amount of G y present in the cell at time

Assuming the presence of a single copy of the gene in the yeast nucleus of volume 1 μm 3 , G y ≈ 1 nM. Therefore, at steady state, C y = α P x K D + P x and G y = αK D K D + P x . The redundancy in the original set of ODEs is eliminated by substituting for C y .

Appendix B. Reaction schemes and rate equations of other regulatory mechanisms Appendix B.1. Elementary Repressor

The reaction scheme for this system is exactly the same as for the case of the two-gene elementary activation X → Y, except that the trancription of G y is now through a basal transcription rate of k + M y . Gene G y is switched to inactive state C y Appendix C. Expression for τ * in the case of a single-gene.

The best way to give an intuition of the time-covariance function is by applying it to the simplest of cases, that of a single gene. Here, we derive the timecovariance function between the mRNA and protein molecules in the case of a single gene. The model includes a single gene which spontaneously switches to the ON or active state G * at a rate k on and back to the OFF or inactive state G at the rate k o f f . Therefore the stationary distribution of G * is of the Binomial type. In the ON state the mRNA M is transcribed out of the gene at a rate of k + M and is further translated to the protein P at a rate k + P . The reactions steps for such a system are:

Let G * represent the average amount of active genes in concentration form and in units of nano-Molar. Then the rate of production of G * follows the deterministic dynamics of the system given by,

Due to the conservation of the gene molecules,

, where G(t 0 ) is the initial amount of the DNA molecule present in the medium. Using this substitution in the rate equation for G * we have the following rate equations of the three variables describing their deterministic behaviour:

Taking the vector of the system variables as X = [G * , M, P], the corresponding Jacobian matrix derived from the rate equations is:

Further, the vector of the deterministic rates and the stoichiometric matrix are also derived from the same rate as:

The diffusion matrix given by ν diag(R( )) ν T is:

The next step is to solve for the stationary covariance matrix:

where Var and Cov represent the stationary auto-covariance and covariances. Each element of the above matrix is evaluated by solving the Lyapunov equation systematically element-wise. For example, the first element of AC + CA T + BB T = 0