Time-dependent propagators for stochastic models of gene expression: an analytical method.

The inherent stochasticity of gene expression in the context of regulatory networks profoundly influences the dynamics of the involved species. Mathematically speaking, the propagators which describe the evolution of such networks in time are typically defined as solutions of the corresponding chemical master equation (CME). However, it is not possible in general to obtain exact solutions to the CME in closed form, which is due largely to its high dimensionality. In the present article, we propose an analytical method for the efficient approximation of these propagators. We illustrate our method on the basis of two categories of stochastic models for gene expression that have been discussed in the literature. The requisite procedure consists of three steps: a probability-generating function is introduced which transforms the CME into (a system of) partial differential equations (PDEs); application of the method of characteristics then yields (a system of) ordinary differential equations (ODEs) which can be solved using dynamical systems techniques, giving closed-form expressions for the generating function; finally, propagator probabilities can be reconstructed numerically from these expressions via the Cauchy integral formula. The resulting 'library' of propagators lends itself naturally to implementation in a Bayesian parameter inference scheme, and can be generalised systematically to related categories of stochastic models beyond the ones considered here.

Introduction

Motivation

Understanding the process of gene expression in the context of gene regulatory networks is indispensable for gaining insight into the fundamentals of numerous biological processes. However, gene expression can be highly stochastic in nature, both in prokaryotic and in eukaryotic organisms; see e.g. the work by Elowitz et al. (2002), Raj and Oudenaarden (2008), Shahrezaei and Swain (2008b), and the references therein. This inherent stochasticity has a profound influence on the dynamics of the involved species, in particular when their abundance is low. Therefore, gene expression is often appropriately described by stochastic models (Bressloff 2014; Karlebach and Shamir 2008; Thattai and Oudenaarden 2001; Wilkinson 2009). A schematic of the canonical model for gene expression is depicted in Fig. 1. Here, the processes of transcription, translation, and degradation are approximated by single rates.

Fig. 1

The canonical model of gene expression. Transcription of mRNA occurs with rate ; mRNA is translated to protein with rate . Both mRNA and protein decay, with rates and , respectively.

Figure courtesy of Shahrezaei and Swain (2008a) (Copyright (2008) National Academy of Sciences, U.S.A.)

To test the validity of such stochastic models, a comparison with experimental data needs to be performed. The development of experimental techniques, such as time-lapse fluorescence microscopy (Coutu and Schroeder 2013; Elowitz et al. 2002; Larson et al. 2011; Muzzey and Oudenaarden 2009; Raj et al. 2006; Young et al. 2011), allows for real-time tracking of gene expression dynamics in single cells, providing mRNA or protein abundance time series, such as those depicted in Fig. 2. To select between competing hypotheses on the underlying regulatory networks given measurement data, as well as to infer the values of the corresponding model parameters, we can apply Bayesian inference theory to calculate the likelihood of a given model, which is constructed as follows.

Fig. 2

Sketch of a potential time series of protein abundance n, sampled at times , with regular, fixed, sampling interval

The abundance of a protein, denoted by n, is sampled at times , see Fig. 2, yielding a list of measurement pairs from which transitions between states can be extracted; here, is the regular, fixed, sampling interval. Next, the underlying stochastic model with parameter set is used to calculate the probabilities of these transitions, which are denoted by . These so-called propagators give the probability of protein being present after time , given an initial protein abundance of . The log-likelihood of the parameter set , given the observed data D, is now defined in terms of the propagators asTo infer the values of parameters in the model, propagators are calculated for a wide range of parameter combinations, resulting in a ‘log-likelihood landscape’; the maximal value of the log-likelihood as a function of the model parameters yields the most likely parameter values, given the experimental data. An example realisation of the above procedure can e.g. be found in the work by Feigelman et al. (2015).

To calculate accurately the log-likelihood in (1.1), it is imperative that the values of the propagators can be extracted from the underlying stochastic model for any desired combination of parameters in . In particular, we need to be able to calculate the propagator as a function of time t for any initial protein number . Unfortunately, the highly complex nature of the stochastic models involved makes it very difficult to obtain explicit expressions for these probabilities. Some analytical progress can be made when a steady-state approximation is performed, i.e. when it is assumed that the system is allowed to evolve for a sufficiently long time, such that it converges to a time-independent state. However, the sampling interval used for obtaining experimental data, as seen in Fig. 2, is often short with respect to the protein life time. As that life time represents a natural time scale for the system dynamics, it follows that the evolution of the probabilities should be studied over short times, in contradiction with the steady-state or long-evolution-time approximations which have previously been employed to derive analytical results (Bokes et al. 2012a; Hornos et al. 2005; Iyer-Biswas and Jayaprakash 2014; Shahrezaei and Swain 2008a).

The complex nature of stochastic models for gene expression has led to the widespread use of stochastic simulation techniques, such as Gillespie’s algorithm (Gillespie 1977), with the aim of predicting values for the associated propagators from these models; see Feigelman et al. (2016) for recent work combining stochastic simulation with a particle filtering approach. However, these approaches can still be very time-consuming, due to the (relatively) high dimensionality of the model parameter space, combined with the fact that, for each combination of parameter values, the stochastic model has to be simulated sufficiently many times to yield a probability distribution that can be used to infer the corresponding propagator. For that reason, it is desirable to be able to obtain explicit expressions for the propagator directly in terms of the model parameters, if necessary in an appropriate approximation.

Analytical method

In the present article, we develop an analytical method for the efficient evaluation of time-dependent propagators in stochastic gene expression models, for arbitrary values of the model parameters. The results of our analysis can be implemented in a straightforward fashion in a Bayesian parameter inference framework, as outlined above.

To demonstrate our approach, we analyse two different stochastic models for gene expression. The first model, henceforth referred to as ‘model A’, is a model that incorporates autoregulation, where transcription and translation are approximated by single rates and protein can either stimulate or inhibit its own production by influencing the activity of DNA; see Fig. 3. That model was first studied by Iyer-Biswas and Jayaprakash (2014) via a steady-state approximation. The second model, henceforth referred to as ‘model B’, models both mRNA and protein explicitly and again incorporates DNA switching between an active and an inactive state; see Fig. 4. That model was first studied by Shahrezaei and Swain (2008a) in a long-evolution-time approximation.

Fig. 3

Schematic of model A, a gene expression model with autoregulation.

Base figure courtesy of Shahrezaei and Swain (2008a). (Copyright (2008) National Academy of Sciences, U.S.A.)

Fig. 4

Schematic of model B, a gene expression model that explicitly incorporates transcription (Shahrezaei and Swain 2008a).

Figure courtesy of Shahrezaei and Swain (2008a). (Copyright (2008) National Academy of Sciences, U.S.A.)

Both model A and model B are formulated in terms of the chemical master equation (CME), which is the accepted mathematical representation of stochastic gene expression in the context of the model categories considered here; cf. Iyer-Biswas and Jayaprakash (2014) and Shahrezaei and Swain (2008a), respectively. Mathematically speaking, the CME is an infinite-dimensional system of linear ordinary differential equations (ODEs) that describes the evolution in time of the probabilities of observing a specific state in the system, given some initial state. Numerous approaches have been suggested for the (approximate) solution of the CME; see e.g. Popović et al. (2016) and the references therein for details. Our method relies on a combination of various techniques from the theory of differential equations and dynamical systems; specifically, we perform three consecutive steps, as follows.We emphasise that the ‘characteristic system’ of ODEs which is obtained in step 2 is low-dimensional, in contrast to the underlying CME system, as well as that it exhibits additional structure, allowing for the derivation of a closed-form analytical approximation for the associated generating function.

CME system PDE system: We introduce a probability-generating function to convert the CME into a (system of) partial differential equations (PDEs).

PDE system ODE system: Applying the method of characteristics—combined, if necessary, with perturbation techniques—we transform the system of PDEs obtained in step 1 into a dynamical system, that is, a system of ODEs.

ODE system Explicit solution: Making use of either special functions (model A) or multiple-time-scale analysis (model B), we obtain explicit solutions to the dynamical system found in step 2.

To convert the results of the above procedure into solutions to the original stochastic model, the three steps involved in our analysis have to be reverted. To that end, we require the following three ingredients:The various steps in our analytical method, as indicated above, are represented in Fig. 5. It is important to mention that the implementation of Bayesian parameter inference, as outlined in Sect. 1.1, is not a topic for the present article; rather, the aim here is to describe our method, and to present analytical results which can readily be implemented in the context of parameter inference. The article hence realises the first stage of our research programme; the natural next stage, which is precisely that implementation, will be the subject of a follow-up article by the same authors.

Fig. 5

Schematic overview of the analytical method

Initial conditions are originally stated in terms of the CME, and first have to be reformulated in terms of the corresponding system of PDEs to ensure well-posedness; then, initial conditions can be extracted for the dynamical system that was obtained via the method of characteristics, reverting step 3.

To transform solutions to the characteristic system into solutions of the underlying PDE system, the associated ‘characteristic transformation’ has to be inverted, reverting step 2.

Lastly, solutions of the CME have to be extracted from solutions to the resulting PDE system, reverting step 1. Although the correspondence between the two sets of solutions is exact, theoretically speaking, the complexity of the expressions involved precludes the efficient analytical reconstruction of propagators from their generating functions. Therefore, we propose a novel hybrid analytical-numerical approach which relies on the Cauchy integral formula.

Outline

The present article is organised as follows. In Sect. 2, we apply the analytical method outlined in Sect. 1.2 to model A, the gene expression model with autoregulation. Here, we use a perturbative approach to incorporate the autoregulatory aspects of the model; the resulting dynamical system can be solved in terms of confluent hypergeometric functions, see §13 in NIST Digital Library of Mathematical Functions . In Sect. 3, the same method is applied to model B, the model that explicitly incorporates transcription. We also indicate how autoregulation can be added to that model, and how the resulting extended model can be analysed on the basis of our treatment of model A. The analysis carried out in Sects. 2 and 3 yields a ‘library’ of explicit asymptotic expressions for the probability-generating functions associated to the underlying stochastic models. To obtain quantifiable expressions for their propagators, we introduce a novel hybrid analytical–numerical approach in Sect. 4, which can be readily implemented in the Bayesian parameter inference framework that provided the motivation for our analysis; see Sect. 1.1. We conclude with a discussion of our results, and an outlook to future work, in Sect. 5.

Model A: gene expression with autoregulation

We first demonstrate our analytical method in the context of an autoregulatory stochastic gene expression model, as presented by Iyer-Biswas and Jayaprakash (2014); see also Fig. 3. In the original article (Iyer-Biswas and Jayaprakash 2014), a Poisson representation was used to obtain analytical descriptions for time-independent solutions to the model. For a visual guide to the upcoming analysis, the reader is referred to Fig. 5.

Stochastic model and CME

The basic stochastic model for gene expression is represented by the reaction schemeThe gene can hence switch between the inactive state D and the active state , with switching rates and , respectively. The active gene produces protein (P) with rate , while protein decays with rate .

The autoregulatory part of the model is implemented as either positive or negative feedback: In the case of autoactivation, viz. (2.2a), protein induces activation of the gene with activation rate a, thereby accelerating its own production; in the case of autorepression, viz. (2.2b), protein deactivates the active gene with repression rate r, impeding its own production.

The CME system that is associated to the reaction scheme in (2.1), with autoactivation as in (2.2a), is given by Here, () represents the probability of n protein being present at time t while the gene is either inactive (0) or active (1). The time variable is nondimensionalised by the protein decay rate ; other model parameters are scaled asAnalogously, the CME system for the case of autorepression, as defined in (2.2b), is given by

Remark 2.1

A priori, it is possible to incorporate both autoactivation and autorepression in a single model, by merging systems (2.3) and (2.5). However, since autoactivation and autorepression precisely counteract each other, a partial cancellation would ensue, resulting in effective activation or repression. It can hence be argued that the simultaneous inclusion of both effects would introduce superfluous terms and parameters, which could be considered as poor modelling practice. Therefore, we choose to model the two autoregulation mechanisms separately.

Generating function PDE

Rather than investigating the dynamics of (2.3) and (2.5) numerically, using stochastic simulation, we aim to employ an analytical approach. To that end, we define the probability-generating functions () as follows; see e.g. Gardiner (2009):In the case of autoactivation, the generating functions can be seen to satisfy if the coefficients in (2.6) obey the CME system (2.3); likewise, in the case of autorepression, (2.3) gives rise to Both (2.7) and (2.8) are systems of coupled, linear, first-order, hyperbolic partial differential equations (PDEs). Systems of this type are typically difficult to analyse; existing techniques only provide general results (Courant and Hilbert 1962; Taylor 2011).

To allow for an explicit analysis of systems (2.7) and (2.8), we make the following assumption:

Assumption 2.2

We assume that the autoactivation rate a in (2.2) is small in comparison with the other model parameters; specifically, we writewhere is sufficiently small. Likewise, we assume that the autorepression rate r is small in comparison with the other model parameters, writing

Previous work on the inclusion of autoregulatory effects in model selection by Feigelman et al. (2016) suggests that, in the context of Nanog expression in mouse embryonic stem cells, autoregulation rates are indeed small compared to other model parameters.

Based on Assumption 2.2, we can expand the generating functions () as power series in :Substitution of (2.11) into (2.7) yields analogously, we substitute (2.11) into (2.8) to find We observe that, in (2.12) and (2.13), the same leading-order differential operator acts on both and , which allows us to apply the method of characteristics to solve the equations for () simultaneously. In particular, we emphasise that, mathematically speaking, the resulting perturbation is regular in the perturbation parameter .

Dynamical systems analysis

In this section, we apply the method of characteristics to derive the ‘characteristic equations’ that are associated to the PDE systems (2.12) and (2.13), respectively; the former are systems of ODEs, which are naturally analysed in the language of dynamical systems.

Autoactivation

We first consider the case of autoactivation; to that end, we rewrite system (2.12) as where we have introduced the new variableThe differential operator in Eq. (2.14) gives rise to characteristics that obey the characteristic equationthese characteristics can thus be expressed asSince the partial differential operators in (2.14) transform into we arrive at the following system: Note that Eq. (2.19) is a recursive (nonhomogeneous) system of ordinary differential equations for (). Henceforth, we will therefore refer to (2.19) as such, while retaining the use of partial derivatives due to the presence of in the corresponding right-hand sides.

To solve system (2.19), we rewrite it as a second-order ODE for : we hence obtainwhich can be solved recursively to determine for any . To simplify (2.20), we introduce the variablewhich transforms the partial derivatives and intoEquation (2.20) hence readsUsing (2.19a), we can express the second component in terms of asAt leading order, i.e. for , (2.23) reduces tothe solutions of which can be expressed in terms of the confluent hypergeometric function , see §13 of NIST Digital Library of Mathematical Functions , to yieldUsing (2.24), we can determinewith as given in (2.26); for an explicit expression, see Eq. (A.1) in Appendix A.

The expression for in (2.26) allows us to determine the first-order correction : substituting in (2.23), we obtainNext, we apply the method of variation of constants to express the solution to (2.28) aswherewith as given in (2.26). Finally, we may again use (2.24) to determinewith and as given in (2.26) and (2.29), respectively.

At this point in the analysis, the constants and in (2.26) and the integration limits and in (2.29) remain undetermined. To fix these constants, and thereby determine a unique solution to (2.19), we have to prescribe appropriate initial conditions.

Autorepression

Given the analysis of autoactivation in the previous subsection, the case of autorepression can be analysed in an analogous manner. Employing the same characteristics as before, recall (2.17), we obtain from system (2.13); cf. Eq. (2.19). Next, we rewrite (2.32) as a second-order ODE for , using again the variable transformation in (2.21):which can be solved recursively to obtain for any . The first component can be expressed in terms of asto leading order, we thus obtainandThe corresponding equation for the first-order correction readswhich can be solved via the method of variation of constants to givehere,with as in (2.35). The first-order correction to the first component is hence given bywith and as in (2.35) and (2.38), respectively. As in the case of autoactivation, the constants and in (2.35) and the integration limits and in (2.38) remain undetermined, and have to be fixed through suitable initial conditions.

Initial conditions

To determine appropriate initial conditions for the dynamical systems (2.19) and (2.32), we consider the original CME systems (2.3) and (2.5), respectively.

At time , we impose an initial protein number , which impliesfor the probabilities (); here, denotes the standard Kronecker symbol, with for and otherwise. Using the definition of the generating functions in (2.6), we find thattaking into account the change of variables in (2.15). Thus, (2.42) provides an initial condition for the PDE systems (2.7) and (2.8). Given the power series expansion in (2.11), we infer that the coefficients () satisfy which holds for both (2.12) and (2.13).

To be able to interpret the initial conditions in (2.43) in the context of the dynamical systems (2.19) and (2.32), we revisit the method of characteristics, which was used to map the PDE systems (2.12) and (2.13) to the former, respectively.

The characteristics of the differential operator in (2.19) and (2.32) are the integral curves of the vector field (1, v). Geometrically speaking, these characteristic curves foliate the (t, v)-plane over the v-axis. Therefore, each characteristic can be identified by its base point, which is the point where the characteristic curve intersects the v-axis, at ; see Eq. (2.17) and Fig. 6.

Fig. 6

The (t, v)-coordinate plane, on which the PDE systems (2.12) and (2.13) are solved. The characteristics are integral curves of the vector field (1, v), indicated in blue. The black characteristic curve intersects the v-axis at (colour figure online)

Equivalently, each characteristic can be represented as a graph over the t-axis. Indeed, the differential equation for the v-component of a characteristic curve (2.16) can be solved to obtain the natural parametric description of a characteristic ‘fibre’ with base point , which is given by . Here, the parameter s along the characteristic is chosen such that its point of intersection with the v-axis lies at . Given that choice, it is natural to identify the parameter along the characteristic (s) with the time variable (t). Hence, the initial conditions in (2.43), which determine a relation between and on the v-axis, can be interpreted on every characteristic as which again holds for both (2.19) and (2.32).

Remark 2.3

For CME systems such as (2.3) and (2.5), it is customary to impose a ‘normalisation’ condition of the formas and represent probabilities. Recast in the framework of generating functions, recall (2.6), the above normalisation condition yields the boundary conditionIt is worthwhile to note that (2.46) is automatically satisfied whenever (2.43) is imposed: by adding the two equations in system (2.12)—or, equivalently, in (2.13)—one sees that satisfiesThe line is represented by the ‘trivial’ characteristic, which is identified by ; therefore, on that characteristic, we impose (2.43) to find which implies that for all s, as well as that vanishes identically for all . Substituting these results into the power series representation for () in (2.11), we obtain (2.46).

At this point, it is important to observe that Eq. (2.43) determines a line of initial conditions in the phase spaces of the dynamical systems (2.19) and (2.32). Therefore, at every order in , we can only fix one of the two free parameters that arise in the solution of the corresponding differential equations. In particular, (2.43) only fixes either or in (2.26), and either or in (2.29), and so forth. That indeterminacy motivates us to introduce a new parameter , which is defined as follows:

Definition 2.4

Consider the CME systems (2.3) and (2.5). We define to be the probability that the gene modelled by the reaction scheme in (2.1) is switched off at time , given an initial protein number .

Definition 2.4 immediately specifies initial conditions for systems (2.3) and (2.5) viathe above expression, in turn, provides us with a complete set of initial conditions for the PDE systems (2.7) and (2.8), to wit Here, we allow for the fact that may depend on other model parameters, and in particular on the autoregulation rates a and r; for a discussion of alternative options, see Sect. 4.1. We therefore expand as a power series in ,The above expansion can be used to infer a complete set of initial conditions for the PDE systems (2.12) and (2.13), yielding By the same reasoning that inferred (2.44) from (2.43), we can conclude that the complete set of initial conditions for the dynamical systems (2.19) and (2.32) is given by

We can now use the conditions in (2.53) to determine the free constants and in (2.26), which yields Analogously, for and in (2.35), we obtainwith and as in (2.54). Using conversion formulas found in §13.3(i) of NIST Digital Library of Mathematical Functions , one can show that the expressions resulting from (2.35) and (2.36) match those in (2.27) and (2.26), respectively, as expected; see also Eq. (A.1). Explicit expressions for the integration limits and in (2.29), as well as for the corresponding limits and in (2.38), can be found in Appendix A.

Inverse transformation

The final step towards providing explicit solutions to the PDE systems (2.7) and (2.8) consists in interpreting the solutions to the dynamical systems (2.19) and (2.32), with initial conditions as in (2.53), as solutions to the PDE systems (2.12) and (2.13), respectively. To that end, we again consider the corresponding characteristics from a geometric viewpoint.

As mentioned in Sect. 2.4, the (t, v)-coordinate plane is foliated by characteristics, which are the integral curves of the vector field (1, v); recall Fig. 6. Hence, any point (t, v) lies on a unique characteristic; flowing backward along that characteristic to its intersection with the v-axis, we can determine the corresponding base point by the inverse transformationsince we have identified the parameter along the characteristic (s) with the time variable (t).

To determine the value of (), interpreted as a solution to the PDE system (2.12) or (2.13), we proceed as follows. We first apply the inverse transformation in (2.56) to establish on which characteristic the coordinate pair (t, v) lies. For that characteristic, identified by its base point , we then find the solution to the dynamical system (2.19) or (2.32), which is a function of s and . Next, we substitute and into that solution to obtain an explicit expression for the solution to the PDE system (2.12) or (2.13):where () on the right-hand side denotes the solution to (2.19) or (2.32), with initial conditions as in (2.53). Lastly, we substitute into the power series in (2.11) to obtain an explicit solution to (2.7) or (2.8), to satisfactory order in .

Remark 2.5

The geometric interpretation of characteristics that was used to motivate the inverse transformation in (2.56) also shows that the introduction of the new system parameter in Definition 2.4 is necessary for the generating functions to be determined uniquely as solutions to (2.7) or (2.8), even if we are only interested in their sum . The crucial point is that any free constants obtained in the process of solving the dynamical systems (2.19) and (2.32)—or, equivalently, their second-order one-dimensional counterparts (2.20) and (2.33), respectively—are constant in s. In other words, they are constant along the particular characteristic on which the dynamical system is solved. These constants—such as e.g.  and in (2.26)—can, and generally will, depend on the base point of the characteristic; see for example (2.54). The inverse transformation in (2.56) that is used to reconstruct the solution to the original PDE from that of the corresponding dynamical system would then yield undetermined functions in the resulting solutions to (2.7) and (2.8), respectively.

Summary of main result

To summarise Sect. 2, we combine the analysis of the previous subsections to state our main result.

Main result: The PDE system (2.7) can be solved for sufficiently small autoactivation rates a; see Assumption 2.2. Its solutions () are expressed as power series in the small parameter ; recall (2.11). The coefficients in these series, written in terms of the shifted variable v defined in (2.15), can be found byLikewise, we can solve the PDE system (2.8) for sufficiently small autorepression rates r; see Assumption 2.2. Its solutions are again expressed as power series in the small parameter ; cf. (2.11). The coefficients in these series, written in terms of the shifted variable v defined in (2.15), can be found byTo illustrate the procedure described above, we state the resulting explicit expressions for the leading-order solution to (2.7)—or, equivalently, to (2.8)—in terms of the original variables z and t: Note that similar expressions were derived by Iyer-Biswas et al. (2009), where an analogous generating function approach was applied under the assumption that the gene is initially inactive—i.e. that , see Definition 2.4—as well as that the initial protein number is zero. With these choices, the expression for from (2.58) can be seen to coincide with that found in Equations (5) through (7) of Iyer-Biswas et al. (2009), using §13.3.4 and §13.2.39 of NIST Digital Library of Mathematical Functions .

solving recursively the second-order ODE (2.20) and using the identity in (2.24), incorporating the initial conditions in (2.53);

and, subsequently, applying the inverse transformation in (2.56) to the resulting solutions.

solving recursively the second-order ODE (2.33) and using the identity in (2.34), incorporating the initial conditions in (2.53);

and, subsequently, applying the inverse transformation in (2.56) to the resulting solutions.

Model B: gene expression with explicit transcription

In this section, we apply our analytical method to model B, a stochastic gene expression model presented by Shahrezaei and Swain (2008a) which explicitly incorporates the transcription stage in the expression process, as well as DNA switching; see also Fig. 4. In the original article by Shahrezaei and Swain (2008a), a generating function approach was used to obtain analytical expressions for the time-independent (‘stationary’) solution to the model. For a visual guide to the upcoming analysis, the reader is again referred to Fig. 5.

The model for stochastic gene expression considered here is given by the reaction schemeThe modelled gene can hence switch between inactive and active states which are denoted by D and , respectively, with corresponding switching rates and . The active gene is transcribed to mRNA (M) with rate ; mRNA is translated to protein (P) with rate . Finally, mRNA decays with rate , while protein decays with rate .

As in model A, autoregulatory terms can be added to the core reaction scheme in (3.1). Since both mRNA and protein are modelled explicitly, we can identify four distinct autoregulatory mechanisms, in analogy to those in (2.2a) and (2.2b): Autoactivation can be achieved either by mRNA or by protein, with rates or , respectively; similarly, autorepression can occur either through mRNA or through protein, with respective rates or .

The CME system associated to the reaction scheme in (3.1) is given by Here, () represents the probability of m mRNA and n protein being present at time t while the gene is either inactive (0) or active (1). As in (2.3) and (2.5), the time variable is nondimensionalised by the protein decay rate ; other model parameters are scaled asWe note that the above scaling was also used by Shahrezaei and Swain (2008a). The effects of incorporating the autoregulatory mechanisms in (3.2) into the CME system (3.3) are specified in Table 1.

Table 1

Contribution to the right-hand sides of (3.3) that is due to incorporation of the autoregulatory mechanisms in (3.2)

Autoregulation type	Contribution to (3.3a)	Contribution to (3.3b)
mRNA autoactivation	-a_Md_1mPm,n(0)	+a_Md_1mPm,n(0)
mRNA autorepression	+r_Md_1mPm,n(1)	-r_Md_1mPm,n(1)
Protein autoactivation	-a_Pd_1nPm,n(0)	+a_Pd_1nPm,n(0)
Protein autorepression	+r_Pd_1nPm,n(1)	-r_Pd_1nPm,n(1)

Since the probabilities () in (3.3) depend on both the mRNA number m and the protein number n, we introduce probability-generating functions that are defined by double asymptotic series:For coefficients that obey the CME system (3.3), the associated generating functions satisfy The effects of incorporating the autoregulatory mechanisms in (3.2) into system (3.6) are specified in Table 2.

Table 2

Contribution to the right-hand sides of (3.6) that is due to incorporation of the autoregulatory mechanisms in (3.2)

Autoregulation type	Contribution to (3.6a)	Contribution to (3.6b)
mRNA autoactivation	-a_Md_1w∂_wF^(0)	+a_Md_1w∂_wF^(0)
mRNA autorepression	+r_Md_1w∂_wF^(1)	-r_Md_1w∂_wF^(1)
Protein autoactivation	-a_Pd_1z∂_zF^(0)	+a_Pd_1z∂_zF^(0)
Protein autorepression	+r_Pd_1z∂_zF^(1)	-r_Pd_1z∂_zF^(1)

As before, the PDE system (3.6) for the generating function can be reformulated as a system of ODEs via the method of characteristics. The differential operator now gives rise to the characteristic system where is the derivative along the characteristic, which is parametrised by s. For simplicity, we have introduced the new variables u and v, which are defined asrespectively. On the resulting characteristics, Eq. (3.6) yields the system of ODEs First, we note that the characteristic system (3.7) can be solved explicitly in terms of the incomplete Gamma function , see §8.2(i) of NIST Digital Library of Mathematical Functions ), yielding However, due to its complex nature, the expression for u(s) given in (3.10a) cannot be used to obtain explicit expressions for the generating functions that solve system (3.9). Therefore, inspired by the analysis by Shahrezaei and Swain (2008a) and Popović et al. (2016), we make the following assumption:

Assumption 3.1

We assume that the decay rate of protein () is smaller than the decay rate of mRNA (); specifically, we writewhere is sufficiently small.

The resulting scale separation between mRNA and protein decay rates is well-documented in many microbial organisms, including in bacteria and yeast (Shahrezaei and Swain 2008a; Yu et al. 2006).

For clarity of presentation, we make an additional assumption here.

Assumption 3.2

We assume that all other model parameters , , , and , as defined in (3.4), are in .

Remark 3.3

Although Assumption 3.2 is not strictly necessary for the upcoming analysis, it is beneficial. It is worthwhile to note that the analytical scheme presented in this section can be applied in a straightforward fashion in cases where Assumption 3.2 fails, which is particularly relevant in relation to previous work (Shahrezaei and Swain 2008a; Feigelman et al. 2015), where the CME system (3.3) is studied for parameter values far beyond the range implied by Assumption 3.2.

Using Assumption 3.1, we can write the characteristic system (3.7) as Since is assumed to be small, we can classify (3.12) as a singularly perturbed slow-fast system in standard form; see Kuehn (2015). A comprehensive slow-fast analysis of Eq. (3.12) was carried out by Popović et al. (2016); we highlight some relevant aspects of that analysis here.

System (3.12) gives rise to a critical manifold . If is asymptotically small, orbits of (3.12) can be separated into slow and fast segments, using Fenichel’s geometric singular perturbation theory (Kuehn 2015). The critical manifold is normally repelling; in other words, orbits converge to in backward time at an exponential rate. For initial conditions asymptotically close to , orbits initially follow closely for some time, after which they move away from under the fast dynamics; see also Fig. 7.

Fig. 7

Phase space dynamics of systems (3.12) and (3.13). The slow flow along is indicated by single arrows; the fast dynamics transverse to are denoted by double arrows

To analyse the fast dynamics of system (3.12), we introduce the fast variable ; in terms of , (3.12) is hence expressed as where . We can solve (3.13b) explicitly and write the result as a power series in , which yieldsExpressing the solution to (3.13a) as a power series in , as well, i.e. writingsubstituting (3.15) into (3.13a), and making use of (3.14), we findwith initial conditionsThe first two terms in (3.15) are thus given by

We now employ the expansion in (3.15) to obtain explicit expressions for the generating functions (). In the fast variable , system (3.9) takes the form As in the analysis of model A, recall Sect. 2.3, we rewrite system (3.19) as a second-order ODE for to findNext, we use (3.19a) to express in terms of asTo incorporate the expansion for u in (3.15), we also expand () in powers of , writingSubstitution of (3.22) into (3.20) then yields andfor . By combining (3.21) with (3.22), we obtainandfor . We can solve Eqs. (3.23) and (3.24) iteratively—taking into account the additional condition on in (3.25)—to find for the first three terms which, upon substitution into (3.26), yields here, and are free constants, to be determined by initial conditions; see Sect. 3.4.

Contrary to common practice in the study of slow-fast systems such as (3.19), we forego a detailed analysis of the slow system (3.9), continuing our discussion with the determination of appropriate initial conditions; cf. Sect. 2.4. For details on why the slow dynamics is disregarded, the reader is referred to Remark 3.5.

Initial conditions and inverse transformation

To complete our analytical method, we discuss the determination of initial conditions and the reconstruction of the solution to the original PDE system (3.6), as was done for model A in Sects. 2.4 and 2.5, respectively.

Initial conditions

We follow the reasoning of Sect. 2.4, and determine appropriate initial conditions by considering the original CME system (3.3).

At time , we prescribe initial mRNA and protein numbers and , respectively. As in Sect. 2.4, we again introduce the parameter , which is defined as follows; compare Definition 2.4:

Definition 3.4

Consider the CME system (3.3). We define to be the probability that the gene modelled by the reaction scheme in (3.1) is switched off at time , given initial mRNA and protein numbers and , respectively.

For the probabilities (), Definition 3.4 implies that Using the definition of the generating functions in (3.5), we find taking into account the change of variables in (3.8).

Next, we can infer e.g. from (3.10) and well-known properties of the incomplete Gamma function—see §8.2 of NIST Digital Library of Mathematical Functions —that solutions to (3.7) exist globally, i.e. for all s. In particular, given an arbitrary triple , we can apply the inverse flow of (3.7) to flow backward for time . We conclude that the characteristics of the operator , which are given by the orbits of system (3.7), foliate the (t, u, v)-coordinate space over the -plane when the parameter along the characteristic (s) is identified with the time variable (t). We can therefore uniquely identify such a characteristic—interpreted as a fibre over the -plane—by its base point. Because orbits of (3.7) provide a parametrisation of the underlying characteristics, the coordinates of that base point are given by , where are the initial values of the corresponding orbit of (3.7). Hence, the conditions in (3.30) yield, on each characteristic, As implies , we can apply the initial conditions in (3.31) to the solutions of the ODE system (3.19) in order to determine the free constants and in (3.27) and (3.28). Combining the power series expansion in (3.22) with (3.31), we obtain which implies in (3.27) and (3.28).

Inverse transformation

Since the (t, u, v)-coordinate space is foliated by the characteristics of the operator , any point (t, u, v) lies on a unique characteristic. Flowing backward along that characteristic to its intersection with the -plane, we can determine the corresponding base point by inverting the relations in (3.10). Since the dynamics of the v-coordinate do not depend on u, we may use (2.56) to express in terms of t and v only. Taking the resulting expression as input for inverting (3.10a), we obtain the inverse characteristic transformationwith Under Assumption 3.1, we can employ the power series expansion in (3.15), in combination with the recursive set of ODEs in (3.16) and initial conditions as in (3.17), as an alternative to (3.35a) to obtain as a function of u, v (or ), and . Rewriting the result as a power series in , we findnaturally, (3.14) gives rise toSince the independent variable in (3.16) is , and since s is naturally identified with t, we replace with in (3.36) to obtain a perturbative expansion for .

The solution to the PDE system (3.6) can now be found to satisfactory order in by applying the inverse transformation in (3.34) to the solutions given in (3.27) and (3.28), taking into account the values of and in (3.33). In other words,where () on the right-hand side of the above expression denotes the solution to Eqs. (3.23a)–(3.26), with initial conditions as in (3.32).

Remark 3.5

The absence of any detailed analysis of the characteristic system in its slow formulation, Eq. (3.12), can be argued as follows, by considering the corresponding phase space, as depicted in Fig. 7.Thus, although the construction of a composite ‘(initially) slow—(ultimately) fast’ expression of as a solution to systems (3.12) and (3.13) certainly makes sense from a dynamical systems perspective, the extreme lack of sensitivity of orbits on their initial conditions may prevent such a composite expansion from being useful for obtaining solutions to the original PDE system (3.6).

For arbitrary initial conditions , the dominant dynamics are fast, since the critical manifold is normally repelling. In other words, solutions are generally repelled away from under the fast dynamics.

All orbits that have their initial conditions on the same Fenichel fibre are exponentially close (in ) to each other near the slow manifold that is associated to the critical manifold . Therefore, flowing backward from (u, v) to —as expressed through the inverse transformation in (3.34) that yields the corresponding PDE solution—may introduce exponentially large terms in the transformation, precluding any sensible series expansion.

Remark 3.6

In Sect. 3.3, explicit expressions are given for the expansion of up to only, cf. (3.27); similarly, is approximated to in (3.28), for the sake of brevity. It is worthwhile to note that a lower bound on the order of the expansion is stipulated by the application; recall Sect. 1.1: the sampling time can be considered as a minimum time interval over which the results of our analysis should be (reasonably) accurate. To that end, we have to compare with , the parameter defining the fast time scale on which the above analytical results have been derived. We can then apply the classical theory of Poincaré expansions (Verhulst 2000) to infer that, if for example —which implies in the fast time variable —the generating functions should at least be expanded up to for the resulting approximation to be accurate to .

Autoregulation

The inclusion of any type of autoregulation into system (3.6)—which is equivalent to the addition of model terms from Table 2 to the right-hand sides of the corresponding equations—precludes the direct application of the method of characteristics, as the resulting partial differential operators in Eqs. (3.6a) and (3.6b) do not coincide anymore. To resolve that complication, we follow the approach of Sect. 2, making the following assumption:

Assumption 3.7

We assume that the autoregulation rates , , , and , as defined in (3.2), are small in comparison with the protein decay rate ; specifically, we writewhere is sufficiently small.

Next, we expand the generating functions () as power series in ; recall (2.11):To demonstrate the procedure, we include mRNA autoactivation in (3.6), see again Table 2; the analysis of the remaining autoregulatory mechanisms can be performed in a similar fashion. Substitution of (3.40) now yields cf. (2.12). System (3.41) is then amenable to the method of characteristics. In fact, employing the same characteristics as in the unperturbed setting, recall (3.7), we findwhile the partial differential operators in Table 2 transform into here, and are as given in (3.10). Thus, the mRNA autoactivation system (3.41) transforms to To obtain explicit solutions to system (3.44), we adopt Assumption 3.1 and revert to the fast time scale to write the dynamical system (3.44) as a second-order ODE for , which yieldsUsing (3.44), we can express in terms of asTo solve (3.45) (recursively), we expand () in powers of :recall Eq. (3.22). Together with the series expansion for u in (3.15), we thus obtain andfor ; compare with Eqs. (3.23a), (3.23b), and (3.24). The coefficients in the above expression are defined from an expansion of the autoregulation term asFrom (3.46), we obtainandfor ; recall (3.25) and (3.26). To solve Eqs. (3.48) through (3.53) iteratively, we fix m—the order of the expansion in —and determine the solution to satisfactory order in n, the order of the expansion in . Then, we increase m to and take the result as input for the dynamics at order . The resulting repeated iteration procedure yields an explicit expression for the generating functions () as double asymptotic series in both and .

The determination of appropriate initial conditions is largely analogous to the non-autoregulated case; see Sect. 3.4.1. However, with the inclusion of autoregulation into model B, we need to incorporate the possibility that depends on the corresponding autoregulation rates. As in the case of model A, we expand as a power series in :recall (2.51). In that case, the initial conditions for can be inferred from (3.32) to give Solutions to Eqs. (3.48)–(3.53) that incorporate the conditions in (3.55) for all types of autoregulation introduced in (3.2) can be found in Appendix B. Finally, our previous results on the inverse characteristic transformation in the non-autoregulated case from Sect. 3.4.2 can now be applied in a straightforward fashion to give solutions to the PDE system (3.6) with added autoregulation.

To summarise Sect. 3, we combine the analysis of the previous subsections to state our main result.

Main result: The PDE system (3.6) can be solved for sufficiently large values of ; see Assumption 3.1. Its solutions () are expressed as power series in the small parameter ; cf. (3.22). The coefficients in these series, written in terms of the shifted variables u and v defined in (3.8), can be found byTo illustrate the procedure described above, we state the resulting explicit expressions for the leading-order solution to (3.6) in terms of the original variables w, z, and t here: Note that the sum corresponds precisely to the leading-order fast expansion found in Equation (21) of Bokes et al. (2012b).

solving recursively the second-order ODEs (3.23a) through (3.24) and using the identities in (3.25) and (3.26), incorporating the initial conditions in (3.32);

subsequently applying the inverse transformations in (3.36) and (3.37) to the resulting solutions;

and, finally, substituting .

If autoregulation as in (3.2) is incorporated into model B, the main result can be formulated as follows.

Main result (autoregulatory extension): The PDE system (3.6) incorporating any one type of autoregulation from Table 2 can be solved as long as is sufficiently large and is sufficiently small; see Assumptions 3.1 and 3.7, respectively. Its solutions () are expressed as double power series in the small parameters and , viz.The coefficients in these series, written in terms of the shifted variables u and v defined in (3.8), can be found by

solving recursively the second-order ODEs (3.48) through (3.50) for fixed m and using the identities in (3.52) and (3.53), incorporating the initial conditions in (3.55);

increasing m to , and repeating step (1a) until a sufficient accuracy in (and ) is attained;

subsequently applying the inverse transformations in (3.36) and (3.37) to the resulting solutions;

and, finally, substituting .

From generating function to propagator

The final step in our analytical method consists in reconstructing the probabilities (model A) and (model B), respectively, from the explicit expressions for the associated generating functions (), which were the main analytical outcome of Sects. 2 and 3.

In principle, the relation between probabilities and probability-generating functions is clear from the definition of the latter, and is given in (2.6) and (3.5), respectively. Specifically, probabilities can be expressed in terms of derivatives of their generating functions as follows: However, explicit expressions for the nth and th order derivatives, respectively, of these generating functions become progressively unwieldy with increasing m and n. Indeed, from the expressions for the generating functions obtained previously, which combine (2.26) for specific initial conditions as in (2.54) with the inverse characteristic transformation in (2.56), it is clear that finding explicit expressions for derivatives of arbitrary order is very difficult indeed, if it is possible at all.1

To complete successfully the final step towards approximating propagators for parameter inference in the present setting, we abandon the requirement of deriving explicit expressions for the probabilities and . Instead, we use the standard Cauchy integral formula for derivatives of holomorphic functions to write here, is a suitably chosen contour around , while is a (double) contour around .

The above expression of probabilities as integrals is well suited for an efficient numerical implementation, which is naturally incorporated into the realisation of the parameter inference scheme discussed in Sect. 1. From a numerical perspective, the integral formula in (4.2) has the additional advantage that the values of on (a discretisation of) the integral contours and , respectively, only have to be determined once to yield propagators for any values of m and n and fixed initial states and . Moreover, we are free to choose the integration contours and , which allows us to accelerate the calculation of these integrals; see Bornemann (2011). Here, we note that the choice of circular integration contours with unit radius, and subsequent discretisation of those contours as M-sided and N-sided polygons, respectively, coincides with the ‘Fourier mode’ approach, as presented by Bokes et al. (2012a).

Remark 4.1

By introducing the Cauchy integral formula for derivatives of holomorphic functions in (4.2), we implicitly assume that the integration contours and are chosen such that they lie completely within the open neighbourhoods of the origin in and , respectively, where the canonical complex extensions of the generating functions —which exist by the Cauchy-Kowalevski theorem—are holomorphic. In other words, and must be chosen such that any poles of lie outside of these integration contours. The expansion for in (3.36) shows that this is not a moot point: the generating functions resulting for model B, as established in Sect. 3, will generically have a pole at , i.e. at . As is positive, by (3.4), choosing the z-contour of to be a circle with at most unit radius allows us to avoid that pole.

Incorporation of

In the course of the analysis presented in Sects. 2 and 3, the introduction of the parameter was necessary to obtain definite, explicit expressions for the generating functions as solutions to the PDE systems (2.7), (2.8), and (3.6); see Definitions 2.4 and 3.4. The successful implementation of these expressions in a parameter inference scheme requires us to decide how to incorporate that new parameter. We identify three options here.Note that step 2 has been anticipated in the analysis of model A, by introducing the series expansion in (2.51). Indeed, by Assumption 2.2, we can expand as Likewise, when autoregulation is added to model B, Assumption 3.7 implies an expansion for of the form

Before implementing parameter inference, we can marginalise over the new parameter to eliminate it altogether, using a predetermined measure , which adds an additional integration step to the requisite numerical scheme.

We can make a choice for that is based on the specifics of the model under consideration. Thus, exploiting the Markov property of the stochastic models underlying (2.3), (2.5), and (3.3), we may use the switching rates and any autoregulation rates to express in model A as the corresponding expressions for model B read

We can determine ‘experimentally’ by including the latter in the parameter set that is to be inferred in the (numerical) process of parameter inference.

Discussion and outlook

In the present article, we have developed an analytical method for obtaining explicit, fully time-dependent expressions for the probability-generating functions that are associated to models for stochastic gene expression. Moreover, we have presented a computationally efficient approach which allows us to derive model predictions (in the form of propagators) from these generating functions, using the Cauchy integral formula. It is important to note that our method does not make any steady-state or long-evolution-time approximations. On the contrary, the perturbative nature of our approach naturally optimises its applicability over relatively short (or intermediate) time scales; see also Remark 3.6. As is argued in Sect. 1.1, such relatively short evolution times naturally occur in the calculation of quantities such as the log-likelihood, as defined in Eq. (1.1). Therefore, our analytical approach is naturally suited to an implementation in a Bayesian parameter inference scheme, such as is outlined in Sect. 1.1.

As mentioned in Sects. 2.2 and 3.3, the introduction of Assumptions 2.2 and 3.7 in our analysis of the systems of PDEs and ODEs that are obtained via the generating function approach is necessary for determining explicit expressions for the generating functions themselves. Therefore, we can only be certain of the validity of our approach if we assume that the autoregulation rates are small in comparison with other model parameters, as is done there. Moreover, in the analysis of model B, we have to assume that the protein decay rate is smaller than the decay rate of mRNA; recall Assumption 3.1. That assumption is valid for a large class of (microbial) organisms (Shahrezaei and Swain 2008a; Yu et al. 2006); however, it is by no means generic, as the two decay rates are often comparable in mammalian cells (Schwanhäusser et al. 2011; Vogel and Marcotte 2012). Since the accuracy of approximation of the explicit expressions for the generating functions derived here is quantified in terms of orders of the perturbation parameter(s), see e.g. Remark 3.6, violation of Assumption 2.2, 3.1, or 3.7 will decrease the predictive power of the results obtained by the application of the analytical method developed in the article.

The method which is described in Sect. 1.2, and outlined visually in Fig. 5, hence provides a generic framework for the analysis of stochastic gene expression models such as model A (Fig. 3 and Sect. 2.1) and model B (Fig. 4 and Sect. 3.1). Note that, for example, the steady-state and long-evolution-time approximations derived by Shahrezaei and Swain (2008a) could be extended to autoregulatory systems via the same approach. However, as is apparent from the (differences between the) analysis presented in Sects. 2 and 3, the ‘path to an explicit solution’ is highly model-dependent. The decision on which analytical techniques to apply, such as the perturbative expansion postulated in (3.22), has to be made on a case-by-case basis. The success of the method presented in the article fully depends on whether the resulting dynamical systems can be solved explicitly. To that end, it is highly beneficial that the systems (2.19) and (3.19) obtained here are linear, which is a direct consequence of the fact that all reactions described in the reaction schemes in (2.1) and (2.2), as well as in (3.1) and (3.2), are of first order. Inclusion of second-order reactions would introduce both nonlinear terms and second-order differential operators in the PDE systems for the corresponding generating functions, which would severely increase the complexity of these systems, thus preventing us from obtaining explicit solutions.

The method presented in this article, and the results thus obtained, can be seen, first and foremost, as the natural extension of previous work by Popović et al. (2016). Analytical results for the classes of models studied here can be found in several earlier articles. We mention the article by Shahrezaei and Swain (2008a), where a leading-order approximation was obtained in a long-evolution-time and steady-state limit. Bokes et al. (2012a) derived analytical expressions for stationary distributions in a model that is equivalent to that considered by Popović et al. (2016). Also for that model, a time-scale separation was exploited by Bokes et al. (2012b), in a manner that is similar to the present article, to obtain leading-order analytical expressions on both time scales. The model that is referred to as Model A in Sect. 2 was analysed in a steady-state setting by Iyer-Biswas and Jayaprakash (2014) via the Poisson integral transform. A similar model was studied by Hornos et al. (2005), were a generating function approach was used; making a steady-state Ansatz, the authors were able to obtain an explicit solution for the generating function in terms of special (Kummer) functions; see also NIST Digital Library of Mathematical Functions . The same model was later solved in a fully time-dependent context by Ramos et al. (2011), after a cleverly chosen variable substitution, in terms of another class of special (Heun) functions; cf. again NIST Digital Library of Mathematical Functions .

Other authors have attempted to solve several classes of CMEs directly, i.e. without resorting to generating function techniques or integral transforms. A noteworthy example is the work of Jahnke and Huisinga (2007) on monomolecular systems. Another, more recent example can be found in the work by Iserles and MacNamara (2017), where exact solutions are determined for explicitly time-dependent isomerisation models.

It is important to emphasise that the ‘time dependence’ referred to in the title of the present article is solely due to the dynamic nature of the underlying stochastic process, and that it hence manifests exclusively through time derivatives in the associated CMEs, such as e.g. in (2.3). In particular, none of the model parameters are time-dependent, as opposed to, for example, the system studied by Iserles and MacNamara (2017). The inclusion of such explicitly time-dependent parameters would be a starting point for incorporating the influence of (extrinsic) noise in the context of the model categories considered in the article.

The availability of analytical expressions for generating functions does, in principle, allow one to try to obtain insight into the underlying processes by studying the explicit form of said expressions, as has been done e.g. by Bokes et al. (Bokes et al. 2012a, b). However, the complex nature of the processes we analyse here seems to preclude such insights. For example, the integrals over confluent hypergeometric functions, which appear in (2.29), cannot themselves be efficiently expressed in terms of (other) special functions. Still, that complication does not necessarily pose an obstacle to the application we ultimately have in mind, i.e. to Bayesian parameter inference. As the last step in our method—the extraction of propagators from generating functions, see Sect. 4—is numerical, the precise functional form of the generating function is not of importance. The mere fact that an explicit expression can be obtained is sufficient for the application of the Cauchy integral formula, where these generating functions enter into the calculation of the appropriate integrals; see again Sect. 4.

The analytical approach explored in the article does not, of course, represent the only feasible way of obtaining numerical values for propagator probabilities, which can, in turn, serve as input for a Bayesian parameter inference scheme. For an example of a direct numerical method in which the Cauchy integral plays a central role, the reader is referred to the work by MacNamara (2015). Our main motivation for pursuing an analytical alternative is reducing the need for potentially lengthy numerical simulations. An efficient implementation of the resulting expressions can result in (significantly) reduced computation times; see, for example, the work by Bornemann (2011). The optimisation of the underlying numerical procedures is, however, beyond the scope of the present article in particular, and of our research programme in general.

The analytical results obtained thus far, as presented in the article, are ready for implementation in a Bayesian parameter inference framework. An analysis of the performance of the resulting approximations to the associated generating functions in the spirit of the article by Feigelman et al. (2015), where parameter inference is tested on simulated data based on specific stochastic models, is ongoing work. Moreover, the successful application of our analytical method to specific model categories, such as are represented by model A and model B, suggests several feasible expansions of the ‘model library’ for which explicit expressions for the corresponding generating functions can be constructed. Thus, stochastic models comprised of multiple proteins represent a natural next stage, bringing the analysis of toggle switch-type models within reach. In addition, one could begin exploring the vast field of gene regulatory networks by considering a simple two-protein system with, for example, activator-inhibitor interaction. Under the assumption of small interaction rates, the resulting PDE system for the associated generating function would be directly amenable to the analytical method described in the article. The analysis of these and similar systems could be a topic for future research.