Accuracy limitations and the measurement of errors in the stochastic simulation of chemically reacting systemsIntroductionBackgroundSSA and the tau -leaping methodSchl ouml gl reactionDistribution functionsKolmogorov distanceDensity distance areaSelf distanceNumerical experimentsConclusion and discussionAcknowledgmentsProof of Lemma 4.2Proof of Lemma 4.3ReferencesAccuracy limitations and the measurement of errors inthe stochastic simulation of chemically reacting systemsYang Cao*, Linda Petzold*,1Department of Computer Science, University of California, Santa Barbara, CA, United StatesReceived 11 November 2004; received in revised form 9 May 2005; accepted 21 June 2005Available online 18 August 2005AbstractThis paper introduces the concept of distribution distance for the measurement of errors in exact and approximatemethods for stochastic simulation of chemically reacting systems. Two types of distance are discussed: the Kolmogorovdistance and the histogram distance. The self-distance, an important property of Monte-Carlo methods that quantifiesthe accuracy limitation at a given resolution for a given number of realizations, is defined and studied. Estimation for-mulas are established for the histogram and the Kolmogorov self-distance. These formulas do not depend on the dis-tribution of the samples, and thus show a property of the Monte-Carlo method itself. Numerical results demonstratethat the formulas are very accurate. Application of these results to two problems of current interest in the simulation ofbiochemical systems is discussed. 2005 Elsevier Inc. All rights reserved.1. IntroductionIn microscopic systems formed by living cells, the small numbers of reactant molecules can result indynamical behavior that is discrete and stochastic rather than continuous and deterministic [1–9]. The sto-chasticity (often called biochemical noise by biologists) in microscopic systems has been implicated in thelysis/lysogeny decision of the bacteria k-phage [3] and the loss of synchrony of Circadian clocks [4].Tostudy the stochasticity in microscopic systems, engineered gene circuits have been designed and imple-mented in the laboratory. The effects of stochasticity have been observed in biological experiments [6–9].0021-9991/$ - see front matter  2005 Elsevier Inc. All rights reserved.doi:10.1016/j.jcp.2005.06.012*Tel.: +1 805 893 5362; fax: +1 805 893 5435.E-mail address: [email protected] (L. Petzold).1Department of Mechanical and Environmental Engineering, University of Santa Barbara, Santa Barbara, CA 93106-5070, UnitedStates.Journal of Computational Physics 212 (2006) 6–24www.elsevier.com/locate/jcpDetailed models [1,2,10,11] for the expression of a single gene and gene networks have been proposed toexplain these experi ments.The comparison between models and experiments is verified through Monte-Carlo simulations. Suchsimulations are based on GillespieÕs stochastic simulation algorithm (SSA) [12,13], which yields an exactstochastic simulation for well-stirred chemically reacting systems. However, it can be prohibitively expen-sive for realistic biochemical simulation. Approximate simulation methods have been proposed, such as thes-leaping method [14], implicit s-method [15] and hybrid methods [16,17]. These methods can achievegreater efficiency and give a close approximation to the SSA method. Two natural questions are concernedhere. One is: How should we measure the difference between the experiments and the Monte-Carlo simu-lations? The other is: How should we measure the accuracy of approximate methods? We seek a quantita-tive measurement. Typically, what is of interest from an experiment or a simulation are the stochasticproperties of the solution variables or of some function of the solution values, as opposed to the valuesfrom one simulation. One possibili ty for measuring the error is to compute the errors in solution momentssuch as the mean and variance. However, often the problems for which stochastic sim ulation makes a bigdifference have a bistable distribution. A simple example of this type is given by the Schlo¨gl [18] reaction.Some well-known problems from biology [3,7,19] also have this property. For such problems, the low-ordermoments such as mean and variance do not have much relevance. Rather, we need to know how well thedetailed model capture s the probability distribution, or how well approximate methods capture the ÔexactÕprobability distributions for the variables and properties of interest. To describe this error, in this paper weadopt the concept of distribution distance. Two types of distribution distance are considered here. One is theKolmogorov distance [20,21], defined to measure the distance between cumulative distribution functions(cdf). The other is the density distance area, defined as the L1distance between the probability density func-tions (pdf).With the distribution distance, we can measure the accuracy of the distributions given by differentapproximation formulas. There are very few systems for which we can analytically solve for the distri-bution. Thus, for most problems, we must collect a large number of samples by experiments or by sim-ulation methods such as SSA. The distribution is then estimated by the empirical distribution function[22] (edf), or the hist ogram of the samples. For a sufficiently large number of samples, the edf is close tothe cdf, while the histogram is close to the pdf. However, in both experiment and computation thenumber of samples is always limited. Thus, such a process is subject to an inherent error due to therandomness of the variables of interest, which is traditionally called ‘‘statistical fluctuation’’ in the liter-ature of Monte-Carlo simulations. In practice, we can measure this statistical fluctuation by the distancebetween two sets of independent samples with the same distribut ion. We call it ‘‘self distance’’. Thisconcept is similar to that of the round-off error in classical numerical analysis. In that context, dueto the limited length of the storage for each variable, there is a ‘‘round-off error’’. In stochastic simu-lation a simple limitation also exists, on the number of samples. For example, if there are two setsof samples for which the distribution distance between them is smaller than their self distance, we can-not tell whether or not these two sets of samples represent different distribut ions unless we increa se thenumber of