1Refining the Weighted Stochastic Simulation Algorithm Dan T. Gillespie a) Dan T Gillespie Consulting, 30504 Cordoba Pl., Castaic, California 91384 Min Roh and Linda R. Petzold Department of Computer Science, University of California Santa Barbara, Santa Barbara California 93106 Abstract: The weighted stochastic simulation algorithm (wSSA) recently introduced by Kuwahara and Mura [J. Chem. Phys. 129,165101 (2008)] is an innovative variation on the stochastic simulation algorithm (SSA). It enables one to estimate, with much less computational effort than was previously thought possible using a Monte Carlo simulation procedure, the probability that a specified event will occur in a chemically reacting system within a specified time when that probability is very small. This paper presents some procedural extensions to the wSSA that enhance its effectiveness in practical applications. The paper also attempts to clarify some theoretical issues connected with the wSSA, including its connection to first-passage time theory and its relation to the SSA. Version of: 20 Mar 2009 - To appear in Journal of Chemical Physics -2 I. I&TRODUCTIO& The weighted stochastic simulation algorithm (wSSA) recently introduced by Kuwahara and Mura1 is an innovative variation on the standard stochastic simulation algorithm (SSA) which enables one to efficiently estimate the probability that a specified event will occur in a chemically reacting system within a specified time when that probability is very small, and the event is therefore “rare”. The difficulty of doing this with the standard SSA has long been recognized as a limitation of the Monte Carlo simulation approach, so the wSSA is a welcomed development. The implementation of the wSSA described in Ref. 1 does not, however, offer a convenient way to assess the accuracy of its probability estimate. In this paper we show how a simple refinement of the original wSSA procedure allows estimating a confidence interval for its estimate of the probability. This in turn, as we will also show, makes it possible to improve the efficiency of the wSSA by adjusting its parameters so as to reduce the estimated confidence interval. As yet, though, a fully automated procedure for optimizing the wSSA is not in hand. We begin in Sec. II by giving a derivation and discussion of the wSSA that we think will help clarify why the procedure is correct. In Sec. III we present our proposed modifications to the original wSSA recipe of Ref. 1, and in Sec. IV we show how these modifications allow easy estimation of the gain in computational efficiency over the SSA. In Sec. V we give some numerical examples that illustrate the benefits of our proposed procedural refinements. In Sec. VI we discuss the relationship between the wSSA and the problem of estimating mean first-passage times, using as an example the problem of spontaneous transitions between the stable states of a bistable system. In Sec. VII we summarize our findings and make an observation on the relationship between the wSSA and the SSA. II. THEORETICAL U&DERPI&&I&GS OF THE wSSA We consider a well-stirred chemical system whose molecular population state at the current time t is x. The next firing of one of the system’s M reaction channels 1, ,MR R… will carry the system from state x to one of the M states j+xν ( 1, , )j M=…, where jν is (by definition) the state change caused by the firing of one jR reaction. The fundamental premise of stochastic chemical kinetics, which underlies both the chemical master equation and the SSA, is that the probability that an jR event will occur in the next infinitesimal time interval dt is ( )ja dtx, where ja is called the propensity function of reaction jR. It follows from this premise that (a) the probability that the system will jump away from state x between times tτ+ and t dτ τ+ + is 0( )0( ) eaa dττ− xx, where 01( ) ( )Miia a=≡∑x x, and (b) the probability that the system, upon jumping away from state x, will jump to state j+xν, is 0( ) ( )ja ax x. Applying the3multiplication law of probability theory, we conclude that the probability that the next reaction will carry the system’s state to j+xν between times tτ+ and t dτ τ+ + is { }0( )00( )Prob in ( , ) ( )e( )jajat t d a daττ τ τ τ−→ + + + + = ×xxx x xxν. (1) In the usual “direct method” implementation of the SSA, the time τ to the next reaction event is chosen by sampling the exponential random variable with mean 01 ( )ax, in consonance with the first factor in Eq. (1), and the index j of the next reaction is chosen with probability 0( ) ( )ja ax x, in consonance with the second factor in Eq. (1). But now let us suppose, with Kuwahara and Mura1, that we modify the direct method SSA procedure so that, while it continues to choose the time τ to the next jump in the same way, it chooses the index j, which determines the destination j+xν of that jump, with probability 0( ) ( )jb bx x, where {}1, ,Mb b… is a possibly different set of functions from {}1, ,Ma a…, and 01( ) ( )Miib b=≡∑x x. If we made that modification, then the probability on the left hand side of Eq. (1) would be ()0( )0 0( )e ( ) ( )aja d b bττ−×xx x x. But we observe that this “incorrect” value can be converted to the “correct” value, on the right hand side of Eq. (1), simply by multiplying by the factor 00( ) ( )( )( ) ( )jjja awb b=x xxx x. (2) So in some sense, we can say that an j→ +x xν jump generated using this modified procedure, and accorded a statistical weight of ( )jwx in Eq. (2), is “equivalent” to an j→ +x xν jump generated using the standard SSA. This statistical weighting of a single reaction jump can be extended to an entire trajectory of the system’s state by reasoning as follows: A true state trajectory is composed of a succession of single reaction jumps. Each jump has a probability (1) that depends on the jump’s starting state, but not on the history of the trajectory that leads up to that starting state. Therefore, the probability