Efficient stochastic sampling of first-passage times with applications to self-assembly simulations

Models of reaction chemistry based on the stochastic simulation algorithm (SSA) have become a crucial tool for simulating complicated biological reaction networks due to their ability to handle extremely complicated reaction networks and to represent noise in small-scale chemistry. These methods can, however, become highly inefficient for stiff reaction systems, those in which different reaction channels operate on widely varying time scales. In this paper, we develop two methods for accelerating sampling in SSA models: an exact method and a scheme allowing for sampling accuracy up to any arbitrary error bound. Both methods depend on analysis of the eigenvalues of continuous time Markov model graphs that define the behavior of the SSA. We demonstrate these methods for the specific application of sampling breakage times for multiply-connected bond networks, a class of stiff system important to models of self-assembly processes. We show theoretically and empirically that our eigenvalue methods provide substantially reduced sampling times for a wide range of network breakage models. These techniques are also likely to have broad use in accelerating SSA models so as to apply them to systems and parameter ranges that are currently computationally intractable.


I. INTRODUCTION
Stochastic simulation methods have become increasingly widespread as a means of simulating and analyzing biochemical reaction kinetics 1 .The chemical master equation, which governs the reaction kinetics for well-mixed systems, forms the basis for the stochastic simulation algorithm (SSA), proposed by Gillespie 2,3 .SSA models a reaction system as a Continuous Time Markov Model (CTMM) in which states of the system are defined by counts of reactants present at a given point in time and transitions between states correspond to individual reaction events.This SSA approach is valuable in part because it provides a model of reaction noise, which can become significant for reaction networks on cellular scales 4 .Furthermore, SSA models can provide significant computational advantages over continuum models for networks characterized by extremely large sets of possible reaction intermediates.
The computational value of the SSA approach lies in the fact that for a large class of networks, the random walk visits only a small fraction of the state space before equilibrium is established.As a result, kinetics on complicated networks can be simulated "on the fly," requiring explicit construction of the CTMM network only in the immediate vicinity of those states visited on a given trajectory.This property is an essential requirement for any feasible simulation algorithm, since the size of the state space describing the master equation is astronomical even for modest system sizes.Successful applications of SSA include gene regulatory networks 4 and self-assembly of complicated structures, such as virus capsids 5,6 .Furthermore, the SSA approach has now been adopted by several approaches for whole-cell modeling 7,8 and modeling generic complex reaction networks 9,10 .
The relaxation time of the SSA can, however, be extremely sensitive to the transition rates controlling the reaction kinetics.A pure SSA model has difficulty with stiff reaction systems, i.e., those where important events occur in parallel on very different time scales.In such cases, a simulation can become bogged down by sampling fast events to the exclusion of the slow events.Hybrid discrete/stochastic models 11,12,13 can resolve this problem in some domains, but not when the fast reactions make use of too many intermediates to allow them to be modeled continuously.One important example of such a stiff reaction system is the breaking of bond networks, where individual bonds may break and repair repeatedly before a sufficiently large bond group is broken to fracture the network.Another form of stiff SSA network occurs near the critical concentration of a self-assembly system, where high-order I n v a r i a n t d e n s i t y nucleation events can be orders of magnitude slower than individual binding reactions.In these stiff systems, an SSA model can become "trapped" for many steps in a small subset of the state space, resulting in negligible simulation progress for long periods of time.
To understand these "trapped" systems, it is useful to consider the graph theoretic representation of the SSA method.An SSA model is represented by a graph in which each node corresponds to one possible state of the full model.Edges connect nodes whose states can be reached from one another by a single reaction event, e.g., two molecules binding to one another.At each simulation step, SSA considers only the immediate neighbors of the current state.As a result, the simulation is prone to traps that can result from irregularities in the invariant density of the embedded Markov chain (EMC) implemented by SSA for a given CTMM.For example, consider a 3-state CTMM represented by a simple path (Fig. 1(a)), where the backward transition rate a is much larger than the forward transition rate b.The average number of SSA steps to reach state 3 from initial state 1 is O(a/b) because once SSA visits state 2, it will jump to state 3 only a b/a fraction of the time.
Nodes 1 and 2 collectively define a trapped subgraph from which the model must escape.In general, for an N-path, where each forward rate b is smaller than the backward rate a, SSA takes O(a/b) N −2 steps to traverse the path (see Theorem III.1 for an analogous problem).
One way such a trapped subgraph can arise in a physical system is through models of the breakage of bond networks.Fig. 1(c) shows the graph arising from a model of the breakage of a three-cycle bond network, which behaves similarly to the 3-state CTMM by establishing a trapped inner graph of four states -the unbroken state and three states with a single broken bond -from which the model must escape to reach any broken network state.We can alternatively understand the trapping problem in terms of a probability landscape view of a reaction system.The SSA is sluggish whenever its equilibrium landscape is irregular, consisting of valleys and hills.The broader and deeper these are, the slower SSA becomes.
To overcome the presence of traps or landscape irregularity, we propose two non-local simulation algorithms that rely on the spectral decomposition of the Kolmogorov matrix (for a CTMM) or the transition matrix (for the Embedded Markov Chain (EMC)).These eigenvalues and their associated eigenvectors describe global modes of relaxation of the full graph or any of its sub-graphs.Since eigenvalues are global properties of a graph, spectral methods are much less sensitive to local landscape traps.These methods can be applied to quickly sample first passage times on small CTMM graphs such as those in Fig. 1 or to sample escape times from trapped subgraphs when the full model is prohibitively large.
Previous attempts at simulating rare events include the Forward Flux Sampling (FFS) technique of Allen et al. 14 and related methods 15,16 .The approach breaks a rare event into a series of relatively more probable stages and uses estimates of waiting times for the successive stages to develop an aggregate transition rate for the full event.This aggregate rate can then be used to approximate the first passage time density as a single exponential random variable.However, while the exponential tail dominates the density for stiff systems and is therefore a highly accurate approximation in many cases, the true probability density has a peak at short times followed by a mixed exponential tail.The methods developed in the present work, by contrast to the FFS-like methods, sample first passage times from the entire density to within arbitrary an error bound.Recently, another method called the slow-scale SSA was proposed by Cao et al. 17,18 , which relies on a technique called the Partial Equilibrium Approximation (PEA).PEA essentially assumes that the set of fast reactions are always in equilibrium and the method approximates transition rates between slowly varying reactant species by their expected value in the partial equilibrium state.While these methods can provide significant benefits for some CTMMs, there are several limitations in using PEA or similar approximations for arbitrary graphs.First, a clear distinction between fast and slow species may not be obvious in a given problem.For example, in rule-based simulation of bond networks, stiffness is built in through the association/dissociation rates of individual bonds rather than being species dependent.Secondly, these methods always need to be supplemented with approximations involved in computing the mean values of the reaction propensities.Furthermore, PEA will be inaccurate whenever fluctuations in the reaction propensities within the partial equilibrium state are comparable to their mean values.
The goal of the present work is to develop efficient methods for some important classes of stiff SSA model for which the above techniques are unsuitable, with a particular emphasis on models important to simulations of self-assembly reactions.The methods proposed in this paper can be applied to Markov processes on arbitrary graphs.Furthermore, they can be made accurate to within arbitrary error bounds.The remainder of this paper is organized as follows: Section II A sets up some basic notation and a description of the sampling problem for general CTMM.In Section II B we introduce a spectral method which relies on the eigen decomposition of the master equation describing the CTMM.We use a complete spectral decomposition of the first passage time density and rejection sampling to return sample first passage times for arbitrary CTMMs.In section II C we introduce another spectral method which works as a hybrid between the purely local SSA and the completely nonlocal Master Equation method.The latter method proceeds by adaptively constructing a basis in which to simulate the Markov chain until the system state has relaxed to its slowest eigenvector.If first-passage out of the trapped subgraphs does not occur by that time, we use the appropriate eigenvalue to sample the time to first-passage as an exponential random variable.In section II D we introduce a method for automated discovery of trapped regions in stiff Markov models.This technique allows efficient implementation of spectral methods for large state spaces by isolating regions repeatedly visited by a given random trajectory and using spectral sampling to escape any such subgraph.In section III we present theoretical results on the time complexity of SSA for bond networks followed by experiments on some special classes of bond networks to compare the simulation efficiency of each method discussed.In Section IV we evaluate the automated discovery variants of the method by applying them to models of a nucleation-limited assembly system with a state space too large to explicitly construct.Section V concludes the paper with a discussion of results and directions for future research.

A. The chemical master equation and the stochastic simulation algorithm
The SSA identifies reaction kinetics for networks of biochemical subunits as a Markov process governed by an appropriate Chapman-Kolmogorov equation or, equivalently, its differential version -the master equation.Let S = {1, 2, . . ., N S } be the state space for the CTMM, each node representing a possible state for the simulated system.The time evolution of probability densities is governed by a Kolmogorov matrix W , which specifies the transition rates W nm from the state m to n.
where, p n (t) denotes the probability to be in state n at time t.The matrix elements W nm satisfy two necessary conditions: 1. W nm ≥ 0 for n = m.
Under these conditions, it is well known that the matrix has a steady state solution |Π = n π n |n that is an eigenvector of W with eigenvalue zero and that all initial distributions relax to |Π in the limit of long times 19 .In addition, we will require W to satisfy the detailed balance condition, which states that at equilibrium, the sum of probability current exchanged between any pair of states (n, m) is zero, i.e., W nm π m = W mn π n .This in turn allows one to define a scalar product on the state space such that W is self-adjoint: This condition ensures that we can construct an orthogonal eigenbasis and compute time evolved versions of any given initial probability distribution using spectral decomposition.

B. Spectral Sampling 1: Master Equation approach
Given a Kolmogorov matrix W on a state space S and an arbitrary initial state i ∈ V ⊂ S, the first-passage time T F (i) is a random variable which gives the time at which the trajectory first reaches any state in some subset of the state space F = S − V .The standard method of solving a first passage problem is to set up the master equation for V with an absorbing boundary over F (zero Dirichlet boundary condition) 19 .Let P V be a projection operator onto the subspace V and let N be the cardinality of V .Then, M = P V W P V is the effective Kolmogorov matrix that governs time evolution over V .From detailed balance, M is selfadjoint over L 2 π −1 .Hence, the eigenvectors of M form a complete basis {|ψ α }.A consequence of the spectral theorem is the completeness relation for the properly normalized eigenbasis, i.e., ψ α |ψ β = δ αβ .Given any vector |η : 1. Spectral decomposition of the first-passage time distribution In terms of the vertex set basis, the completeness relation over L 2 π −1 is I = n∈V P n , where P n = π n |n n| is the projector onto vertex state |n .Given an initial probability density p n (t = 0) = δ ni the probability for state n ∈ V evolves as: The transition to an element f ∈ F outside of V , is governed by the following equation: The probability for a first passage to the state f between t and t + dt is hence given by ρ(T f = t)dt = N α=1 c α,f e −λαt dt.

Exact sampling for the first-passage time distribution
In this section we describe a method for returning a sample time from the computed firstpassage density ρ(t) = N i=1 c i e −λ i t to any state f ∈ F .A general method for sampling from complicated distributions is to use the method of rejection sampling, which first chooses a random variable from a convenient envelope density and accepts or rejects the sample based on a second random sample that depends on the tightness of the envelope fit.The rejection rate is low if the envelope curve closely approximates the given curve.A simple envelope curve is provided by a pure exponential of the most slowly decaying eigenvalue, with a coefficient equal to the sum of all positive terms c i >0 c i in the computed density ρ(t).However, there is no guarantee that the rejected part is small.Since each eigen mode encloses an area c i /λ i , cancellations between near-degenerate eigenvalues can in principle lead to a high rejection ratio.We therefore present a method for choosing an envelope curve g(t) which eliminates these cancellations.Furthermore, in section III B we show that the envelope curve is exact for bond networks generated by cycle graphs C N .We sample from g(t) using a decomposition into a discrete mixture of densities f α (t) = d α (e −λαt −e −λ α+1 t ) and an efficient rejection step.Here d α are constants, one for each component f α of the envelope curve g(t).The next theorem proves that the density f α (t) can be sampled efficiently using a rejection method.
Theorem II.1.The expected rejection ratio for f α (t) is bounded from above by 1.5.
Proof.We will use a simple exponential h α (t) = Cλ exp[−λt] as the envelope function.In order to minimize C we choose h α (t) such that and where t * is defined implicitly by the condition These constraints yield a unique solution t * = 2 The rejection ratio is given by: where x = λ α /λ α+1 ∈ (0, 1).To upper-bound C, note that the exponent increases monotonically with x and its maximum is As a final comment, we note that for general graphs the average time complexity of this algorithm is dominated by the computation of the eigenvectors and eigenvalues, which gives us the following theorem: Theorem II.2.The average time complexity for spectral decomposition of the master equation is O(N 3 ) for a graph of N vertices 20 .

C. Spectral Sampling 2: Modified embedded Markov chain method
The efficiency of the SSA is dependent on the relaxation time of the embedded Markov chain (EMC).We use this observation to modify the basis in which the EMC is simulated.
The standard method of executing a random walk is to consider the transition between adjacent states, each of which is localized at a vertex of the CTMM.However, correct simulation only requires that these states form a basis, not that they are orthogonal.If we can choose a set of states which are increasingly likely to appear during the simulation of the Markov chain, we are unlikely to make repeated visits to the same state.In order to identify such a basis starting from an initial state |i , we first identify the transition matrix for an embedded Markov chain that correctly describes the given CTMM.Consider the vertex set V = {1, 2, . . ., N} and the basis constructed from V , B = {|1 , |2 , . . ., |N }.
At any given time t, let the state of the time-evolved Markov chain be |ψ(t) = N i=1 ψ i |i .Let V t = {i ∈ V |ψ i = 0} be the vertex subset populated by the current state vector.We construct the EMC for the subgraph induced by V t at each step of the algorithm.Given the projection of the Kolmogorov matrix M over the vertex set V , choose r = max(−M ii |ψ i = 0) to be the effective rate of transition to the next state and choose an exponentially distributed . . .random time step τ with mean waiting time 1/r.Then, L t = −(1/r) * M is the Laplacian governing the EMC and Q t = I − L t is the effective transition matrix at that time step.
The next state vector is chosen to be |φ = Q t |ψ(t) .The reason for choosing this particular value of r is to ensure that no term in Q t becomes negative, a necessary condition for a transition matrix.
Theorem II.3.The choice of next state is consistent with the master equation governing the CTMM.
Proof.Rewrite the master equation in terms of the {|ψ , |φ } basis (where the other N − 2 linearly independent basis vectors can be chosen arbitrarily): Since there is a unique decomposition for any vector in terms of a linearly independent basis set, Eq. 10 proves that starting from |ψ the next state is uniquely determined to be |φ .
In general, the next state |φ will have a total probability  the simulation.This sequence will continue until the state has relaxed to its slowest eigen vector |λ min , such that M|λ min = −λ min |λ min to within a user-defined relative error ǫ.
Once that state is achieved, we just need one more exponentially distributed random sample time τ with mean 1/λ min to escape the network.
SSA chooses a stochastic trajectory by sampling both the next neighbor and the time for the next step at random.The EMC method, on the other hand, evolves deterministically in our modified basis and only the time between transitions is stochastic.At each time step, transition to the absorbing boundary states is governed by the matrix elements connecting each of the transient states to the absorbing boundary.The advantage of such an approach is that it allows us to automatically compute the most slowly decaying eigenvector during the simulation.For completeness, we note the following result: Theorem II.4.For a graph of degree bounded by d and V of cardinality N, each step of this algorithm takes O(N * d) time.

D. Automated discovery of trapped subgraphs
As previously mentioned, stiffness in Markov model graphs results from repeated visits by a typical random trajectory to a small subset of vertices of the entire graph.Since the performance of spectral methods is sensitive to the size of the vertex set, it would prove useful if we could somehow identify these "trapped" subgraphs for stiff Markov models and apply spectral methods directly to those.In this section we present one such method, which we call "Automated Discovery" (AD) and which we show to be formally applicable to arbitrary bounded-degree graphs.
Let there be a state space S over which a CTMM is defined and consider a subgraph G(V, E) with vertex set V ⊂ S and edge set E. Starting from an initial state i ∈ V , we are interested in the time T F (i) to first passage out of V .Consider the subgraph H i induced by the vertex set U i ⊆ V visited by a random trajectory executing the SSA random walk before it escapes V and let N i = |U i | be the cardinality of U i .If T F i is the number of steps a SSA random walk takes to escape V , then a Markov model will show simulation stiffness whenever the expected values satisfy, SSA until S(N) ≤ C * N 3 for some constant C, to discover the trapped subgraph K and then use spectral sampling to escape the discovered graph.This alternative approach could be less efficient in some circumstances, but would guarantee that the overhead for spectral sampling is no more than a constant factor beyond that of the standard SSA.Fig. 5 shows the pseudocode for implementing AD for a given graph by this method.The algorithm generates a sample trajectory using SSA till such time that the trajectory spends O(N 3 ) steps within a trapped graph K of vertex set cardinality N = |K|.Then either of the spectral methods described in section II B or II C are used to sample the first passage outside K, to a vertex i.In general the state of the system at the time of first passage outside K will be a discrete probability mixture of more than one vertices.In such cases, the vertex i is randomly selected in accordance with the appropriate probability weight.The algorithm then resumes SSA execution over the enlarged graph K {i}.*E*Further investigation is, however, required to search for algorithms that may further improve the performance of AD.
In section IV B we prove that for at least one important class of graphs, namely models of chemically reacting species, we can indeed reduce the time complexity to its optimal value to within a constant factor, i.e., O(f (N i )).In order to validate the methods, we instantiate them for some specific challenging systems.We begin by demonstrating the non-AD variants of the methods for the problem of sampling the time required to break a network of bonds.This problem is an example of a stiff SSA on a generally small graph.It is also of independent interest because of its importance in modeling self-assembly processes on long time scales.Given such a system, we are interested here in the first passage time to the subset of states corresponding to disconnected graphs V b ⊂ S. Since each bond can occur in two states, intact or broken, a network of d bonds can be represented as a vertex on a unit hypercube in d dimensions.
The state space generated by the bond network before it becomes disconnected will usually be a truncated unit hypercube.An N-cycle C N generates the simplest non-trivial example, where the absorbing boundary is placed at all points on the hypercube at distance 2 from the fully-connected state.Fig. 1(c) illustrates this absorbing boundary for C 3 .Given a dbond network, we will represent the µ th bond-breaking rate by b µ and association or binding rate by a µ .It is convenient to represent a vertex on this hypercube by a binary d-tuple i = {i d , . . .i µ , . . .i 1 }, where i µ = 0 implies that the µ th bond is intact (see Fig. 1(c) for the graph corresponding to a trimer).From here on, we will use the notation μ = {δ dµ , . . ., δ 1µ } for the vector describing a state of the model with only the µ th bond broken.For such a graph, the time complexity of each SSA step is O(d).In the rest of this paper we will use this model of truncated hypercubes to represent bond networks.Morris and Sinclair 21 have proven that in the case of unweighted graphs, a random walk on a hypercube truncated by a hyperplane relaxes to equilibrium in polynomial time bounded by O(d) 9/2+ǫ for any ǫ > 0. However, as we have argued in the introduction, the mean hitting time, i.e., the number of random walk steps between a pair of vertices, can be extremely sensitive to the parameters governing the walk.We formalize this observation in the Theorem III.1 below, which bounds the expected number of SSA steps before the network is disconnected.Let Theorem III.1.The expected number of SSA steps required to break a k-connected network with k > 1 and r > 1 is Ω(r k−1 ).
A detailed proof of the theorem is provided in the appendix.Figs. 6 and 7 provide an empirical demonstration of the theorem.Fig. 6 analyzes the number of steps required in 100 trials of the SSA algorithm for simulating the breakage of a set of cycle graphs C N ranging in size from three to seven.Each model was examined using ratios of forward to backward rate from 1 to 20 in increments of 1. Breakage times for the cycle graphs increase linearly with rate ratio, although they also fall monotonically with cycle size (Fig. 6(a)).Fig. 7 analyzes the number of steps required to break k-connected hypercube graphs of dimensions k = {2, 3, 4, 5}.Fig. 7(a) also shows that the slope of a log-log plot approaches the predicted exponent k − 1. Fig. 6(b) and 7(b) suggest why a spectral approach might be effectiveas the reaction rate increases, steps to first passage behave more like a geometric random variable (as mean → ∞, standard deviation → mean), as expected for a slowly decaying eigen mode of the transition matrix.More detailed explanations of the simulation protocol for these figures is provided in section III C. We illustrate the Master Equation spectral method using the cycle graph C N as an example.This is a graph of N vertices and N edges, connected together in a loop such that exactly two edges need to be removed to disconnect the graph (called a separation pair ).The state space is S = {0} N µ=1 { μ} N µ=2 ν<µ { μ + ν}.In this case, the subspace V b = N µ=2 ν<µ { μ+ ν} defines the absorbing boundary and the subspace V c = S −V b defines the space of transient states.We begin with the most general form for M, the projection of In what follows, we assume that all the eigenvalues of M are negative (as they must be over the subset of transient states since n M nm ≤ 0 ensures any positive probability density decays to zero) and that the set of rates {a i } and {b j } are positive (ensured by property 1 of W ). For economy of notation, let us define k n = a n + m =n b m .Also, in what follows we assume that the bond indices have been labeled such that In the case of a C N network, the T f (i) distribution can be efficiently sampled due to certain properties of the eigenvalue distribution and the form of the eigenvectors.Since the sampling technique for a general CTMM will be an extension of this special case, it will be helpful to illustrate the method by investigating the spectral properties of C N .The next few results establish bounds on the eigenvalues of M as a special case of the interlacing eigenvalue theorem 22 .
12 satisfy the following: If n such diagonal elements are identical then the eigenvalue is (n − 1)-fold degenerate.

There is at least one eigenvalue of M in the interval
Proof.The eigenvalue condition Det|M − λI| = 0 implies that the eigenvalues λ are the zeroes of an (N + 1) th order polynomial: We establish bounds on the roots by calculating the sign of f (λ) over the set of points 1.Each term inside the summation sign in f (λ) contains n − 1 factors of (k i + λ).Hence −k i is an (n − 1)-fold degenerate eigenvalue.In what follows we assume that the remaining k j are all distinct.
2. The sign of the function f (λ) at λ = −k i is (−1) i .Hence ǫ i encloses at least one root of f (λ).
The eigenvectors of M {|ψ α } are mutually orthogonal for the set of non-degenerate eigenvalues.In the case of non-degenerate eigenvalues (−λ α = −k m ), these eigenvectors are: where N α is a normalization constant.For degenerate eigenvalues, an orthogonal basis can always be chosen using the Gram-Schmidt procedure.As will become apparent later, however, these eigenvalues do not contribute to the sampling in the case of a first-passage problem beginning with the unbroken loop i = 0.
Theorem III.3.The envelope curve g(t) defined by our method is identical to the first passage density for a C N network.
Proof.Beginning with an unbroken state at t = 0, the probability the model occupies a given state n at time t is given by: where ψ α,i is an eigenvector of M with eigenvalue −λ α .Note that only those λ α = k i contribute, for otherwise ψ α,0 = 0. Assuming λ α < λ α+1 , the coefficients satisfy c α,n < 0 for α > n.Since the partial sum S N,n = N α=0 c α,n = N α=0 π n ψ α,n ψ α,0 = π n n|0 = 0 , all other partial sums satisfy S β,n = β α=0 c α,n ≥ 0. These observations provide a means of decomposing the probability density into the following discrete mixture with positive coefficients: Since b n > 0, the combined rate of decay to any one of the broken states is given by: where, and C. Simulation models used for bond networks Although our methods can in principle sample escape times from any subnetwork of a CTMM graph, we have validated them here for the specific case of breaking networks of bonds due to the importance of this problem for self-assembly modeling.In rule-based models of self-assembly, a simulation is initialized with a set of assembly subunits, each with a complement of pre-specified binding sites.As the simulation progresses, the system evolves into a state with an assembly of disjoint networks.The binding interactions between two disconnected pieces of the network usually occurs on a slower scale than individual bond breaking reactions 6 .For bi-connected networks, however, the association rate within a connected network is much larger than the bond breaking rate since there is no entropy penalty in associating bonds between constituent subunits.Such models allow for a natural partitioning of the state space into subgraphs corresponding to the bi-connected components of the entire network.The first set of experiments that we performed were on such biconnected networks.The simplest non-trivial example of a bi-connected bond network is the graph generated by an N-cycle (C N ).More complicated networks of N bonds can be viewed as special cases of a truncated unit hypercube in N dimensions.We therefore carried out simulations for the network generated by C N as well as the full hypercube (Z N ).
Theorem III.1 guarantees that the expected number of SSA steps for a k-connected network of d bonds is P (d, k)Ω(r k−1 ), where P (d, k) is some combinatorial function dependent on the topology of the network.
Each model is parameterized by a rate of bond formation, a, and a rate of bond breaking, b.These values were varied in different simulations.Each of the bonds had different binding/breaking rates but the ratio was maintained at the same order of magnitude for each simulation.Specifically, for a d bond network b µ = b(1.0+ 0.05µ/d) and a µ = a.These slight variations in rates from bond to bond were used to avoid giving our methods an unfair advantage, as they will generally be more efficient when the transition matrix has degenerate eigenvalues.

D. Experiments
We conducted a series of simulations to determine the performance of the SSA, Master Equation, and EMC methods for bond network first-passage times.All simulations were implemented in Mathematica.Run time simulations were executed on a Macintosh machine with a 1.8GHz G5 processor and 512 MB RAM.For the EMC based spectral method, we allowed each component of the state vector to converge within a relative error of ǫ = 0.01.
Each data point reported was the average over 500 simulations except for run time data, which were averaged over 100 simulations.
We first examined the efficiency of the Master Equation method by assessing the number of rejection steps needed to sample each first-passage time.We carried out simulations for cycle graphs (C N ) varying the cycle length from 3 to 7 and the rate ratio a/b from 1 to 20 in increments of 1.These experiments were then repeated for unit hypercubes (Z N ) with dimension varied from 2 to 5 and rate ratio a/b from 1 to 10 in increments of 1.For each condition, we recorded the number of rejection steps required for each of 500 simulations and computed the mean and standard deviation across the 500 trials.
We next examined the number of steps required by the EMC method for sampling times to network breakage.We examined the same models as those used to validate the Master Equation method: cycles of length 3 to 7 with rate ratios from 1 to 20 in increments of 1 and hypercubes of dimension 2 to 5 with rate ratios from 1 to 10 in increments of 1.We similarly recorded the number of EMC steps required for each of 500 simulations and computed the mean and standard deviation across the 500 trials.We also computed the fraction of models that reached the first passage time before relaxing to the slowest decay mode.
We next tested the total run time of each of the three methods on a broader set of parameter ranges.We evaluated run times for each method for cycle networks of sizes 3 through 7. We performed two sets of evaluations for each.The first set varied the rate ratio a/b from 500 to 5000 in increments of 500 to provide a broad view of the relative run times of the three methods.These numbers span ranges of values likely for protein assemblies.For example, Zlotnick et al. 23 have estimated a binding free energy of ∆G = 4.2 kcal/mole for ODE based simulation of the kinetics of the Hepatitis B virus, which yields a/b = exp (∆G/RT ) ∼ 1200.We then examined ratios of SSA to Master Equation and SSA to EMC run times for each data point based on averages over 100 simulations per parameter set.In a second set of experiments, designed to give a finer view of where each method is dominant in parameter space, we varied the rate ratio a/b from 30 to 300 in increments of 30.We then identified the most efficient of the three methods for each point, again using averaged run times over 100 trials per data point.
We then performed analogous experiments for hypercube graphs in order to test performance on networks with higher connectivity.For each graph Z 2 to Z 5 , we carried out simulations for rate ratio a/b from 3 to 30 in increments of 3. We were limited to small ratios because the SSA method becomes prohibitively costly for high-connectivity networks at higher ratios.Each simulation was repeated 100 times to yield average run times for each parameter set and for each of the three methods.For each parameter set, we computed the ratio of run times for SSA versus Master Equation and SSA versus EMC.We further evaluated which of the three methods produced the shortest average run time for each parameter set.

E. Results
We first present results on the efficiency of the rejection sampling scheme for the Master Equation method.The expected run time of the method is proportional to the expected number of trials needed to produce a successful sample.A low number of steps is therefore preferable, with a value of one being ideal.Fig. 8(a) shows the rejection ratio for cycle graphs C 3 through C 7 .The mean number of rejection steps is consistently below 1.5, as expected from theorem III.3 and II.1.The number of rejection steps drops with increasing rate ratio but increases with increasing cycle length.These results together establish the efficiency of the method.Fig. 8(b) shows that the method is also robust, with standard deviation consistently below 0.9 for the experiments shown here.The standard deviation also decreases with increasing rate ratio but increases with cycle size.and (c) provide the explanation for this feature.For small rate ratio, multiple eigen modes are responsible for the decay (see part (c)), which corresponds to increasing EMC steps before first passage, similar to SSA.However, as rate ratio increases further, relaxation time to the Fig. 11 shows comparable results for hypercube graphs.Fig. 11(a) shows that mean numbers of steps drop substantially between ratios 1 and 2 but quickly level off to an apparent constant for each graph.The number of steps increases with increasing hypercube dimension.Figs.11(b) and (c) again show that the method has high variability for low rate ratios, where multiple eigen modes contribute significantly to the time distribution and the method must behave similarly to the standard SSA.At higher ratios, though, the slowest mode quickly dominates and the number of steps required becomes highly reproducible.
We next examined total run times of the three methods, beginning with the cycle graphs method.The ratio grows rapidly with increasing rate ratio, although it falls with increasing cycle size.Fig. 12(b) shows the comparison of SSA to the EMC method.The SSA:EMC ratio likewise peaks for large rate ratios and small cycle sizes.The EMC method appears generally superior to the Master Equation method, beginning to dominate at a lower rate ratio and reaching a higher peak.Fig. 12(c) shows for a narrower rate range where each of the three methods dominates.The EMC method is the fastest for most of the range examined, with the standard SSA superior at the extreme of low ratios and large cycle sizes.
We then examined run times on the hypercube graphs Z 2 to Z 5 .Fig. 13(a) shows run time ratios for SSA versus the Master Equation method and Fig. 13(b) for SSA versus the EMC method.Both spectral methods show sizable improvements over the pure SSA method for larger rate ratios and higher hypercube dimensions.SSA appears much more sensitive to rate ratio as compared to the spectral methods.Even for a rate ratio of 30, the spectral methods were more than three orders of magnitude more efficient than SSA for Z 5 .
For hypercubes, unlike cycle graphs, the Master Equation method appears generally more efficient than the EMC-based method, even for small rate ratios.Fig. 13(c) shows where each method is dominant.The Master Equation method is dominant for most of the parameter range examined, with the SSA method superior at the limit of lowest degree and smallest rate ratios and the EMC dominant for low degree and higher rate ratios.This result is expected from Fig. 7(a), since the average number of steps seems to increase monotonically with the connectivity of the graph for the EMC-based method.The efficiency of the Master Equation method, on the other hand, depends primarily upon the size of the complete graph, since matrix diagonalization is the eventual efficiency bottleneck.In this section, we apply the AD variants of the methods to a different system type also motivated by self-assembly modeling.The rate of a self-assembly processes is often limited by the time required to build the first stable multi-subunit complex, called a nucleus, which then acts as a seed for assembly of the rest of a larger structure.Because partially formed nuclei are unstable, considerable trial-and-error may be needed before one reaches completion.The time to complete a single nucleus can thus be orders of magnitude longer than the inter-subunit binding rate.These nucleation-limited assembly systems are one example of the broader class of stiff models for chemically reacting species.The state space of any such a system can be represented as a lattice corresponding to the populations individual species.These models are similar to those treated in earlier studies of accelerated SSA methods 14,18 .We apply one such model, representing the formation of simple trimeric nuclei, to demonstrate and evaluate the AD variants of the spectral methods.

A. Integer lattice models
The second model we consider is again an assembly of bond networks where monomers (m) with two identical binding sites combine to form dimers (d) and trimers (t).In order to show stiffness with respect to a single parameter, the trimers were assumed to be completely where 1/v is an entropy penalty due to the finite volume of the system.We initialize the system at the state (0, 0) for a given monomer count N and sample the first passage time until the trimer count reaches a given value.This system will show stiffness if the parameter ρ ≡ N/v is small.For small ρ, which corresponds to low concentration and/or small binding energy, trimer formation will be much slower than dimer breaking/binding reactions.

B. Automated Discovery for integer lattice
Efficient simulation over an integer lattice, where one pair of species react on a much faster timescale than the others requires a partitioning of the entire lattice into subgraphs with fixed trimer count (since trimer formation occurs on a much longer timescale than monomer-dimer reactions).These subgraphs are simple paths with vertex set V (N t ) = FIG. 14: Pseudocode for Automated Discovery for a simple path {(N t , 0), (N t , 1), . . .(N t , [(N − 3N t )/2])}, where N t represents the fixed trimer count and square brackets represent the largest integer smaller than the enclosed expression.Fig. 14 presents a procedure for implementing automated discovery on such graphs which works in optimal time, to within a small constant factor (the vertices are represented by dimer count for simplicity).At each step of automated discovery the method enlarges the graph by a factor proportional to its present length.The scale factor m can be optimized for any given sampling algorithm to optimize run time.For example, if a given spectral sampling algorithm works in time f (N) = N α , where N is the cardinality of the vertex set; total number of steps used in automated discovery of a graph sized N i can be bounded from above by the following quantity S N i , where, C = 2 2α/α+1 /(2 α/α+1 − 1).In general the Master Equation based method works in O(N 3 ), however for simple paths the Kolmogorov matrix is sparse and effective power is expected to be more like α = 2. Since S N i is more sensitive to deviations for smaller values of m, we chose m = 1.3 > 2 1/3 in the experiments reported here.
The method reported here can in principle be generalized for arbitrary lattice graphs in d dimensions.The size of the discovered graph in such cases would overestimate the actual trapped graph by a factor of m d , for a scaling factor m.For small dimensions, this may still be more efficient than the method discussed in section II D which exactly samples the trapped subgraph.

C. Experiments
We performed two sets of experiments to compare the performance of spectral methods with SSA.The first set of experiments compared the Master Equation method implemented in conjunction with Automated Discovery for the trimer model with SSA.Each experiment compared the ratio of run times for sampling first passage times to reach a trimer count N t = 100, starting from an initial monomer count N. The state space was partitioned into subgraphs corresponding to fixed trimer counts and AD was used to identify the trapped regions for spectral decomposition.We then performed a total of 50 comparative run time simulations varying N from 1000 to 9000 in steps of 2000 and varying ρ from 10 −5 to 1.9 × 10 −4 in steps of 2.0 × 10 −5 .All run times were averaged over 50 samples.The scale factor for AD was set at 1.3.The second set of experiments compared the run time ratio for the EMC based method and SSA for first passage time to reach a trimer count N t = 100, starting from an initial monomer count N. The state space was again partitioned into subgraphs of fixed trimer counts.AD was not required for these simulations since the method automatically selects the trapped region of the subgraph according to the evolving probability distribution.We then performed a total of 25 comparative run time simulations varying N from 1000 to 9000 in steps of 2000 and varying ρ from 10 −4 to 0.9 × 10 −3 in steps the behavior for 5 different initial monomer counts N. The EMC method is effective at substantially larger values of ρ than is the Master Equation method.As is to be expected, the spectral methods do not scale as well as the usual SSA method with increasing N.Even for relatively large networks, though, the performance gain obtained by spectral sampling is appreciable.The reason for this is that for most cases the slowest eigen mode is reasonably well approximated by a vector populating only a small fraction of the subgraph vertices.As a result we can look at the EMC method as a generalization of other accelerated sampling schemes which only use one vertex, the mean value of the slowest decay mode, as in the PEA based methods.

V. DISCUSSION
We have investigated the problem of efficiently simulating stochastic reaction models and introduced two methods for accelerating sampling on problems characterized by multiple time scales.Both methods are based on spectral analysis of CTMMs equivalent to the SSA model.We have applied these methods in the present work to two special cases of these models that are important to simulations of molecular self-assembly: sampling times to break multiply-connected bond networks and simulating growth in nucleation-limited assembly systems.Collectively, these two applications demonstrate the use of the proposed spectral methods on small CTMM graphs known a priori and on automatically discovered subgraphs of large CTMMs.We have shown theoretically and empirically that the new methods are substantially more robust to variations in the ratios of reaction rates than is the basic SSA method for these problems.
While we have applied these methods here to models used in self-assembly simulations, the basic methods can be expected to have much broader application.Both methods can be applied to sample first passage times for any arbitrary subset of states of any SSA CTMM graph.Both can also be applied to sample escape times from any subgraph of such a graph, using automated discovery to identify "trapped" regions of the CTMM graph.The latter distinction is important because CTMM graphs for complicated biological systems are generally far too large to represent explicitly.These spectral methods might be extended to incorporate "on the fly" graph construction techniques, like those used by rule-based methods widely used for SSA simulations 9,10 .The EMC method, especially, would seem to be a candidate for such an extension.For example, if at each iteration, instead of adding all the possible next neighbors to the system state, we add only a subset of them depending upon their transition probabilities then we will get a natural, non-local generalization of the SSA.Such an approach could provide a precise and general method for pruning full SSA graphs to achieve more efficient pathway sampling in extremely large state spaces.
graph will then have the same hitting time T b0 as the original graph.
We next need three auxiliary results about properties of the resulting graph in order to prove our main theorem.
Lemma A.1.For an ergodic Markov chain, the cover time C ij ≡ T ij + T ji between any two states i and j satisfies 24 Proof.
1.The transition probability corresponding to the matrix element connecting i to j is : 2. First consider any state μ with one bond broken: Proof.In order to apply lemma A.1 to bound the hitting time we need to look at graphs with k odd or even separately.If k is even we can apply lemma A.1 directly to C 0b for the 2-step chain Q 2 even .However, if k is odd, we need to consider the cover time between b and each

FIG. 1 :
FIG. 1: Illustration of trapped subgraphs in SSA models.(a) A simple CTMM on a 3-path with transition rates a and b.(b) The probability landscape for the model.SSA is slow whenever the invariant density for the corresponding Markov chain is irregular.Here, the SSA takes O(a/b) steps to reach vertex 3. (c) CTMM model of a trimer assembly system with three subunits.Graph of possible configurations joined by reaction rates.States in which the trimer is broken are surrounded by solid lines and others by dashed lines.

FIG. 3 :
FIG.3: Schematic of the EMC-based spectral method.Vertices in black are the currently occupied nodes.Simulation advances the system state as a discrete mixture until such time as the state has relaxed to its slowest eigenstate |λ min .At each step, direct transitions to the absorbing vertex (grey) are computed according to the Kolmogorov matrix.

FIG. 6 : 2 sFIG. 7 :Log 10 − 2 sB. Master equation for bond networks 1 .
FIG. 6: Number of SSA steps until first passage for the network generated by an N-cycle C N .(a) Average number of steps s (b) Relative deviation √ δs 2 s

Fig. 9 (FIG. 8 :FIG. 9 :
Fig.9(a) shows mean numbers of rejection steps for hypercube graphs.Since the envelope curve for hypercubes is not exact, these experiments provide information about how well FIG. 10: Number of EMC steps until first passage for the network generated by C N (a) Average number of steps s (b) Standard deviation √ δs 2 (c) Fraction of times the trajectory escapes before relaxing to the slowest decay mode.
FIG. 11: Number of EMC steps until first passage on Z N (a) Average number of steps s (b) Standard deviation √ δs 2 (c) Fraction of times the trajectory escapes before relaxing to the slowest decay mode.

FIG. 12 :
FIG. 12: Comparative run times for the network generated by C N (a) Ratio of SSA to Master Equation run times (b) Ratio of SSA to EMC run times (c) Region in 2D parameter space where each method is optimal

FIG. 13 :
FIG. 13: Comparative run times for first passage on Z N (a) Ratio of SSA to Master Equation run times (b) Ratio of SSA to EMC run times (c) Region in 2D parameter space where each method is optimal

E
[C ij ] = E[T ij ] + E[T ji ] = 1/(π j P r[T jj > T ij ]) (A1)Lemma A.2.The transition matrix Q satisfies the following conditions:1.If i and j = i+ μ are two neighboring states with n and n+1 bonds broken, respectively, then Q j,i ≤ (n * r) −1 for any n > 0.2.For any initial state i containing n broken bonds, the n-step transition probability to0 is bounded from below by Q n 0,i ≥ (1 + d/r) −n .3.Let T be any stopping time for the transition matrix Q with expectation valueE[T ] = n=1 nP r[T = n] = n=0 P r[T > n].For any integer l > 1 consider the expectation value of T for the l-step transition matrix Q l defined as E (l) [T ] = n nP r[(n − 1) * l < T ≤ n * l] = n=0 P r[T > n * l].Then, l(E (l) [T ] − 1) ≤ E[T ] ≤ lE (l) [T ].

( 1 2 ( 1 +Theorem A. 1 .
μi ≥ (1 + d/r) −n for all n broken bond states n i=1 μi , then μi Q P µ i =ν μi , P n+1 i=1 μi≥ (1 + d/r) −n ν∈{µ 1 ,...,µ n+1 } a ν n+1 i=1 a µ i + η ={µ i } b η ≥ (1 + d/r) −n−1 (A3) Since Q 0, μ > (1 + d/r) −1, the assertion holds for all n > 1 by induction.wherewe have used Lemma A.2 part (2).Let us define p = (1 + d/r) −k .In terms of p, the previous inequality and Lemma A.2 part (3) implyP r[T 0b > n * k] = 1 − P r[T 0b ≤ n * k] ≤ (1 − p) n ⇒ E (k) [T 0b ] ≤ ∞ n=1 − p) n = 1/p E (k) [T 0b ] ≤ (1 + d/r) k ⇒ E[T 0b ] ≤ k(1 + d/r) k (A9)An immediate consequence of the previous lemma is that for k evenk/2 ≤ E (2) [T 0b ] ≤ (k/2)(1 + d/r) k + 1 (A10)Similarly, if k is odd, using the fact that P r[T 0b < T μb ] = 0 we getk − 1 2 ≤ E (2) [T μb ] ≤ E[T 0b ]/2 + 1 ≤ k d/r) k + 1 (A11)Let us define the equilibrium probability for b as πb , thenπb = π b i∈Veven π i if k is even = π b i∈V odd π i if k is odd (A12)We can finally compute upper (U) and lower (L) bounds on the hitting time T b0 which are asymptotically equivalent in the limit r → ∞.The following theorem implies that ∆(r) ≡ U(r) − L(r) is monotonically decreasing in r and lim r→∞ ∆(r)/L(r) = 025 .The expected number of SSA steps before first passage on a k-connected graph is bounded within 2 (1 − d/r) −1 − ((k − 1)/2) πb πb ≥ E[T b0 ] ≥ 2 1 − (k/2(1 + d/r) k + 2)π b πb (A13) ln(λ α+1 /λα) ∆λ .Since d 2 fα dt 2 > 0 for t > t * , the slope of ln[f α ] monotonically increases to −λ α as t → ∞.The corresponding envelope rate Algorithm:Prepare Discrete Mixture Input: First passage time probability density ρ Output: Envelope g , r}; Algorithm:Check Convergence Input: Next state |φ , present state |ψ , rate r Output: The first passage time t for i ∈ V do P ← P + φ i ; if |φ i − ψ i | > ǫψ i thenEigenMode ← No; end end Generate a uniform [0, 1] random variable U; if EigenMode = Yes then λ min ← r * (1 − P ) ; Return t ← t − 1 λ min ln U; ) since this would imply certain vertices in U i are being visited repeatedly.AD works by progressively sampling larger regions of V until it identifies a subgraph K i induced by a vertex set W i ⊆ V such that U i ⊆ W i .Once K i is identified either of the spectral methods can be used directly over K i .The method will be efficient as long as |W i | ∼ |U i | and the number of steps taken to identify K i is comparable to the computational cost of using spectral sampling over K i .Formally, H i can be exactly discovered by repeatedly enlarging the discovered graph to include the last vertex outside K i visited by the trajectory.If spectral sampling for a graph of vertex set size N works in time f (N), this procedure would ensure that implementing spectral sampling in conjunction with automated discovery takes O(N i * f (N i )) steps.The stiffness condition (Eqn.11) would usually ensure that this procedure is still efficient.
*B* Another method for discovering the trapped subgraph would be to implement the SSA random walk for a specified number of steps S(N)(depending on the size of the vertex set N). Since eigenvalue methods are in general O(N 3 ), we can implement Algorithm:Automated Discovery Input: Discovered vertex set K, current state i, time elapsed t Output: Final state vertex i / ∈ V , first passage time t stable.If the total number of monomer subunits is N, the state space is the intersection of the plane N m + N d + N t = N with the positive octant of the 3 dimensional lattice formed by integer counts of the monomer (N m ), dimer (N d ) and trimer (N t ) populations.Let us represent each vertex of this graph by the pair (N t , N d ).The reaction propensities α t , N d ) are (to within an overall constant) Algorithm:Automated Discovery Input: Discovered vertex set K, current state i, time elapsed t Output: Final state vertex i / ∈ V , first passage time t