Eugenio Cinquemani
email: eugenio.cinquemani@inria.fr

On observability and reconstruction of promoter activity statistics from reporter protein mean and variance profiles

Keywords: Gene regulation, Doubly stochastic process, Spectral analysis

Reporter protein systems are widely used in biology for the indirect quantitative monitoring of gene expression activity over time. At the level of population averages, the relationship between the observed reporter concentration profile and gene promoter activity is established, and effective methods have been introduced to reconstruct this information from the data. At single-cell level, the relationship between population distribution time profiles and the statistics of promoter activation is still not fully investigated, and adequate reconstruction methods are lacking. This paper develops new results for the reconstruction of promoter activity statistics from mean and variance profiles of a reporter protein. Based on stochastic modelling of gene expression dynamics, it discusses the observability of mean and autocovariance function of an arbitrary random binary promoter activity process. Mathematical relationships developed are explicit and nonparametric, i.e. free of a priori assumptions on the laws governing the promoter process, thus allowing for the decoupled analysis of the switching dynamics in a subsequent step. The results of this work constitute the essential tools for the development of promoter statistics and regulatory mechanism inference algorithms.

Introduction

A common experimental technique to monitor gene expression is the use of reporter proteins [START_REF] De Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF], i.e. fluorescent or luminescent proteins that are synthesized upon expression of the gene of interest. Light intensity measurements collected at different points in time are proportional to the amount of reporter molecules. This provides a quantitative, however indirect, readout of the activity of the gene, since reporter abundance depends on gene activation via its own transcription and translation dynamics.

When cellular populations are observed as a whole, such as in automated microplate readers, an average reporter profile is obtained. An estimate of the average gene activation over the population of cells may thus be obtained by regularized inversion of the reporter synthesis dynamics [START_REF] Zulkower | Robust reconstruction of gene expression profiles from reporter gene data using linear inversion[END_REF]. Provided accurate knowledge of the latter, reconstruction of the promoter activity allows one to investigate gene expression regulatory mechanisms, a crucial step toward inference of gene regulatory networks [START_REF] Stefan | Inference of quantitative models of bacterial promoters from time-series reporter gene data[END_REF].

When individual cells are observed, for instance via flow-cytometry or fluorescence videomicroscopy, a statistical distribution of gene expression levels over a sample of the population (often called a population snapshot [START_REF] Hasenauer | Identification of models of heterogeneous cell populations from population snapshot data[END_REF]) is obtained at several points in time (reporter traces for individual cells can also be obtained by suitable experimental setups and image processing techniques [START_REF] Zechner | Scalable inference of heterogeneous reaction kinetics from pooled single-cell recordings[END_REF], but we will not analyze this case here). In many cases of interest, this crucially reveals variability of gene expression levels across cells that can be explained in terms of the stochasticity of the gene regulation and expression process [START_REF] Paulsson | Models of stochastic gene expression[END_REF][START_REF] Zechner | Moment-based inference predicts bimodality in transient gene expression[END_REF][START_REF] Munsky | Listening to the noise: Random fluctuations reveal gene network parameters[END_REF]. Reporter statistics thus contains information about the stochastic laws governing gene activation. However, recovering the relevant information from the data is less trivial than in the population average case, and no satisfactory methods exist to date.

With reference to population snapshot data, in [START_REF] Cinquemani | Reconstruction of promoter activity statistics from reporter protein population snapshot data[END_REF][START_REF] Cinquemani | Statistics of Promoter Switching from Reporter Protein Population Snapshot Data[END_REF], we have started addressing the problem of estimating promoter activity statistics (the biological information of interest in gene expression reporting) from reporter mean and variance profiles. In [START_REF] Cinquemani | Reconstruction of promoter activity statistics from reporter protein population snapshot data[END_REF], parametric models of stochastic gene activation have been considered, and the identifiability of promoter switching rates that are fixed over time and across cells has been analyzed. However, due to a priori unknown regulatory mechanisms, switching rates may fluctuate over time and/or across cells (extrinsic noise). To cope with this, in [START_REF] Cinquemani | Statistics of Promoter Switching from Reporter Protein Population Snapshot Data[END_REF], a nonparametric method, i.e. avoiding assumptions on the regulatory mechanisms behind the expression of the gene of interest, has been proposed for the special case of irreversible activation. A rather extensive account of relevant research literature is also contained in these works.

Following up from the developments in [START_REF] Cinquemani | Reconstruction of promoter activity statistics from reporter protein population snapshot data[END_REF][START_REF] Cinquemani | Statistics of Promoter Switching from Reporter Protein Population Snapshot Data[END_REF], for the general case of unmodelled stochastic (possibly time varying) gene expression regulation, we address here the problem of reconstructing second-order statistics of the promoter activity process from reporter mean and variance profiles. The importance of this problem lies in the fact that, in analogy with linear stochastic processes [START_REF] Lindquist | Linear stochastic systems -A geometric approach to modeling, estimation and identification[END_REF], cross-correlation of promoter activity at different points in time (i.e. the autocorrelation function) contains information about the time dynamics of activation and deactivation. Reconstruction of these statistics from data is thus the crucial step for the understanding of the gene regulatory mechanisms at the level of single cells, where stochastic variability offers more to discover than traditional population analysis [START_REF] Munsky | Listening to the noise: Random fluctuations reveal gene network parameters[END_REF][START_REF] Neuert | Systematic identification of signal-activated stochastic gene regulation[END_REF].

The contribution of this paper is the development of explicit relationships between the unknown (first-and) second-order promoter activity statistics and the experimentally measurable reporter mean and variance profiles. Crucially, these relationships rely on nonparametric models of gene activation, i.e. no a priori assumption is made except the absence of stochastic feedback from re-porter abundance to the regulation of the gene itself, a hypothesis that agrees well with the biochemistry of reporter systems. Based on analytic investigation and examples, we show that these relationships are essentially linear, whence invertible in a tractable manner, and allow for the discrimination among different promoter activity regulatory laws. On these basis, the implementation of algorithms for the actual estimation of the statistics of interest is left for future work. For ease of reading, all mathematical proofs are deferred to Appendix A. Appendix B, instead, summarizes results from [START_REF] Cinquemani | Reconstruction of promoter activity statistics from reporter protein population snapshot data[END_REF] that are used in this work.

Background material

Gene expression monitoring over time is commonly operated by the use of fluorescent or luminescent reporter proteins (see [START_REF] De Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF] and references therein). In essence, synthesis of a reporter protein is placed under the control of the promoter of the gene of interest by engineering its coding sequence onto the DNA at an appropriate place. When the gene is expressed, transcription and subsequent translation leads to the formation of new reporter protein molecules. Whether luminescent or fluorescent, reporter protein molecules can be quantified at any time by measuring light intensity at the relevant wavelength, thus providing a dynamical readout of the activity of the gene. To do so, time-lapse microscopy, flow cytometry, microplate reading, or other experimental techniques are used, depending whether single-cell measurements, population histograms, or population-average profiles are sought. Synthesis of reporter proteins is often completed by a maturation step, that takes immature proteins into their mature, visible form.

Stochastic gene expression modelling

Gene expression is commonly described in terms of the synthesis and degradation reactions for mRNA and protein molecules

R 1 : F k M --→ F + M R 2 : M d M --→ ∅ (1)
R 3 : M k P --→ M + P R 4 : P d P --→ ∅ (2)
[5, 10] where M and P denote mRNA and protein species, respectively, and F represents the active promoter species. In the context of this paper, P is the fluorescent or luminescent reporter protein. We will not distinguish between immature (invisible) and mature (visible) protein molecules. If necessary (e.g. for slow, stochastic maturation), an additional first-order reaction P → P mature can be included in the model (along with P mature → ∅) to account for protein maturation (and mature protein degradation). Denote with X 1 ∈ N and X 2 ∈ N the number of copies of M and P , in the same order, and with X 3 ∈ {0, 1} the state of the promoter, i.e. X 3 = 0 when the promoter is inactive (absence of F ) and X 3 = 1 when it is active (presence of F ). Switching promoter dynamics (responsible of mRNA synthesis bursts in single cells) are formally captured by two additional reactions,

R 5 : ∅ λ+•(1-X3) -------→ F, R 6 : F λ- --→ ∅, (3)
representing in the order activation with propensity λ + • (1 -X 3 ) (only enabled if X 3 = 0), and deactivation with propensity λ -• X 3 (only enabled if X 3 = 1). Overall, this is a system of m = 6 chemical reactions over n = 3 different species.

The kinetics of this biochemical reaction system can be expressed in terms of stoichiometry matrix S and reaction rate vector a(x) given by

S =   1 -1 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 1 -1   , a(x) =         k M x 3 d M x 1 k P x 1 d P x 2 λ + (1 -x 3 ) λ -x 3        
, where, for i = 1, . . . , n and j = 1, . . . , m, S i,j denotes the net change in molecule number of species i when reaction R j occurs. At the level of a single cell, X = [X 1 X 2 X 3 ] T is a stochastic process and, for every j, a j x is interpreted as the infinitesimal probability that reaction R j occurs in an infinitesimal time period when X = x molecules of the different species are present in the reaction volume [START_REF] Paulsson | Models of stochastic gene expression[END_REF]. For constant rates λ + and λ -, Eq.( 1)-( 3) together constitute the so-called random telegraph model [START_REF] Paulsson | Models of stochastic gene expression[END_REF]. In general, however, these rates might themselves depend upon the amount of transcription factors regulating the expression of the gene, which one may write as λ + (X ? ) and λ -(X ? ), with X ? denoting the amount of some unspecified species.

The question we are going to investigate is what can be said about the statistics of F , given mean and variance profiles of the amounts of protein P across a population of cells. In practice, fluorescence or luminescence measurements proportional to the actual amount of protein are measured and are possibly affected by error. In this paper, however, we are not concerned with the details of the measurement model, and assume that mean and variance of X 2 are observed directly.

Propagation of moments

Consider an arbitrary biochemical reaction system with n reactants, m reactions, stoichiometry matrix S and reaction rates a(x, u) possibly depending on a deterministic input u. Let X(t) be the corresponding random state vector at time t, and define

µ(t) = E[X(t)] and Σ(t) = Cov X(t) = E X(t) -µ(t) X(t) - µ(t)
T . It can be shown (see e.g. [START_REF] Hespanha | Modelling and analysis of stochastic hybrid systems[END_REF]) that µ and Σ obey the so-called moment

equations μ = SE[a(X, u)], (4)
Σ = SE[a(X, u)(X -µ) T ] + E[(X -µ)a T (X, u)]S T + Sdiag(E[a(X, u)])S T . ( 5
)
Above and in the sequel, time t is omitted from notation where no confusion may arise. If rates are affine in the state, i.e. a(x, u) = W (u)x + w 0 (u) for some W (u) and w 0 (u), these equations simplify to

μ = SW (u)µ + Sw 0 (u), (6)
Σ = SW (u)Σ + ΣW T (u)S T + Sdiag W (u)µ + w 0 (u) S T . (7)
This system of differential equations is closed in the sense that it does not depend on unmodelled moments. If in addition W does not depend on u, then the system is linear in the input (and the initial conditions).

For the system (1)-( 3), Eqs. ( 6)-( 7) apply in the case of constant rates λ + and λ -. In the general case of regulated switching rates λ + (X ? ) and λ -(X ? ), one may instead interpret ( 4)-( 5) as the moment equations for the augmented state composed of X and X ? . Since the laws regulating X ? are unspecified, the full system cannot be spelled out, but one may still work out the equations for the evolution of the moments of X 1 , X 2 and X 3 . Define

z T M P | z T × | z T F = µ M µ P σ M M σ P P σ M P | σ M F σ P F | µ F σ F F ,
(vertical bars denoting vector blocks) where of course µ • and σ •• are the mean and covariance of the states corresponding to the species in subscript (identical subscripts denoting variance). From an engineering viewpoint, z M P is the state of the dynamical sensor for the statistics of F , with sensor output given by the elements [µ P σ P P ] T of z M P . Then one gets

żMP = A M P • z M P + A M P,× • z × + A M P,F • z F , (8) ż×
= A ⊗ • z × + z ⊗ + A ×,F • z F . (9)
for matrices A M P , A M P,× , A M P,F , A ⊗ and A ×,F depending solely on θ M P = (k M , d M , k P , d P ) (see Appendix A), i.e. the parameters of the sensing system. Note that, for every fixed t, F (t) is a Bernoulli random variable. Then σ F F (t) = µ F (t) 1 -µ F (t) for all t (as a consequence, ( 8)-( 9) are somewhat redundant). From ( 8)-( 9) one observes that mean and variance of X 2 , the observed elements of z M P , are thus a dynamical transformation of those of F , i.e. z F , plus a contribution from

z ⊗ = Cov X 1 X 2 , λ + (X ? )(1 -X 3 ) λ -(X ? )X 3 • 1 -1 .
As it will become clear, z ⊗ implicitly brings about a contribution from the correlation structure of F (see later Remark 1).

Marginalization of moments

From now on, abusing notation in favor of simplicity, we will refer to X 1 , X 2 and X 3 by the symbols for the corresponding species, i.e. M , P and F , in the same order. Let f be any possible outcome of F , and let

µ P (t) = E[P (t)], µ f P (t) = E[P (t)|F = f ], M P (t) = E[P (t) 2 ], M f P (t) = E[P (t) 2 |F = f ],
where, unlike the approach in [START_REF] Hasenauer | Method of conditional moments (MCM) for the Chemical Master Equation[END_REF], conditioning is intended over the whole history of F . By marginalization,

µ P = E E[P |F ] = µ f P dP F (f ), M P = E E[P 2 |F ] = M f Y dP F (f ), (10)
with P F the probability distribution of F over all possible binary switching sequences. Let us now state the following assuption.

Assumption 1 (Granger causality [START_REF] Lindquist | Linear stochastic systems -A geometric approach to modeling, estimation and identification[END_REF]) There is no feedback from M and P to F , i.e., at any time t, the future of F is conditionally independent on the past of M and P given the past of F .

This captures the idea that species M and P do not participate in the regulation of the promoter [START_REF] Cinquemani | Statistics of Promoter Switching from Reporter Protein Population Snapshot Data[END_REF][START_REF] Bowsher | The fidelity of dynamic signaling by noisy biomolecular networks[END_REF], and corresponds well to all the reporter systems where reporter and regulatory proteins are physically different molecules. In the light of Assumption 1, the conditional moments µ f P and M f P are those of the reduced system (1)-( 2) with f defining the state of species F at all times. Let

z f M P = µ f M µ f P σ f M M σ f P P σ f M P
be the vector of conditional moments of M and P . Working out the moment equations ( 6)- [START_REF] Hasenauer | Method of conditional moments (MCM) for the Chemical Master Equation[END_REF] for X = M P T and input u = f , one gets that

żf M P = A M P • z f M P + (A M P,F ) 1 • f, (11)
where (A M P,F ) 1 denotes the first column of A M P,F . Then µ f P and σ f P P follow from the solution of this system and M f P = σ f P P + (µ f P ) 2 , while marginalization [START_REF] Kaern | Stochasticity in gene expression: From theories to phenotypes[END_REF] completes the computation of µ P and M P . Note that, because of the relationship M P = σ P P + (µ P ) 2 , we can equivalently consider (µ P , σ P P ) or (µ P , M P ) to be the observed output quantities. We will often exploit this equivalence in the sequel without further notice.

Incidentally, notice that (11) represents a linear switching system with two alternating operational modes, f = 0 and f = 1.

The fixed rate promoter process

In order to investigate how statistics of F reflect into the observed profiles µ P and M P , and how they may possibly be reconstructed from the output, we first focus on the fundamental case where switching rates λ + and λ -are constant.

Define α = λ + + λ -. Proposition 1. Mean µ F (t) = E[F (t)] and autocovariance function ρ F (t, s) = cov F (t), F (s) obey the equations µ F (t) = µ F (0)e -αt + λ + α 1 -e -αt , t ≥ 0, ( 12
)
ρ F (t, τ ) = λ + α + (α -λ + ) α e -α(t-τ ) • µ F (τ ) -µ F (t) • µ F (τ ), t ≥ τ. ( 13
)
In stationary conditions, with an abuse of notation for the arguments of ρ F ,

µ F = λ + α , ρ F (t -τ ) = λ + (α -λ + ) α 2 • e -α(t-τ ) . ( 14
)
Incidentally, the autocovariance in ( 14) is the same as that of an Ornstein-Uhlenbeck process [START_REF] Papoulis | Probability, random variables, and stochastic processes[END_REF]. It can be appreciated that, in transient conditions, the mean profile µ F contains all the information about the statistics of F . Indeed, in this simple case, rates λ + and λ -(or equivalently α), together with the initial condition µ F (0), fully determine the laws of F . In turn, these three quantities have distinct effects on µ F , i.e. they are distinguishable from a transient mean profile. In [START_REF] Cinquemani | Reconstruction of promoter activity statistics from reporter protein population snapshot data[END_REF], it was shown that these and other model parameters, notably θ M P , are also jointly distinguishable from the measured profiles µ P and M P . The result is based on the specialization of ( 8)-( 9) for the case of the fixed rate process, given by [

2]   żMP ż× żF   =   A M P A M P,× A M P,F 0 A × A ×,F 0 0 A F   •   z M P z × z F   +   0 0 u F   , (15)
where A × = -αI + A ⊗ , u F = λ + λ + T and A F depends only on λ + and α as detailed in Appendix A. For known parameters θ M P , we may easily show that λ -, λ + and µ F (0) are also distinguishable from the sole mean µ P . For simplicity, we consider the case where M and P are identically 0 at time 0. By inspection of ( 15), μF = -αµ

F + λ + , μM = -d M µ M + k M µ F , μP = -d P µ P + k P µ M (16)
(the expression of μF above coincides with the differential form of ( 12)). Thus, in terms of Laplace transform,

µ P (s) = λ + k M k P s(α + s)(d M + s)(d P + s) + µ F (0)k M k P (α + s)(d M + s)(d P + s) ,
and one may apply the method of [START_REF] Cinquemani | Reconstruction of promoter activity statistics from reporter protein population snapshot data[END_REF] (also reported in Appendix B) to prove sensitivity of this solution (equivalently, the solution over time) to any change in the three unknown parameters, almost everywhere in the space of the remaining parameters. In practical terms, parameter values can be reconstructed either from µ F as obtained by deconvolution from µ P , or by direct fit of (16) to an observed µ P profile. Now assume that F has reached stationarity. In this case, all relevant statistics of F are determined by λ + and λ -. However, from Proposition 1, mean µ F only conveys information about the ratio λ + /α, and, because σ 2 F = µ F (1 -µ F ) at any point in time, no more information is contained in the variance. Specific contributions of the two parameters can instead be traced in the autocovariance function ρ F . Indeed, multiplicative factor λ + (α -λ + )/α 2 and decay rate α have distinguishable effects on ρ F (different choices of the two lead to different profiles ρ F (•)) and uniquely determine λ + and λ -. The question arises whether ρ F is observable from the measured profiles µ P and M P (i.e. whether λ + and λ -are also distinguishable from the experimental output). In this section we provide a positive answer in terms of identifiability of λ + and λ -, i.e. for processes F with fixed rates. A more general answer will be provided in the next section.

From Proposition 1, stationary conditions are achieved when µ F is in steady state (i.e. when the factors of ρ F (t, τ ) involving µ F no longer depend on τ ). It then suffices to check identifiability of λ + and λ -from the solution of ( 15) with stationary initial conditions µ F (0) = λ + /α and σ 2 F (0) = λ + /α(1 -λ + /α). Using again the method of [START_REF] Cinquemani | Reconstruction of promoter activity statistics from reporter protein population snapshot data[END_REF], one computes the Laplace response function of this system. The resulting equations are lengthy and not reported here. Then, it can be checked that the Laplace sensitivity condition also reported in Appendix B is verified, i.e. the time profiles of µ P and M P are sensitive to all possible changes of λ + and α, almost everywhere in the space of the parameters θ M P . Statistics for F (dashed lines and circles) and F (solid lines and dots). Lines visualize analytic solutions, markers are for empirical statistics from Gillespie simulations. Gillespie simulations are performed using Stochkit [START_REF] Sanft | Stochkit2: Software for discrete stochastic simulation of biochemical systems with events[END_REF] for the generation of 10 4 sample paths (i.e. simulated cells). Numerical calculations are performed in Matlab.

processes, F and F , are considered. F has λ + = λ -= 0.05, while F has the faster switching dynamics λ + = λ -= 0.5. Starting from the non-stationary conditions F = F = 0 at time 0, means µ F and µ F converge both at 0.5 at different rates (Figure 1, left), thus resulting into different output profiles µ P (not shown). In other words, the two processes are distinguishable from the mean. In stationary conditions, instead, the means for F and F are the same. Yet the stationary autocovariance functions ρ F and ρ F differ in the two cases (Figure 1, center). This results in different observed profiles of σ P P (Figure 1, right). In other words, in stationary conditions, F and F are distinguishable from the output variance.

Remark 1. Equations ( 15) are obtained from ( 8)-( 9) by developing the expression of z ⊗ . This results in expressions depending on matrices A × and A F , which bring in the role of α, the decay rate of ρ F , into the propagation of second-order moments from z F to z M P . This fact is indeed in agreement with the discussion of z ⊗ at the end of Section 2.2.

To summarize, we have shown that constant switching rates, whence all statistics, of a promoter activity process F can be reconstructed from the output mean if F is not in stationary conditions. In stationary conditions, the promoter statistics cannot be determined from the output mean, but rather from the output variance since this reflects differences in the autocovariance function of F . Analytic expressions and a case study have been developed to support our arguments.

General promoter switching processes

We now wish to study how first-and second-order moments of switching process F reflect into outputs µ P and M P , and how to possibly reconstruct the former from the latter, without a priori knowledge on F . In particular, we do not assume that switching rates λ + and λ -are fixed. We only assume that F has continuous (mean and) autocovariance ρ F (t, s), and that Assumption 1 holds. For simplicity, we focus on the case where z M P (0) is null (M and P equal to zero at time zero).

From Equation [START_REF] Koopmans | The Spectral Analysis of Time Series[END_REF], for some final time T > 0, the conditional moments µ f P (t) and σ f P P (t) over [0, T ) are the output of a linear dynamical system with (zero initial conditions and) input f . We may then introduce linear operators, L 1 and L 2 , and abstract the transformation from function f to µ f P and σ f P P as µ f P = L 1 f and σ f P P = L 2 f . When necessary, for any t ∈ [0, T ), we will write

µ f P (t) = (L 1 f )(t) as L t 1 f and σ f P P (t) = (L 2 f )(t) as L t 2 f . Of course, for k = 1 and k = 2, L t k f = t 0 dτ k (t, τ )f (τ ), k (t, τ ) = C k e A M P (t-τ ) (A M P,F ) 1 ,
with C 1 = 0 1 0 0 0 (mean readout) and C 2 = 0 0 0 1 0 (variance readout).

Observability and reconstruction of the process mean

From the first equality in [START_REF] Kaern | Stochasticity in gene expression: From theories to phenotypes[END_REF], one has that

µ P = (L 1 f )dP F (f ) = L 1 f dP F (f ) = L 1 µ F .
Not surprisingly at this point, µ P thus follows from the linear dynamical transformation of µ F already found in [START_REF] Hespanha | Modelling and analysis of stochastic hybrid systems[END_REF]. Observability of µ F from µ P essentially depends on the spectrum of L 1 . Since

µ P (s) = k M k P (d M + s)(d P + s) µ F (s),
for strictly positive parameters θ M P , the transformation is invertible over the whole spectrum, i.e µ F can be perfectly reconstructed from µ P . In practice, this amounts to a deconvolution problem of rather easy solution [START_REF] Cinquemani | Reconstruction of promoter activity statistics from reporter protein population snapshot data[END_REF].

Observability and reconstruction of the process covariance

We begin with the following result.

Proposition 2. For any time t ∈ [0, T ), it holds that

M P (t) = L t 2 µ F + E[(L t 1 F ) 2 ]. ( 17
)
Clearly the autocovariance function of F plays a role in the term E[(L t 1 F ) 2 ]. To study this term further, consider the Karhunen-Loève decomposition [START_REF] Papoulis | Probability, random variables, and stochastic processes[END_REF] of process F , given by

F -µ F = ∞ i=1 a i φ i ,
where the φ i are the mutually orthogonal, unit norm eigenfunctions of the operator K : φ → dτ ρ F (•, τ )φ(τ ), i.e. Kφ i = σ 2 i φ i , and the a i are mutually uncorrelated, zero-mean random variables with variance equal to the eigenvalues σ 2 i (function norm is in L 2 and the decomposition holds in the mean-square sense). Then

ρ F (t, τ ) = ∞ i=1 σ 2 i φ i (t)φ i (τ ), σ F F (t) = ∞ i=1 σ 2 i φ 2 i (t). Proposition 3. It holds that E[(L t 1 F ) 2 ] = (L t 1 µ F ) 2 + ∞ i=1 σ 2 i (L t 1 φ i ) 2 .
In sums, from Propositions 2-3 and using the fact that (µ

P ) 2 = (L 1 µ F ) 2 , σ P P (t) = M P (t) -µ 2 P (t) = L t 2 µ F + ∞ i=1 σ 2 i (L t 1 φ i ) 2 . ( 18
)
Comparing the expressions of σ P P and σ F F one notices that, besides term L t 2 µ F , the functions composing F and characterizing its autocovariance structure are transformed by L t 1 into contributions that make up the variance of P at time t. Were L t 1 an evaluation operator, i.e. L t 1 φ i = φ i (t), then σ P P (t) would degenerate to L t 2 µ F +σ F F (t), i.e. information about the autocovariance structure of F would be lost. For every t, it is the integral nature of L t 1 that channels information about the whole ρ F (•, •) into σ P P (t). Another viewpoint on this is given in what follows.

Equation [START_REF] Stefan | Inference of quantitative models of bacterial promoters from time-series reporter gene data[END_REF] explains the nature of the information transfer from ρ F to σ P P . For reconstruction purposes, however, we seek a more explicit relationship between σ P P and ρ F . The following result relies on the convolutional form of L 1 .

Proposition 4. It holds that

σ P P (t) = L t 2 µ F + dτ dv 1 (t, τ ) 1 (t, v)ρ F (τ, v). (19)
Hence ρ F undergoes itself a linear transformation H defined by

H t ρ = (Hρ)(t) = t 0 dτ t 0 dv 1 (t, τ ) 1 (t, v)ρ(τ, v).
In particular, suppose that F is stationary. Then, by a change of variables,

H t ρ F = t 0 dτ t 0 dv 1 (t, τ ) 1 (t, v)ρ F (τ -v) = t -t dδ h(t, δ)ρ F (δ) with h(t, δ) = min{t,t-δ} max{-δ,0} dv 1 (t, v + δ) 1 (t, v).
In the light of these results, the problem of the observability of ρ F , or better the joint observability of ρ F and µ F from µ P and σ P P , is thus equivalent to that of the invertibility of the linear operator

(µ F , ρ F ) → (L 1 µ F , L 2 µ F + Hρ F ) (20)
(with relevant simplifications if stationarity of F is hypothesized). We note that, besides term L 2 µ F , the relationship between ρ F and σ P P is analogous to that pertaining linear transformations of second-order processes. In particular, using the fact that 1 (t, •) is the impulse response of a time-invariant dynamical system, the second term of ( 19) can be seen as the autocovariance of the output of a linear filter with response 1 fed with an input process with autocovariance ρ F . It is then natural to frame observability analysis of ρ F in the context of spectral analysis [START_REF] Koopmans | The Spectral Analysis of Time Series[END_REF][START_REF] Lindquist | Linear stochastic systems -A geometric approach to modeling, estimation and identification[END_REF]. This analysis is left for future work. Here we limit ourselves to the discussion of an illustrative example.

Example 2. Refer to Figure 2. We consider a promoter activity process F with randomly regulated rates, and compare its statistics with those of relevant fixedrate processes F and F . All processes are analyzed in stationary regime and have rate λ -identically set to 0.5, i.e. their definition only differs in the activation rate. The activation rate of

F is λ + (R) = 1 • R.
Regulator R is another random binary process with switch-off rate equal to 0.1 and switch-on rate, equal to 0.2217, chosen so as to guarantee that the stationary mean of F is µ F = 0.5. Process F is defined as in Example 1, i.e. it has λ + = 0.5, again resulting in µ F = 0.5. Finally, process F has activation rate λ + = E[λ + (R)] = 0.6892, i.e. Statistics for a random-rate promoter process F (dash-dotted lines) and relevant fixed-rate promoter F (same as in Figure 1, dashed lines) and F (dotted lines). Left: autocovariance functions ρF (dots: estimates from Gillespie simulations; line: interpolation), ρ F (from ( 14)) and ρ F (from ( 14)); Center: Observed output mean µP for F (Gillespie simulation), F (solution of ( 15)) and F (solution of ( 15)) -curves for F and F are superimposed; Right: The observed output variance of P for F (diamonds: numerical computation of [START_REF] Zechner | Moment-based inference predicts bimodality in transient gene expression[END_REF], based on the profile of ρF from Gillespie simulations; line: estimate from Gillespie simulation -diamonds and line are superimposed), F (solution of ( 15)) and F (solution of ( 15)) -curves for F and F are superimposed. Gillespie simulations are performed using Stochkit [START_REF] Sanft | Stochkit2: Software for discrete stochastic simulation of biochemical systems with events[END_REF] for the generation of 10 5 sample paths (i.e. simulated cells). Numerical calculations are performed in Matlab.

a switch-on rate equivalent on average to that of F . This results in a different mean, µ F = 0.5795.

The autocovariance function of F (Figure 2, left) is markedly different from those of F and F , which are similar. Because µ F = µ F = µ F , F can be distinguished from F , but not from F , from the output mean µ P (Figure 2, center). However, because of the different autocovariance function, F can be distinguished from F from the output variance σ P P . Interestingly, the output variance profiles for F and F are quite similar, a sign that the differences between F and F in mean and autocovariance compensate each other in this case. This is possible since output variance depends not only on the autocovariance but also on the mean of the promoter activity process.

Finally, in the light of the linearity of [START_REF] Zechner | Scalable inference of heterogeneous reaction kinetics from pooled single-cell recordings[END_REF], joint estimation of µ F and ρ F from possibly noisy and sampled measurements of µ P and σ P P can be seen as a linear inversion problem. Regularized solutions for both the stationary and the nonstationary case may be developed in accordance with the vast literature on the subject (see e.g. [START_REF] De Nicolao | Nonparametric input estimation in physiological systems: Problems, methods, and case studies[END_REF] and references therein). Note that, because µ F can be reconstructed from the sole mean µ P , the problem may also be reduced to that of estimating ρ F from σ P P -L t 2 µ F via (regularized) inversion of H. In summary, we have analyzed the relationship between second-order promoter activity statistics and mean and variance profiles of the reporter protein P in the case of promoter processes with randomly regulated rates. In particular, we have developed explicit relationships between the autocovariance function of F and the readout variance profile of P , showing that this integral relationship is essentially linear. By this we provided the basis for a full spectral analysis of observability of ρ F and its linear reconstruction from reporter protein mean and variance statistics. We also illustrated the relevance of our results by investigating the distinguishability of a random-rate and relevant fixed-rate processes on an example.

Conclusions

We have studied the relationships between second-order statistics of random promoter activity and the mean and variance profiles of gene expression reporter proteins typically observed in biological experiments. For both fixed and randomly regulated (thus also possibly time-varying) switching rates, we developed explicit mathematical formulas showing that these relationships are linear, and provided first results about the observability of the promoter process statistics from gene reporter data. Based on analytic considerations as well as on example case studies, we showed when and how analysis of second-order moments allows for discrimination of different promoter activation statistics.

This work provides the basis for an extensive observability analysis of promoter processes from gene reporter data at a single-cell level, and the development of promoter statistics reconstruction algorithms that are fully nonparametric, i.e. independent of a priori knowledge about the promoter activity laws. Our results show that both observability and estimation can be framed in the well-studied context of linear operators. Subsequent research work will henceforth focus on the application of the relevant spectral analysis and regularized linear inversion techniques. On these bases, we will then address a key challenge of this research effort, namely the identification and discrimination among alternative promoter activity regulatory mechanisms on the basis of the reconstructed promoter activation statistics and data from candidate regulators.

A Definitions and proofs

Matrix definitions. A

M P A M P,× A M P,F are given by       -d M 0 0 0 0 k P -d P 0 0 0 d M 0 -2 d M 0 0 k P d P 0 -2 d P 2 k P 0 0 k P 0 -d M -d P       ,       0 0 0 0 2 k M 0 0 0 0 k M       ,       k M 0 0 0 k M 0 0 0 0 0      
, in the same order, while

A ⊗ = -d M 0 k P -d P , A ×,F = 0 k M 0 0 , A F = -α 0 α -2λ + -2α .
Proof of Proposition 1. Process F is a homogeneous continuous-time binary Markov chain. Letting p(t) = Prob{F (t) = 0} Prob{F (t) = 1} T , for any t and τ it holds that

p(t) = e Q(t-τ ) p(τ ), Q = -λ + λ - λ + -λ - .
Mean µ F = Prob{F (t) = 1}. Using the fact that ṗ = Qp, the differential equation for µ F , the second element of p,

is μF = λ + (1 -µ F ) -λ -µ F = -αµ F + λ + .
The solution of this equation relative to µ F (0) yields the expression in the statement. Covariance ρ

F (t, τ ) = Prob{F (t) = 1, F (τ ) = 1} -µ F (t)µ F (τ ). By Bayes'law, Prob{F (t) = 1, F (τ ) = 1} = Prob{F (t) = 1|F (τ ) = 1}•Prob{F (τ ) = 1}.
Second factor is equal to µ F (τ ), while the first factor is given by the entry of row 2 and column 1 of e Q(t-τ ) . Computing the matrix exponential thus yields the result. Stationary versions of µ F and ρ F are found simply by taking the limit of µ F (t) as t → +∞ and replacing the result for µ F (τ ) and µ F (t) in the expression of ρ F (t, τ ).

Proof of Proposition 2. Starting from the second relation in [START_REF] Kaern | Stochasticity in gene expression: From theories to phenotypes[END_REF], where the last row vanishes since E P -E[P |F ]|F = 0. Then, using the definitions of µ F P and σ F P P , the chain of equalities continues with = E σ F P P (t) + E µ F P (t)

M P (t) = E E[P 2 |F ] = E E (P -E[P |F ]) + E[P |F ] 2 |F = = E E (P -E[P |F ]) 2 |F + E E E[P |F ] 2 |F + + 2 • E E (P -E[P |F
2 = E L t 2 F + E (L t 1 F ) 2 ] = L t 2 µ F + E[(L t 1 F ) 2 ].
Proof of Proposition 3. The following chain of inequalities hold:

E[(L t 1 F ) 2 ] = E L t 1 i a i φ i 2 = E i,j L t 1 (a i φ i )L t 1 (a j φ j ) = = i,j E[a i a j ](L t 1 φ i )(L t 1 φ j ) = i σ 2 i (L t 1 φ i ) 2 , ( 21
)
where the latter equality follows from the mutual uncorrelation of the a i .

Proof of Proposition 4. Expanding the last term of (17) one gets

E[(L t 1 F ) 2 ] = dP F (f )(L t 1 f ) 2 = dP F (f ) t 0 dτ 1 (t, τ )f (τ ) t 0 dv 1 (t, v)f (v) = t 0 dτ t 0 dv 1 (t, τ ) 1 (t, v) dP F (f )f (τ )f (v) = t 0 dτ t 0 dv 1 (t, τ ) 1 (t, v) ρ F (τ, v) + µ F (τ )µ F (v)
where the last integrand is of course the autocorrelation of F at τ and v. Therefore

σ P P (t) = M P (t) -µ 2 P (t) = L t 2 µ F + t 0 dτ t 0 dv 1 (t, τ ) 1 (t, v) ρ F (τ, v) + µ F (τ )µ F (v) - t 0 dτ 1 (t, τ )µ F (τ ) t 0 dv 1 (t, v)µ F (v) ,
and the result follows by collecting integrals and simplifying.

B Laplace sensitivity method for the analysis of parameter identifiability

This section reports the identifiability analysis method of [START_REF] Cinquemani | Reconstruction of promoter activity statistics from reporter protein population snapshot data[END_REF]. Let Y θ (t) be a vector function of t ∈ R depending on parameters θ. Typically Y θ (•) is an observed response of a dynamical system defined in terms of θ. In the present paper, the Laplace transforms that are used to discuss identifiability belong to this last class (see [START_REF] Cinquemani | Reconstruction of promoter activity statistics from reporter protein population snapshot data[END_REF]), whence Corollary 1 applies. In practice, these conditions can be easily checked by the use of the Matlab Symbolic Math Toolbox and evaluation of the rank conditions based on a finite set of heuristically chosen points S L (see again [START_REF] Cinquemani | Reconstruction of promoter activity statistics from reporter protein population snapshot data[END_REF]).

Example 1 .

 1 Fig.1. Statistics for F (dashed lines and circles) and F (solid lines and dots). Lines visualize analytic solutions, markers are for empirical statistics from Gillespie simulations. Gillespie simulations are performed using Stochkit[START_REF] Sanft | Stochkit2: Software for discrete stochastic simulation of biochemical systems with events[END_REF] for the generation of 10 4 sample paths (i.e. simulated cells). Numerical calculations are performed in Matlab.

  Fig.2. Statistics for a random-rate promoter process F (dash-dotted lines) and relevant fixed-rate promoter F (same as in Figure1, dashed lines) and F (dotted lines). Left: autocovariance functions ρF (dots: estimates from Gillespie simulations; line: interpolation), ρ F (from (14)) and ρ F (from (14)); Center: Observed output mean µP for F (Gillespie simulation), F (solution of (15)) and F (solution of (15)) -curves for F and F are superimposed; Right: The observed output variance of P for F (diamonds: numerical computation of (19), based on the profile of ρF from Gillespie simulations; line: estimate from Gillespie simulation -diamonds and line are superimposed), F (solution of (15)) and F (solution of (15)) -curves for F and F are superimposed. Gillespie simulations are performed using Stochkit[START_REF] Sanft | Stochkit2: Software for discrete stochastic simulation of biochemical systems with events[END_REF] for the generation of 10 5 sample paths (i.e. simulated cells). Numerical calculations are performed in Matlab.

  ]) • E[P |F ]|F = = E E (P -E[P |F ]) 2 |F + E E[P |F ] 2 + + 2 • E E P -E[P |F ]|F • E[P |F ] ,

Definition 1 .Proposition 5 .Corollary 1 .

 151 The parametric family (of functions) {Y θ : θ ∈ Θ}, with Θ ⊆ R N , N ∈ N, is (a) locally identifiable at θ * if a neighborhood B θ * ⊆ Θ of θ * exists such that the implication holds ∀θ ∈ B θ * ; (b) locally identifiable if (a) holds for almost every (a.e.) θ * ∈ Θ.For any givenθ let Y (s, θ) be the Laplace transform of Y θ (•). Let ∇Y (s, θ) = ∂Y ∂θ (s, θ) = ∂Y ∂θ1 • • • ∂Y ∂θ N (s, θ). If, for some L ∈ N, a set of points S L = {s 1 , . . . , s L } ⊆ R (or C) exists such that the matrix ∆(S L , θ * ) = ∇Y (s 1 , θ * ) T • • • ∇Y (s L , θ * ) T Thas full column rank, then {Y θ : θ ∈ Θ} is locally identifiable at θ * (in the sense of Definition 1(a)). Now assume that the elements of Y (s, θ) are ratios of polynomials in the entries of θ. If, for a given set of points S L and a given θ * , matrix ∆(S L , θ * ) is full column rank, then {Y θ : θ ∈ Θ} is locally identifiable (a.e. in the sense of Definition 1(b)).