Approximation and inference methods for stochastic biochemical kinetics - a tutorial review

Stochastic fluctuations of molecule numbers are ubiquitous in biological systems. Important examples include gene expression and enzymatic processes in living cells. Such systems are typically modelled as chemical reaction networks whose dynamics are governed by the Chemical Master Equation. Despite its simple structure, no analytic solutions to the Chemical Master Equation are known for most systems. Moreover, stochastic simulations are computationally expensive, making systematic analysis and statistical inference a challenging task. Consequently, significant effort has been spent in recent decades on the development of efficient approximation and inference methods. This article gives an introduction to basic modelling concepts as well as an overview of state of the art methods. First, we motivate and introduce deterministic and stochastic methods for modelling chemical networks, and give an overview of simulation and exact solution methods. Next, we discuss several approximation methods, including the chemical Langevin equation, the system size expansion, moment closure approximations, time-scale separation approximations and hybrid methods. We discuss their various properties and review recent advances and remaining challenges for these methods. We present a comparison of several of these methods by means of a numerical case study and highlight some of their respective advantages and disadvantages. Finally, we discuss the problem of inference from experimental data in the Bayesian framework and review recent methods developed the literature. In summary, this review gives a self-contained introduction to modelling, approximations and inference methods for stochastic chemical kinetics.


Introduction
Understanding the functioning of living cells and biological organisms at the system level has gained increasing attention in recent years and defines a key research programme for the next decades. Experimental techniques are developing at breathtaking speed producing a wealth of data at finer and finer resolutions. However, such experimental data does not by itself reveal the function of such biological systems. The underlying processes typically involve large numbers of interacting components giving rise to highly complex behaviour [1]. Moreover, experimental data are generally corrupted by measurement noise and incomplete, thus posing the mathematical and statistical challenge to infer the relevant biological information from such measurements.
We focus here on mathematical and statistical modelling of chemical reaction networks in biological systems in which random fluctuations of molecule numbers play an important role. Recent experiments have shown this to be the case in many biological processes, gene expression being a prominent example [2]. Such random fluctuations or stochastic effects in gene expression have been found to lead to dramatically differing behaviours of genetically identical cells. From a modelling point of view, stochastic systems are considerably harder to analyse than their deterministic counterparts.
The rate equations give a macroscopic, deterministic description of the dynamics of chemical reaction networks. They consist of a set of ordinary differential equations governing the dynamics of the mean concentrations of the species in a system, thereby ignoring stochastic fluctuations. The rate equations have been successfully applied to various problems and have the advantage of being relatively straightforward to analyse. They typically give an accurate description of systems with large numbers of molecules, which is generally the case in in vitro experiments. However, the rate equations do no longer provide a valid description in cases where the effects of stochastic fluctuations become significant. This is typically the case when some species in a system occur at low molecule numbers, a common feature of chemical networks in cells. In this case, the chemical master equation (CME) constitutes the accepted probabilistic description of the resulting stochastic process [3]. In this framework, the state of a system is given by the molecule numbers of different species, and the CME governs the single-time probability of being in a certain state at a given time.
Despite its simple structure, no analytic solutions to the CME are known for all but the simplest systems. It is however possible to simulate exact sample paths of the underlying stochastic process, by means of the stochastic simulation algorithm [4,5]. The latter allows us to draw exact samples from the process. However, since the stochastic simulation algorithm simulates each and every chemical reaction event in a system, it is computationally expensive and quickly becomes infeasible for larger systems. Accordingly, significant effort has been spent in recent decades on the development of approximation methods of the CME, and a large variety of different methods has emerged. Similarly, due to the growing amount of experimental data becoming available, a considerable amount of work has been devoted to the development of statistical inference methods for such data, i.e. methods that allow to calibrate a model to observational data.
There exists a significant amount of literature on modelling of stochastic chemical kinetics. However, most reviews are mainly concerned with simulation based methods [6][7][8][9][10][11][12]. Others are rather technical and require a high level of pre-existing and/or mathematical knowledge [13][14][15][16][17]. None of these references gives a thorough introduction into modelling of stochastic chemical kinetics in the biological context that is accessible for non-experts, neither do they give an overview or comparison of state of the art approximation and inference methods.
With this article, we aim at filling these gaps in several ways: 1. We give a self-contained introduction to deterministic and stochastic modelling techniques for chemical reaction networks. First, we introduce and give a historic motivation to the deterministic rate equations. Next, we derive and give a detailed discussion of the CME. We review exact solution methods for the CME as well as simulation methods. 2. We give a detailed derivation and discussion of the following approximation methods of the CME: the chemical Langevin equation and its associated chemical Fokker-Planck equation; the system size expansion; moment closure approximations. Moreover, we give an introduction and overview of other types of approximations, including time-scale separation based methods and hybrid approximations. 3. We perform a numerical case study comparing the various approximation methods mentioned before. 4. We give an introduction to inference methods in the Bayesian framework and review existing methods from the literature.
The presentation is written to be accessible for non-experts that are new to the field of stochastic modelling. Even though this review is motivated by stochastic effects in systems biology, it is important to stress that many systems in other scientific fields are frequently modelled by means of Master Equations. The methods discussed in this article can therefore readily be applied to such systems. Examples include ecology [18][19][20], epidemiology [21][22][23][24][25][26], social sciences [27][28][29] and neuroscience [30][31][32][33][34][35][36].
This article is structured as follows. We start by discussing the importance of stochasticity in biological systems in section 2. We describe the underlying mechanisms giving rise to stochastic fluctuations in cells (section 2.1) and discuss key experimental studies that measured such fluctuations (section 2.2). Next in section 3 we discuss deterministic and stochastic modelling methods for chemical kinetics. We start by introducing the concept of chemical reaction networks and deterministic descriptions in terms of macroscopic rate equations in section 3.1. Next, we introduce stochastic modelling techniques in terms of the CME and stochastic simulation algorithms in sections 3.2 and 3.3, respectively. We discuss analytic solution methods for certain classes of reaction systems in section 3.4. Section 4 is devoted to approx imation methods of the CME. We give detailed introductions and discussions of the chemical Langevin equation (section 4.1), the system size expansion (section 4.2) and moment closure approximations (section 4.3). Next, we discuss how approximate moment values obtained from these methods can be used to construct probability distributions by means of the maximum entropy principle in section 4.4. In section 4.5 we review existing software packages implementing the approximation methods discussed in sections 4.1-4.3. We give an introduction to other approximation methods in section 4.6, including the finite state projection algorithm, time-scale separation based approximations and hybrid methods.
In section 5 we perform a numerical case study and compare the chemical Langevin equation, the system size expansion and moment closure approximations. Section 6 deals with the problem of inference for CME type systems from observational data. We introduce the Bayesian approach to this problem and review existing methods from the literature. Finally, we summarise and conclude in section 7.

Stochasticity in biological systems
Stochastic effects play an important role in many chemical reaction networks in living cells. Examples are enzymatic/catalytic processes, transduction of external signals to the interior of cells or the process of gene expression, to name just a few. Here we discuss different sources of stochasticity and illustrate their emergence in biochemical networks using gene expression as an example in section 2.1. Subsequently in section 2.2 we explain how the different sources of stochasticity can be measured and distinguished experimentally and highlight the importance of stochasticity for living cells.

Emergence of stochasticity
The term 'gene expression' denotes the process of synthesis of functional gene products such as proteins. The mechanism is illustrated in figure 1 for prokaryotic cells, namely cells lacking in-membrane-bound organelles (e.g. nucleus, mitochondria), such as bacteria. The process includes two main steps: transcription during which mRNA molecules are produced, and translation during which protein molecules are synthesised [37]. The transcription process involves the enzyme RNA polymerase. For the mechanism to initiate, an RNA polymerase enzyme must bind to the initial sequence of a gene. It then slides along the gene and produces an mRNA molecule that reproduces the DNA code of the gene. The RNA polymerase molecules move around randomly in the cell, a process which can be approximately described as Brownian motion [15]. This means that the RNA polymerase binding to the gene is a stochastic event that happens randomly in time. As it turns out, not only the binding of the RNA polymerase to the gene, but also the sliding along the gene happens stochastically. Therefore, the production of mRNA molecules is a stochastic process.
The production of protein molecules from mRNA during translation is conducted by ribosomes, which are RNA and protein complexes. The ribosomes and mRNA diffuse in the cell and hence meet randomly before translation can occur. Translation is thus also a stochastic process. Similarly, the degradation of mRNA molecules and proteins is conducted by certain enzymes and hence happens stochastically.
The transcriptional regulation of gene expression is frequently modulated by regulatory proteins named transcription factors. Transcription factors are gene specific and bind to the promoter region, which is located upstream of the gene's encoding region. Upon binding, transcription factors tune the affinity of the RNA polymerase molecules for the promoter, thereby modulating the rate at which transcription initiation events occur and hence the overall transcription rate. A transcription factor can either increase or decrease the binding rate of RNA polymerase and thus either enhance or suppress gene expression.
In eukaryotes, i.e. cells possessing in-membrane-bound organelles such as a nucleus, gene expression happens in a similar but somewhat more complicated way involving more steps. For instance, genes are generally located in the nucleus and transcription hence happens in the nucleus. Processes such as mRNA maturation and diffusion of mRNA to the cytoplasm need to happen before the mature mRNA can be translated into protein.
The inherent stochasticity of chemical processes leads to fluctuations of molecule numbers in time. As an example figure 2(a) shows fluorescence intensity measurements of a fluorescent protein in individual Escherichia coli cells, indicating strong fluctuations of protein numbers in time. Temporal fluctuations that are due to the stochasticity of chemical processes is what we call intrinsic noise [2]. Differences in molecule numbers of a certain species, say proteins, between different cells can originate from this type of fluctuations. However, such differences  The time trajectories correspond to the difference of the measured intensities in individual cells and the population average. We observe that the protein numbers fluctuate strongly over time [39]. From [2]. Reprinted with permission from AAAS. (b) Fluorescence intensities for the dual reporter technique measured in [2] for the two E. coli strains M22 and D22. Permission kindly granted by the publisher. Each triangle corresponds to the measured fluorescence intensity in one cell. Variations along and perpendicular to the diagonal correspond to extrinsic and intrinsic noise, respectively. We observe that in this system both intrinsic and extrinsic noise contribute significantly to the overall fluctuations. Reprinted by permission from Macmillan Publishers Ltd: Nature [39], copyright 2006. can also stem from other effects, such as physiological differences between cells or differing environmental conditions. For example, the numbers of RNA polymerase or ribosomes may differ between different cells, or different cells may be exposed to varying nutrient concentrations due to environmental fluctuations. Such differences that are not due to the stochasticity of chemical reactions are referred to as extrinsic noise [38].

Experimental evidence
As explained above, stochastic fluctuations are inherent to biochemical processes such as gene expression. The question hence arises: what effects do these fluctuations have on the functioning of cells and which processes dominate the emergence of these fluctuations? To answer such questions it is crucial to be able to empirically distinguish and quantify intrinsic and extrinsic noise.
One of the first experiments that aimed at explicitly separating the effects of intrinsic and extrinsic noise on gene expression was conducted by Elowitz et al in 2002 [2] by means of the dual reporter technique. In this study, the authors integrated two distinguishable but identically regulated fluorescent reporter genes into the chromosome of E. coli cells, one expressing cyan fluorescent protein (CFP) and one yellow fluorescent protein (YFP). The two genes were controlled by identical promoters, which means that they experienced the same external effects, i.e. the same extrinsic noise. Assuming that the two promoters are independent, fluctuations of the two proteins originating from the stochasticity of chemical reactions, i.e. intrinsic noise, should be completely uncorrelated. On the other hand, extrinsic effects should influence both expression systems equally, and the corresponding fluctuations should therefore be strongly correlated. This can be visualised as in figure 2(b) which is taken from [2]. The figure shows the measured intensities of individual cells for two different strains, with the YFP intensity on the y-axis and the CFP on the x-axis. The width of the cloud along the diagonal corresponds to correlated fluctuations and hence extrinsic noise, while the width orthogonal to the diagonal corresponds to uncorrelated fluctuations, i.e. intrinsic noise. The figure indicates that for this particular system both intrinsic and extrinsic noise significantly contribute to the overall noise.
Many other experimental studies have been conducted investigating the origins and role of intrinsic and extrinsic noise in cells (for example [39][40][41][42][43][44]). For an overview of experimental techniques see [45,46]. In [43], for example, the authors measured the expression levels of 43 genes in E. coli under various different environmental conditions. They observed that the variance of protein numbers scales roughly like the mean protein numbers and that for intermediate abundance proteins the intrinsic fluctuations are comparable or larger than extrinsic fluctuations. In another large scale study on E. coli in [44], the authors found that fluctuations are dominated by extrinsic noise at high expression levels.
The reported fluctuations can have significant influence on the functional behaviour of cells [1]. A particularly important example are stochastic cell fate decisions [47], where genetically identical cells under the same environmental conditions stochastically differentiate into functionally different states. Such differing cell fate decisions of genetically identical cells are believed to be beneficial for a population of cells experiencing fluctuating environments [47].
Most of the mentioned experimental studies measure certain statistics of gene expression systems, such as the equilibrium mean and variance of protein numbers. While experimental methods such as the dual reporter technique allow to estimate the magnitude of intrinsic and extrinsic noise, they do not explain where the stochasticity originates from and what the underlying processes are. In order to answer such questions, mathematical methods are needed that allow researchers to model the dynamics and stochasticity of such biochemical reaction systems.

Stochastic chemical kinetics
After having given a biological motivation for the importance of stochasticity in chemical reaction networks in the previous section, we consider next the question of how such systems can be described mathematically. First, we introduce the concept of chemical reaction networks in section 3.1 and describe the classical mathematical description of such systems in terms of deterministic rate equations. Next, we discuss the stochastic description of such systems in section 3.2. We derive the chemical master equation and discuss its validity and solutions. Next, we discuss the stochastic simulation algorithm in section 3.3. Finally, we review methods for exact analytic solutions of the CME for certain classes of systems in section 3.4 and give an overview of available analytic solutions of specific systems that do not fall under any of these categories.

Chemical reaction networks and deterministic rate equations
Biological processes such as gene expression generally consist of complicated mechanisms involving several different types of molecules and physical operations. For a mathematical description of certain processes, one typically does not model all these mechanisms explicitly, but rather replaces them by an effective single chemical reaction event. In the context of gene expression, for example, transcription or translation may be modelled as single chemical reactions. A finite set of chemical species that interact via a set of such chemical reactions constitutes what we call a chemical reaction network. Given a set of chemical species X i , i = 1, ..., N, we define R chemical reactions by the notation where the stoichiometric coefficients s ir and ′ s ir are non-negative integer numbers denoting numbers of reactant and product molecules, respectively. We say that the rth reaction is , i.e. if it involves m reactant molecules. We further call a reaction 'unimolecular' if m = 1, 'bimolecular' if m = 2 and a system 'linear' if ⩽ m 1 for all reactions occurring in the system. We further call a system 'open' if it contains a chemical process that generates molecules, and 'closed' otherwise. The quantity k r in equation (1) is called the reaction rate constant of the rth reaction.
Classically, the dynamics of a chemical reaction system as in equation (1) has been modelled by the law of mass action which was developed by Guldberg and Waage in its first version in the 1860s in the context of macroscopic in vitro experiments, i.e. macroscopic amounts of chemical molecules in solutions [48][49][50][51]. The law of mass action states that the rate of a reaction is proportional to the product of the concentrations of reactant molecules. Specifically, if we define , where φ i denotes the concentration of species X i , the rate g r of a reaction as in equation (1) is given by We call g r the macroscopic rate function of the rth reaction. Let us further define the stoichiometric matrix S as The entry S ir corresponds to the net change of X i molecules when the rth reaction occurs. Consequently, ( ) φ S g ir r is the rate at which the concentration φ i of species X i changes due to the rth reaction. Summing up the contributions from all reactions in a given system one obtains [52] ( ) This set of ordinary differential equations is called the rate equations. The law of mass action and the associated rate equations assume continuous concentrations, which means that they ignore the discreteness of molecule numbers. Moreover, they constitute a deterministic method, i.e. for a fixed initial condition the state of the system is exactly determined for all times. However, in the previous section we have seen that chemical reactions occur stochastically leading to fluctuations of molecule numbers in time. The reason why the rate equations give an accurate description in experiments such as the ones of Guldberg and Waage is that they studied chemical systems in solutions which typically contain a large number of substrate molecules. For large molecule numbers, it has been found experimentally that the relative fluctuations, i.e. the standard deviation divided by the mean value of molecule concentrations, scales like the inverse square root of the mean concentration and hence becomes small for systems with large molecule numbers [43]. One can therefore expect the rate equations to be an accurate description whenever a system has large molecule numbers for all species.
Example. As an example, consider the gene system in figure 3. This system can be viewed as a simplified version of the gene expression process described in figure 1: we replace the process of transcription (gene produces an mRNA molecule) and translation (mRNA molecule produces a protein) by a single effective reaction in which the gene directly produces the protein.
The example in figure 3 depicts negative autoregulation. The latter is one of the most common regulatory mechanisms: the gene's product downregulates its own expression. In the bound state G off , the gene does not produce any protein. The protein thus suppresses its own production, which means the system is an example of a negative feedback loop. The corresponding reactions in the notation of equation (1) read  Illustration of a gene expression system with negative feedback loop. When the gene is in the 'on' state it produces proteins of a certain type. The protein can decay or bind to the gene's promoter. In the bound state the gene is in the 'off' state, i.e. no protein is produced. The corresponding reactions are given in equation (5). The protein hence suppresses its own production.
where we call the gene in the bound and unbound state G on and G off , respectively and the protein P. In our nomenclature, → + G P G on off is a second-order or bimolecular reaction, while the other three reactions are first order or linear. By ' → ∅ P ' we indicate that P leaves the system under consideration. This could for example mean that P becomes converted into different types of chemical species or degraded in its constitutive elements that are not included in our model.
Let us order reactions according to the rate constants in equation (5), and species as G P , on and G off , respectively. Consider the stoichiometric coefficients s ir and ′ s ir defined in equation (1), which correspond to the number of reactant and product molecules, respectively, of species X i in the r th reaction. For the reaction → + G G P on on , for example, we have = = s s 1, 0 11 21 and s 31 = 0, since there is one G on on the lhs of the reaction but no P or G off , and = = ′ ′ s s 1, 1 11 21 and = ′ s 0 31 since there is one G on , one P and no G off molecules on the rhs. Proceeding similarly for the other reactions in equation (5), we find Accordingly, the stoichiometric matrix = − ′ S s s reads 1 2 and φ 3 denote the concentrations of G P , on and G off , respectively. The macroscopic rate vector , with the g r defined in equation (2), is obtained by using s in equation (6) and reads .
Using equations (4), (7) and (8) it is easy to write down the corresponding rate equations. However, note that the system has a conservation law in particle numbers which we can use to find a simplified description by reducing the number of variables: the total number of genes in the 'on' state and genes in the 'off' state is constant, i.e. φ φ + = ≡ g const. 1 3 0 . We can thus reduce the system to a two species system by using φ φ = − g 3 0 1 . The matrices s and ′ s for the reduced system are obtained from equation (6) by dropping the last row, and the stoichiometric matrix and propensity vector of the reduced system read accordingly Using equation (4) we hence obtain the rate equations The gene system in figure 3 will be used to showcase several different methods in this article, and we will throughout use the reduced system.

Stochastic methods
The rate equations discussed in the previous section constitute a deterministic description of chemical reaction systems in terms of continuous molecule concentrations. As we have seen in section 2, however, fluctuations become important in many biological processes such as gene expression.
The reason for this is that often some species occur at low molecule counts. A more general, stochastic description would keep the discrete nature of the X i molecules in equation (1) intact. One general approach is to explicitly model the spatial positions of molecules and to model their movement as Brownian motion, with chemical reactions happening stochastically under certain rules. Exact analytic results for such systems are generally not known. Simulations, on the other hand, have to keep track of every single particle and quickly become computationally unfeasible. However, under certain conditions, which we discuss in the following section, simplified descriptions can be employed making spatial models and the simulation of single particles unnecessary.

The chemical master equation.
We now consider a chemical reaction system as in equation (1) in a closed compartment of volume Ω. We seek a stochastic description of the system under well-mixed and dilute conditions. By 'well-mixed' we mean that the diffusion of particles in the compartment constitutes the fastest time scale of the system, in the sense that the expected distance travelled by each particle between successive reactive collisions is much larger then the length scale of the compartment. This implies that the spatial positions of molecules can be ignored and the dynamics of the system only depends on the total molecule numbers. By 'dilute' we mean that the combined volume of all the considered molecules is much smaller than the total volume, which means that the molecules can be considered as point particles.
If these two conditions are met, it can be shown [3] that the state of the system at any time is fully determined by the state vector where n i is the molecule number of species X i in the compartment. In particular, the spatial locations and diffusion of molecules does not have to be modelled, and the system corresponds to a continuous-time Markov jump process. It can further be shown that the probability for the rth reaction to happen in an infinitesimal time step . The number of pairs with one A and one B molecule is n n A B , where n A and n B are the molecule numbers of A and B, respectively. The corresponding propensity function is hence given by /Ω k n n r A B . The scaling with the volume Ω stems from the fact that the probability for two molecules to collide is proportional to /Ω 1 . Generalising these arguments to reactions as in equation (1) (12) Propensity functions of this form are called mass-action kinetics type [14]. Although equation (2) is often, due to historical reasons, stated to be the law of mass action, from a microscopic perspective it is more accurate to state equation (12) as the law of mass action. Equation (2) can be viewed as the macroscopic version of (12) obtained in the limit of large molecule numbers and small fluctuations. Specifically, for reactions of up to order two equation (12) becomes Let us now consider how we can mathematically describe the dynamics of such a system. To this end, consider the probability distribution ( ) | P t t n n , , 0 0 for the system to be in state n at time t given that it was in state n 0 at time t 0 . We will use the shorthand ( ) ( ) = | P t P t t n n n , , , 0 0 in the following and implicitly assume conditioning on an initial state. The probability ( ) + P t t n, d after an infinitesimal time step t d is given by ( ) P t n, plus the probability to transition into state n from a different state * n minus the probability to leave state n, which leads us to where we used the fact that probability of the rth reaction to happen in an infinitesimal time interval t d is given by ( ) f t n d r and where S r is the rth column of the stoichiometric matrix S. Subtracting ( ) P t n, , dividing by t d and taking the limit → t d 0 gives the chemical master equation (CME) [3], Since n is a discrete-valued vector, equation (15) is a coupled system of linear ordinary differential equations. Note that this system is infinite whenever n is unbounded. Despite its simple structure, there are generally no analytic solutions known to the CME. We call a distribution ( ) P t n, a steady-state solution of the CME if it solves the CME and fulfils ( ) ∂ = P t n, 0 t . In the context of chemical reactions, one of the first applications of the CME was done by Delbrück in 1940 for describing the dynamics of an autocatalytic reaction [53]. It has later been applied to linear reactions by Bartholomay [54] and McQuarrie [55], and bimolecular reactions by Ishida [56] and McQuarrie et al [57]. Gillespie derived the CME from molecular physics of a dilute and well-mixed gas in 1992 [3], which he extended to liquids in 2009 [58].
Example. Consider again the gene system in figure 3 with reactions given in equation (5). The corresponding stoichiometric matrix S is given in equation (10). Let n 1 and n 2 be the molecule numbers of the gene in the unbound state G on and the protein P, respectively. Let us assume that there is only one gene in the system, which implies that the number of bound genes G off is 1 − n 1 . Using equation (12) we find the propensity vector ( ) ( ( ) ) = Ω − k n k n n k n k n f n , , 1 , .
The corresponding CME becomes (see equation (15)) P n n t k n P n n t k n n P n n t k n P n n t k n P n n t k n k n n k n k n P n n t , , , Despite having a relatively simple system here with effectively only two species, no timedependent solution for its CME in equation (17) is known to our knowledge. A solution in steady state has been derived in [59], but for most other systems not even a steady state solution is available. Therefore, one generally needs to rely on stochastic simulations or approximations to study the behaviour of such systems.

Moment equations.
Suppose we are not interested in the whole distribution solution of the CME but only in its first few moments, say the mean and variance. Starting from the CME in equation (15) Here, ⟨ ⟩ ⋅ denotes the expectation with respect to the solution ( ) P t n, of the CME in equation (15). For moments of up to order two equation (18) becomes t i j r R jr i r ir r j i r jr r 1 (20) We see that if all ( ) f n r are zeroth or first-order polynomials in n, i.e. the system is linear without any bimolecular or higher order reactions, the equation of a moment of order m depends only on moments of order m or lower, i.e. the equations are not coupled to higher order equations. The equations up to a certain order hence constitute a finite set of linear ordinary differential equations which can be readily solved numerically or by matrix exponentiation as will be described in the context of the CME in section 3.4. Note that no approx imation has been made here which means that the equations describe the exact moments of the CME's solution. This means that for linear systems, the exact moments up to a finite order of the process can be obtained by solving a finite set of linear ordinary differential equations.
If the systems is non-linear, i.e. contains bimolecular or higher order reactions, the equation of a certain moment depends on higher order moments. This means that the moment equations of different orders are coupled to each other, leading to an infinite hierarchy of coupled equations. This can obviously not be solved directly but gives the basis of approximation methods such as moment closure approximations which we will discuss in section 4.3.
Example. Let us again consider the gene system in figure 3 with reactions given in equation (5). The corresponding stoichiometric matrix and propensity vector are given in equation (10). Using these in equations (19) and (20) one obtains where we introduced the shorthand Note that the equations for the first order moments in equations (21) and (22) depend on the second moment y 1,2 , and that the equations for the second moments (23)-(25) depend on the third moments y 1,1,2 and y 1,2,2 . Similarly, it is easy to see that moment equations of any order depend on higher order moments, which means that we have an infinite system of coupled equations. Note that all terms in equations (21)-(25) depending on higher order moments are proportional to the rate constant k 2 of the bimolecular reaction in equation (5), illustrating that the moment equations decouple in the absence of bimolecular (and higher order) reactions. This could be achieved here by setting k 2 = 0 for which the moment equations would decouple and could thus be solved numerically.

Non-mass-action propensity functions.
So far we only considered propensity functions of mass-action type defined in equation (12). However, in the literature propensity functions that are not of mass-action kinetics type are frequently used, such as Michaelis-Menten or Hill functions to model the dependence of the production rate of a product on the substrate concentration in an enzymatic catalysis or the dependence of the expression level of a gene on its transcription factor. Such propensity functions typically arise in reduced models where an effective reaction replaces several microscopic reactions. For a gene that becomes regulated by a transcription factor, for example, the binding of the transcription factor to the promoter of the gene is not modelled explicitly, but the effect of its concentration included in the modified propensity function of the expression reaction. Such non-mass-action type reactions should thus normally be seen as an effective reaction replacing a set of mass-action type reactions. A possible reduced model of the gene expression system in figure 3, for example, eliminates the gene from the system and combines the binding, unbinding and protein production reactions into one effective reaction → ∅ P with a Michaelis-Menten type propensity function /( ) + k k n k n k P P 1 3 2 3 , which means that the reaction producing protein now depends on the protein number n P . Such reductions often emerge from separations of time-scales in a system and become exact in certain limits, see section 4.6 for more details.
Independently of the origin of a non-mass-action propensity function, as long as it is interpreted as the firing rate of the corresponding reaction, the CME description remains valid. Unless explicitly stated otherwise, we will assume mass-action kinetics in this article.

Stochastic simulations
As mentioned above, there are generally no analytic solutions known to the CME. However, it is possible to directly simulate the underlying process. The stochastic simulation algorithm (SSA) is a popular Monte Carlo method that allows one to simulate exact sample paths of the stochastic process described by the CME.
The stochastic simulation algorithm was first proposed in the context of chemical kinetics by Gillespie [4,5], and several variants have been proposed in the literature, see [9][10][11] for reviews. The basic idea is to simulate reaction events explicitly in time and to update the time and state vector accordingly. The stochastic process described by the CME is a continuoustime Markov jump process, which has the important property that waiting times, i.e. time intervals between successive reaction events, are exponentially distributed [15]. Since it is easy to sample from exponential distributions, it is straightforward to simulate the occurrences of chemical reactions.
One example is the so-called direct method [4], which samples the time step τ for the next reaction to happen and subsequently which of the different reactions occurs. Specifically, let ( ) τ| p t n, be the probability for the next reaction to happen in an infinitesimal time interval t d around τ + t, given that the state of the system is n at time t, and ( ) | p r t n, the probability that the next reaction is a reaction of type r. Using that ( ) f t n d r is the probability for the rth reaction to happen in t d , it can be shown that [4] ( Samples from equations (27) and (28) can be respectively obtained as where u 1 and u 2 are uniform random numbers between 0 and 1. The direct SSA method iteratively updates the state vector and time of the process by first sampling the time point of the next reaction event according to equation (29) and subsequently sampling which reaction happens according to equation (30).
Unfortunately, the applicability of the SSA is severely limited due to its computational cost. Since each and every reaction is simulated explicitly, the SSA becomes computationally expensive even for systems with few species. This is particularly the case if the molecule numbers have large fluctuations or if many reactions happen per unit time. In the first case a large number of samples have to be simulated to obtain statistically accurate results, whereas in the second case single simulations become expensive since the time between reaction events becomes small.

Extrinsic noise.
So far we have only discussed stochastic simulations of systems with propensity functions that depend solely on the system's molecule numbers and the system's volume. This implies that the propensity functions are constant between reactions events, leading to inter-reaction intervals being exponentially distributed. However, as already mentioned in section 3.2.1, it is frequently of interest to consider stochastically varying propensity functions. This could for example represent the effect of a randomly fluctuating environment on a cell.
To model the effect of extrinsic noise, one typically includes a stochastic variable in one or several propensity functions, whose dynamics may for example be governed by a stochastic differential equation. This means that the corresponding propensities become explicit (stochastic) functions of time. While the CME in equation (15) is still valid in this case, stochastic simulation algorithms discussed above are not since the inter-reaction times are now not exponentially distributed anymore.
The simplest approach is to assume the propensity functions to be constant between reaction events, which allows one to use the standard simulation algorithm [4,60]. However, this is of course an approximation and can lead to inaccuracies if the extrinsic process (and hence the corresponding propensity functions) fluctuates strongly between consecutive reaction times. The most straightforward way to simulate such processes exactly is to integrate the propensity functions step-wise over time until a certain target value is reached [61][62][63]. Apart from numerical integration errors this class of algorithms is exact. However, due to the numerical integration it becomes computationally highly expensive. More recently an alternative exact method has been proposed [64]. This method introduces an additional reaction channel in such a way that the total propensity (i.e. λ in equation (27)) is constant between successive reactions. This in turn allows one to use a standard SSA algorithm on the augmented system. It is hence generally more efficient than the aforementioned integral methods.

Exact results
As pointed out before, for most systems no analytic solutions of the CME are known. However, for certain classes of systems analytic solutions do exist. For some classes the general timedependent case can be solved while others can be solved only in steady state. Moreover, analytic solutions have been derived for various simple example systems that do not fall under any of these classes. We give here first an overview of general classes of systems that can be solved analytically, and subsequently list some additional systems for which analytic solutions are known.
In section 3.2 we have noted that in this case the CME in equation (15) is a finite system of coupled ordinary differential equations. If , the CME can be written in the matrix form as where E is a × M M matrix whose elements can be easily derived for a given system using equation (15). The solution of equation (31) is simply given by Therefore, for systems with finite state space a solution of the CME for all times can be in principle computed using equation (32). However, even for systems with finite state space, the dimension of equation (32) is often quite large in practice, making matrix exponentiation computationally expensive. Efficient numerical methods for matrix exponentiation have been developed in recent years [65,66], but for many chemical reaction systems of practical interest it remains an intractable task. For systems with large or infinite state space equation (32) is hence not of direct practical use. It forms the basis for certain approximation methods, however, see section 4.6.

Linear systems.
In section 3.2.2 we found that the moment equations decouple from higher order equations for linear reaction systems, i.e. systems without bimolecular or higher order reactions. The moment equations up to a finite order can hence be directly solved numerically for such systems. For non-linear systems, in contrast, this is not possible since the moment equations couple to higher orders. Similarly, exact solutions for the whole distribution solution of the CME are easier to obtain for linear systems. Often this can be done by means of the generating function method which transforms the CME into a partial differ ential equation [15]. The form of the latter becomes particularly simple for linear systems and can often be solved analytically.
In 1966 Darvey and Staff used the generating function approach to derive a multinomial solution for ( ) P t n, for multinomial initial conditions and for systems with for all reactions r, i.e. exactly one reactant and exactly one product molecule for all reactions [67]. Note that these are closed systems. Later Gardiner used the so-called 'Poisson representation' to show that for systems with for all reactions r and Poisson product initial conditions, the solution of the CME is a Poisson product for all times [68]. The auto-and cross-correlation functions for such systems have more recently been computed in [69]. These results have been generalised to systems with for all reactions r, with arbitrary initial conditions in [70]. In this case the solution can be written as a convolution of multinomial and product Poisson distributions. In all these cases, the solutions can be constructed from the solution of the deterministic rate equations discussed in section 3.1 and can hence be obtained efficiently by (numerically) solving a finite set of ordinary differential equations.
We would like to stress that the methods described here do not apply to all linear reactions, but only to those with less than two product molecules ( ). This excludes reactions of the type → + A A B which are common in biology, such as translation of a protein B from its mRNA, A.

Non-linear systems.
As discussed in the previous section, for linear systems where each reaction has at most one product molecule, an analytic solution of the CME can be derived. In contrast, for non-linear systems or linear systems including reactions with two product molecules no analytic solutions are known in general. However, for certain subclasses of systems it is possible to derive steady-state solutions, i.e. solutions satisfying ( ) ∂ = P t n, 0 t . One example are reversible systems that obey detailed balance, which is also called thermodynamic equilibrium. By 'reversible' we mean that each reaction possesses a corresponding reversing reaction. Note that the steady-state condition ( ) ∂ = P t n, 0 t in the CME in equation (15) merely means that the probability flow out of a state (second term in the CME) equals the probability flow into a state (first term in the CME). These probabilities consist of sums over all reactions. Detailed balance is a stronger condition, as it requires the flow of each single reaction in each state to be balanced by its corresponding reversible reaction. Specifically, let r and ′ r be the indices of two reactions that revert each other and let ( ) P n be a steady-state solution. The detailed balance condition reads [15] which needs to be fulfilled for all reactions and all n. Note that = − ′ S S r r . The left side of equation (33) is the flow out of state n due to the rth reaction, while the right side is the flow into state n due to the backward reaction ′ r . An analogue detailed balance condition can be formulated for the deterministic rate equations. It can be shown that for a reversible reaction system with mass-action kinetics, the stochastic system is in detailed balance if and only if the deterministic system is in detailed balance [71]. If the detailed balance condition is fulfilled, the solution of the CME is a product of Poisson distributions times a function accounting for conservation laws in molecule numbers, where the mean values of the Poisson distributions are simply given by the detailed balance solutions of the rate equations [15,72].
As an example, consider a closed, bimolecular system with reactions and n B be the molecule numbers of A and B respectively. The quantity is conserved. The corresponding steady-state solution of the CME in detailed balance reads where α φ = Ω 1 1 , φ 1 0 and φ 2 0 are the steady-state solutions of the rate equations of A and B, respectively, Ω is the system volume, and δ is the delta function accounting for the conservation law. These results have more recently been generalised to weakly reversible reaction systems [73]. A system is called weakly reversible when for each reaction the change in state space can be reversed by a chain of other reactions. For example, a system with the reactions → → → X X X X

Exact results for special cases.
Most reaction systems of practical interest are neither linear nor do they satisfy the detailed/complex balance conditions. However, the CME has been solved for several specific systems that do not fall under any of these categories. In particular for systems with only a single reaction it is often straightforward to obtain a time-dependent solution of the CME. Some examples are given below. Most of these solutions are obtained either by means of the aforementioned generating function method or the Poisson representation. While the former transforms the CME into a partial differential equation, the latter maps the discrete process underlying the CME to a continuous diffusion process governed by a Langevin equation defined in a complex-valued space [68]. We would like to stress that both approaches are exact. However the Poisson representation has not been used much in the literature, probably because of its relative complexity compared to other methods and hence we do not cover it in detail here. We refer the reader to the following references for specific applications [76][77][78][79][80][81][82]. For more details on the Poisson representation and the generating function approach we refer the reader to [15].

Single reaction systems.
Time-dependent solutions have been derived for Gene regulatory networks.
: a steady-state solution has been derived if there is one gene in the system [59]. Note that if we ignore the last two reactions this system corresponds to the gene system in figure 3.
• A system similar to the previous one but with geometrically distributed bursts of protein production has been solved in steady state in [84].

Enzyme systems.
• ↔ → + + E S C E P: A time-dependent solution has been derived for the case of a single enzyme in the system [85].
• Adding the reaction → + E P C to the previous system makes it reversible and a corresponding steady-state solution has been derived in [86].
• The system with added substrate input reaction → ∅ S has been solved in steady state in [87]. • A similar system but with multiple substrate species competing for the enzyme has been solved in steady state in [88].
Most of these examples only consider steady-state solutions and all of them consider systems with few species. The lack of solutions for more complicated systems makes the development of approximation methods necessary.

Approximation methods
As discusses in the previous section, analytic solutions of the CME are known only for very restrictive classes of systems and few simple special cases. For most systems of practical interest, no analytic solutions are known to date. Stochastic simulation algorithms introduced in section 3.3 allow to simulate exact samples of the underlying stochastic process, but they quickly become computationally infeasible for larger systems. For these reasons significant effort has been spent in the literature on the development of approximation methods. We give here an introduction to a wide range of such methods. First, we give a detailed introduction to a few approximation methods that can be applied to (almost) arbitrary systems without any pre-knowledge of the system necessary. We derive the approximations, give examples and discuss their properties. Subsequently, we give an overview of other types of approximation methods developed in the literature. The first method aims at approximating the solution of the CME, namely the chemical Fokker-Planck equation (CFPE) and the associated chemical Langevin equation (CLE), which we introduce in section 4.1. The CFPE/CLE define an approximating stochastic process. An alternative method for approximating the solution of the CME is given by the system size expansion, which constitutes a systematic expansion of the CME in the inverse volume size. The system size expansion includes the popular linear noise approximation (LNA) and also provide an efficient way of approximating the moments of a process. We discuss the system size expansion in section 4.2. Next, in section 4.3 we introduce a certain class of moment closure approximations, which approximate the moments of a process. In section 4.4 we show how such approximate moments can be used to construct distributions using the maximum entropy principle. Next, in section 4.5 we review software packages implementing the discussed approximation methods. Finally, in section 4.6 we give a brief overview of other approximation methods found in the literature. As we shall see, many of these methods use the CLE, the system size expansion or moment closure approximations as building blocks.

The chemical Langevin equation
The chemical Langevin equation (CLE) and the corresponding chemical Fokker-Planck equation (CFPE) constitute a popular diffusion approximation of the CME. Kramers and Moyal derived the CFPE by applying a Taylor expansion to the CME which upon truncation leads to a partial differential equation approximation of the CME [89,90]: Suppose we let the variables in the CME in equation (15) become continuous and let , where x i is the continuous variable denoting the molecule number of species X i . Performing a Taylor expansion to second order around x in the first term of the rhs of equation (15) gives (with n replaced by x) , .
Inserting this into the CME in equation (15), we see that the first term on the rhs of equation (35) cancels the last term on the rhs of equation (15) leading to the CFPE: where the drift vector A and diffusion matrix B are respectively given by Note that the drift vector ( ) A x and diffusion matrix ( ) B x do not depend on time. Note also that whereas the state variables denote discrete molecule numbers in the CME, they denote continuous real numbers in the CFPE. The CFPE in equation (36) is equivalent to the CLE which is an Ito stochastic differential equation [15]. W in equation (39) is a multi-dimensional Wiener process. It can be shown [15] that the distribution of a process described by equation (39) agrees exactly with the solution of the corresponding Fokker-Planck equation in (36). One can thus interpret the CLE in equation (39) as a generator of realisations of the stochastic process described by the corresponding CFPE. In this sense, the CLE and CFPE are considered to be equivalent to each other. An alternative derivation, which leads directly to the CLE in equation (39), was given by Gillespie in [91]. Generally there exist different choices for ( ) C x in equation (39) corresponding to different factorisations of the matrix ( ) B x ; these lead to as many different representations of the CLE. One possibility is given by This representation of the CLE is the one most commonly used in the literature [91,92].
Example. Let us come back to the gene expression system in figure 3 with reactions in equation (5) and consider the corresponding CFPE and CLE. Using the stoichiometric matrix in equation (10) and propensity vector in equation (16) we obtain for the drift vector and diffusion matrix defined in equation (37) and equation (38), respectively, where x 1 and x 2 are the (continuous) particle numbers of G on and P, respectively. To obtain the CLE in equation (39), we have to compute C, which is the square root of B and thus generally not uniquely defined. One possibility as given in equation (40) reads This gives rise to the CLE where the W i are independent Wiener processes. Note that it does not make a difference if one changes the signs in front of the square roots in equation (43), or equivalently in front of the noise terms in equation (44), as long as one does so simultaneously for each occurrence of a specific term, i.e. changes the sign of whole columns in equation (43). To see that such changes are equivalent note that the diffusion matrix = B CC T of the CFPE is invariant under such changes. This can also be seen directly from the CLE in equations (44) and (45) since the Wiener processes are symmetric.

Stochastic simulations.
As for the CME, there are no analytic solutions known for the CLE for most systems. However, the computational cost of CLE simulations scales with the number of species, rather than with the rate of reaction events as CME simulations. This means that CLE simulations are often more efficient than simulations of the CME, typically if the molecule counts in a system are not too small.
A popular method to simulate the CLE is the Euler-Maruyama algorithm which discretises time into intervals dt and simulates the process iteratively as [93] where ( ) N 0, 1 is a normal distribution with mean 0 and variance 1. The smaller dt, the better the true process is approximated by equation (46) but also the slower the algorithm becomes. The right choice of dt is therefore a tradeoff between accuracy and efficiency. It is important to point out that more efficient simulation methods exist, see for example [93][94][95][96].

Properties and recent developments.
The CLE is typically a good approximation whenever the molecule numbers of the system are not too small and is particularly useful if one is interested in time trajectories or distributions of a process. It can be shown that the differences between the CLE and the CME tend to zero in the limit of large molecule numbers [97]. In other words, the CLE becomes exact in the thermodynamic limit.
By multiplying the CFPE in equation (36) by … x x i l and integrating over the whole state space, one obtains ordinary differential equation for the moment ⟨ ⟩ … x x i l of the process described by the CLE. Importantly, it turns out that the equations for moments of up to order two are exactly the same as the corresponding equations derived from the CME. These are given in equations (19) and (20). Note however that since they are generally coupled to higher order moments for which the evolution equations derived from the CLE and CME do not agree, the first two moments (and higher order moments) of the CLE do not generally agree with the ones of the CME. However, since the moment equations decouple for linear systems as shown in section 3.2.2, we obtain the important result that the moments up to order two of the processes described by the CLE and CME agree exactly for linear reaction systems.
Note that the CME has a natural boundary at zero molecule numbers, i.e. for a sensible initial condition with zero probability to have a negative number of molecules, this probability remains zero for all times. Until recently it has not been clear how this boundary condition behaves when approximating the discrete process underlying the CME by a continuous process using the CLE. As it turns out, this boundary issue leads to an ill-definedness of the CLE due to occurrences of square roots of negative expressions in finite time with finite probability [99,100]. Since traditionally the domain of the CLE is (implicitly) assumed to be that of real numbers, the CLE is not well-defined in this case. More recently, it has been shown that this problem is independent of the chosen factorisation of the CFPE's noise matrix B (see equation (39)), which means that the CLE is not well-defined for real variables [87]. In this study it has been shown that the same can be expected for the majority of reaction systems.
Several modified versions of the CLE have been proposed that try to keep the state space real, for example [98,99]. However, these are ad hoc modifications and have been found to introduce high inaccuracies for some non-linear reaction systems [87]. Importantly, they also have been found to violate the CLE's exactness for the moments up to order two for linear systems. Alternatively, the ill-defined problem can be solved by extending the state space of the CLE to complex variables. It has been proven that this leads to real-valued moments, real-valued autocorrelation functions and real-valued power spectra, and to restore the CLE's exactness for the moments up to order two for linear systems [87]. This complex CLE was found to be highly accurate for some non-linear systems in comparison to other modified versions. For one example the results from [87] are shown in figure 4. The complex-valued CLE has other drawbacks, however. For example, it does not directly give approximations of the process or of distributions. Rather, the results have to be projected to real space.
In recent years the CLE has been frequently used for approximating the dynamics of intermediate or high abundance species in so-called hybrid methods, see section 4.6.4. In this case, the probability of negative concentrations becomes very small. If the partitioning into low and high abundance species is done adaptively, the problem of square roots of negative expressions is avoided completely. In this case the complex-valued CLE automatically reduces to the real-valued CLE. The CLE is particularly useful in simulation based hybrid methods since it gives an approximation of the whole process (rather than for example only the process's moments).
Recently, a tensor-based method has been proposed in [101] for the direct numerical solution of the CFPE. This is a promising approach since it avoids computationally expensive ensemble averaging and has been used for sensitivity and bifurcation analysis [101]. Recall that the CLE becomes exact in the limit of large system sizes. As one may expect, the CLE therefore generally captures the multimodality of the CME's solution if the multimodality persists for large volumes, since this is the case if the deterministic rate equations are multistable. Surprisingly, however, it has recently been found that the CLE is also able to reproduce noise-induced multimodality of the CME's solution, i.e. multimodality that occurs only if the system volume is decreased below a certain critical value [102]. However, it was found that this is not the case for all reaction systems [103]. In [104], the CLE as a diffusion approximation was applied to Petri nets.

The system size expansion
Suppose we are not interested in approximating the whole process but only in its distribution or first few moments. Running stochastic simulations of the CLE seem like unnecessary computational effort in this case and the question arises if there exist more efficient approximation methods to achieve this. We next discuss the system size expansion which aims at The results are normalised by the exact CME result, which means that the horizontal dashed line corresponds to the exact value. CLE-C corresponds to the complex-valued CLE, CLE-R to a real-valued implementation with rejecting boundaries at zero molecules, and CLE-DR to the real-valued version proposed in [98]. LNA and 2MA correspond to the linear noise approximation and the second-order normal moment closure, which we will introduce in sections 4.2 and 4.3, respectively. We observe that both real-valued implementations give large deviations from the CME result, significantly more than the LNA and 2MA. The complex CLE on the other hand is significantly more accurate than all the other methods. Reprinted from [87], with the permission of AIP Publishing.
approximating the distribution or the first few moments of a process. The system size expansion is a perturbative expansion of the CME in the inverse system size originally developed by van Kampen [14,105]. The idea is to separate concentrations into a deterministic part, given by the solution of the deterministic rate equations, and a part describing the fluctuations about the deterministic mean.
4.2.1. Derivation. The system size expansion splits the instantaneous particle numbers n i of the CME into a deterministic part and a fluctuating part as where Ω is the volume of the system, φ i is the solution of the deterministic rate equations in (4) and we introduced the new variables ε i representing fluctuations about the deterministic mean. Following [106], the system size expansion proceeds by transforming the master equation to these new variables ε i . To this end, we rewrite the CME in equation (15) as where we introduced the step operators − E i S ir , which are defined in terms of their action on a general function of the state space, ( ) The system size expansion now proceeds in three steps: 1. Time derivative. Due to the change of variables in equation (47) we need to transform the distribution ( ) P t n, into a distribution ( ) Π ε t , of the new variables. Note that the time derivate in the CME in equation (48) is taken at constant n, which implies . This in turn leads to 2.
Step operator. Using equation (47) in the definition of the step operator in equation (49), we obtain Expanding the rhs in powers of / Ω −1 2 we can expand the term of the CME involving the step operators as 3. Propensity functions. We assume that the propensity functions in the CME can be expanded as If we assume mass-action kinetics for which the propensity functions take the form in equation (12), the propensity functions are polynomials and it is easy to see that such an expansion always exists. The same is also true for many propensity functions that are not of mass-action kinetics type, such as Michaelis-Menten and Hill-type propensity functions. Using equation (47) and expanding the rhs of equation (53) in where we have identified is the macroscopic rate function of the rth reaction introduced in equation (2) in context of the deterministic rate equations.
Applying equations (50), (52) and (54) to equation (48) we obtain , , where we defined Here we use the shorthand / ∂ = ∂ ∂ε i i and assume summations over doubly occurring indices for notational simplicity. We note that in equation (55) also two terms of order / Ω 1 2 occur. However, these terms cancel each other because they just correspond to the deterministic rate equations. In equations (56)-(58) we defined Note that both the ( ) J i and ( ) D i matrices depend on the solution φ of the rate equations in (4) and are thus generally time dependent. Comparing with equation (4) we find that ( ) J 0 in equation (60) is just the Jacobian of the rate equations.

The linear noise approximation.
The popular linear noise approximation (LNA) [14,105] is obtained if we truncate the expansion in equation (55) to zeroth order, leading to Equation (61) is a Fokker-Planck equation with drift and diffusion linear and constant in ε , respectively, and hence has a multivariate normal solution under appropriate initial conditions. By multiplying equation (61) with ε i and ε ε i j and integrating over all ε , one obtains ordinary differential equations for the first and second-order moments, ⟨ ⟩ ε i and ⟨ ⟩ ε ε i j , respectively. By doing so one finds that if the mean is initially zero, , it remains zero for all times. This is normally the case since otherwise the initial condition of the rate equations would not agree with the initial mean value of the stochastic system. If we additionally assume deterministic initial conditions, we also have ⟨ ⟩| = = ε ε 0 . The solution of equation (61) is thus a multivariate normal distribution with zero mean. Since n and ε are related by the linear transformation given in equation (47), the distribution of n is also given by a multivariate normal distribution. The mean of the latter satisfies the rate equations in equation (4) and the covari- The LNA describes the lowest order fluctuations of the system size expansion about the deterministic mean and is valid in the limit of large volumes, i.e. the thermodynamic limit. An alternative derivation of the LNA is given in [15], where it is derived from the chemical Langevin equation given in equation (39).
Example. Let us consider the rate equations and LNA for the gene system in figure 3. The corresponding stoichiometric matrix, the macroscopic rate vector and the rate equations were already derived in section 3.1 and are given in equations (10) and (11), respectively. The LNA defined in equation (61) (59) and (60), respectively, for which we obtain Note that ( ) D 0 and ( ) J 0 are functions of the time-dependent solutions φ 1 and φ 2 of the rate equations in (11). The solution of the LNA is a normal distribution in ( ) = n n n , 1 2 . Its mean is obtained by (numerically) solving the rate equations in (11), and the covariance by subsequently solving equation (62) using equations (63) and (64).

Higher order corrections. Suppose we want to include higher orders in
/ Ω −1 2 in equation (55) beyond the LNA. In contrast to the LNA, now the PDE in equation (55) can generally not be solved analytically. However, we can derive ordinary differential equations for the moments of the system accurate to the corresponding order in / Ω −1 2 , as follows. Assuming an expansion of the distribution ( ) Π ε t , as ,, one obtains an expansion of the moments Multiplying equation (55) with … ε ε j k , using equations (65) and (66) and integrating over all ε , we obtain ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε Here, ( ↔ ↔ ) i j k denotes a sum over the previous term over all cyclic permutations of the indices i,j and k, and similarly for four indices.
holds for all times. The moments in molecule numbers n i can now be obtained from the moments of the ε i variables by using the ansatz in equation (47). For the mean and covariance this leads to To order Ω 0 , equation (76) corresponds to the macroscopic rate equations, while equation (77) to order Ω −1 corresponds to the LNA estimate for the covariance. Including terms of order Ω −1 in equation (76) gives the leading order corrections to the mean given by the rate equations.
To this order, the system size expansion equations have been called 'effective mesoscopic rate equations' in the literature [107]. The next leading order corrections to the covariance, which have been called 'inverse omega square' in the literature [106], are obtained by keeping terms of order Ω −2 in equation (77).

Properties and recent developments.
As long as one is only interested in the LNA or higher order corrections to the moments, rather than the distributions of higher order truncations, the system size expansion amounts to the solution of finite sets of ordinary differential equations. It is thus generally significantly more efficient than stochastic simulations of the CME or of the CLE.
Since the system size expansion is an expansion around the deterministic mean, it cannot be used for deterministically multistable systems, i.e. systems whose rate equations have two or more positive stable steady states, unless one is only interested in the short-time behaviour of a process (for example for the purpose of inference from time-series data, see section 6).
In many cases, the lowest order of the system size expansion, i.e. the LNA, already gives remarkably accurate results and has been used succesfully in various applications [109][110][111][112][113]. However, for some systems it gives rise to significant deviations from the exact result [114]. Higher order approximations beyond the LNA have frequently been found to be highly accurate for such systems. Examples include oscillations in networks of coupled autocatalytic reactions [115,116], predator-prey systems [117] and circadian oscillators [118].
We note that for truncations of higher orders than zero, the resulting PDEs involve higher order derivatives than two and hence have no probabilistic interpretation due to non positivedefinite solutions [119]. Moreover, the solutions to the PDEs are generally not longer known in this case. However, for systems with only one species a general solution to all orders has recently been derived [108]. A similar approach has been developed in [120] for discrete-time models in neuroscience. Note that the variables described by the system size expansion are typically assumed to be continuous. In [108] it has been found that if the support is assumed to be continuous, the first few leading orders of the system size expansion lead to oscillations in the tail of the Poisson distribution for a linear birth-death process. In contrast, if one assumes a discrete support, the distribution is captured accurately. In the same work, a modified system size expansion has been employed. Instead of expanding the variables of the CME around the deterministic mean value given by the rate equations, one expands the variables around the mean value given by the system size expansion itself to the considered order (which is not known a priori). The authors call this method 'renormalised approximation'. This approach has been found to give an improved approximation of the distribution for a non-linear process as shown in figure 5.

Moment closure approximations
Another popular class of methods that approximate the first few moments of a process are socalled moment closure approximations. We present a popular class of such methods here that are particularly easy to derive and implement.

General formulation.
In section 3.2.2 we showed that the CME gives rise to ordinary differential equations for the moments of a process. For linear reaction systems, these equations can be solved numerically up to a certain desired order. For non-linear systems, however, we found that the moment equations of a certain order couple to higher order equations, leading to an infinite hierarchy of equations which can hence not be solved directly. Moment closure approximations truncate this infinite set of equations at a certain order M in some approximate way.
One popular class of moment closure approximations close the moment equations by expressing all moments above a certain order M as functions of lower-order moments. One way to achieve this is for example by assuming the distribution of the system to have a particular functional form, for example a normal distribution. This decouples the equations of the moments up to order M from higher-order moments, and hence leads to a finite set of coupled ordinary differential equations which can then be solved (numerically). We refer to such a moment closure as 'Mth order moment closure'. Let denote raw moments, central moments and cumulants of order k, respectively. We call … y i i 1, , , and a 'mixed moment' otherwise, and similarly for central moments and cumulants. The function g(s) in equation (80) is the cumulant generating function and is defined as We note that all three types of moments are respectively invariant under permutations of their indices. Therefore, only one representative combination of each permutation class has to be considered. Taking this symmetry into account significantly reduces the number of variables and moment equations. Some popular moment closure methods can be defined as Note that the cumulants of order higher than two are zero for a normal distribution, hence the name 'normal moment closure'. We refer to the normal moment closure approximations to second and third order as '2MA' and '3MA', respectively. • 'Poisson moment closure' [124]: the cumulants of a one-dimensional Poisson distribution are all equal to the mean value. We assume here a multivariate distribution to be a product of uni-variate Poisson distributions. Accordingly, for the Poisson moment closure approximation of order M we set all diagonal cumulants to the corresponding mean and all mixed cumulants to zero, i.e.
and , , , for some 1, , , • 'Log-normal moment closure' [125]: let m and S be the mean vector and covariance matrix of a multi-dimensional normal random variable. Then the logarithm of this random variable has a multivariate log-normal distribution whose moments can be expressed in terms of m and S as [126] ⎜ ⎟ ⎛ ⎝ , and g m is the number of i j 's having the value m. This allows one to express m and S in terms of the first two moments y i and y i,j which then in turn allows one to express higher-order moments in terms of y i and y i,j , too.
• 'Central-moment-neglect moment closure' (also called 'low dispersion moment closure' in the literature) [127]: all central moments above order M are set to zero: • 'Derivative matching': the idea of this method is to express moments above order M by lower order moments in such a way that the time derivatives of the moments of the closed system approximate the time derivatives of the moments up to order M of the exact system at some initial time point. In [128] a method is derived that allows one to produce the corresponding expressions.
Each of the closure methods allows one to express all raw moments above a specified order M in terms of lower order moments and hence to close the moment equations. Note that third order central moments and third order cumulants are identical. Therefore, whenever the moment equation of up to order two do not depend on moments of order higher than three, the normal and central-moment-neglect moment closure are identical. Similarly, the normal and Poisson moment closure approximations can be equivalent for certain systems.

Example.
As an example, consider again the gene system in figure 3. The moment equations up to order two for this system are given in equations (21)- (25). To close these we need to express the third order moments y 1,1,2 and y 1,2,2 in terms of lower order moments. We do so by means of the normal moment closure defined in equation (82) to second order, i.e. we set the third order cumulants c 1,1,2 and c 1,2,2 to zero: Using these expressions in equations (21)-(25) the equations decouple from higher order moments. We give here the resulting equations in terms of central rather than raw moments (see equation (79)) Note that the equations do not depend on third or higher order moments and hence are closed and can be integrated.

Properties and recent developments.
Moment closures are a popular class of approximations of the moments of the CME with many useful properties. First of all, they are easy to derive and implement. Moreover, since they amount to solving a finite set of ordinary differential equations no ensemble averaging is needed, which means they are computationally significantly more efficient than stochastic simulations of the CME or the CLE, and comparable to the system size expansion. However, this computational efficiency comes at a cost: moment closure methods only give approximate moments and not an approximation of the process or distributions. One advantage of moment closures over the system size expansion is that they can be applied to deterministically multistable systems. However, care must be taken for such system since moment closures can lead to unphysical results, see below.
In contrast to the system size expansion, moment closures are not an expansion in any small parameter, and one can therefore generally not expect that increasing the closure order leads to a higher accuracy. Numerical case studies suggest that this is often indeed the case, however [129][130][131]. One explanation for this was given in [129] where it has been shown that increasing orders of the normal moment closure agree with increasing orders of the system size expansion in the limit of large volumes. These results can to some extent be generalised to other moment closure methods. For monostable systems, one may therefore expect higher order closures to become more accurate for large enough system sizes. For small system sizes on the other hand, the closures can generally not be expected to converge.
Moreover, since moment closure approximations are an ad hoc approximation, it is not even clear if they always give physically meaningful results. As it turns out, this is not always the case, but they sometimes give rise to unphysical behaviour, such as negative mean values, negative variances or negative higher central moments, as well as diverging trajectories [124,131]. In numerical case studies, it has been found for several non-linear reaction systems that the normal, log-normal, Poisson and central-moment neglect moment closure show such unphysical behaviour for system volumes below a certain critical volume [131,132]. An example of such unphysical behaviour is shown in figure 6. This may be expected, since smaller volumes correspond to stronger non-linearity and hence stronger fluctuations. However, it has been found for deterministically oscillatory and deterministically multistable systems, that the moment closure methods give rise to non-physical oscillations and non-physical multistability, respectively, for system volumes above a certain critical volume [131,132]. In [132], the normal moment closure was found to give physically meaningful results for larger ranges of system volumes than the other three methods. In terms of accuracy, the four methods performed similar to each other. We emphasise that these were numerical case studies, and it remains open to what extent these results can be generalised to arbitrary reaction systems. In [133] and [134], for example, the log-normal moment closure was found to be significantly more accurate than the normal and central-moment-neglect moment closure for several reaction systems.

Construction of distributions from moments
The system size expansion and moment closure approximations generally only provide approximations of the moments of a process. Suppose we are however interested in approximating distributions of a process. This can be achieved by running Monte Carlo simulations using the SSA or simulations of the chemical Langevin equation. These methods are computationally quite expensive, however. An alternative approach is to first compute approximations of the moments up to a certain order by means of the system size expansion or moment closure methods, and subsequently constructing a distribution from these moments. If the second step can be done efficiently, this can provide a significantly more efficient method than stochastic simulations since one avoids ensemble averaging. One popular method to construct distributions from moments relies on the maximum entropy principle.
0, , where the burst size m is a geometrically distributed random number. The top panel shows results obtained from exact stochastic simulations (SSA) and the lower panel results obtained using the normal moment closure to third order (3MA). The parameter k corresponds to the system size. We find that the results from stochastic simulations converge to a positive fixed point for all system sizes and all initial conditions, and that the mean and variance remain non-negative for all times. The 3MA, in contrast, converges to a fixed point with negative variance for the smallest system size k = 20. For k = 60, the fixed point is now positive, but some initial conditions still lead to trajectories that have a negative variance for some time and are hence physically not meaningful. Only for the largest system size k = 150 do all trajectories lead to nonnegative trajectories. Reprinted from [131], with the permission of AIP Publishing.
For simplicity, assume we have a one-dimensional problem with discrete-valued variable n and that we have approximate values for the first K moments µ µ … , , K 1 , for example obtained by means of the system size expansion. The goal is to construct a distribution P(n) matching these moments. The entropy H of a distribution P(n) is defined as [135] The maximum entropy method aims at finding a distribution P(n) in a certain family of distributions that maximises the entropy. In our case, the family of distributions is given by the constraint that the first K moments of the distribution should match µ µ … , , K 1 . This is a non-linear constrained optimisation problem, which can be solved by means of Lagrange multipliers as follows. We define the Lagrangian Optimising L with respect to P and the λ i corresponds to solving the constrained optimisation problem. The variation of L with respect to P(n) reads We observe that the approximations accurately capture the skewed distribution, with increasing accuracy at higher order. The figure is taken from [137].
with normalisation constant Z. Inserting equation (99) into equation (97) one obtains an unconstrained optimisation problem which can be efficiently solved using standard numerical methods. The maximum entropy method for the construction of the marginal distributions of single species has recently been used in combination with moment closure methods and the system size expansion in [136,137]. Figure 7 shows the result for one example system. A related approach has been developed in [138], where the maximum entropy principle has been used directly to close the moment equations. For time-dependent approximations this method is computationally expensive because the moment equations have to be solved iteratively in small time steps, and in each time step a multi-variate optimisation problem has to be solved. For steady-state approximations, however, the method was found to be efficient and accurate.

Software
Several software packages for using the discussed approximation methods, as well as for deterministic rate equations and stochastic simulation algorithms, are freely available. For exact stochastic simulations, for example, available packages include the Java-based program Dizzy [139], stand-alone COPASI [140], stand-alone StochKit [141] and the python package StochPy [142]. The system size expansion for orders beyond the LNA is implemented in the stand-alone package iNA [143], which also allows one to perform exact stochastic simulations. Various moment closure approximations are available in the Matlab toolbox StochDynTools [144], and the normal moment closure in Python package MomentClosure [145]. These packages require programming knowledge and are only applicable to mass-action propensity functions. The mathematica package MOCA [132] extends the applicability to non-polynomial and time-dependent propensity functions and does not require any programming skills. Similarly, the python package MEANS [146] extends moment closure methods to non-mass-action propensity functions. The Matlab package CERENA [147] implements several of the mentioned methods, including exact stochastic simulations, the system size expansion and moment closure methods. For a detailed overview of available software packages and their capabilities see figure 1 in [147].

Other approximations
One reason why the CLE, the system size expansion and moment closure approximations discussed above are so popular is that they are easy to implement, do not require any preknowledge of the system (except monostability for the system size expansion), are generally efficient computationally, and often give accurate approximations. The three methods have been frequently applied successfully in the literature. However, there are many scenarios where the three methods give quite inaccurate results. In particular if one or several of the species occur in very low copy numbers, the three methods often perform poorly. A large number of other approximation methods have been developed in the literature, and we give an overview of these methods here. Many of these methods are more sophisticated but only apply to certain classes of systems and require pre-knowledge and/or fine-tuning. As we shall see, the CLE, the system size expansion or moment closure approximations form building blocks of many of these methods. 4.6.1. State space truncation. In section 3.2 we saw that the CME can be solved exactly by matrix exponentiation whenever the state space is finite (see equation (32)). However, for many chemical systems of interest the state space is infinite and the solution in equation (32) can not be computed. The idea of the finite state projection algorithm is to truncate the state space to a finite subspace and use matrix exponentiation to obtain an approximation of the distribution on this subspace [148]. It also provides an estimate of the error, which can be systematically reduced by increasing the truncated space.
While this is an efficient method in some cases, often a large truncated space has to be chosen to achieve a reasonable accuracy, making matrix exponentiation intractable, despite recent progress on numerical algorithms, see [16,66] for overviews. Similarly, sometimes the state space of a system is finite but too large to compute the matrix exponential. Accordingly, several modified versions of the finite state projection algorithm have been developed, see [149] for a review.

Tau-leaping.
Tau-leaping methods are approximate ways of simulating a chemical network with the goal of being more efficient than exact stochastic simulations [150]. The basic idea is to 'leap' along time in certain steps τ during which several reaction events occur, thereby avoiding the simulation of each single reaction event as exact stochastic simulations do. The time step has to be chosen small enough such that the propensity functions do not change significantly during one time step. In that case different reactions become independent of each other, and the number of reaction events ( ) τ N t n ; , r of the rth reaction during the time interval τ given that the state at time t is n becomes a Poisson random variable ( ( ) ) τ P f n r with mean ( )τ f n r . Accordingly, the tau-leaping algorithm updates the state ( ) t n of the system iteratively as where S r is the rth row of the stiochiometric matrix. If many reactions are happening per time step the algorithm is more efficient than exact simulations. Increasing τ increases the algorithm's efficiency but obviously decreases its accuracy. If a system contains species with very low particle numbers τ may have to be chosen smaller than the average inter-reaction time to ensure approximately constant propensity functions, which would lead to the algorithm becoming less efficient than exact simulations. Roughly speaking, tau-leaping therefore works best for systems with not too small average molecule numbers.
Despite its simplicitly, implementing the tau-leaping method bears several difficulties. Most importantly, the trade-off between accuracy and efficiency makes the choice of the step size τ non-trivial. A bad choice can lead to highly inefficient or inaccurate results and also to other problems, such as negative molecule numbers. Recent years have therefore seen a wide range of studies developing modified versions of the tau-leaping method that aim at solving these problems, including implicit methods [151,152], binomial methods [153][154][155][156], multinomial methods [157], K/R-leaping [156,158] and a post-hoc correcting method [159]. See [8,9,160] for overviews of these methods.

Time-scale separation.
Many biochemical systems involve processes with highly varying time scales. Given a certain separation of time scales in a system, it is often possible to derive a reduced model of the system which either allows more efficient simulations or more accurate approximations. Such methods have been well-known for deterministic ordinary differential equation models of biochemical kinetics for several decades [161] but have only more recently been developed for stochastic methods [162]. We start here by describing the deterministic setting. The deterministic case. In the following we describe two popular methods for reducing deterministic rate equations under time-scale separation. The first method assumes that the chemical reactions in a given system can be divided into 'slow' and 'fast' reactions, i.e. reactions that happen very infrequently or very frequently on a certain time scale of interest. In this case it is sometimes possible to derive a reduced model by eliminating the fast reactions, typically by assuming that certain reactions balance each other. Such methods are often called quasiequilibrium approximations (QEA) [163].
The second method separates a system into 'slow' and 'fast' species, rather than reactions. The fast species are assumed to be asymptotically in steady state on the time scale of the slow species. This is known as the quasi-steady-state approximation (QSSA). The idea is to eliminate the fast species from the system and to include their steady-state effect on the slow species.
Example. Let us illustrate the QEA and QSSA in the deterministic setting by means of a simple example, the well-studied Michaelis-Menten system: A substrate S reversibly binds to an enzyme E to form the substrate-enzyme complex C, from which a product molecule P becomes catalysed. The corresponding rate equations are where φ φ φ , , S E C and φ P denote the concentrations of S, E, C and P, respectively. Let us first make the assumption that the binding and unbinding of the substrate and the enzyme in equation (101) are fast reactions and that the catalysis reaction is slow. This is the case whenever k k This corresponds to the quasi-equilibrium approximations (QEA). Note that the total concentration E 0 of enzyme molecules is conserved, i.e. φ φ . Using this to eliminate φ E in equation (103) be the production rate of the product P. According to equation (102) this becomes which is the well-known Michaelis-Menten equation [164]. Next, we illustrate the quasi-steady-state approximation (QSSA). To this end, we assume that the complex is approximately in steady state: It can be shown that this is a reasonable approximation whenever φ φ [165]. This corresponds to the complex C and free enzyme E being fast species, and the substrate S being a slow species.
Using equation (105) in equation (102) and the conservation law φ φ This is also called Briggs-Haldane kinetics [166]. Note that equation (106) has the same form as equation (104) but with a different constant ′ K M , showing that the QEA and QSSA are generally not equivalent. Depending on the parameters, only one of the two conditions (or none) may be satisfied. This is illustrated in figure 8.
The stochastic case. For stochastic systems described by the CME, reductions based on time-scale separations are not as straightforward as in the deterministic case. Ideally, one would like to derive a reduced CME and/or a stochastic simulation algorithm for a reduced system. As it turns out, under the quasi-equilibrium assumption of fast and slow reactions it is indeed possible to derive a reduced CME [167]. In contrast, the quasi-steady-state assumption of fast and slow species does not necessarily allow to derive a reduced master equation since it leads to a non-Markovian process for the fast species [168]. One can use a reduced master equation with non mass-action propensity functions that corresponds to the reduced deterministic system under the quasi-steady-state assumption. This heuristic master equation has been found to be accurate for some examples [169][170][171]. However, even if the quasi-steady-state assumption is valid for the deterministic system, the reduced master equation can be highly inaccurate [172,173]. The validity of a reduction in the deterministic case does not generally imply a valid reduction in the stochastic case. The relationship between the two descriptions with respect to time-scale separation approximations has for example been studied in [173]. In [174] a reduced description of stochastic models under the quasi-steady-state assumption was derived by means of the linear noise approximation.
It is however possible to derive reduced descriptions based on quasi-steady-state assumptions for a stochastic system if one requires stronger conditions on a system than the conditions needed to reduce the deterministic system. In [175], for example, reactions are split into fast and slow reactions, and slow species are defined as those involved in slow reactions only and fast species as those participating in at least one fast reaction and any number of slow reactions. As mentioned before, the fast variables conditioned on the slow ones are not Markovian. The authors in [175] solve this problem by approximating the fast process by a virtual fast process that is affected by fast reactions only and therefore Markovian.
A large number of other approximation methods for CME type systems based on timescale separations have been developed in recent years, see for example [167,[176][177][178][179][180][181][182][183][184][185]. Most of these methods split up the full CME into a CME for the fast variables conditioned on the slow reactions/species, and a CME for the slow reactions/variables with marginalised fast variables. On the time scale of the slow dynamics, the conditional CME of the fast dynamics is assumed to quickly reach steady state. Therefore, the CME of the slow dynamics only depends on the steady-state distribution of the fast dynamics. Most methods rely on SSA simulations of the slow dynamics and assume that the fast dynamics is in steady state [175][176][177][178][179]183].
For some systems, the reduced CMEs under time-scale separation have been solved exactly, for example for gene expression [186][187][188]. In [189] a perturbative expansion based on time-scale separation for a two-stage gene expression model is derived that allows one to systematically include higher order corrections to the solution.
The studies mentioned so far mainly rely on stochastic simulations of the reduced equations or on analytic solutions for special cases. An alternative approach is to combine timescale separation reductions with other approximations, such as approximating the fast variables by the chemical Langevin equation [183,190,191], deterministic rate equations [190], or moment closure methods [180]. Other methods combine time-scale separation with the finite state projection algorithm [192] (see section 4.6.1), tau-leaping [193], or the linear noise approximation [174,194,195]. An adiabatic approximation derived form a stochastic path integral description has been developed in [196]. Other methods derive reduced models for systems that allow a partition in more than two typical time scales, see for example [179,197]. In [198] a conditional linear noise approximation for gene expression systems with finite, slow promoter states has been derived. 4.6.4. Hybrid methods. In many biochemical reaction networks of practical interest no timescale separation assumptions apply. However, often some species occur in low and others in high copy numbers, which motivates the combination of different simulation and/or approx imation methods for these two groups of species, similar to the time-scale separation case. One possibility would be to partition species into a discretely modelled group and a continuously modelled group, and accordingly reactions into discrete reactions involving discrete species and continuous reactions that do not involve discrete species. In this case species that are involved in both continuous and discrete reactions describe diffusion-jump processes. Such methods are broadly referred to as hybrid methods. Note that many of the methods discussed in section 4.6.3 fall under this definition. We give here an overview of hybrid methods that do not rely on separation of time scales, although a clear distinction between the two cases is often not possible.
Consider the case where we split the species of a system into a low abundance group which we model discretely using the SSA and a high abundance group which we model continuously using deterministic rate equations or chemical Langevin equations. Accordingly, we group reactions into discrete reactions involving discrete species and continuous reactions that do not involve discrete species. Between two discrete reaction events, the continuous variables follow rate equations or chemical Langevin equations which can be solved numerically or simulated in a standard way. The propensity functions of the discrete reactions, however, may depend on the continuous variables and hence depend on time in such a hybrid approach. This is akin to the case of extrinsic noise discussed in section 3.3.1, where a discrete system's propensities where assumed to depend on some external stochastic process. Correspondingly, the discrete system of our hybrid system can not be simulated by a standard SSA. As in the extrinsic noise case, one possibility is to numerically integrate the time-dependent propensity functions over time until a target value is reached and the corresponding reaction occurs. The discrete system is then updated accordingly which may also change the propensity functions of the continuous reactions.
A key challenge of such an approach is to decide which species should be modelled discretely and which continuously, in particular since some species may fluctuate strongly during a simulation. More sophisticated algorithms therefore partition species and reactions adaptively during simulation. Many different methods addressing these and other issues have been developed in the literature, see for example [199][200][201][202][203][204][205][206][207][208][209][210][211]. They differ mainly in how the partitioning is conducted and in the simulation methods for the resulting reduced system. Simulation-based approaches include combination of the SSA for the low abundance species with tau-leaping [204], chemical Langevin equations [202,206,208,209,212], the chemical Fokker-Planck equation [210] and deterministic rate equation approximations [199-201, 205, 211] for the abundance species. Other approaches split species into more than two sets [203].
Since many of these methods are heuristic, it is not straightforward to assess their performance or to prove their convergence to the exact system in some limit. In some cases, this is possible, however. For example, error bounds for some hybrid methods have been derived [208,213,214]. In [215] the convergence of different types of hybrid methods to the exact system have been studied, and in [184] criteria have been developed for the convergence of a CME to discrete-deterministic hybrid approximations. An error analysis of various methods has recently been conducted in [216].
While simulation-based hybrid methods are often orders of magnitude more efficient than the standard SSA, they all still rely on stochastic simulations for the low abundance species and some of them also for the high abundance species. They therefore can still become computationally expensive depending on the studied system. Some methods aim at circumventing expensive simulations by applying additional approximations to the reduced systems. In [217], for example, the dynamics of the large abundance species is formulated in terms of moment equations conditional on the low abundance (discrete) species. For non-linear systems these can be closed by means of moment closure approximations (see section 4.3). This method amounts to the solution of a differential algebraic equation which is difficult to solve. The authors in [217] propose several simulation based algorithms. This method has shown to give accurate approximations to the distributions and moments of some gene systems, but its implementation is non-trivial and requires additional approximations and/or simulations. In [218] a simpler conditional moment closure has been proposed for two-state gene systems. Here, the conditional moment equations are closed by means of the second-order normal closure or derivative matching (see section 4.3). This method amounts to solving a coupled set of ordinary differential equations. In [219] a hybrid method is developed where the abundance species is described by the (non-conditional) rate equations. This is a simplification of the previously mentioned methods that model the abundance species by rate equations conditional on the low abundance species. However, this simplification allows one to derive analytic solutions of the CME for the low abundance species for certain systems [219].

Comparison of approximation methods
In the previous section we gave a detailed introduction to three popular approximation methods of the CME: the chemical Langevin equation (CLE, section 4.1), the system size expansion (section 4.2) and moment closure approximations (section 4.3). These methods have been successfully used in many applications in the literature, but there exist only very few studies comparing their accuracy. It thus remains unclear how the different methods compare to each other.
In sections 4.1-4.3 we discussed advantages and disadvantages of the different methods. Here, we perform a numerical case study to enable the reader to understand the differences between the methods. First, we study an enzyme reaction system of the Michaelis-Menten type in section 5.1 which can be viewed as a catalysed degradation of a spontaneously produced protein. We then extend the model in section 5.2 by including transcription of mRNA and translation of the protein from mRNA, which allows for bursts in protein production. Finally in section 5.3 we summarise the results and give an overview of advantages and disadvantages of the different methods.
Implementation. Implementation details for the different methods are: • CLE: we simulate the CLE using the Euler-Maruyama algorithm introduced in section 4.1.1. As pointed out in section 4.1, the CLE is traditionally defined for real variables, but suffers from the occurrences of square roots of negative expressions for which it is not well-defined. We therefore implement two versions of the CLE here: (i) a real-valued implementation which keeps the variables positive by rejecting Euler-Maruyama steps that lead to negative variables. We term this version 'CLE-R'; (ii) the recently proposed complex CLE which is defined for complex variables and has been shown to give realvalued moments in the ensemble average [87]. We term this version 'CLE-C'.
The Euler-Maruyama algorithm requires a step size dt for time discretisation. For the simulation of steady-state moments or distributions, we take a total number of M samples from a single long trajectory at time steps separated by ∆t. The chosen values for these parameters will be given for each of the examples. • System size expansion and stochastic simulation algorithm: we use the stand-alone software package iNA [143] for both the system size expansion and the stochastic simulation algorithm (section 3.3). We study the zeroth order system size expansion, i.e. the linear noise approximation (termed 'LNA', section 4.2.2), as well as the first order corrections to the mean and variance (both termed 'SSE-1') given in equations (76) and (77), respectively. • Moment closure approximations: we study here the second-order normal moment closure (termed '2MA') introduced in section 4.3, which is probably the most commonly used one in the literature. For its implementation we use the Mathematica software package MOCA [132].

Enzymatic protein degradation
Let us consider the well-known Michaelis-Menten system with reactions We studied this system without the first reaction in section 4.6.3 in the context of time-scale separations in a deterministic setting. A substrate molecule S is created spontaneously and binds reversibly to a free enzyme molecule E to form a complex C, which catalyses the substrate into a product molecule P. We consider S to be a protein in the following.
Here, we are interested in the accuracy of the different approximation methods as compared to exact stochastic simulations of the corresponding stochastic system. Note that the total number of enzymes E 0 is conserved. The system has a steady state in the protein numbers if and only if , which means that the input rate must be smaller than the maximum turnover rate. The parameter α can hence be viewed as a saturation factor. Figure 9 shows the relative mean and variance of the protein S as a function of α for a system with a large number of total enzymes, E 0 = 2500. We find that both the real and complex CLE implementations (CLE-R and CLE-C), the first order system size expansion (SSE-1) and the second-order normal moment closure (2MA) give good approximations for the mean value. Only the LNA, which corresponds to the deterministic rate equations for mean values, shows significant deviations from the exact result. For the variance we observe larger deviations for all methods, with the LNA again being the least accurate. The two CLE implementations show the best performance, being more accurate than the 2MA, which in turn is more accurate than the SSE-1. Note that the CLE-R and the CLE-C give very similar results, meaning that the inclusion of a rejecting boundary has a negligible effect in this case.
Overall, we find that the approximation methods give rise to larger errors for larger values of α, which can be explained as follows. As mentioned before, the system only possesses a steady state for α < 1 and becomes unstable for α > 1. We therefore expect larger fluctuations for values of α close to unity. Moreover, since most of the enzymes are in the complex state C in this limit, we expect a skewed distribution for the substrate. In the other limit, → α 0, most enzymes are in the free state E, which reduces the non-linear effect of the bimolecular reaction → + S E C. We therefore expect the approximation methods to perform well in this limit, and to lead to larger errors for → α 1. This is exactly what we observe in figure 9. Next, we study the steady-state probability distribution of the protein as predicted by the SSA, CLE-R, CLE-C and the LNA. Note that the CLE-C gives complex-valued samples. To obtain a distribution in real variables, we take the real parts of these samples. We reduce the total number of enzymes to E 0 = 60 to study small molecule number effects. The results are shown in figure 10. We find that the LNA strongly underestimates the true mean value (obtained using the SSA) and does not reproduce the distribution very accurately. However, surprisingly, the real-valued CLE (CLE-R), does even worse than the LNA, both in terms of mean value and distribution. The complex valued CLE (CLE-C) on the other hand, predicts the true mean value with a negligible error, and even captures the highly skewed distribution very well (right panel of figure 10). This demonstrates that a naive fixing of the boundary problem of the CLE (CLE-R) can lead to highly inaccurate results.

Bursty protein production
We now extend the system in equation (107) by considering the protein S to be produced in a gene expression motif, i.e. via transcription of mRNA M followed by translation to protein: Depending on the parameters k dM and k s this system can give rise to bursts in the production of protein S, namely whenever each produced mRNA molecule produces on average several proteins during its lifetime. The average number of proteins per mRNA molecule is given by / = b k k s dM and also called 'burst size' [220]. Large b correspond to large bursts which in turn lead to large fluctuations.
Here, we are interested in the accuracy of the different approximation methods as a function of the burst size b, since we expect larger errors for larger fluctuations. To isolate the effect of the burst size b, we vary k 0 and k dM such that the average production rate of protein / = k k k 115 s dM 0 in steady-state conditions is held constant and equal to the production rate used in the previous section in figure 10. Figure 11 shows the mean and variance of the protein S as a function of b as predicted by the different approximation methods. Similar to the previous section in figure 9 we find that the LNA performs worse than the other methods for the mean and the variance, with the exception of the real-valued CLE (CLE-R) which performs even worse. Similarly to figure 10, we find here that the complex-valued CLE (CLE-C) performs surprisingly well, demonstrating again the significance of the boundary problem of the CLE. In fact, the CLE-C performs significantly better than all the other methods, including the 2MA and SSE-1. The latter two perform similarly and both significantly better than the LNA. One should keep in mind, however, that the 2MA, SSE-1 and LNA are computationally significantly more efficient than the CLE-C since they do not rely on sampling.
Next, we consider the steady-state distribution of the protein for a large burst size of b = 7.5. Figure 12 shows the results predicted by the stochastic simulation algorithm, the two CLE implementations, and the LNA. We find that the CLE-R and LNA give highly inaccurate results, both for the mean value and the distribution, with even larger deviations than for the simpler system in the previous section in figure 10. This is to be expected, since the system now leads to large bursts in protein production, while the average production rate of protein and the rate constants of the enzymatic reactions are the same as in figure 10. Similar to the previous section, we find here that the complex-valued CLE gives a highly accurate approximation for the distribution, despite the large skewness of the later. Once again we find that the boundary problem of the CLE can lead to highly inaccurate results if not treated carefully.

Numerical results.
In this section we gave a numerical comparison of the CLE, the LNA, the normal moment closure to second order (2MA) as well as the next leading order corrections to mean and variance of the system size expansion beyond the LNA (SSE-1). We implemented two versions of the CLE: a real-valued version with rejecting boundary and a complex-valued version. We considered an enzymatic protein degradation system (section 5.1) and an extension thereof including bursty protein production (section 5.2). One important observation we made is the large discrepancy between the two CLE implementations: while we found a negligible difference in the case of large total enzyme numbers (figure 9), the complex version was significantly more accurate in the case of smaller total enzyme numbers (figures [10][11][12]. Crucially, in the latter cases, the real-valued CLE performed worse than all the other methods, while the complex-valued CLE performed better than all the other methods. This illustrates the boundary problem of the CLE, and how significant inaccuracies can arise when fixing it in a naive way. Another important observation is that the LNA was throughout found to be less accurate than all the other methods (except the CLE-R). This result is not very surprising, since the LNA corresponds to the deterministic rate equations on the mean level, and gives the zeroth order fluctuations (in terms of the system size expansion) about the deterministic mean for the variance. The SSE-1 includes the leading order corrections to both the mean and variance. Similarly, the 2MA and CLE capture effects of fluctuations on the mean. It is hence not surprising that these methods perform better than the LNA. In terms of the mean and variance (figures 9 and 11), we found that the SSE-1 and 2MA performed similar to each other, except for the variance in figure 9 where 2MA was more accurate. Surprisingly however, the CLE-C was found to be significantly more accurate than the 2MA and SSE-1 in most cases.
In terms of distributions (figures 10 and 12) we found the CLE-C to be highly accurate, even though the distributions were highly skewed. The LNA, which predicts a Gaussian distribution, was obviously not able to capture these skewed distributions.

5.3.2.
Advantages and disadvantages. The advantage of the CLE over the other methods is that it gives approximations of the process and distributions in contrast to moment closure methods and the system size expansion. Moment closure approximations give only approximations to the first few moments of a process. The system size expansion in principle predicts distributions. However, closed-form solutions for higher orders beyond the LNA have so far only been derived for one-dimensional systems [108]. It is not clear if the same is possible for multi-species systems. The higher orders of the system size expansion have therefore mainly been used to approximate the moments of a process. If one is interested in approximating the whole process or its distributions, the CLE is therefore a useful method.
Suppose now that we are only interested in the moments of a process. In the numerical case study performed before, we found the CLE to be more accurate than the LNA, the SSE-1 and the 2MA. However, the CLE is computationally significantly more expensive than the other methods. While the CLE requires a large number of stochastic simulations to obtain the moments of a process, the other methods only require the numerical solution of a finite set of ODEs and are hence typically orders of magnitude faster. Moreover, if defined for real-valued variables, the CLE suffers from a boundary problem at zero molecule numbers, and realvalued modifications lead to inaccurate results. The boundary problem is solved by extending the CLE to complex-valued variables, which is however less efficient to simulate [87]. Due to these reasons, it seems preferable to use the system size expansion or moment closure approximations if one is only interested in the moments of a process.
Next, the question arises if the system size expansion or moment closure approximations are preferable. While the system size expansion is a systematic expansion in a small parameter, , corresponding to a burst size of moment closure approximations are an ad hoc approximation. This fact makes the system size expansion more appealing, since it is guaranteed to be accurate for large system volumes. The same cannot generally be expected to be true for moment closure approximations. On the other hand, the system size expansion has the disadvantage that it is not applicable to systems that are deterministically multi-stable, a limitation not shared by moment closure methods. Moreover, higher order corrections of the system size expansion are significantly harder to derive and implement than higher order moment closure methods. Current software packages implementing the system size expansion only implement two orders beyond the LNA for the mean, and one order beyond the LNA for the covariance [143,147]. Moment closure approximations, on the other hand, are implemented to various orders [132,146,147]. Due to these reasons, it depends on the problem at hand as to decide which method is preferable.

Inference
So far, we have focussed on the forward problem of approximating marginal distributions of a fully specified process. Such distributions depend naturally on the parametrisation of the process: it is not uncommon for e.g. steady-state distributions to exhibit qualitatively different behaviours depending on the specific value of reaction propensities. In many concrete applications, the model parameters may only be known approximately: direct measurements of kinetic reaction parameters are difficult to obtain, and, even in cases when good estimates are available, in vivo parameters of a concrete system embedded in a cell may be influenced by a plethora of additional factors, leading to significant uncertainty. It is therefore of considerable interest to also address the inverse problem: using (noisy) observations of a system to constrain the uncertainty over model parameters and/ or predictions. This is a well-studied problem in statistics and machine learning: we give here a brief review of recent developments within the Bayesian approach to solving this inverse problem, with particular attention to methodologies which have employed the approximation methods described earlier.

General problem formulation
The setup we will consider is the following: we consider a stochastic process ( ) θ | p x T 0: as a measure over the space of trajectories x T 0: of the system in the time interval [0, T ], with x t representing the value of the state variable at time t. In the case of a CME system, such trajectories will be piecewise constant functions from [0,T ] onto a discrete space, while for a continuous approximation (e.g. CLE or LNA) x t will be real-valued. The stochastic process is assumed to depend on a set of kinetic parameters θ, whose a priori uncertainty is captured by a prior distribution ( ) θ p . Additionally, we assume the existence of a measurement process which associates each trajectory x T 0: with an observed random variable y; in the simplest case, the observed variable y may just be thought of as the state of the system at a particular set of time points, corrupted by random observation noise. We account for such experimental errors through an observation model encoded in a probability distribution ( ) | p y x T 0: , i.e. the likelihood of a measurement given the true state of the system, which may depend on additional parameters (here omitted for notational conciseness). We restrict our interest here to the case where the observations are in the form of a time series of state variable observations with independent and identically distributed noise, i.e. measurements of all or of a subgroup of the species in the system. Hence, the observation vector will take the form ( ) = … y y y , , T 0 for some discrete time points 0,...,T. More general cases where the observations take other forms, such as continuous-time constraints or penalties over particular areas of the state space, are treated for example in [221,222].
Also, we consider all observations to come from a single trajectory, or from replicate trajectories with the same (unknown) kinetic parameter values; in other words, we exclude the case where parameters can also be variable, e.g. to accommodate extrinsic noise between cells [223].
The general inference problem is then the problem of computing the joint posterior meas- The second term in equation (110) can be obtained by the Backward algorithm, a recursion procedure similar to filtering (but running backward in time) [225]. Furthermore, a modification of the algorithm, the so called forward filtering/backward sampling algorithm [227], can be used to draw sample trajectories from the posterior process, which can then be used in MCMC approaches for joint state and parameter inference.

Parameter inference
Suppose we are not interested in state inference but only in the parameters of the system. In this case we do not need to compute the posterior marginal in equation (110). Rather, it is sufficient to either run the forward algorithm (filtering) or the backward algorithm, since each of them independently deliver the marginal likelihood ( ) p y . The backward algorithm computes the likelihood ( ) | > p y x i t t of future data conditioned on the current state, and ( ) p y is simply the end result of the recursion algorithm at time 0.
To see that the forward algorithm also allows to compute the likelihood ( ) p y , note that due to the Markov property the latter can be written as We find the factors on the rhs of equation (112) are just the normalisations factors of the Bayesian updates of the filtering procedure in equation (112). Optimising ( ) ( ) θ = | p p y y with respect to the parameters yields asymptotically consistent parameter estimates, also called maximum likelihood estimate. In a Bayesian framework, one would combine the likelihood with a parameter prior ( ) θ p to give the posterior over the parameters according to Bayes' law which also allows to quantify uncertainty of the inferred parameter values.

Computational methods for Bayesian inference in stochastic chemical reaction networks
6.4.1. Methods for general networks. The primary difficulty in applying the forward-backward approach to inference in chemical reaction networks is the requirement for forward integrability of the system dynamics: calculating the transition probabilities ( ) | − p x x i i 1 requires solving the CME, which is generally not possible analytically. Some approaches resort to numerical integration of the equations, including the variational approach of [228] and the uniformisation sampler of [229]. These approaches can be effective for closed systems with low molecular numbers; however their application to open systems invariably requires an artificial truncation of the state space, introducing a bias which is hard to quantify. Truncations are also used in the recent work of [230]; however, here a random truncation scheme guarantees unbiasedness of the results, as well as leading to substantial computational savings. Other approaches that can handle open systems either introduce additional latent auxiliary variables (such as the number of reactions in the time interval as in [231]), or resort to a sequential Monte-Carlo scheme which relies on multiple simulations from different initial conditions (particle filtering, [232]). Both such schemes incur potentially large computational overheads.
The computationally intensive nature of inference methodologies adopting a microscopic system description constitutes a formidable obstacle to inference in large-scale reaction networks. This has justified a considerable interest in the use of mesoscopic approximations for inference. One of the earliest attempts [233] relied on the chemical Langevin equation approximation to the CME (section 4.1); this provides a more efficient inference scheme compared to the auxiliary variable approach of [231], however the computational costs remain high due to the need to compute transition probabilities for non-linear diffusion processes. In this light, the LNA provides a more promising avenue, since the sufficient statistics of the (Gaussian) single-time marginals can be efficiently computed by integrating a system of ordinary differential equations. Several authors have therefore proposed inference schemes which integrate the LNA approximation [234][235][236]. Moment closure approximations provide an alternative approximation scheme with similar computational complexity to the LNA, however they do not generally compute a marginal distribution, rather only a few moments of a generally unknown distribution. Their use for time series inference is therefore limited to second-order normal moment closure schemes, where a Gaussian approximation is taken. This approach has been proposed in [237], where it has been shown to yield accurate results with modest computational overheads. In [238] the second-order normal moment closure has been combined with the chemical Langevin equation and integrated into an expectation-propagation algorithm for intractable likelihoods and continuous-time constraints. Another interesting opportunity for moment-based inference is offered by flow-cytometry data: here, simultaneous measurements of millions of cells enable an empirical characterisation of the marginal moments directly (albeit potentially corrupted by extrinsic noise), which can then be fitted to a moment-approximation of the CME [239]. 6.4.2. Inference for gene expression data. The inference methods mentioned above do not assume any knowledge about a given system. While this makes them in principle applicable to any type of reaction network, more efficient and/or accurate methods can often be employed by including a priori knowledge in a model. Gene expression systems constitute a particularly important example where this is often the case. Cells typically possess only one or very few copies of a gene. Proteins and mRNA molecules on the other hand often occur at copy numbers that are orders of magnitude larger. Such systems are therefore often suitable for hybrid methods (see, section 4.6.4) that model some of the species as discrete variables (e.g. the genes) and others as continuous variables (e.g. the mRNA and proteins).
While it is often not straightforward to integrate hybrid methods into inference schemes, significant progress has recently been made in this respect. In [240], for instance, the different promoter states of a gene are modelled as a change-point process which drives a linear SDE representing the protein dynamics, and an efficient MCMC inference method is developed. The same model has been integrated into a variational inference method which additionally allows transcriptional feedback in [241,242]. In [243], a particle MCMC scheme is developed based on a hybrid method that approximates the continuous variables by the linear noise approximation (see section 4.2.2). More recently a different hybrid approach combining several types of approximations has been integrated into an efficient Bayesian inference scheme in [244].

Summary
This bird eye survey of inference for stochastic chemical reaction networks highlights the diversity of statistical research in the area. Such diversity can appear baffling to the outsider, and a major problem for the greater diffusion of these ideas is the lack of standard software tools. Inference tools for stochastic chemical reaction networks often form a small subsection of software tools for parameter estimation of deterministic methods [245,246], and there haven't been systematic comparisons of various inference schemes that investigate the relative merits of the different algorithms on a number of relevant examples. Furthermore, inference approaches inevitably construct an approximation of a posterior distribution; while much statistics research investigates the convergence properties of these approximations, it remains entirely unclear how inference errors combine with approximation errors when inference schemes are deployed on approximate dynamics. Due to these reasons, we have chosen not to include an explicit numerical comparison between inference methods in this tutorial, as we feel this would deserve a separate review on its own.

Conclusions
Recent years have seen an explosion of experimental studies revealing the crucial role that stochastic fluctuations in chemical reaction networks play for living cells. Driven by these discoveries a plethora of methods for the mathematical and statistical analysis of such systems has evolved. The goal of this review is to give a self-contained introduction to the field of modelling for stochastic chemical kinetics. Moreover, it introduces key approximation and inference methods for this field and gives an overview of recent developments.
The chemical master equation (CME) constitutes the accepted non-spatial description of stochastic chemical networks. Recent years have seen a burst of analysis and approximation methods based on the CME. We gave here an introduction to the CME modelling framework and discussed stochastic simulation and analytic solution methods. Next, we introduced various approximation methods with particular focus on the chemical Langevin equation, the system size expansion and moment closure methods. We also gave an introduction to time-scale separation based approximations, as well as hybrid methods and reviewed the existing literature. Finally, we gave an introduction to the problem of statistical inference from experimental data in a Bayesian framework, and reviewed existing methods. The presentation is aimed to be a self-contained introduction for scientists from different disciplines.
In a numerical case study we compared the chemical Langevin equation, the zeroth order system size expansion (linear noise approximation, LNA), the first-order corrections of the system size expansion to mean and covariance (SSE-1) and the second-order normal moment closure (2MA) with exact results obtained using stochastic simulations. In terms of moments, we found that a naive real-valued implementation of the CLE (CLE-R) enforcing positive concentrations was less accurate than all the other methods. A complex version of the CLE (CLE-C), in contrast, was found to be the most accurate of all methods. The SSE-1 and 2MA performed similar two each other and significantly better than the LNA. In terms of steadystate distributions, we compared the CLE-R, the CLE-C and the LNA with exact results obtained using the SSA. We found that the CLE-R and the LNA were not able to accurately capture the distributions, the LNA being more accurate than the CLE-R, however. The CLE-C, in contrast, gave accurate approximations even for highly skewed distributions.
The CME is a valid description of systems that are well-mixed and sufficiently dilute, i.e. the diffusion of particles constitutes the fastest time scale of the system and the total volume of all molecules in the model is much smaller than the system volume. While these assumptions are valid in some cases, it turns out that they are not met by many biological systems. Whenever this is the case, models need to be employed that take spatial positions and diffusion of particles into account. The main family of such models goes under the name of stochastic reaction-diffusion processes (SRDPs), Markovian models where independent particles diffuse in space and react whenever they come in contact (or sufficiently close). The evolution equation for the marginal probabilities of an SRDP, the spatial analogue of the CME, is complicated by the fact that the number of particles varies in time, and needs to be defined on an infinite-dimensional Fock space [247,248]. Solving, or even approximating, such an equation is essentially impossible, and hence SRDPs are mainly analysed in an algorithmic way: each particle is simulated performing Brownian motion in continuous space and chemical reactions between particles happen stochastically under certain rules [249]. This is computationally extremely expensive and significant effort has been spent in improved simulation methods [250,251]. An alternative modelling framework is given by the reaction-diffusion master equation (RDME) which coarse-grains an SRDP by assuming a compartmentalisation of space and locally homogeneous conditions within each compartment [252]. While simulations in this framework are generally more efficient than in the continuous case, they are typically still expensive [253], generally significantly more expensive than in the non-spatial CME case. More importantly, the RDME is not a systematic discretisation of an SRDP, in the sense that its continuum limit generally does not lead to the original SRDP in two or more spatial dimensions [252], because bimolecular reactions become vanishingly infrequent in the continuum limit.
Due to these reasons, analysis methodologies for spatial models are much less developed than in the non-spatial CME case. Some studies investigating spatial stochastic phenomena and comparing the SRDP and RDME approaches include [249,[254][255][256][257][258][259]. In contrast to the CME case, very few studies have attempted analytical approximations for SRDPs. However, some progress has been made in recent years in this respect. In [260] and [261], for example, the linear noise approximation was extended to spatial systems, and in [262] higher orders of the system size expansion for effective one-species systems. Even fewer studies have addressed inference for SRDPs. Only a handful of studies have approached the issue of inference working directly with the SRDP or RDME framework, using likelihood-free sampling methods [263] or variational approximations [264]. In [265], the linear noise approximation has been integrated into an inference scheme for SRDPs. This method is however limited to certain classes of one-dimensional systems. More recently, in [82] it was shown that SRDPs can be approximated by spatio-temporal point processes, a popular class of models from statistics [266]. This was used to derive an efficent inference algorithm for general SRDPs.
In summary, we have presented an introduction to modelling, approximation and inference methods for stochastic chemical kinetics based on the chemical master equation. We hope that this review will help scientists from other disciplines to dive into this exciting field, and that it will stimulate research in the presented areas.