Efficient simulation of intrinsic, extrinsic and external noise in biochemical systems.

MOTIVATION: Biological cells operate in a noisy regime influenced by intrinsic, extrinsic and external noise, which leads to large differences of individual cell states. Stochastic effects must be taken into account to characterize biochemical kinetics accurately. Since the exact solution of the chemical master equation, which governs the underlying stochastic process, cannot be derived for most biochemical systems, approximate methods are used to obtain a solution.
RESULTS: In this study, a method to efficiently simulate the various sources of noise simultaneously is proposed and benchmarked on several examples. The method relies on the combination of the sigma point approach to describe extrinsic and external variability and the τ -leaping algorithm to account for the stochasticity due to probabilistic reactions. The comparison of our method to extensive Monte Carlo calculations demonstrates an immense computational advantage while losing an acceptable amount of accuracy. Additionally, the application to parameter optimization problems in stochastic biochemical reaction networks is shown, which is rarely applied due to its huge computational burden. To give further insight, a MATLAB script is provided including the proposed method applied to a simple toy example of gene expression.
AVAILABILITY AND IMPLEMENTATION: MATLAB code is available at Bioinformatics online. CONTACT: flassig@mpi-magdeburg.mpg.de. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

1 Introduction

Variability and heterogeneity are fundamental properties of biological systems. Cells differ in all kinds of attributes including cell size, protein abundances and morphology (Spiller ), which is caused by various sources of noise. In this study we refer to intrinsic noise as an inherent stochastic biochemical process, extrinsic noise as cell-to-cell variability and external noise as environmental fluctuations. Intrinsic noise is very dominant for small biochemical reaction systems involving low copy numbers of chemical species (e.g. gene networks), whereas extrinsic and external noise increase with system size (Patnaik, 2006). Since the cellular abundance of numerous chemical species span all scales from just a few (e.g. genes) to several millions (e.g. proteins) all sources of noise contribute to biological heterogeneity. The interplay between intrinsic, extrinsic and external noise and their effects on system dynamics is hardly exploited due to experimental and numerical challenges (Spiller ; Fernandes ; Delvigne ). Monte Carlo (MC) techniques are the standard approach to tackle stochastic biochemical reaction networks, but they suffer from an immense computational burden especially in optimization problems, where a system has to be simulated numerous times. Approximations have to be used in order to make computations feasible.

Figure 1A illustrates a cell population that is corrupted by intrinsic, extrinsic and external noise. In this scenario, a heterogeneous population of cells under realistic conditions, which differ in size, shape and number of organelles is situated in an inhomogeneous medium containing concentration gradients (gray background). In addition to that, intrinsic noise caused by switching of a single gene between different states contributes to the overall variability. To clarify our understanding of intrinsic, extrinsic and external noise their impact on a simple decay process describing the degradation of the protein P is shown in Figure 1B–D. Here, intrinsic noise is modeled by the Gillespie algorithm (SSA) (Gillespie, 1977) that treats reactions as stochastic events, whereas extrinsic noise is modeled by distributed initial conditions, which accounts for cell-to-cell variability. Both effects rely on different mechanisms, but result in a probability density function (PDF) characterizing the abundance of the protein (red). Adding both effects yields a further spread of the resulting PDF indicated by an increase of the standard deviation, see Figure 1D. Note that external noise was left out in the investigation of the decay process, because in the case of a stationary random process its mathematical treatment is identical to extrinsic noise. In case the external noise is a stochastic dynamic process, its mathematical treatment is identical to intrinsic noise. For simplicity, we refer from now on only to the terms extrinsic and intrinsic noise, representing a stationary random variable or a stochastic dynamic process, respectively.

Fig. 1

Noise in biochemical systems. (A) Various sources of noise corrupt biochemical reaction systems. The influence of intrinsic noise (B), extrinsic noise (C) and intrinsic noise combined with extrinsic noise (D) on a decay process results in differences of the probability density function (PDF), the mean and the standard deviation (std)

In this study, we propose an approximate method to model distributed stochastic processes by means of stochastic differential equations. The method relies on the combination of the sigma point (SP) approach (Julier ) accounting for extrinsic noise and the τ-leaping algorithm (Cao ) capturing probabilistic reactions. The article is organized as follows. In Section 2 a brief introduction to modeling of stochastic biochemical systems and a short overview of existing methods is given. In Section 3 the idea and advantages of the proposed method are described. In Section 4 a detailed benchmark on an example of gene expression is performed followed by the application of our method on several parameter optimization examples. Section 5 summarizes and discusses our results. We provide Supplementary Material with additional information concerning our methodology and the benchmark systems. The gene expression example is illustrated with a MATLAB script using the proposed method.

2 Theoretical background

The simplest and most common approach to model biochemical systems is by means of ordinary differential equations (ODEs)

They describe the temporal evolution of the continuous state vector by reaction rate equations. In this context denotes the stoichiometric matrix, the reaction rates, a set of parameters and external deterministic forcing. ODEs in this form neglect stochasticity introduced by intrinsic noise and fail to model the underlying process correctly. A more detailed description can be achieved with the chemical master equation (CME) (Gillespie ) taking into account the inherent stochastic nature of biochemical reactions. The CME governs the temporal evolution of the probability P to find the system in the discrete state . a denotes the propensity of reaction k and the k’s column of the stoichiometric matrix. For simplicity and are left out in Equation (2), but they can be easily introduced by interpretation as additional reaction channels (Sanft ). From the solution P for all reachable states a PDF describing the stochastic variable can be reconstructed for every time point. can be interpreted as a vector, whose entries represent the probability of a certain abundance interval of a chemical species. For most biochemical systems it is not possible to find an exact solution of the CME, so approximate methods have to be used. The SSA and its derivatives (Gillespie ), such as leaping (Cao ; Fu ) and time-scale separation approaches (Marchetti ) are powerful methods, which rely on a statistical mechanics ansatz treating chemical reactions as discrete molecular collision events. Although these algorithms are very popular they fail to capture variability introduced by extrinsic noise. The finite state projection algorithm (Munsky and Khammash, 2006) obtains a solution by integration of the CME. This algorithm accounts for intrinsic and extrinsic noise, but due to the curse of dimensionality it is not applicable to systems with a large state space. Another approach is the method of moments (Lakatos ), which relies on the integration of coupled ODEs for the statistical moments. This method incorporates intrinsic and extrinsic noise, but does not reconstruct the state vector’s PDF and numerical instabilities make it difficult to handle (Lee ; Azunre ). A promising approach consists of solving the CME directly via tensor trains, but the derivation of the tensor trains is nontrivial (Kazeev ). To overcome these drawbacks, we propose a novel method paving the way to further understanding of variability and heterogeneity in biochemical systems.

3 Efficient modeling of intrinsic and extrinsic noise

A straight forward approach to model intrinsic and extrinsic noise simultaneously is to perform MC sampling for extrinsic noise e.g. distributed initial conditions or parameters and to compute intrinsic noise i.e. the temporal evolution of a stochastic process with the SSA. In the limit of infinite function evaluations the obtained histogram is equal to the solution of the CME. Since both methods rely on the generation of random numbers they are extremely time consuming and computationally intense. In order to speed up the algorithm, it is possible to lower the resolution of the histograms and use fewer, broader bins or to perform kernel smoothing, which results in the summation of kernel distributions while keeping a high resolution. This is illustrated for the Schlögl model (Schlögl, 1972) in Figure 2A and B. In this study, only Gaussian kernels are considered with the bandwidth denoting the kernel’s standard deviation. The idea of the here proposed method is to further accelerate the computation by approximating the PDF of the extrinsic noise and propagating it through time by a stochastic process, see Figure 2C. We suggest using the SP approach to account for the extrinsic noise. In contrast to MC sampling, the SP approach chooses only sigma points deterministically and estimates mean E and covariance C of a random variable given by a nonlinear transformation h (see Supplementary Material for further information). In this context refers to the abundance of chemical species and h governs their temporal evolution due to chemical reactions. denotes the number of distributed initial conditions , but it can also include kinetic parameters of the right hand side of Equation (2). The samples are chosen according to where denotes the ith row of the covariance matrix in the original domain. We chose the free parameter κ according to (see Supplementary Material for more information). The transformed samples are used to compute the mean and covariance with the weights Assuming normality, log-normality or any other appropriate PDF characterized by mean and variance-covariance, the underlying distribution can be estimated from the mean and covariance (Julier ).

Fig. 2

Approximation of MC methods. Histogram binning (A) and kernel smoothing (B) represent common methods to derive PDFs from MC trajectories (gray) sampled from distributed initial conditions (red, left). The proposed SP approach (C) provides an efficient simulation technique sampling from starting points (red, left). The red distribution on the right indicates the pseudo exact solution obtained by MC sampling combined with the SSA. The black histogram illustrates the histogram binning in (A) and the black distributions in (B) and (C) are summed up to obtain an approximate solution (not shown) to the exact solution. For better visualization all distributions are scaled. The accuracy (D) and convergence (E) elucidate the computational advantages of the proposed method

To incorporate intrinsic noise due to stochastic reactions on top of the extrinsic variability, the solution of CME at each point in time can be attributed a corresponding transformation of a probability function. Alternatively, when approximating the CME by sampling the underlying stochastic process with the SSA, single realizations of the SSA yield realizations of a transformation of realizations of in . Note that is a mapping between functional spaces with elements , whereas maps from i.e. . To overcome the SSA’s huge computational load we used an approximate version, which is the τ-leaping algorithm implementation of StochKit2 (Sanft ). Note that by combining intrinsic and extrinsic noise we have to deal with a distributed CME. As for real valued stochastic variables, the corresponding ensemble of probability functions or PDFs of the CME may be characterized by statistical moments. To obtain an estimate of the true average PDF , the scheme is repeated n times and the resulting distributions are used as where we assumed that . A scheme of the algorithm is given in Figure 3 and the corresponding pseudo code can be found in the Supplementary Material.

Fig. 3

Workflow of the proposed algorithm. Note that in 4, one may assume any PDF that is characterized by mean and covariance and also in 5 one may use different weights ω for each sample k. However, the choice is a priori not clear, and we therefore suggest to use a Gaussian PDF and an equal weighting

The sigma points have been applied frequently to deterministic ODEs e.g. (Flassig and Sundmacher, 2012; Schenkendorf ; Toni and Tidor, 2013), whereas (Toni and Tidor, 2013) have proposed a combination of sigma points and Ω-expansion to describe extrinsic and intrinsic variability. We are not aware of applications of sigma points to stochastic ODEs, which we apply here for sampling realizations of the CME to ultimately get an approximate solution to a distributed CME. However, the accuracy of deterministic functions given in (Julier ) should also apply for each incremental time step τ in the τ-leaping realization, since for each realization we have formally a deterministic mapping, which we can expand into a Taylor series as done in (Julier ) for the accuracy analysis. We note, however, that the convergency of is in general not guaranteed. This depends on the ensemble of and the choice of ω and . Even though our results show that the proposed approach works well, a thorough analysis is required at this point which is, however, out of scope of this contribution.

4 Benchmarking

4.1 Comparison to extensive Monte Carlo sampling

To compare all methods illustrated in Figure 2 the Euclidean distance is used as similarity measure (Cha, 2007) between the approximate solution and the pseudo exact solution obtained by MC sampling combined with the SSA. Since the final distribution of all methods is a random variable due to the underlying stochastic process, which is used for temporal evolution of the system, the computation was repeated times to derive a statistical statement. The mean and standard deviation of the Euclidean distance of the nstat distributions are exploited as measures for accuracy and convergence. In Figure 2D and E the dependency of accuracy and convergence on the number of function evaluations for the Schlögl model introduced in the previous section is illustrated. For the histogram and kernel method, the number of function evaluations is given by the number of MC samples and for the proposed method by . As can be seen, our proposed method outperforms the others in accuracy as well as convergence for up to 310 function evaluations, which highlights its computational benefit. For this illustration we used a kernel bandwidth of 10 and a histogram binwidth of 75.

To further demonstrate the computational efficiency of the proposed method an example of gene expression resulting in multimodal probability distributions is investigated. The gene expression system consists of the following reactions A single gene is considered which is able to flip between an active and inactive state. In the active state the gene produces protein A, which is degraded by protein B. This example involves intrinsic noise due to the low abundance of the gene, but also extrinsic noise due to the distributed initial conditions of protein B. In Figure 4A the proposed method is illustrated for protein A. Several distributions of (gray, dashed) are averaged to obtain the resulting approximate distribution (gray, solid). The approximate solution mimics the distributed character of the pseudo reference obtained by MC sampling combined with the SSA, which constitutes a bimodal distribution with—regarding any approximation approach—a challenging sharp peak for low abundances. The temporal evolution of the distributions of protein A and B are shown in Figure 4B and E. The approximation is very similar to the solution obtained by MC sampling combined with the SSA demonstrating the capability to qualitatively model intrinsic and extrinsic noise simultaneously. Furthermore, accuracy and convergence were investigated in comparison to the kernel approach. Since a priori the optimal selection of the kernel bandwidth is not clear several bandwidths were tested systematically. In Figure 4C and F the difference of the mean Euclidean distance of the proposed method and the kernel approach is shown for the time point . If this term takes positive values (shades of red) the proposed method outperforms the kernel approach and for negative values (white) the kernel approach excels. The same applies to the difference of the standard deviation of the Euclidean distance . For protein A, it can be seen that for less than function evaluations the accuracy of the proposed method is always higher than the accuracy of the kernel method. For more function evaluations an optimal bandwidth yields better results at the price of higher computational costs, see Figure 4C. For protein B the proposed method is superior for all bandwidths and function evaluations, see Figure 4F. Concerning the convergence it was found that the proposed method outperforms the kernel approach for both proteins, except for kernel bandwidths much larger than the spread of the underlying distribution resulting in very low accuracy, see Figure 4D and G. With this example including multimodality, low and high copy numbers of chemical species the strengths of the proposed method has been clearly demonstrated. Note that beforehand the optimal kernel bandwidth is not known and changes with time meaning that this parameter has to be estimated for the optimal representation of every PDF. The proposed method avoids this computational demanding task making it easy to handle and applicable.

Fig. 4

Comparison of performance. (A) Averaging of (gray, dashed) yields an approximate solution (gray, solid) of the CME for protein A at the time point 250 s. The pseudo exact solution obtained by MC sampling combined with the SSA is shown in red. For better visualization, the are scaled and hence smaller than the approximate distribution. The corresponding temporal evolution of the probability distributions for protein A and B are illustrated in (B) and (E). In (C) and (F), the difference in accuracy of the proposed method and the kernel approach is shown for protein A and B. Shades of red indicate superior and white the inferior accuracy of the proposed method (all negative values are marked white). The difference in convergence is illustrated in (D) and (G) for protein A and B with a similar color code

4.2 Application to parameter optimization

Having shown that the proposed method yields convincing results, we utilize it for optimization of several biochemical reaction networks. In biology it is not possible to measure all parameters directly, which are necessary for computational modeling, leading to parameter estimation problems. The exact simultaneous simulation of intrinsic and extrinsic noise is computational extremely intense, wherefore approximate methods are needed. In this subsection, we use the proposed method to estimate the unknown rate constants of five different example systems described in Table 1. Therefore, six reference measurements computed with the SSA and MC sampling with equal time spacing were used for comparison with the results of our approximate method. For the example of gene expression only protein A and B are utilized for optimization, whereas for the other examples all chemical species were used. The objective function was the sum of the eucledian distance between the pseudo reference and the approximate solution for all chemical species and time points. The optimization was performed with a genetic algorithm, since gradient based algorithms fail for stochastic systems and get easily stuck in local extrema (Poovathingal and Gunawan, 2010). For all example systems 500 generations with a population size of 40 were used. In Table 1 the true parameters k and the optimized ones k are provided. For all example systems, the optimized parameters are very close to the true parameters indicating that the proposed method yields excellent results and depicts a promising technique for optimization of stochastic biochemical reaction systems. In the Supplementary Material further comments regarding the benchmarking can be found.

Table 1.

Benchmark model description

Model	Initial conditions			ktrue	kopti
Gene Model:
gene_off⇄k_2k_1gene_on	NON	0	k1	10^−2	8.8·10^−3
gene_on→k_3∅gene_on+A	NOFF	1	k2	10^−3	6.6·10^−4
A+B→k_4∅	NA	0	k3	5·10^−1	5.2·10^−1
	NB	e^N(μ=10^3,σ=1.5·10^2)	k4	5·10^−7	5.2·10^−7
Schlögl Model (Schlögl, 1972):
2X+A⇄k_2k_13X	NX	e^N(μ=2.75·10^2,σ=5)	k1	3·10^−7	3.1·10^−7
B⇄k_4k_3X	NA	105 (const.)	k2	10^−4	10^−4
	NB	2·10^5 (const.)	k3	10^−3	10^−3
			k4	3.5	3.5
Michaelis–Menten-Kinetics (Michaelis and Menten, 1913):
E+S⇄k_2k_1ES	NE	2.5·10^2	k1	10^−4	10^−4
ES→k_3P	NS	N(μ=10^4,σ=10^3)	k2	5·10^−3	2.4·10^−3
	NES	0	k3	10^−1	10^−1
	NP	0
Virus Model (Gupta and Rawlings, 2014):
V→k_1G	NV	e^N(μ=50,σ=5)	k1	1.5·10^−1	1.4·10^−1
G→k_2G+M	NG	0	k2	2·10^−2	2.1·10^−2
G→k_32G	NM	0	k3	5·10^−2	5.1·10^−2
M→k_4M+P	NP	0	k4	1	9.9·10^−1
Yeast Model (Poovathingal and Gunawan, 2010):
PC_3→k_1PC_3+mRNA_G	NPC	e^N(μ=50,σ=10)	k1	5·10^−1	5.1·10^−1
mRNA_G→k_2∅	N_mRNA_G	0	k2	5	5.1
mRNA_G→k_3mRNA_G+yEGFP	NyEGFP	0	k3	1	1
yEGFP→k_4∅	N_mRNA_R	0	k4	2·10^−2	2.1·10^−2
∅→k_5mRNA_R	NTetR	0	k5	5·10^−1	5.3·10^−1
mRNA_R→k_6∅			k6	2·10^−1	2·10^−1
mRNA_R→k_7mRNA_R+TetR			k7	1	1
TetR→k_8∅			k8	2·10^−2	2.1·10^−2

5 Discussion

In this study, the problem of efficiently simulating biochemical reaction systems containing intrinsic and extrinsic noise is addressed. We propose a novel algorithm relying on the SP approach and the τ-leaping algorithm, which computes approximate solutions of the distributed CME. The method is benchmarked on several examples illustrating its computational benefit compared to others. Choosing only SPs deterministically the proposed method converges very fast to an approximate solution. The accuracy of our approximative approach is barely dependent on system size but rather on the complexity of the stochastic mapping itself. This can be seen when the estimated PDFs of the Schlögl model are compared to the Virus model (one versus four states). Therefore, our results are very likely to also apply for large reaction systems. The strength of the method i.e. fast convergency at reduced sample size, also introduces some of its weaknesses. Since only a few samples from the initial distribution are chosen it might be sampled too sparsely. Additionally the choice of PDF kernels is a priori not clear. Although the proposed method might not converge to the exact solution of the CME (see saturation behavior in Fig. 2D) it qualitatively describes important characteristics of the underlying distribution (see Fig. 4A, B and E, and also additional illustrations of approximate versus pseudo exact solution in the Supplementary Material).

The method may serve as a tool for rapid analysis and optimization of a stochastic system. For precise predictions, we suggest to use our method as a starting point for refinement by means of numerically demanding, but asymptotically exact methods. Although there is a lot of development making stochastic analysis tools with state of the art algorithms and fast implementations (Drawert ; Fan ; Somogyi ) available for a broad community the optimization of stochastic biochemical systems is hardly performed due to its challenging computational effort (Poovathingal and Gunawan, 2010). With this study, we demonstrated that the proposed method is capable of describing biochemical systems containing intrinsic and extrinsic noise and it represents a promising tool suited for optimization and analysis of distributed, stochastic biochemical reaction systems.

Click here for additional data file.