Population Continuous-Time Markov Chains (PCTMCs) provide a widely used framework to capture stochastic interactions between groups of identical agents. This subclass of Continuous-Time Markov Chains (CTMCs) is used to describe the stochastic dynamics of systems in various domains. Prominent applications are chemical reaction networks in quantitative biology , epidemic spreading , performance analysis of technical and information systems [9, 20] as well as the behavior of collective adaptive systems .
For the quantitative analysis of CTMCs, many approaches have been developed, where properties of interest are often expressed in terms of temporal logics such as CSL [1, 4, 3], MTL , and timed-automata specifications [13, 37]. In addition, there exist efficient software tools [28, 34, 15]
. A central problem in this context is the computation of reachability probabilities.
Popular exact methods for CMTCs rely on numerical approaches that explicitly consider each system state individually. A major problem is that these methods cannot scale in the context of population models with large copy numbers of agents. A popular alternative to tackle this problem is statistical model checking, which is based on stochastic simulation . For PCTMCs arising in the context of chemical reaction networks, trajectories of the process are usually generated using the Stochastic Simulation Algorithm (SSA) 
. However, since the number of possible interactions grows with the number of agents, stochastic simulations of PCTMCs are time-consuming. Moreover, they are subject to inherent statistical uncertainty and give only statistically estimated bounds.
Recent work concentrates on numerical methods for PCTMCs that approximate the statistical moments of the system without the need to consider the probability of each state. For groups of identically behaving agents, it is possible to derive systems of differential equations for the evolution of the statistical population moments[8, 47, 10, 19, 46, 20]
. However, as the system of exact moment equations is not closed, approximation schemes typically rely on certain assumptions about the underlying probability distribution. For example, one might employ a “low dispersion closure” which assumes that higher-order moments are the same as those of a normal distribution. Such approximations are, by nature, ad-hoc and do typically not come with any guarantees. Here, we do not need such closure schemes and retain guaranteed results up to the numerical accuracy of the computations.
Moment-based methods often scale well in terms of population sizes. However, it is not possible to control the effects of such approximations, which in some cases can lead to large errors . This issue reverberates on the application of these methods to compute reachability probabilities and mean first passage times [25, 10, 11]. Moreover, they can suffer from numerical instabilities, in particular, when the maximum order of the considered moments has to be increased to more appropriately describe the underlying distribution.
Here, we build on recent work on moment bounds [43, 18] to propose a method to compute bounds for Mean First Passage Times (MFPTs) in PCTMCs. For a set of states, the MFPT within a fixed time horizon directly characterizes the probability of reaching that set within time units. Thus, safe upper and lower bounds on MFPTs can constitute a core component for the verification of properties in PCTMCs. Our approach is based on a martingale formulation of the stopped process that we derive from the exact moment equations. From this formalization, we deduce a set of linear moment constraints from which we derive upper and lower moment bounds using semi-definite programming (SDP). Monotone sequences of both upper and lower bounds can be obtained by increasing the order of the relaxation. Crucially, no closure approximations are introduced. Therefore the bounds are strict up to the numerical accuracy of the SDP solver. For our numerical investigations, we concentrate on examples from biology and find encouraging results already for a small number of moments. For instance, in one of our case studies
SSA runs are necessary to achieve a relative width of 0.9% for the MFPT confidence interval. The SDP solver, however, returns a guaranteed interval with a relative width of 0.3%.
In summary, this paper presents the following novel contributions:
the derivation of moment constraints for bounding first passage times and reachability probabilities using a convex programming scheme
the extension of this scheme to stochastic hybrid systems exhibiting multimodal behavior
a scaling strategy for improved robustness during optimization
The paper is structured as follows: Section 2 covers work related to the analysis of first passage times in PCTMCs and recent work on moment bounds. Section 3 introduces the PCTMC framework and its semantics. In Section 4 we derive a martingale from the moment dynamics of a PCTMC. Based on this process, in Section 5 we formulate linear and semi-definite constraints to state a semi-definite program to compute bounds on the MFPT and reachability probabilities. In Section 6, we discuss the practical considerations of the SDP implementation and provide results on a set of case studies. Finally, in Section 7 we provide concluding remarks and directions of future work.
2 Related Work
Considerable effort has been directed at the analysis of first passage time distributions in PCTMCs. Most works can either focus on an explicit state-space analysis [5, 39, 33, 32] or employ approximation techniques for which, in general, no error bounds can be given [45, 25, 11]. For some model classes such as kinetic proofreading, analytic solutions are possible [39, 6, 29].
Barzel and Biham  propose a recursive scheme that consists of one equation for each state, expressing the average time the system needs to transition from that state to the target state. Kuntz et al. 
propose to employ moment bounds in a linear programming approach to compute exit time distribution using state-space truncation schemes. In Ref. the authors propose a finite state-space projection scheme to bound first passage time distributions
Hayden et al.  use moment closure approximations and Chebychev’s inequality to gain an understanding of first passage time dynamics. Schnoerr et al.  also employ a moment closure approximation and further approximate threshold functions to derive an approximate first passage time distribution. Bortolussi and Lanciani  use a mean-field approximation which is required to reach the target region.
Recently, several groups independently suggested the use of semi-definite optimization for the computation of moment bounds for the limiting distribution [21, 17, 31, 43]. In this approach, the differential equations describing the moment dynamics are set to zero and form linear constraints. Alongside, semi-definite constraints can be placed on the moment matrices. These give a semi-definite program that can be solved efficiently.
This approach has been extended to the transient case [18, 44]. The approach is similar in both works and is a cornerstone of the MFPT analysis presented here. They differ mainly by the fact that Sakurai and Hori apply a polynomial time-weighting , while Dowdy and Barton use an exponential one . We adopt the former approach because it can be naturally adapted to the description of densities over time. The resulting forms can also be adapted to statistical estimation problems .
Semi-definite programming has been applied to a wide range of problems, including stochastic processes in the context of financial mathematics [36, 30]. For good introductions and overviews of application areas, we refer the reader to Parrilo  and, more recently, Lasserre .
Particularly relevant for this work is the application of convex optimization to first passage times. Helmes et al.  formulated a linear program using the Hausdorff moment conditions to bound moments of the first passage time distribution in Markovian processes. Semi-definite optimization has been successfully applied in financial mathematics by Kashima and Kawai , as well as Lasserre et al.  to bound prices of exotic options. Here, the approach by Lasserre is adapted to PCTMCs.
A Population Continuous-Time Markov Chain (PCTMC) describes the interactions among a set of species in a well-stirred reactor111In the sequel, we will also use other letters than as species names.. Since we assume that all reactant molecules are equally distributed in space, we only keep track of the overall copy number of molecules of each species. Therefore the state-space is . The interactions are expressed as reactions
with a certain gain and loss of molecules, given by the non-negative integer vectorsand for some reaction , respectively. Such a reaction is denoted as
The reaction rate constant determines the propensity function of the reaction. If just a constant is given, mass-action propensities are assumed, where for we define
The system’s behavior is described by a stochastic process . We denote the abundance of a given species in by . The propensity gives the infinitesimal probability of a reaction occurring, given a state . That is, for and a small time step ,
Therefore, given a system of reactions, the semantics of is given by a continuous-time Markov chain (CTMC) on with infinitesimal generator matrix with entries
Accordingly, given an initial distribution on , the time-evolution of the process’ distribution is given by the Kolmogorov forward equation. For a single state, it is commonly referred to as the chemical master equation (CME)
where and .
In this work, we are interested in first passage times of such processes. That is the time, the process first enters a set of target states . Naturally, the analysis of first passage times is equivalent to the analysis of times at which the process exits the complement . More formally, the first passage time for some target set
is defined as the random variable
Consider the following simple non-linear PCTMC as an example.
Model 3.1 (Dimerization)
We first examine a simple dimerization model on an unbounded state-space with reactions
and initial condition . The semantics is given by a CTMC , where . The reaction propensities according to (2) are and . The change vectors , , , and . Consequently, and .
In this example, we are interested in the time at which exceeds some threshold . With the framework presented in the sequel, one can bound the expected value of this time. Further, it is possible to impose a time-horizon , and find bounds on the probability of for some . The employed framework is centered around semi-definite relaxations of the generalized moment problem . These require linear constraints on the moments of measures. In the following section, we derive such constraints.
4 Martingale Formulation
Next, we will discuss equations for the evolution of the statistical moments of the process and a related martingale formulation. This is later used to derive linear constraints on the moments of appropriate measures that can be used to bound MFPTs.
In particular, we consider the dynamics of raw moments for and a fixed probability measure. The order of a moment is given by its exponent sum, i.e. . We can derive the time evolution of the raw moments directly from the CME in (5). Note, that the notion of the expected value can be generalized to any measure on a Borel-measurable space . There the -th raw moment is .
Let be a polynomial function,
. We can easily derive ordinary differential equations (ODEs) to describe the dynamics of. Specifically,
When choosing and and , for Model 3.1, for example, we get the following system of ODEs for the change of the first and second statistical moment of species
where we let for ease of notation. These ODEs cannot be integrated because the system is not closed. The right-hand side for moment always contains . To solve an initial value problem, one typically resorts to ad-hoc approximations of the highest order moments to close the system. Here we do not need such approximations because we do not numerically integrate such equations.
If we now assume that and remain finite for all , we can interchange summation and integral of a monomial and pull all expectation operators outside, i.e. for a polynomial
Hence, for (10) pulling the expectation operator outside yields a martingale , where
with respect to . When choosing with and it takes the form
where , , and are finite sequences resulting from the substitution of and and expansion of (11). We will use this martingale in the following section to derive linear constraints for the semi-definite program used to bound MFPTs.
5 Bounds for Mean First Passage Times
We now turn to the analysis of first passage times within some time-bound . Given some set the first passage time is given by the random variable
For this work, we only look at threshold hitting times, i.e. we set a threshold for species and thus .222Note, that this framework allows for a more general class of target sets, which are discussed in Section 5.4. In the sequel, we will use as a stopping time in our martingale formulation and consider instead of . Since (12) defines a martingale, remains a martingale by Doob’s optional sampling theorem . In particular, this implies that .
5.1 Linear Moment Constraints
To simplify our presentation, we fix an initial state , i.e. . Using and the form (12) for yields the following linear constraint on expected values.
where . For the ease of exposition, we now turn to first passage times of one-dimensional processes w.r.t. an upper threshold . In particular, we will consider moments of a one-dimensional process for . The approach proposed in the sequel, however, can be straightforwardly extended to multi-dimensional processes and more complex target sets .
Consider again Model 3.1 and assume that we are interested in the time at which species exceeds threshold . Since the abundance of does not influence , we can ignore species and treat the process as one-dimensional. Figure 1 shows three example trajectories: Two reach an upper threshold , while one reaches the final time horizon .
We notice, that (14) expresses a relationship between the process dynamics up to the hitting time via expected values of the time-integrals and the final process state at the hitting time via . In particular, we can distinguish between the following two positive measures [35, Chapter 9.2]:
Expected Occupation Measure supported on :
Exit Location Probability supported on :
where is a measurable set, i.e. and are elements of the Borel -algebras on and , respectively.
Using Figure 1, one can gain an intuition for these two measures. The expected occupation measure is shaded in blue. As the name implies tells us how much time the process spends in up to restricting to the time instants belonging to . In particular, . The exit location probability , while being a two-dimensional distribution, can be viewed as a composition of a density describing the time at which the process reaches (if it does) and a probability mass function on the states of the process if the time-horizon is reached without exceeding . We split the measure into and by conditioning on . Thus, and . To refer to the moments of these measures, we define partial moments
for some polynomial and some indicator function . Then
Therefore the linear moment constraints have the form
Next, we consider infinite sequences of partial moments by letting and range over the natural numbers. Let , , and denote the moment sequences of , , and , respectively.
Thus, the variable corresponding to becomes the objective of the optimization problem that we describe in the sequel.
5.2 Semi-Definite Constraints
It is a necessary condition for a positive measure that the moment matrices are positive semi-definite. A matrix is positive semi-definite, denoted by if and only if
As an example, let us consider a one-dimensional random variable with moment sequence . For moment order , the entries of the moment matrix are given by the raw moments. In particular, for where and the maximum order in the matrix is . For instance,
needs to be positive semi-definite. By Sylvester’s criterion this means and . We can easily see, that this entails
This restriction is natural since the variance is always non-negative.
It is crucial to restrict the measures , , and to their supports. This can be done, by defining polynomials that are non-negative on the intended support of the measure. For example, has support . We can now define
with the constraint , where the application of a polynomial such as to a moment matrix is formally defined for the multidimensional case in Section 5.4. Similarly, let to restrict to .
5.3 A semi-definite program to bound MFPTs
With the linear constraints on the measures (15) and (16) and the semi-definite constraints discussed in the previous sections, we can now formulate a semi-definite program (SDP). An SDP is an optimization over the cone of positive semi-definite -matrices under linear constraints:
with constant matrices , and constants , . Such a problem is convex and can be solved efficiently .
The derived linear equations and linear matrix inequalities can now be used to formulate an SDP. The full optimization problem has infinitely many constraints because there are infinitely many moments. We relax this problem by constructing the SDP using by choosing a finite order for the moment matrices . With each moment sequence we associate a sequence proxy variables used in the optimization problem. Now we can state the SDP relaxation to the MFPT problem for any order
5.4 Multi-Dimensional Generalization
For a general multi-dimensional moment sequence , the moment matrix is 
where row and column indices, and , are ordered according to the canonical basis
Equivalently, . For a moment sequence the semi-definite restriction must hold.
Measures can be restricted to semi-algebraic sets , where , are polynomials . This is done by placing restrictions on the localizing matrices. For each polynomial with coefficient vector , i.e. , the localizing matrix is
Requiring that this matrix is positive semi-definite restricts the measure to . This way we can, for example, restrict the moment sequence to measures that are positive w.r.t. dimension . Simply letting and requiring for gives us this restriction.
6 Implementation and Evaluation
The main challenge of finding a solution to the SDP problem in (20) is numerical stability. Usually, the moment sequences vary by many orders of magnitude. For an SDP solver to work, the moment matrices need to be re-scaled  such that moments only vary by few orders of magnitude. In other scenarios such as the bounding of general transient or steady-state moments, the scaling can be particularly difficult, because the magnitude of moments is generally not known a priori. However, for the MFPT problem, we propose the following moment scaling.
6.1 Moment Scaling
Using the fact that is often finite, it is possible to derive trivial bounds, which can be used to scale moments. If, for example, we have a one-dimensional process with a.s. and are interested in the hitting time of an upper threshold until time for
Thus, we fix a scaling vector with entries in the same order as the canonical base vector (21). Using this scaling vector, we can define a scaling matrix . Clearly, . Now we can formulate the optimization (20) over a scaled version instead of . The moment matrices of the exit location probabilities are scaled in the same way. Alternatively, one can use approximations such as moment closures or bounds obtained by lower-order relaxations.
6.2 Case Studies
As a first case study, we use Model 3.1 with parameters and . In this model, we are interested in the time at which the number of agents of type surpasses a threshold of 25 before some time-horizon , i.e. . First, we set no finite time horizon , i.e. . This is achieved by dropping the moments of measure in the linear constraints (20). The empirical FPT distribution based on 100,000 SSA simulations is given in Figure 2a and the bounds, given different moment orders, are given in Figure 2b. As we can see in Figure 2b, the bounds capture the MFPT precisely for orders 5, 6. The difference between upper and lower bound decreases roughly exponentially with increasing relaxation order . We found that this trend was consistent among the case studies presented here (cf. Figure 4).
Next, we look at first passage times within a finite time-horizon . In Figure 3a we summarize the bounds obtained for the MFPT over . While low-order relaxations (light) give rather loose bounds, the bounds are already fairly tight when using . In many cases, hitting probabilities, that is, the probability of reaching the threshold before time , are of particular interest. This is done by switching the optimization objective in (20) from the mass of the expected occupation measure to the mass of . In terms of moments, the objective changes from to . The need for such a scenario often arises in the context of model checking, where one might be interested in the probability of a population exceeding a critical threshold. By varying the time horizon, we are able to recover bounds on the cumulative density of the first passage time (Fig. 3b).
As a second study, we consider a 2-dimensional model by combining two independent dimerizations.
Model 6.1 (Parallel independent dimerizations)
As a FPT we consider the time at which either or surpasses a threshold of 200 or a time horizon of is reached, i.e.
As before, we ignore the products and since they do not influence . Still, the possible state-space reaches a size of . The SSA (using runs) gives the estimate which is captured tightly by the SDP bounds (cf. Table 1). For higher relaxation orders numerical issues prevented the solution of the corresponding SDPs.
6.3 Hybrid Models and Multi-Modal Behaviour
The analysis of switching times is a particularly interesting case of FPTs that arises in many contexts. Often mode switching in such systems can be described a modulating Markov process whose switching rates may depend on the system state (e.g. the population sizes). In biological applications, mode switching often describes a change of the DNA state [24, 49] and the analysis of switching time distribution is of particular interest [48, 5]. In the context of PCTMCs, the state-space of such models can be given as
This state is modeled by a population variables with binary domains. Therefore, at each time point, the state of these modulator variables is given by a set of Bernoulli random variables. When considering the moments of such a variable , clearly for all .
We apply a split of into the high count part and the binary part to the expectations in (7). Similarly, we split and with a case distinction over the mode variable, we arrive at a similar result as in :
Similarly to the general moment case, we can derive a constraint, by multiplying with a time-weighting factor and integrating.
|Double Dim. (Model 6.1)||lower||0.0010||0.0250||0.0275||0.0280||—|
|Gene Expression (Model 6.2)||lower||4.0000||6.0028||6.2207||6.3377||6.3772|
For simplicity, here we assume . Fixing appropriate sequences , , , and the constraint has the following form.
This way we can decompose the moment matrices such that for each mode , we have moment matrices composed of the respective partial moments. To this end, let be the partial moment w.r.t. . The moment constraint over the partial moments has a linear structure:
As an instance of a multi-modal system, we consider a simple gene expression with self-regulating negative feedback which is a common pattern in many genetic circuits .
Model 6.2 (Negative self-regulated gene expression)
This model consists of a gene state that is either on or off, i.e. , . Therefore the system has two modes.
The model parameters are and , a.s.
As a first passage time we consider
The results are summarized in Table 1. The estimated MFPT based on SSA samples is at confidence level. Note that our SDP solution for yields tighter moment bounds than the statistical estimation.
In Fig. 4 we summarize our results about the decrease of the interval widths for increasing relaxation order by plotting them on a log-scale. We see an approximately exponential decrease in . The semi-definite programs above were all solved within at most a few seconds.
State-based methods to compute reachability probabilities and first passage times for continuous-time Markov chains are not scalable due to state-space explosion, an issue exacerbated in population models. Moment-based methods offer an alternative for PCTMCs, which scales with the number of different populations in the system, but are approximate methods with little or no control of the error. In this paper, we bridge this gap by proposing a rigorous approach to derive bounds on first passage times and reachability probabilities, leveraging a semi-definite programming formulation based on appropriate moment constraints.
Our proposed scaling mitigates numerical instabilities of the SDP solvers, which are caused by the fact that moments typically span several orders of magnitude. However, the scaling only addresses the moment matrices but not the linear constraints which still contain values with varying orders of magnitudes. We, therefore, plan as future work to introduce an appropriate scaling for the linear constraints or to redefine the moment constraints (e.g. using an exponential time weighting ). Based on this investigation, we expect to make this approach applicable to more problems including, for example, the computation of bounds of rare event probabilities. Numerical instabilities due to moment values of largely differing orders of magnitudes are a current limitation of all moment-based methods. We expect that the development of more sophisticated scaling techniques will improve approximate moment-based methods, as well.
We would like to thank Andreas Karrenbauer for helpful comments on the usage of SDP solvers and Gerrit Großmann for the valuable comments on this manuscript. This work is supported by the DFG project “MULTIMODE”, and partially supported by the italian PRIN project “SEDUCE” n. 2017TWRCNB.
-  Aziz, A., Sanwal, K., Singhal, V., Brayton, R.: Verifying continuous time Markov chains. In: International Conference on Computer Aided Verification. pp. 269–276. Springer (1996)
-  Backenköhler, M., Bortolussi, L., Wolf, V.: Control variates for stochastic simulation of chemical reaction networks. In: Bortolussi, L., Sanguinetti, G. (eds.) Computational Methods in Systems Biology. pp. 42–59. Springer, Cham (2019)
-  Baier, C., Haverkort, B., Hermanns, H., Katoen, J.P.: Model-checking algorithms for continuous-time Markov chains. IEEE Transactions on software engineering 29(6), 524–541 (2003)
-  Baier, C., Haverkort, B., Hermanns, H., Katoen, J.P.: Model checking continuous-time Markov chains by transient analysis. In: International Conference on Computer Aided Verification. pp. 358–372. Springer (2000)
-  Barzel, B., Biham, O.: Calculation of switching times in the genetic toggle switch and other bistable systems. Physical Review E 78(4), 041919 (2008)
-  Bel, G., Munsky, B., Nemenman, I.: The simplicity of completion time distributions for common complex biochemical processes. Physical biology 7(1), 016003 (2009)
-  Bernardo, M., De Nicola, R., Hillston, J. (eds.): Formal Methods for the Quantitative Evaluation of Collective Adaptive Systems, Lecture Notes in Computer Science, vol. 9700. Springer International Publishing, Cham (2016)
-  Bogomolov, S., Henzinger, T.A., Podelski, A., Ruess, J., Schilling, C.: Adaptive moment closure for parameter inference of biochemical reaction networks. In: International Conference on Computational Methods in Systems Biology. pp. 77–89. Springer (2015)
-  Bortolussi, L., Hillston, J., Latella, D., Massink, M.: Continuous approximation of collective system behaviour: A tutorial. Performance Evaluation 70(5), 317–349 (May 2013)
-  Bortolussi, L., Lanciani, R.: Model checking Markov population models by central limit approximation. In: International Conference on Quantitative Evaluation of Systems. pp. 123–138. Springer (2013)
-  Bortolussi, L., Lanciani, R.: Stochastic approximation of global reachability probabilities of Markov population models. In: Computer Performance Engineering - 11th European Workshop, EPEW 2014, Florence, Italy, September 11-12, 2014. Proceedings. pp. 224–239 (2014)
-  Chen, T., Diciolla, M., Kwiatkowska, M., Mereacre, A.: Time-bounded verification of CTMCs against real-time specifications. In: International Conference on Formal Modeling and Analysis of Timed Systems. pp. 26–42. Springer (2011)
-  Chen, T., Han, T., Katoen, J.P., Mereacre, A.: Quantitative model checking of continuous-time Markov chains against timed automata specifications. In: 2009 24th Annual IEEE Symposium on Logic In Computer Science. pp. 309–318. IEEE (2009)
-  David, A., Larsen, K.G., Legay, A., Mikučionis, M., Poulsen, D.B., Sedwards, S.: Statistical model checking for biological systems. International Journal on Software Tools for Technology Transfer 17(3), 351–367 (2015)
-  Dehnert, C., Junges, S., Katoen, J.P., Volk, M.: A storm is coming: A modern probabilistic model checker. In: International Conference on Computer Aided Verification. pp. 592–600. Springer (2017)
Diamond, S., Boyd, S.: CVXPY: A Python-embedded modeling language for convex optimization. Journal of Machine Learning Research17(83), 1–5 (2016)
-  Dowdy, G.R., Barton, P.I.: Bounds on stochastic chemical kinetic systems at steady state. The Journal of chemical physics 148(8), 084106 (2018)
-  Dowdy, G.R., Barton, P.I.: Dynamic bounds on stochastic chemical kinetic systems using semidefinite programming. The Journal of chemical physics 149(7), 074103 (2018)
-  Engblom, S.: Computing the moments of high dimensional solutions of the master equation. Applied Mathematics and Computation 180(2), 498–515 (2006)
-  Gast, N., Bortolussi, L., Tribastone, M.: Size expansions of mean field approximation: Transient and steady-state analysis. Performance Evaluation 129, 60 – 80 (2019). https://doi.org/https://doi.org/10.1016/j.peva.2018.09.005
-  Ghusinga, K.R., Vargas-Garcia, C.A., Lamperski, A., Singh, A.: Exact lower and upper bounds on stationary moments in stochastic biochemical systems. Physical biology 14(4), 04LT01 (2017)
-  Gihman, I., Skorohod, A.: The theory of stochastic processes ii. 1975
-  Gillespie, D.: Exact stochastic simulation of coupled chemical reactions. The Journal of Physical Chemistry 81(25), 2340–2361 (1977)
-  Hasenauer, J., Wolf, V., Kazeroonian, A., Theis, F.J.: Method of conditional moments (MCM) for the chemical master equation. Journal of mathematical biology 69(3), 687–735 (2014)
Hayden, R.A., Stefanek, A., Bradley, J.T.: Fluid computation of passage-time distributions in large Markov models. Theoretical Computer Science413(1), 106–141 (2012)
-  Helmes, K., Röhl, S., Stockbridge, R.H.: Computing moments of the exit time distribution for Markov processes by linear programming. Operations Research 49(4), 516–530 (2001)
-  Hespanha, J.: Moment closure for biochemical networks. In: 2008 3rd International Symposium on Communications, Control and Signal Processing. pp. 142–147. IEEE (2008)
-  Hinton, A., Kwiatkowska, M., Norman, G., Parker, D.: Prism: A tool for automatic verification of probabilistic systems. In: International Conference on Tools and Algorithms for the Construction and Analysis of Systems. pp. 441–444. Springer (2006)
-  Iyer-Biswas, S., Zilman, A.: First-passage processes in cellular biology. Advances in Chemical Physics 160, 261–306 (2016)
-  Kashima, K., Kawai, R.: Polynomial programming approach to weak approximation of lévy-driven stochastic differential equations with application to option pricing. In: 2009 ICCAS-SICE. pp. 3902–3907. IEEE (2009)
-  Kuntz, J., Thomas, P., Stan, G.B., Barahona, M.: Rigorous bounds on the stationary distributions of the chemical master equation via mathematical programming. arXiv preprint arXiv:1702.05468 (2017)
-  Kuntz, J., Thomas, P., Stan, G.B., Barahona, M.: Approximation schemes for countably-infinite linear programs with moment bounds. arXiv preprint arXiv:1810.03658 (2018)
-  Kuntz, J., Thomas, P., Stan, G.B., Barahona, M.: The exit time finite state projection scheme: bounding exit distributions and occupation measures of continuous-time Markov chains. SIAM Journal on Scientific Computing 41(2), A748–A769 (2019)
-  Kwiatkowska, M., Norman, G., Parker, D.: Prism 4.0: Verification of probabilistic real-time systems. In: International conference on computer aided verification. pp. 585–591. Springer (2011)
-  Lasserre, J.B.: Moments, positive polynomials and their applications, vol. 1. World Scientific (2010)
-  Lasserre, J.B., Prieto-Rumeau, T., Zervos, M.: Pricing a class of exotic options via moments and sdp relaxations. Mathematical Finance 16(3), 469–494 (2006)
-  Mikeev, L., Neuhäußer, M.R., Spieler, D., Wolf, V.: On-the-fly verification and optimization of dta-properties for large Markov chains. Formal Methods in System Design 43(2), 313–337 (2013)
-  MOSEK ApS: MOSEK Optimizer API for C 188.8.131.52 (2018), https://docs.mosek.com/8.1/capi/index.html
-  Munsky, B., Nemenman, I., Bel, G.: Specificity and completion time distributions of biochemical processes. The Journal of chemical physics 131(23), 12B616 (2009)
-  O’Donoghue, B., Chu, E., Parikh, N., Boyd, S.: SCS: Splitting conic solver, version 2.1.0. https://github.com/cvxgrp/scs (Nov 2017)
-  Parrilo, P.A.: Semidefinite programming relaxations for semialgebraic problems. Mathematical programming 96(2), 293–320 (2003)
-  Porter, M.A., Gleeson, J.P.: Dynamical systems on networks. Frontiers in Applied Dynamical Systems: Reviews and Tutorials 4 (2016)
-  Sakurai, Y., Hori, Y.: A convex approach to steady state moment analysis for stochastic chemical reactions. In: Decision and Control (CDC), 2017 IEEE 56th Annual Conference on. pp. 1206–1211. IEEE (2017)
-  Sakurai, Y., Hori, Y.: Bounding transient moments of stochastic chemical reactions. IEEE Control Systems Letters 3(2), 290–295 (2019)
-  Schnoerr, D., Cseke, B., Grima, R., Sanguinetti, G.: Efficient low-order approximation of first-passage time distributions. Phys. Rev. Lett. 119, 210601 (Nov 2017). https://doi.org/10.1103/PhysRevLett.119.210601
-  Schnoerr, D., Sanguinetti, G., Grima, R.: Comparison of different moment-closure approximations for stochastic chemical kinetics. The Journal of Chemical Physics 143(18), 185101 (Nov 2015). https://doi.org/10.1063/1.4934990
-  Schnoerr, D., Sanguinetti, G., Grima, R.: Approximation and inference methods for stochastic biochemical kinetics—a tutorial review. Journal of Physics A: Mathematical and Theoretical 50(9), 093001 (Mar 2017). https://doi.org/10.1088/1751-8121/aa54d9
-  Spieler, D., Hahn, E.M., Zhang, L.: Model checking csl for Markov population models. arXiv preprint arXiv:1111.4385 (2011)
-  Stekel, D.J., Jenkins, D.J.: Strong negative self regulation of prokaryotic transcription factors increases the intrinsic noise of protein expression. BMC systems biology 2(1), 6 (2008)
-  Ullah, M., Wolkenhauer, O.: Stochastic approaches for systems biology. Wiley interdisciplinary reviews. Systems biology and medicine 2, 385–97 (07 2009). https://doi.org/10.1002/wsbm.78
-  Vandenberghe, L.: The cvxopt linear and quadratic cone program solvers. Online: http://cvxopt. org/documentation/coneprog. pdf (2010)