Parameter inference with analytical propagators for stochastic models of autoregulated gene expression

1 Introduction and background

Gene expression in regulatory networks is an inherently stochastic process [1]. Mathematical models typically take the form of a chemical master equation (CME), which describes the temporal evolution of the probabilities of observing specific states in the network [2]. Recent advances in single-cell fluorescence microscopy [3–7] have allowed for the generation of experimental longitudinal data, whereby the fluorescence intensity of mRNA or protein abundances in single cells is measured. While most common models assume the availability of protein abundance data, abundances of mRNA may equally be of interest, depending on the model [8]. Here, we focus on abundances of protein, which we assume to be measured at regular sampling intervals Δt. A typical data set, denoted by Q, thus consists of protein abundances n _i at N + 1 different points in time; see Figure 1A. We can group these abundances into transitions n _i → n _i+1; cf. Figure 1B. A model-derived propagator P n i + 1 | n i ( Δ t , Θ ) allows for the calculation of the probabilities of such transitions for some set of model parameters Θ. Summing over all these probabilities for all N observed transitions, we can calculate the log-likelihood L(Θ) of that particular parameter set as

which can be evaluated over a range of values for the model parameters to yield a ‘log-likelihood landscape’, the maximum of which corresponds to the most likely parameter set Θ subject to the measured data set Q. Assessing the log-likelihood of a model by evaluating the associated propagator for the observed transitions in a time-lapse experiment is a feasible, established approach that has been successfully applied previously [6], [9], [10]. Due to the complex nature of the underlying regulatory networks, explicit expressions for P n i + 1 | n i are difficult to obtain in general. Hence, a variety of approximations have been proposed, which can be either numerical [9], [11] or analytical [12], [13], to cite but a few examples. Here, we apply the analytical method recently developed by the current authors [14], which was based on ideas presented by Popović et al. [12], to obtain fully time-dependent approximate propagators; an outline of the method is given in Section 2.

Figure 1:

(A) Simulated time series of protein abundance n, with measurements at times t _i and sampling interval Δt. (B) Histogram of the frequency k _j of transitions ( n 0 → n ) j , inferred from a longer time series with 100 transitions.

Our aim in the present article is to demonstrate the applicability of these propagators, as well as to evaluate their performance in the context of parameter inference for synthetic data. Specifically, we showcase the resulting inference procedure for a family of stochastic gene expression models. First, in Section 3, we consider a model that incorporates DNA ‘on’/‘off’ states (‘telegraph model’); see also the work of Raj et al. [15] and Shahrezaei and Swain [10]. Subsequently, in Section 4, that model is extended with an autoregulatory mechanism, whereby protein influences its own production through an autocatalytic reaction. In Section 5, we summarize our results and present an outlook to future research; finally, in Appendix A, we collate the analytical formulae that underly our inference procedure for the family of models showcased here.

2 Method

3 Showcase 1: the telegraph model

To demonstrate our parameter inference procedure, we consider a stochastic gene expression model that incorporates DNA ‘on’/‘off’ states (‘telegraph model’) [10], [15]:

In recent work [14], we presented an analytical method for obtaining explicit, general, time-dependent expressions for the generating function associated to the CME that arises from the model in (3.1). A pivotal element of the application of that method to (3.1) is the assumption that the protein decay rate d ₁ is notably smaller than the decay rate d ₀ of mRNA, which implies that the parameter ε ≔ d 1 d 0 is small; hence, the associated generating function is approximated to a certain order O = k, corresponding with a theoretical accuracy that is proportional to ɛ ^k . For more details on the resulting approximation, we refer to Appendix A.

To obtain synthetic data, we simulate the model in (3.1) using Gillespie’s stochastic simulation algorithm (SSA) [17], for fixed values of the (rescaled) parameters

on the time interval 0 ≤ t ≤ 10, and we measure the protein abundance n with a fixed sampling interval Δt. As our method assumes that Δt is of order ɛ, cf. again Appendix A, we set Δt = ɛ = 0.1, which yields N = 100 transitions. Finally, we consider random initial states, with mRNA and protein numbers chosen uniformly between 0 and 10, and we assume DNA to be in the ‘on’ state with probability κ 0 κ 0 + κ 1 . Based on the simulated measurement data, we perform the parameter inference procedure described in Section 2. As the data consists of protein abundances only, and as propagators for the model in (3.1) depend on abundances of both mRNA and protein, we marginalize over mRNA, assuming a steady-state distribution reported by Raj et al. [15, Supporting Information, Protocol S1, Eq. (1)]; that distribution coincides with the steady-state limit of the associated fully time-dependent distribution considered by Veerman et al. [14]. We assume that the values of κ ₀, κ ₁, ɛ, and d ₁ are known, and calculate the log-likelihood in (2.7) for varying λ and μ. We scan these two parameters in the range 1 0 − 3 ≤ λ ≤ 1 0 3 , 1 0 − 3 ≤ μ ≤ 1 0 2 , using a logarithmically spaced grid of 50 × 40 grid points. Figure 2 shows the resulting log-likelihood landscapes and, in particular, a comparison of the performance of the leading (zeroth) order approximation for the generating function, see Figure 2A, with that of the first order approximation in Figure 2B. We emphasize that our choice of λ and μ as the parameters to be inferred is not guided in any way by our analytical method – indeed, any other choice would have served our purpose, under the assumption that ɛ is sufficiently small. Rather, we chose λ and μ for illustrative purposes.

Figure 2:

Log-likelihood landscapes inferred from a simulation of the telegraph model in (3.1) with N = 100 transitions and parameter values as in (3.2): true value (cross) versus maximum log-likelihood estimate (MLE; dot). (A) Leading (zeroth) order approximation. (B) First order approximation.

Moreover, it is important to note that, in the analytical derivation of the propagators used to produce Figure 2, all model parameters are assumed to be of order ɛ ⁰ = 1 [14, Assumption 3.2], which is reflected in the numerical values in (3.2). In particular, it follows that both the ratio of the transcription and translation rates ν 0 ν 1 = ε λ μ and the ratio of the transcription and mRNA decay rates ν 0 d 0 = ε λ are assumed to be small. However, in Figure 2, the ranges over which λ and μ are scanned significantly exceed these estimates without causing apparent issues for our inference procedure. That could be taken as a sign that, while our analytical propagators were derived under certain assumptions on the order of model parameters, the resulting expressions may in fact be valid over a wider parameter range. On the other hand, it is worthwhile to note that noise-induced bistability is observed in a similar gene expression model if the above assumption on λ is abandoned, as considered with an alternative analytical approach in [18].

To quantify the performance of the method developed by Veerman et al. [14] for parameter inference, we compare four different scenarios:

1. Parameter values as in (3.2), with sampling interval Δt = ɛ = 0.1 on the time interval 0 ≤ t ≤ 10, corresponding to N = 100 transitions, which is the original setup that yields the results shown in Figure 2.

2. As in (a), with the time interval increased to 0 ≤ t ≤ 100, which yields N = 1000 transitions.

3. As in (a), with ɛ = 0.01: the sampling interval is decreased accordingly to Δt = ɛ = 0.01; measurements are taken on the time interval 0 ≤ t ≤ 1, which yields N = 100 transitions.

4. As in (a), with μ = 28.5.

For each scenario, we infer the most likely values of the parameters λ and μ, for increasing approximation order O. The inferred values of λ and μ are compared to the ‘true’ values λ _true and μ _true, where we consider relative errors to quantify the performance of our inference procedure. The results of that comparison are shown in Figure 3. The accuracy of inference for λ clearly increases when the approximation order O is increased from 0 to 1; the increase in accuracy from O = 1 to O = 2 is obfuscated by grid size effects. A ten-fold increase in the number of transitions (b) increases the accuracy of the leading order approximation, while a ten-fold increase in the value of μ _true (d) decreases the accuracy of the leading order approximation. For μ, there is no noticeable increase in accuracy with the approximation order O within the parameter grid used, which could indicate that the number of transitions remains the dominant limiting factor. However, the accuracy of inference for μ increases overall when N is increased (b), the small parameter ɛ is decreased (c), or the value of μ _true is increased (d). An alternative explanation for the apparent insensitivity of the error in μ to the approximation order can be found in the underlying mRNA marginalization procedure; however, one would expect that marginalization to negatively influence the sensitivity of transcription (λ), rather than that of translation (μ).

Figure 3:

Relative error Δ r ( x ) ≔ x − x true x true of the inferred parameters λ (A) and μ (B), for increasing approximation order O; for O = k, the propagators P n i + 1 | n i are approximated up to and including terms of order ɛ ^k . (a) N = 100 transitions, ɛ = 0.1, λ _true = 3.3, and μ _true = 2.85. (b) N = 1000 transitions, other parameters as in (a). (c) ɛ = 0.01, number of transitions N and other parameters as in (a). (d) μ _true = 28.5, number of transitions N and other parameters as in (a).

4 Showcase 2: an autoregulated telegraph model

We extend the telegraph model in (3.1) with an autoregulatory mechanism, where the DNA activation rates are influenced by the presence of protein. Autoregulation is modelled in a catalytic manner, via one of the two following reactions:

The above pair of autoregulation mechanisms was introduced by Hornos et al. [19], and implemented, e.g. by Iyer-Biswas and Jayaprakash [20]; see Section 5 for a discussion of the physical validity of these mechanisms. Propagators for the telegraph model incorporating both autoregulation through protein, as in (4.1), and autoregulation via mRNA were determined in [14]. The protein-autoregulated telegraph model is also discussed in [21], where it is called the ‘full model’; cf. [21, Figure 1] for a comparison and discussion of several reductions of that model.

To assess the performance of our parameter inference procedure, we fix the parameter values as in (3.2). Again, we marginalize over mRNA, assuming a steady-state distribution; see Section 3. Note that the degree to which that steady-state assumption is valid directly depends on the magnitude of the autoregulation rates, as protein levels influence levels of mRNA through DNA activation rates. We generate six data sets, as follows:

1. Simulate the model in (3.1) without autoregulation (‘null model’; a _P = 0 = r _P ) on the time interval 0 ≤ t ≤ 10, which yields N = 100 transitions.

2. As in (A), with the time interval increased to 0 ≤ t ≤ 100, which yields N = 1000 transitions.

3. Simulate the extended model {(3.1),(4.1a)} with autoactivation for a _P δ = 0.3 on the time interval 0 ≤ t ≤ 10, which yields N = 100 transitions.

4. As in (C), with the time interval increased to 0 ≤ t ≤ 100, which yields N = 1000 transitions.

5. Simulate the extended model {(3.1),(4.1b)} with autorepression for r _P δ = 0.3 on the time interval 0 ≤ t ≤ 10, which yields N = 100 transitions.

6. As in (E), with the time interval increased to 0 ≤ t ≤ 100, which yields N = 1000 transitions.

Every data set consists of 10 runs of equal length.

Generating functions for the autoregulated extension, by (4.1), of the telegraph model in (3.1) have been derived in the theoretical companion article [14] to the current work, under the assumption that the autoregulation rate a _P or r _P is small compared to the protein decay rate d ₁. That assumption implies that the ratios a P d 1 ≔ α P δ and r P d 1 ≔ ρ P δ are small.

Parameter inference now proceeds as follows. We fix a data set, and take a single run from that set. For that run, we determine the log-likelihood of the autoactivated model in {(3.1),(4.1a)}, varying 0 ≤ α _P δ ≤ 1; likewise, we determine the log-likelihood of the autorepressed model in {(3.1),(4.1b)}, varying 0 ≤ ρ _P δ ≤ 1. The log-likelihood L of the autoregulated extension is then compared with the log-likelihood L ₀ of the non-regulated model in (3.1); as before, N denotes the number of data points, where L is defined as in (2.7), that is, we simply evaluate the probability of the observed transitions given the model under consideration, after marginalization over mRNA. The log-likelihood difference L − L ₀, which is equal to the logarithm of the likelihood ratio, quantifies the evidence for that model. We repeat the above procedure for all 10 runs in the data set, and we determine the mean and standard deviation; the outcome is illustrated in Figure 4. We observe that 1000 transitions suffice to correctly refute autorepression in (B, D), and to correctly refute autoactivation in (F). In the case of 100 transitions, no conclusion can be drawn from (A) and (C), while (E) correctly refutes autoactivation; however, the small difference between L and L ₀, with |L − L ₀|< 1, indicates low significance.

Figure 4:

Parameter inference for the extended autoregulated model in {(3.1),(4.1)} on the basis of various types of synthetic data, where the performance of the model is quantified via the log-likelihood difference L − L ₀. On the vertical axis, L − L ₀ indicates the difference between the log-likelihood of the autoactivated or autorepressed model and that of the non-regulated model in (3.1); the higher the value of L − L ₀, the more likely the associated model is. In each panel, the solid curve indicates the mean values based on 10 model runs; dashed curves indicate the uncertainty (one standard deviation). On the horizontal axis, the strength of autoregulation is measured by α _P δ (increasing to the right) or ρ _P δ (increasing to the left). (A) Data generated from the null (non-regulated) model in (3.1), with N = 100 transitions. (C) Data generated from the model in (3.1), with autoactivation as in (4.1a), for α _P δ = 0.3 and N = 100 transitions. (E) Data generated from the model in (3.1), with autorepression as in (4.1b), for ρ _P δ = 0.3 and N = 100 transitions. (B, D, F) as (A, C, E), but with N = 1000 transitions; note that the vertical axis has a different scaling. All other model parameters were assumed to be known.

It is important to reiterate that the validity of our assumption that mRNA obeys a steady-state distribution, which was used in the marginalization step, is coupled to the magnitude of the autoregulation rates. While a steady-state assumption can hence be argued to be approximately valid for small δ, an increase in autoregulation strength increases the protein level feedback on DNA activation rates, thereby indirectly, but dynamically, influencing levels of mRNA. We propose that the ensuing breakdown of the steady-state assumption is reflected in the behaviour of the log-likelihood difference shown in Figure 4. The assumption can be argued to induce a bias as α _P δ, ρ _P δ → 1, accordingly skewing the log-likelihood difference, which could also explain the absence of a clear maximum in Figure 4D and F.

5 Discussion

In the present article, we showcase a parameter inference procedure that is based on a recently developed analytical method [14] which allows for the efficient numerical approximation of propagators via the Cauchy integral formula on the basis of asymptotic series for the underlying generating functions. The resulting hybrid analytical–numerical approach reduces the need for computationally expensive simulations; moreover, due to its perturbative nature, it is highly applicable over relatively short time scales, such as those occurring naturally in the calculation of the log-likelihood in (1.1).

We simulate protein expression with a simple stochastic simulation algorithm (SSA) [17]; the resulting protein abundances are in the lower range of experimentally measured values [22], [23], which does not, however, limit our proof of principle. In addition, in the presented application of our approach, we assume that some of the model parameters are known, that protein abundances can be measured at regular sampling intervals, and that there is no noise. We then present results for synthetic data in a family of models for stochastic gene expression from the literature under the central assumption that lifetimes of protein are significantly longer than those of mRNA, which introduces a small parameter ɛ and, hence, a separation of scales. An extensive discussion of the validity of our assumption that ɛ is small can be found elsewhere [10], [14]; see also Section 3. The assumed smallness of ɛ is crucial to the underlying analytical method, as introduced by Veerman et al. [14]. In addition, our approach is specifically attuned to the time variability of the expression process, in the sense that we assume the sampling interval Δt to be small, as well; cf. again Section 3. It is thus ideally suited to describing transient dynamics far from steady state, as is evident in the bursting behaviour seen in Figure 1A.

In Section 3, we discuss a simple (‘telegraph’) model for gene expression without autoregulation, showing that our approach can successfully infer relevant model parameters. Unlike in previous work [24], the underlying implementation avoids potential bias due to zero propagator values and large initial protein numbers through the use of ‘implicit’ series expansions in ɛ; see Appendix A for an in-depth argument. Note that, in the simple telegraph model, transcription factors are assumed to act as single molecules. Our analytical approach would have to be extended to incorporate transcription factor dimer regulation, e.g. via the recently developed linear mapping approximation [25].

In Section 4, we perform a model comparison in an autoregulated extension of the telegraph model. We consider three types of gene regulation: autoactivation, autorepression, and no regulation of DNA activity (null model). For each type, we simulate protein abundance data with 100 and 1000 transitions, respectively. Throughout, we find that 100 data points are insufficient to reject model hypotheses with our approach. With 1000 data points, however, we can successfully reject the non-regulated and the autorepressed model for simulated data from an autoactivated model, in which case we can even infer the correct order of the autoactivation parameter. For simulated autorepression, we can reject the model with autoactivation, but not the non-regulated model. Our approach fails to identify the correct model for data from a non-regulated model for 1000 transitions, where the autoactivated model is clearly, but wrongly, favoured. We believe that more research is needed into the sources of these discrepancies in dependence on both model parameters and the order of our approximation.

In both showcases, we observe a trade-off between the accuracy of inference versus the required computation time. Computation times seem to increase exponentially with the approximation order, at least for the setup realized in this article. For practical purposes, we hence propose an algorithm whereby the fastest, leading order approximation is used to obtain an initial estimate for the underlying model parameters; that estimate can then be improved by including higher order corrections, resulting in a much more computationally efficient procedure.

It is insightful to compare our results with other recent work on parameter inference in regulated gene expression models. In work by Feigelman et al. [9], three models for regulated gene expression with a slightly different structure compared to the models studied in the present article were simulated and inferred via a particle filtering-based inference procedure that employs genealogical information of dividing cells. Interestingly, positive and negative autoregulation could be successfully rejected there for data that was simulated from a no-feedback model. However, the no-feedback model could not be rejected for data originating from the corresponding models with positive or negative feedback. From that comparison with Feigelman et al. [9], we conclude that the structure of the data, the intensity of regulatory feedback, and the chosen inference procedure together will influence the extent of insight which can be obtained from an approach that is based on stochastic models of gene expression.

We emphasize that our analytical method is not restricted to the specific models showcased here. Our aim in the present article is to demonstrate the applicability of the method, and to investigate its performance, rather than to assess the biological validity of a given model. Minimally, our approach can be extended to recent, physiologically more relevant modifications of the telegraph model with autoregulation [19] by Grima et al. [26] and Congxin et al. [27]; another feasible alternative model can be obtained by introduction of a refractory state [28]. Ideally, we would like to test a variety of stochastic gene expression models against a given set of measurement data. Current computational approaches struggle to provide propagators for models with more than one regulated species, which can often only be approximated even in that simple scenario. The principal advantage of our hybrid approach is that propagators can be evaluated in a computationally efficient manner, via a combination of mathematical analysis and numerical integration [14]; other approaches rely either on the calculation of propagators based on direct numerical simulation of the underlying model – which is computationally demanding – or on the assumption that symbolic derivatives of the generating function are explicitly known, which only holds for specific models of relatively low complexity [12].

The input for our propagator-based approach is the abundance of the involved species, viz. of protein. Thus, we assume that absolute protein numbers are measured, which is in practice hampered by an unknown scaling factor between the observed fluorescence and the corresponding abundances, and by noise. While various suggestions for inferring that factor have been made [5], [29], [30], and while a linear scaling is customarily assumed [6], [11], an accurate experimental determination remains extremely challenging.

It is instructive to compare our propagator-based approach to alternative approaches to parameter inference in the context of stochastic gene expression, such as the linear noise approximation [31], a system size expansion with moment closure [32], or tensor-based methods for the corresponding Fokker–Planck equation [33]. Both the linear noise approximation and a system size expansion assume large system sizes and are only valid for large copy numbers, see also [31]; moreover, the impact of moment closure on the relevance of nonlinearities is potentially non-negligible [34]. The relatively recently developed tensor-based methods focus on steady state distributions for large system sizes; their validity away from the thermodynamic limit is as yet unclear, nor is it clear how results are influenced by the particular tensor decomposition method that is chosen. In contrast to these approaches, the propagator-based approach we have employed in this article not only performs well with low copy numbers, i.e. for small system sizes; our analytical propagators are also explicitly time-dependent, in contrast to the usual steady state assumption. Details can be found in our theoretical companion article [14].

Finally, the showcases presented in this article are based on synthetic data that was generated in silico; in the future, we plan to consider experimental data, such as can be found in the work by Suter et al. [6].