The pathway elaboration method for mean first passage time estimation in large continuous-time Markov chains with applications to nucleic acid kinetics
Sedigheh Zolaktaf, Frits Dannenberg, Mark Schmidt, Anne Condon, Erik Winfree
SSubmitted
THE PATHWAY ELABORATION METHOD FOR MEAN FIRST PASSAGETIME ESTIMATION IN LARGE CONTINUOUS-TIME MARKOV CHAINSWITH APPLICATIONS TO NUCLEIC ACID KINETICS B Y S EDIGHEH Z OLAKTAF , F
RITS D ANNENBERG , M
ARK S CHMIDT , A
NNE C ONDON , AND E RIK W INFREE University of British Columbia, * [email protected]; † [email protected]; ‡ [email protected] California Institute of Technology, § [email protected]; ¶ [email protected] Alberta Machine Intelligence Institute
Continuous-time Markov chains (CTMCs) are widely used in many ap-plications, including in modeling nucleic acid kinetics (non-equilibrium dy-namics). A typical issue in CTMCs is that the number of states could belarge, making mean first passage time (MFPT) estimation challenging, par-ticularly for events that happen on a long time scale (rare events). We proposethe pathway elaboration method, a time-efficient probabilistic truncation-based approach for detailed-balance CTMCs. It can be used for estimatingthe MFPT for rare events in addition to rapidly evaluating perturbed pa-rameters without expensive recomputations. We demonstrate that pathwayelaboration is suitable for predicting nucleic acid kinetics by conductingcomputational experiments on 267 measurements that cover a wide rangeof rates for different types of reactions. We utilize pathway elaboration togain insight on the kinetics of two contrasting reactions, one being a rareevent. We compare the performance of pathway elaboration with the stochas-tic simulation algorithm (SSA) for MFPT estimation on 237 of the reactionsfor which SSA is feasible. We further use pathway elaboration to rapidlyevaluate perturbed model parameters during optimization with respect toexperimentally measured rates for these 237 reactions. The testing erroron the remaining 30 reactions, which involved rare events and were notfeasible to simulate with SSA, improved comparably with the training er-ror. Our framework and dataset are available at https://github.com/DNA-and-Natural-Algorithms-Group/PathwayElaboration .
1. Introduction.
A continuous-time Markov chain (CTMC) is a stochastic process ona discrete set of states that have the Markov property, so that future possible states are in-dependent of past states given the current state. The time in a state before transitioning toanother state, the holding time, is continuous; to retain the Markov property, holding timesfollow an exponential distribution with a single rate parameter for each state-to-state transi-tion. CTMCs are widely used in natural and physical sciences, such as for modeling nucleicacid reactions (Schaeffer, Thachuk and Winfree, 2015; Flamm et al., 2000; Dykeman, 2015),protein folding (McGibbon and Pande, 2015), chemical reaction networks (Anderson andKurtz, 2011), and molecular evolution (Liò and Goldman, 1998). Particularly, in elemen-tary step models of nucleic acid kinetics (non-equilibrium dynamics) (Schaeffer, Thachukand Winfree, 2015; Flamm et al., 2000; Dykeman, 2015), which has inspired this work,states correspond to secondary structures and a transition between two states corresponds tothe breaking or forming of a base pair. The transition rates are specified with kinetic mod-els (Metropolis et al., 1953) along with thermodynamic models (Hofacker, 2003; Zadeh et al., * Joint first authorship for first two authors: Sedigheh Zolaktaf and Frits Dannenberg contributed equally tothis work.
Keywords and phrases: continuous-time Markov chain, mean first passage time, nucleic acid kinetics. a r X i v : . [ q - b i o . Q M ] J a n I F (a) Pathway construc-tion
I F (b) State elaboration
I F (c) Transition con-struction
I D (d) δ -pruning I D (e) Updating per-turbed parameters (f) Parameter estimation (g) Obtain functionality
Fig 1: The pathway elaboration method and its applications. Pathway elaboration makespossible MFPT estimation of rare events and the rapid evaluation of perturbed parameters.Here, in the underling detailed-balance CTMC, boxes in a square grid represent states of theCTMC, with transitions between adjacent boxes, initial state I at bottom left and target stateF at top right. ( a ) From state I, sample paths that are biased towards the target state F. Threesampled paths are shown with blue, pink and purple dotted lines. ( b ) From each sampledstate found in the previous step, run short unbiased simulations to fill in the neighborhood.Simulations from two states are shown with green dashed lines. The green states and transi-tions are sampled. ( c ) Include all missing transitions between the states that were sampled insteps a and b. The red transitions are included. ( d ) Prune states that are expected to reach thetarget state quickly by redirecting their transitions into a new state. ( e ) For perturbed modelparameters, keep the topology of the truncated CTMC, but update the transition rates. ( f ) Wecan use truncated CTMCs for perturbed parameters, such as to estimate model parameters or( g ) to predict forward ( k + ) and reverse ( k − ) reaction rate constants as temperature changes.2011) that make equilibrium predictions. Predicting nucleic acid kinetics is desirable forbuilding nanoscale devices whose nucleic acid sequences and experimental conditions needto be carefully designed to control their behaviour, such as RNA toehold switches (Angenent-Mari et al., 2020) and DNA robots (Thubagere et al., 2017). It would facilitate the design ofcomplex molecular devices by reducing, though not eliminating, the need for debugging de-ficiencies with wet-lab experiments.A central quantity of interest in CTMCs is the mean first passage time (MFPT) to reacha set of target states starting from a set of initial states. The MFPT is commonly used toestimate the rate of a process (Schaeffer, 2013; Reimann, Schmid and Hänggi, 1999; Singhal,Snow and Pande, 2004). For a CTMC with a reasonable state space size, a matrix equationcan provide an exact solution to the MFPT (Suhov and Kelbert, 2008). Direct application ofmatrix methods is not feasible for CTMCs that have large state spaces.For such cases researchers may resort to stochastic simulation (Gillespie, 2007; Ripley,2009; Asmussen and Glynn, 2007; Doob, 1942; Gillespie, 1977). For CTMCs in particu-lar, the stochastic simulation algorithm (SSA) (Doob, 1942; Gillespie, 1977) is a widelyused Monte Carlo procedure that can numerically generate statistically correct trajectories. HE PATHWAY ELABORATION METHOD In brief, the algorithm iteratively samples states based on transition rates and holding timeof states. By sampling enough trajectories that reach a target state, an estimate of the MFPTcan be obtained. In turn, direct simulation is inefficient for estimating MFPTs of events thathappen on a long time scale (rare events), such as reactions that involve high-energy bar-rier states. A number of techniques have been developed for efficient sampling of simulationtrajectories relevant to the event of interest (Shahabuddin, 1994; Bolhuis et al., 2002; Allen,Valeriani and ten Wolde, 2009; Escobedo, Borrero and Araque, 2009; Rubino and Tuffin,2009). Such sampling techniques, unfortunately, must generally be re-run if model parame-ters change, which makes it costly to perform parameter scans or to optimize a model.CTMC systems may be treated by methods that truncate the state space to just a subsetof “most relevant” states, so long as those states can be identified and they are few enoughthat matrix methods can compute the MFPT or other properties of interest (Munsky andKhammash, 2006; Kuntz et al., 2019). In such truncation-based methods, after the initialcost of enumerating states, the truncated CTMC can be reused to compute MFPTs for mildyperturbed parameters with accuracy relying on the truncated state space still containing themost relevant states. For example, we can reuse truncated CTMCs to speed up parameterinference or to optimize experimental conditions, such as the temperature, to obtain a desiredfunctionality. The key challenge here is to efficiently enumerate a suitable subset of statesthat are sufficient for accurate estimation and also few enough that the matrix methods aretractable.Here we are interested in a method that successfully addresses all three challenges forCTMCs: large state spaces, rare events, and efficient recomputation for perturbed model pa-rameters. We explore this possibility by developing a method which uses biased and localstochastic simulations to build truncated state spaces relevant to a (possibly rare) event ofinterest. We demonstrate that our method is suitable for predicting nucleic acid kinetics.
Our contributions.
We propose the pathway elaboration method for estimating MFPTs indetailed-balance CTMCs. Pathway elaboration is a time-efficient probabilistic truncation-based approach which can be used for MFPT estimation of rare events and also enables therapid evaluation of perturbed parameters. In pathway elaboration, we first construct a path-way by biasing SSA simulations from the initial states to the targets states. The biased simu-lations are guaranteed to reach the target states in expected time that is linear in the distancefrom initial to target states. Then, we expand the pathway by running SSA simulations for alimited time from every state of the pathway, with the intention of increasing accuracy by in-creasing representation throughout the pathway. Finally, we compute all possible transitionsbetween the sampled states that were not encountered in the previous two steps. For the result-ing truncated CTMC, we solve a matrix equation to compute the MFPT to the target state (orstates). Since solving matrix equations could be slow for large CTMCs, pathway elaborationincludes a δ -pruning step to efficiently prune CTMCs while keeping MFPT estimates withinpredetermined upper bounds. In this way, solving the system for other parameter settings be-comes faster. Figure 1 illustrates a conceptual figure of the pathway elaboration method andits applications.To evaluate pathway elaboration, we focus on prediction of nucleic acid kinetics. Weimplement the method using the Multistrand kinetic simulator (Schaeffer, 2013; Schaeffer,Thachuk and Winfree, 2015). Multistrand provides a secondary-structure level model of thefolding kinetics of multiple interacting nucleic acid strands, with thermodynamic energiesconsistent with NUPACK (Zadeh et al., 2011) and a stochastic simulation method based onSSA. The challenges of large state spaces, rare events, and handling perturbed parameters allarise for nucleic acid kinetics. Since the number of secondary structures may be exponen-tially large in the length of the strands, applying matrix equations is infeasible. Also, SSA Complex A Complex C (a) Hairpin closing Complex A Complex C Complex D (b) Helix dissociation Complex F Complex B Complex H Complex G (c) Toehold-mediated three-way strand displacement Fig 2: Examples of unimolecular and bimolecular interacting nucleic acid strand reactions.( a ) Hairpin closing is a unimolecular reaction. It has one reactant complex ( A ) and one prod-uct complex ( C ). The reverse reaction, hairpin opening, is also a unimolecular reaction. ( b )Helix dissociation is a unimolecular reaction. It has one reactant complex ( A ) and two prod-uct complexes ( C and D ). The reverse reaction, helix association, is a bimolecular reaction.( c ) Toehold-mediated three-way strand displacement is a bimolecular reaction. It has tworeactant complexes ( B and F ) and two product complexes ( G and H ).often takes a long time to complete for rare nucleic acid reactions. Moreover, the rapid eval-uation of mildly perturbed parameters is required, for example to calibrate the underlyingkinetic model or to obtain a desired functionality (see Figures 1f and 1g). We conduct com-putational experiments on a dataset of 267 nucleic acid kinetics (Bonnet, Krichevsky andLibchaber, 1998; Cisse, Kim and Ha, 2012; Hata, Kitajima and Suyama, 2018; Zhang et al.,2018; Machinek et al., 2014) (see Figure 2 and Table 1). The dataset consists of various typesof reactions, such as helix association and toehold-mediated three-way strand displacement,for which experimentally measured reaction rate constants vary over 8.6 orders of magnitude.We partition the 267 reactions into two sets, 237 where SSA is feasible for MPFT estimation,i.e., completes within two weeks, and the remaining 30 for which SSA is not feasible.In our experiments, first we conduct a case study and use pathway elaboration to gaininsight on the kinetics of two contrasting reactions, one being a rare event. Then, to evaluatethe estimations of pathway elaboration, we compare them with estimations obtained fromSSA for the 237 feasible reactions that were feasible with SSA. We use SSA since obtainingMFPTs with matrix equations is not possible for many of these reactions and SSA providesstatistically correct trajectories. We find that for the settings we use, the mean absolute error(MAE) of the log reaction rate constant (or equivalently the MAE of the log MFPT)is . . This is a reasonable accuracy since the log reaction constant predictions of SSAvary over . orders of magnitude (see Figure 6). In our experiments, pathway elaboration ison average 5 times faster than SSA on these reactions. We further use pathway elaborationto rapidly evaluate perturbed model parameters during optimization of Multistrand kineticparameters. We use the same 237 reactions for training the optimizer and the remaining 30as our testing set. Using the optimized parameters, pathway elaboration estimates of reactionrate constants on our dataset are greatly improved over the estimates using non-optimizedparameters. For the training set, the MAE of the log reaction rate constants of pathwayelaboration with experimental measurements reduces from . to . , that is, a . -folderror in the reaction rate constant reduces to a . -fold error on average. The MAE over the30 remaining reactions – which involve rare events and have large state spaces – reducesfrom . to . , that is, a . -fold error in the reaction rate constant reduces to a . -folderror on average. On average for these 30 reactions, pathway elaboration takes less than twodays, whereas SSA is not feasible within two weeks. The entire optimization and evaluationtakes less than five days.
2. Related Work.
There exist numerous Monte Carlo techniques (Madras, 2002; Ru-bino and Tuffin, 2009) for driving simulations towards the target states or to reduce the
HE PATHWAY ELABORATION METHOD variance of estimators. For example, importance sampling techniques (Kuwahara and Mura,2008; Shahabuddin, 1994; Hajiaghayi et al., 2014; Andrieu et al., 2003; Doucet and Johansen,2009) use an auxiliary sampler to bias the simulation, after which estimates are corrected withan importance weight. Moreover, many accelerated variants of SSA have been developed forCTMC models of chemically reacting systems (Gillespie, 2007; Cao, Gillespie and Petzold,2007; Gillespie, 2001), which can be adapted to simulate arbitrary CTMCs. There also existsa proliferation of rare event simulation methods for molecular dynamics (Weinan, Ren andVanden-Eijnden, 2002; Cabriolu et al., 2017; Zuckerman and Chong, 2017; Allen, Frenkeland ten Wolde, 2006; Allen, Valeriani and ten Wolde, 2009; Elber, 2017; Bolhuis et al.,2002). The ideas behind these methods can more or less be adapted for CTMCs and can beused along with SSA for more efficient computations. However, stochastic simulations areusually not immediately reusable for the rapid evaluation of perturbed parameters and haveto be adapted. With the pathway elaboration method we can rapidly evaluate perturbed pa-rameters whilst successfully building the truncated state spaces from initial states to targetstates within a bounded running time. Stochastic simulation methods have been to some ex-tent adapted for the rapid evaluation of perturbed parameters. SSA has been adapted in thefixed path ensemble inference (FPEI) approach (Zolaktaf et al., 2019). In this approach, anensemble of paths are generated using SSA and are then compacted and reused for mildlyperturbed parameters. To estimate MFPTs, a Monte Carlo approach is used based on expectedholding times of states. Despite being useful for parameter estimation in general, this methodis not suitable for rare events, because the paths are generated according to SSA.Coarse-graining using Markovian state models have been effective (separately) for exam-ining rare events (Sarich et al., 2014), but are mainly developed for molecular dynamics mod-els. For example, Singhal, Snow and Pande (2004) use transition path sampling (TPS) (Bol-huis et al., 2002) to build Markov state models for protein dynamics and to estimate MFPTsat perturbed temperatures using matrix equations. In brief, in TPS an ensemble of paths aregenerated using a Monte Carlo procedure. First, a single path is generated that connects theinitial and target states. New paths are then generated by picking random points along thecurrent paths and running time-limited simulations from the points. Paths that reach eitherthe initial or target state define possible new paths and the ones that do not are rejected. Eventhough we could use TPS along with SSA to simulate rare events for CTMCs (Eidelson andPeters, 2012), it is likely that many of the simulated paths could be rejected. For example, ifthe energy landscape has more than one local maximum between the initial and target states,then paths simulated from in between these local maxima could require a long simulationtime to reach either the initial or the target states.Another important problem in CTMCs is computing transient probabilities, that is theprobability distribution of the states over time. Transient probabilities can be computed ex-actly with the master equation (Van Kampen, 1992) for CTMCs that have a feasible statespace size. An important tool that has been developed to quantify the error of transient prob-ability estimations for truncated CTMCs is the finite state projection (FSP) method (Munskyand Khammash, 2006). The FSP method tells us that as the size of the state space of the trun-cated CTMC grows, the approximation monotonically improves. Also, it guarantees that theapproximate solution never exceeds the actual solution and provides bounds on the solution.As the authors of the FSP method mention, there are many ways to grow the state space, forexample by iteratively adding states that are reachable from the state space within a fixednumber of steps. However, applying matrix computations for very large state spaces couldbe time consuming. There have been many attempts to enumerate a suitable set of of states(Dinh and Sidje, 2016). In the Krylov-FSP-SSA approach (Sidje and Vo, 2015) an SSA ap-proach is used to drive the FSP and adaptive Krylov methods are used to efficiently evaluatethe matrix exponential. In brief, the method starts from an initial state space and proceeds iteratively in three steps. First, it drops low-probability states. Second, it runs SSA from eachstate of the remaining state space to incorporate probable states. Third, it adds states that arereachable within a fixed number of steps. Despite its great potential, this way of building thestate space may not be suitable for rare events. However, our pathway elaboration methoduses biased simulations to reach target states efficiently.The idea of optimizing parameter sets by using truncated CTMCs has also been used withthe Krylov-FSP-SSA method (Dinh and Sidje, 2017). Moreover, in related work (Georgoulas,Hillston and Sanguinetti, 2017), an ensemble of truncated CTMCs is used to obtain an un-biased estimator of transient probabilities, which are further used for Bayesian inference.Any success in building more efficient truncated CTMCs will also be useful in ensembleapproaches.Other types of truncation-based methods that are related to our work are probabilisticroadmap planning (Kavraki et al., 1996; Tang et al., 2005; Amato and Song, 2002; Tang,2010) methods. These methods first sample a set of states according to some criteria, suchas stability, to capture potentially important features. The states are then connected to nearbystates to form a roadmap. To generate a truncated CTMC for MFPT estimation, one couldenumerate all states that satisfy a certain criteria. For example, for nucleic acid reactions, onecould enumerate all states below a certain free energy bound. However, this approach hastwo drawbacks. First, setting the boundary too low would mean the reaction pathway is notincluded in the state space, while setting the barrier too high could make the method ineffi-cient as too many states are included. Second, this method would sample states irrespectiveof the transition rates. Instead, we rely on stochastic sampling from the initial states to thetarget states.
3. Background.
In this section, we first describe the continuous-time Markov chain(CTMC) model to which our pathway elaboration method applies and also provide relateddefinitions. Then we provide background for our experiments in Section 5 in which we ex-plain how nucleic acid kinetics can be modeled using CTMCs with the Multistrand kineticmodel (Schaeffer, 2013; Schaeffer, Thachuk and Winfree, 2015).
Continuous-time Markov chain (CTMC).
We indicate a CTMC as a tuple C =( S , K , π , S target ), where S is a countable set of states, K : S × S → R ≥ is the rate ma-trix and K ( s, s ) = 0 for s ∈ S , π : S → [0 , is the initial state distribution in which (cid:80) s ∈S π ( s ) = 1 , and S target is the set of target states. We define the set of initial statesas S init = { s ∈ S | π ( s ) (cid:54) = 0 } . For CTMCs considered here, S target ∩ S init = ∅ . A transitionbetween states s, s (cid:48) ∈ S can occur only if K ( s, s (cid:48) ) > . The probability of moving from state s to state s (cid:48) is defined by the transition probability matrix P : S × S → [0 , where(1) P ( s, s (cid:48) ) = K ( s, s (cid:48) ) E ( s, s ) . Here E : S × S → R ≥ is a diagonal matrix in which E ( s, s ) = (cid:80) s (cid:48) ∈S K ( s, s (cid:48) ) is the exitrate. The time spent in state s before a transition is triggered is exponentially distributed withexit rate E ( s, s ) . The generating matrix Q : S × S → R is Q = K − E . Detailed-balance CTMC.
In a detailed-balance CTMC C R = ( S , K , π , S target , π ) , alsoknown as a reversible CTMC, a probability distribution π : S → [0 , over the states existsthat satisfies the detailed balance condition π ( s ) K ( s, s (cid:48) ) = π ( s (cid:48) ) K ( s (cid:48) , s ) for all s, s (cid:48) ∈ S .The detailed balance condition is a sufficient condition for ensuring that π is a stationarydistribution ( π P = π ). For a detailed-balance finite-state CTMC, π is the unique stationarydistribution of the chain and is also the unique equilibrium distribution (Whitt, 2006). Boltzmann distribution.
In many Markov models of physical systems, eventually thepopulation of states will stabilize and reach a Boltzmann distribution (Schaeffer, Thachuk
HE PATHWAY ELABORATION METHOD and Winfree, 2015; Flamm et al., 2000; Tang, 2010) at equilibrium. With this distribution,the probability that a system is in a state s is(2) π ( s ) = 1 Z e − E ( s ) kBT , where E ( s ) is the energy of the system at state s , T is the temperature, k B is the Boltzmannconstant, and Z = (cid:80) s ∈S e − E ( s ) kBT is the partition function. To ensure that at equilibrium statesare Boltzmann distributed, the detailed balance conditions are(3) K ( s, s (cid:48) ) K ( s (cid:48) , s ) = e − E ( s (cid:48) ) − E ( s ) KBT . Reversible transition.
In this work, a reversible transition between states s and s (cid:48) means K ( s, s (cid:48) ) > if and only if K ( s (cid:48) , s ) > . Trajectories and paths.
A trajectory ( s , t ) , ( s , t , ) , ..., ( s m , t m , ) with m transitionsover a CTMC C = ( S , K , π , S target ) is a sequence of states s i and holding times t i for which K ( s i , s i +1 ) > and t i ∈ R > for i ≥ . We define a path s , s , ..., s m with m transitions overa CTMC C = ( S , K , π , S target ) as a sequence of states s i for which K ( s i , s i +1 ) > . The stochastic simulation algorithm (SSA).
SSA (Gillespie, 1977; Doob, 1942) simu-lates statistically correct trajectories over a CTMC C = ( S , K , π , S target ) . At state s i , theprobability of sampling s i +1 is P ( s i , s i +1 ) . At a jump from state s i , it samples the holdingtime T i from an exponential distribution with exit rate E ( s, s ) = (cid:80) s (cid:48) ∈S K ( s, s (cid:48) ) . Mean first passage time (MFPT).
In a CTMC C = ( S , K , π , S target ) , for a state s ∈ S and a target state s f ∈ S target , the MFPT τ s is the expected time to first reach s f starting fromstate s . For state s , the MFPT from s to s f equals the expected holding time in state s plusthe MFPT to s f from the next visited state (Suhov and Kelbert, 2008), so τ s = 1 E ( s, s ) + (cid:88) s (cid:48) ∈S K ( s, s (cid:48) ) E ( s, s ) τ s (cid:48) . (4)Multiplying the equation by the exit rate E ( s, s ) = (cid:80) s (cid:48) ∈S K ( s, s (cid:48) ) then yields (cid:88) s (cid:48) ∈S K ( s, s (cid:48) )( τ s (cid:48) − τ s ) = − . (5)Now writing t : S \ s f → R ≥ to be the vector of MFPTs for each state, such that t [ s ] = τ s ,we find a matrix equation as(6) ˜ Qt = − , where ˜ Q is obtained from Q by eliminating the row and column corresponding to the targetstate, and is a vector of ones. If there exists a path from every state to the final state s f ,then ˜ Q is a weakly chained diagonally dominant matrix and is non-singular (Azimzadeh andForsyth, 2016). The MFPT from the initial states to the target state s f is found as τ π = (cid:88) s ∈S π ( s ) τ s . (7)If instead of a single target state s f we have a set of target states S target , then tocompute the MFPT to S target we convert all target states into one state s f so that S ∗ = S \ S target ∪ { s f } . For s, s (cid:48) ∈ S ∗ \ { s f } , we update the rate matrix K ∗ : S ∗ → R ≥ by K ∗ ( s, s f ) = (cid:80) s (cid:48)(cid:48) ∈S target K ( s, s (cid:48)(cid:48) ) , K ∗ ( s, s (cid:48) ) = K ( s, s (cid:48) ) , and K ∗ ( s f , s ) is not used in the com-putation of the MFPT (see Eq. 6). Truncated CTMC.
Let ˆ S ⊆ S be a subset of the states over the CTMC C = ( S , K , π , S target ) or detailed-balance CTMC C R = ( S , K , π , S target , π ) and let ˆ S target ⊆ ˆ S . We construct therate matrix ˆ K : ˆ S × ˆ S → R ≥ as(8) ˆ K ( s, s (cid:48) ) = K ( s, s (cid:48) ) . We construct the initial probability distribution ˆ π : ˆ S → [0 , as(9) ˆ π ( s ) = π ( s ) (cid:80) s ∈ ˆ S π ( s ) . We define the truncated CTMC as ˆ C = ( ˆ S , ˆ K , ˆ π , ˆ S target ) and ˆ C R = ( ˆ S , ˆ K , ˆ π , ˆ S target , ˆ π ) for C and C R , respectively. For a detailed-balance ˆ C R , ˆ π : ˆ S → [0 , defined as(10) ˆ π ( s ) = π ( s ) (cid:80) s ∈ ˆ S π ( s ) , satisfies the detailed balance conditions in ˆ C R and is the unique equilibrium distribution of ˆ S in ˆ C R (Whitt, 2006).3.1. The Multistrand Kinetic Model of Interacting Nucleic Acid Strands.
Here we pro-vide background for our experiments in Section 5. We describe how the Multistrand kineticsimulator (Schaeffer, 2013; Schaeffer, Thachuk and Winfree, 2015) models the kinetics ofmultiple interacting nucleic acid strands as CTMCs and how it estimates reaction rate con-stants from MFPT estimates for these reactions.
Interacting nucleic acid strands (reactions) . Following Multistrand (Schaeffer, 2013;Schaeffer, Thachuk and Winfree, 2015), we are interested in modeling the interactions of nu-cleic acid strands in a stochastic regime. In this regime, we have a discrete number of nucleicacid strands (a set called Ψ ∗ ) in a fixed volume V (the “box”) and under fixed conditions,such as the temperature T and the concentration of Na + and Mg cations. This regimecan be found in systems that have a small volume with a fixed count of each molecule, andcan also be applied to larger volumes when the system is well mixed. Moreover, it can beused to derive reaction rate constants of reactions in a chemical reaction network that followsmass-action kinetics (Schaeffer, 2013; Schaeffer, Thachuk and Winfree, 2015).Following Multistrand (Schaeffer, 2013; Schaeffer, Thachuk and Winfree, 2015), a com-plex is a subset of strands of Ψ ∗ that are connected through base pairing (see Figure 2).A complex microstate is the complex base pairs, that is secondary structure. A system mi-crostate is a set of complex microstates, such that each strand ψ ∈ Ψ ∗ is part of exactly onecomplex. A unimolecular reaction with reaction rate constant k (units s − ) has the form A k −→ C + D, (11)and a bimolecular reaction with reaction rate constant k (units M − s − ) has the form B + F k −→ G + H. (12)Each reactant and product is a complex; A , B , C and G are nonempty but D and H may beempty complexes. For example, hairpin closing (Figure 2a) is a unimolecular reaction involv-ing one strand, where complexes A and C are comprised of this one strand, while D is empty.Helix dissociation (Figure 2b) is an example of a unimolecular reaction where complex A hastwo strands while C and D are each of one of these strands. An example of a bimolecular re-action with two reactants and two non-empty products is toehold-mediated three-way stranddisplacement (Figure 2c). We discuss these type of reactions further in Section 5.1. We areinterested in computing the reaction rate constants of such reactions. HE PATHWAY ELABORATION METHOD The Multistrand kinetic model.
Multistrand (Schaeffer, 2013; Schaeffer, Thachuk andWinfree, 2015) is a kinetic simulator for analyzing the folding kinetics of multiple in-teracting nucleic acid strands. It can handle both a system of DNA strands and a sys-tem of RNA strands . The Multistrand kinetic model is a detailed-balance CTMC C R =( S , K , π , S target , π ) for a set of interacting nucleic acid strands Ψ ∗ in a fixed volume V (the “box”) and under fixed conditions, such as the temperature T and the concentration ofNa + and Mg cations. The state space S of the CTMC is the set of all non-pseudoknottedsystem microstates of the set Ψ ∗ of interacting strands. The transition rate K ( s, s (cid:48) ) is non-zero if and only if s and s (cid:48) differ by a single base pair . Multistrand distinguishes betweenunimolecular transitions, in which the number of strands in each complex remains constant,and bimolecular transitions where this is not the case. There are bimolecular join moves,where two complexes merge, and bimolecular break moves, where a complex falls apart andreleases two separate complexes.The transition rates in the Multistrand kinetic model obey detailed balance as(13) K ( s, s (cid:48) ) K ( s (cid:48) , s ) = e − ∆ G ◦ box( s (cid:48) ) − ∆ G ◦ box( s ) RT , where ∆ G ◦ box ( s ) is the free energy of state s (units: kcal mol − ) and depends on the temper-ature T (units: K ) as ∆ G = ∆ H − T ∆ S , and R ≈ . × − kcal K − mol − is the gasconstant. The enthalpy ∆ H and entropy ∆ S are fixed and calculated in the model using ther-modynamic models that depend on the concentration of Na + and Mg cations and also on avolume-dependent entropy term. The detailed balance condition determines the ratio of ratesfor reversible transitions. A standard kinetic model that is used in Multistrand to determinethe transition rates is the Metropolis kinetic model (Metropolis et al., 1953), where all ener-getically favourable transitions occur at the same fixed rate and energetically unfavourabletransitions scale with the difference in free energy. Unimolecular transition rates are given as(14) K ( s, s (cid:48) ) = (cid:40) k uni if ∆ G ◦ box ( s ) < ∆ G ◦ box ( s (cid:48) ) ,k uni e − ∆ G ◦ box( s (cid:48) ) − ∆ G ◦ box( s ) RT otherwise,and bimolecular transition rates are given as(15) K ( s, s (cid:48) ) = (cid:40) k bi u join move, k bi e − ∆ G ◦ box( s (cid:48) ) − ∆ G ◦ box( s )+∆ G ◦ volume RT × M break move,where u is the concentration of the strands (units: M), ∆ G ◦ volume = − RT ln u , k uni > is theunimolecular rate constant (units: s − ), and k bi > is the bimolecular rate constant (units:M − s − ). The kinetic parameters θ = { k uni , k bi } are calibrated from experimental measure-ments (Wetmur and Davidson, 1968; Morrison and Stols, 1993).The distribution π is an initial distribution over the microstates of the reactant complexes,and the set S target is a subset of the microstates of the product complexes, which we deter-mine based on the type of the reaction (see Section 5.1). To set π for unimolecular reactions, Currently, Multistrand does not handle a system of mixed DNA and RNA strands, though it can be extendedto handle such systems using good thermodynamic parameters. A pseudoknotted secondary structure has at least two base pairs in which one nucleotide of a base pair isintercalated between the two nucleotides of the other base pair. A non-pseudoknotted system microstates doesnot contain any pseudoknotted secondary structures. Currently, Multistrand excludes pseudoknotted secondarystructures due to computationally difficult energy model calculations. Multistrand allows Watson-Crick base pairs to form, that is A-T and G-C in DNA and A-U and G-C in RNA.Additionally, it provides an option to allow G-T in DNA and G-U in RNA. we use particular complex microstates. One illustrative example is the (unimolecular) hair-pin closing reaction, where we set π ( h ) = 1 for the system microstate that has no base pairsand π ( s ) = 0 for all other structures, and S target is the system microstate where the Fora bimolecular reaction, when the bimolecular transitions are slow enough between the twocomplexes, it is valid to assume the complexes each reach equilibrium before bimoleculartransitions occur and therefore are Boltzmann distributed (Schaeffer, 2013). Let CM be theset of all possible complex microstates of a complex B in a volume. A distribution π b isBoltzmann distributed with respect to complex B if and only if π b ( c (cid:48) ) = e − ∆ G ( c (cid:48) ) /RT (cid:80) c ∈CM e − ∆ G ( c ) /RT (16)for all complex microstates c (cid:48) ∈ CM . In a bimolecular reaction of the form in Eq. 12, for asystem microstate s that has complex microstates c and c (cid:48) corresponding to complexes B and F , we define the initial distribution as π ( s ) = π b ( c ) × π b ( c (cid:48) ) . For all other states, we define π ( s ) = 0 .Following the conventions of Multistrand (Schaeffer, 2013), we estimate the reaction rateconstant for a reaction from its MFPT τ π (Eq. 7). For a reaction in the form of Eq. 11,(17) k = 1 τ π . In the limit of low concentrations for a reaction in the form of Eq. 12,(18) k = 1 u τ π .
4. The Pathway Elaboration Method.
We are interested in efficiently estimatingMFPT of rare events in detailed-balance CTMCs and also the rapid evaluation of mildly per-turbed parameters. Our approach is to create a reusable in-memory representation of CTMCs,which we call a truncated CTMC, and to compute the MFPTs through matrix equations(Eqs. 6 and 7).We propose the pathway elaboration method for building a truncated detailed-balanceCTMC ˆ C R for a detailed-balance CTMC C R . We call this approach the pathway elaborationmethod as we build a truncated CTMC by elaborating an ensemble of prominent paths in thesystem. The method has three main steps to build a truncated CTMC, and an additional stepfor the rapid evaluation of perturbed parameters.1. The “pathway construction” step uses biased simulations to find an ensemble of shortpaths from the initial states to the target states. This step is inspired by importance sam-pling (Madras, 2002; Rubino and Tuffin, 2009; Andrieu et al., 2003; Hajiaghayi et al.,2014) and exploration-exploitation trade-offs (Sutton and Barto, 2018).2. The “state elaboration” step uses SSA from every state in the pathway to add additionalstates to the pathway, with the intention of increasing accuracy. This step is inspired bythe string method (Weinan, Ren and Vanden-Eijnden, 2002).3. The “transition construction” step creates a matrix of transitions between every pair ofstates obtained from the first and second steps.4. The “ δ -pruning” step prunes the CTMC obtained from the previous steps to facilitate therapid evaluation of perturbed parameters.These steps result in a truncated detailed-balance CTMC ˆ C R = ( ˆ S , ˆ K , ˆ π , ˆ S target , ˆ π ) . Fig-ure 1, parts ( a ) to ( d ), illustrates the key steps of the pathway elaboration method, and Algo-rithm 1 provides high-level pseudocode. We next describe these steps in detail. HE PATHWAY ELABORATION METHOD Algorithm 1:
The pathway elaboration method.
Function
PathwayElaboration( C R , N , β , K , κ, π (cid:48) ) ( S , K , π , S target , π ) = C R S ← ConstructPathway( C R , N , β , π (cid:48) ) ˆ S ← ∅ for s ∈ S do S (cid:48) ← ElaborateState( s, C R , K , κ ) // Run SSA K times from s with atime limit of κ and return the visited states. ˆ S ← ˆ S ∪ S (cid:48) ˆ K ← Construct rate matrix from ˆ S and K // Eq. 8. return ˆ C R = ( ˆ S , ˆ K , ˆ π , ˆ S target , ˆ π ) // For ˆ π and ˆ π , see Eq. 9 and Eq. 10, respectively. Function
ConstructPathway( C , N , β , π (cid:48) ) ( S , K , π , S target ) = CS ← ∅ for n = 1 to N do Sample s ∼ π S ← S ∪ { s } Sample s b ∼ π (cid:48) for t =1,2, ... doif s = s b then breakSample z ∼ Uniform (0 , if z < β then // Bias simulations towards s b using Eq. 19. Sample s (cid:48) | s ∼ P ( ·| X t − = s ) else Sample s (cid:48) | s ∼ ˘ P s b ( ·| X t − = s ) S ← S ∪ s (cid:48) s ← s (cid:48) return S Pathway construction.
We construct a pathway by biasing N SSA simulations towardsthe target states. We bias a simulation by using the shortest-path distance function d : S × S target → R ≥ from every state s ∈ S to a fixed target state s b ∈ S target (Kuehlmann,McMillan and Brayton, 1999; Hajiaghayi et al., 2014). For every biased path, we can use adifferent s b . Therefore, in general, we can sample s b from a probability distribution π (cid:48) overthe target states. Given s b , we use an exploitation-exploration trade-off approach. At eachtransition, the process randomly chooses to either decrease the distance to s b or to explorethe region based on the actual probability matrix of the transitions.Let D s b ( s ) be the set of all neighbors of s whose distance with s b is one less than thedistance of s with s b , and let P ( s, s (cid:48) ) be as in Eq. 1. Instead of sampling states according to P , we use ˜ P : S × S → R ≥ where(19) ˜ P ( s, s (cid:48) ) = P ( s, s (cid:48) ) = K ( s,s (cid:48) ) (cid:80) s (cid:48)(cid:48)∈S K ( s,s (cid:48)(cid:48) ) ≤ z ≤ β, ˘ P s b ( s, s (cid:48) ) = K ( s,s (cid:48) ) { s (cid:48) ∈D sb ( s ) } (cid:80) s (cid:48)(cid:48)∈S K ( s,s (cid:48)(cid:48) ) { s (cid:48)(cid:48) ∈D sb ( s ) } β < z ≤ . Here z is chosen uniformly at random from [0 , , β is a threshold, and { . } is an indicatorfunction that is equal to 1 if the condition is met and 0 otherwise. When β = 1 , then ˜ P ( s, s (cid:48) ) = P ( s, s (cid:48) ) . P ROPOSITION
Let d max be the maximum distance from a state in a CTMC to targetstate s b . Then when ≤ β < / , the expected length of a pathway that is sampled accordingto Eq. 19 is at most d max − β . P ROOF . Based on the distance of states with s b , we can project a biased path that is gener-ated with Eq. 19 to a 1-dimensional random walk R , where coordinate x = 0 corresponds to s b and coordinate x > corresponds to all states s (cid:54) = s b with d ( s, s b ) = x . From the definitionof ˜ P and since all states have a path to s b by a transition to a neighbor state that decreasesthe distance by one, at each step, the random walk either takes one step closer to x = 0 withprobability at least − β or one step further from x = 0 with probability at most β . If we let E ( R, k ) denote the expected time for random walk R to reach from k , then we have thatwhen ≤ β < / ,(20) E ( R, k ) ≤ k − β , which follows from classical results on biased random walks—see Feller XIV.2 (Feller,1968). Therefore, if ≤ β < / , the proposition holds, and the state space built with N biased paths from the initial state s to a target state s b has expected size E [ | ˆ S| ] ≤ N · d ( s , s b )1 − β ≤ N · d max − β . (21)If for each biased path, the initial state is sampled from π and the target state is sampled from π (cid:48) , then we sum over the N sampled (initial state, target state) pairs, and the total expectedstate space size is still bounded by N · d max − β .For efficient computations, we should compute the shortest-path distance efficiently. Forelementary step models of interacting nucleic acid strands, we can compute d ( s, s b ) by com-puting the minimum number of base pairs that need to be deleted or formed to convert s to s b . Multistrand provides a list of base pairings for every complex microstate in a systemmicrostate (state) and we can calculate the distance between two states in a running time ofO( b ), where b is the number of bases in the strands. State elaboration.
By using Eq. 19, a biased path could have a low probability of reachinga state that has a high probability of being visited with SSA. For example, in some helixassociation reactions (Zhang et al., 2018), intra-strand base pairs are likely to form beforecompleting hybridization. However, the corresponding states do not lie on the shortest pathsfrom the initial states to the target states. Let c be the minimum number of transitions from s that are required to reach s but which increase the distance to s b . Let the random walk R bedefined as the previous step. Let P denote the probability of reaching s b before reaching s for this random walk. Following classical results on biased random walks (Feller, 1968), for β (cid:54) = 1 / , P ≥ ( β − β ) c − β − β ) dsb ( s )+ c − . In the extreme case if β = 0 , then P = 1 and the probabilityof reaching s will be 0.Therefore, for detailed-balance CTMCs, we elaborate the pathway to possibly includestates that have a high probability of being visited with SSA but were not included with ourbiased sampling. Here, we use SSA to elaborate the pathway; we run K simulations fromeach state of the pathway for a maximum simulation time of κ , meaning that a simulationstops as soon as the simulation time becomes greater than κ . By simulation time we mean thetime of a SSA trajectory, not the wall-clock time. K and κ are tuning parameters that affectthe quality of predictions. The running time of elaborating the states in the pathway with thisapproach is O( | ˆ S| Kκk max ), where ˆ S is the state space of the pathway and k max is the fastest HE PATHWAY ELABORATION METHOD Fig 3: In the elaboration step, the simulation finds s and s (cid:48) but not s (cid:48)(cid:48) . Without detailedbalance, a slow transition from s (cid:48) to s could make the MFPT to F large. However, in the fullstate space, s (cid:48) might quickly reach F via a fast transition to s (cid:48)(cid:48) .rate in the pathway. Alternatively, we could use a fixed number of transitions instead of afixed simulation time. Another approach is to add all states that are within distance r of everystate of the pathway. However, with this approach, the size of the state space could explode,whereas by using SSA the most probable states will be chosen.Note that any elaboration which stops before hitting the target state might be problematicfor non-detailed-balance CTMCs. Trajectories that stop while visiting a state for the first timemight effectively be introducing a spurious sink into the enumerated state space. Without re-versibility that last transition of the elaboration might be irreversible. Sink states that are nota target state make the MFPT to the target states infinite. For example in Figure 3, assumein the elaboration step, the simulation finds s and s (cid:48) but not s (cid:48)(cid:48) (or any other neighbor of s (cid:48) ). Then without the reversible transition, s (cid:48) will be a sink state and the MFPT to the targetstate F will be infinite. Moreover, having reversible transitions that do not obey the detailedbalance condition may make MFPT estimations large. For example, in Figure 3 assume thatthe reversible transitions between s and s (cid:48) do not obey detailed balance. Also, assume π ( s ) and π ( s (cid:48) ) are both high, and that K ( s, s (cid:48) ) is large whereas K ( s (cid:48) , s ) is small. Therefore, if theelaboration stops at s (cid:48) it will make the MFPT large. However, in the full state space, s (cid:48) mightquickly reach F through a fast transition to s (cid:48)(cid:48) . Thus, the state elaboration step may not besuitable for non-detailed-balance CTMCs. Transition construction.
After the previous two steps, fast transitions between the states ofthe pathway could still be missing. To make computations more accurate, we further com-pute all possible transitions in ˆ S that were not identified in the previous two steps. In relatedroadmap planning work (Tang et al., 2005; Thomas et al., 2013; Kavraki et al., 1996), statesare connected to their nearest neighbors as identified by a distance metric. We can include allmissing transitions by checking whether every two states in ˆ S are neighbors in O ( | ˆ S| ) or bychecking for every state in ˆ S whether its neighbors are also in ˆ S in O( | ˆ S| m } ), where m isthe maximum number of neighbors of the states in the original CTMC. δ -pruning. Given a (truncated) CTMC in which we can compute the MFPT from every stateto the target state, one question is: which states and transitions can be removed from theMarkov chain without changing the MFPT from the initial states significantly? This questionis especially relevant for the rapid evaluation of perturbed parameters, where MFPTs need tobe recomputed often.Given a CTMC C = ( S , K , π , S target ) and a pruning bound δ , let the MFPT from anystate s to S target be τ s and let the MFPT from the initial states to S target be τ π . Let S δp = { s ∈ S | τ s < δτ π and π ( s ) = 0 } be the set of states that are δ -close to S target andthat are not an initial state. We construct the δ -pruned CTMC C δ = ( S δ , π , K δ , { s d } ) over the pruned set of states S δ = S \ S δp ∪ { s d } , where s d is the new target state. For s, s (cid:48) ∈ S δ \ { s d } , we update the rate matrix K δ : S δ → R ≥ by K δ ( s, s d ) = (cid:80) s (cid:48) ∈S δp K ( s, s (cid:48) ) and K δ ( s, s (cid:48) ) = K ( s, s (cid:48) ) . Note that K δ ( s d , s ) is not used in the computation of the MFPT(Eq. 6), so we can simply assume K δ ( s d , s ) = 0 . Alternatively, to retain detailed-balance con-ditions, we can define the energy of s d as E ( s d ) = − RT log (cid:80) s (cid:48)(cid:48) ∈S δp e − E ( s (cid:48)(cid:48) ) RT (see Eqs.7.1and 7.2 from Schaeffer (Schaeffer, 2013)) and define K δ ( s d , s ) = e − E ( s ) − E ( sd ) RT K δ ( s, s d ) . Forthe pruned CTMC C δ = ( S δ , π , K δ , { s d } ) , let the MFPT τ δπ be given as usual (Eq. 7). Thenby construction(22) τ δπ ≤ τ π δ . We can calculate the MFPT from every state to the target states by solving Eq. 6 oncefor CTMC C . Therefore, the running time of δ -pruning depends on the running time of thematrix equation solver that is used. For a CTMC with state space S , the running time of adirect solver is at most O( |S| ). For iterative solvers the running time is generally less thanO( |S| ). After the equation is solved, the CTMC can be pruned in O( |S| ) for any δ . Note thatfor a given bound δ , the running time for solving Eq. 6 for the pruned CTMC C δ might stillbe high. In that case, a larger value of δ is required. To set δ in practice, it could be useful toconsider the number of states that will be pruned for a given δ , that is |S δp | . Updating perturbed parameters.
We are interested in rapidly estimating the MFPT to targetstates given mildly perturbed parameters. Our approach is to reuse a truncated CTMC formild perturbations. The MFPT estimates will be biased in this way. However, we could havesignificant savings in running time by avoiding the cost of sampling and building truncatedCTMCs from scratch. We would still have to solve Eq. 6, but it could be negligible comparedto the other costs. For example in Table 2, on average, solving the matrix equation is fasterthan SSA by a factor of 47 and is faster than building the truncated CTMC by a factor of 10.A perturbed thermodynamic model parameter affects the energy of the states. Therefore,to update the transition rates, we would also have to recompute the energy of the states. Aperturbed kinetic model only affects the transition rates. A perturbed experimental conditioncould affect both the energy of the states and the transition rates. Therefore, assuming theenergy of a state can be updated in a constant time, the truncated CTMC can be updatedin O( | ˆ S| + | ˆ E| ), where ˆ E is the set of transitions of the truncated CTMC. For nucleic acidkinetics the energy of a state can be computed from scratch in O( b ) time, or in O( ) timeusing the energy calculations of a neighbor state which differs in one base pair (Schaeffer,2013). Quantifying the error.
After we build truncated CTMCs, we need to quantify the error ofMFPT estimates when experimental measurements are not available. It would help us setvalues for N , β , K and κ for fixed model parameters, and also evaluate when a truncatedCTMC has a high error for perturbed model parameters. For exponential decay processes, onepossible approach is to adapt the finite state projection FSP (Munsky and Khammash, 2006)method that is developed to quantify the error of truncated CTMCs for transient probabilities.We adapt it as follows. We combine all target states into one single absorbing state s f . Weproject all states that are not in the truncated CTMC into an absorbing state s o and we redirectall transitions from the truncated CTMC to states out of the CTMC into s o . Then we use thestandard matrix exponential equations to compute the full distribution on the state space ata given time. However, we only care about the probabilities that s f and s o are occupied. Wesearch to compute the half-completion time t / with bounds by(23) (cid:40) t min s.t. p ( s f ; t min ) + p ( s o ; t min ) = ,t max s.t. p ( s f ; t max ) = , HE PATHWAY ELABORATION METHOD where p ( s ; t ) is the probability that the process will be at state s at time t starting fromthe set of initial states. Since s f and s o are the only absorbing states, then t min exists andclearly t min ≤ t / . Based on FSP, p ( s f ; t max ) is an underestimate of the actual probabilityat time t max , if it exists. A possible way to determine if a solutions exists is to determinethe probability of reaching state s f compared to state s o from the initial states, which can becalculated by solving a system of linear equations (see Eq. 2.13 from Metzner, Schütte andVanden-Eijnden (2009)). If the probability is greater or equal to then a solutions exists. If asolution does not exist for the given statespace, then based on FSP the error is guaranteed todecrease by adding more states and we can eventually find a solution to Eq. 23. The searchfor t max can be completed with binary search. Thus, the true t / is guaranteed to satisfy t min ≤ t / ≤ t max . For exponential decay processes, the relation between the half-completiontime and the MFPT is (Cohen-Tannoudji et al., 1977; Simmons, 1972)(24) t / = ln λ and τ = 1 λ → τ = t / ln , where λ is the rate of the process. Thus, t min ln ≤ τ ≤ t max ln .A drawback of this approach is that we might need a large number of states to find a solu-tion to Eq. 23, which might make the master equation or the linear system solver infeasiblein practice. Efficiently quantifying the error of MFPT estimates in truncated CTMCs for ex-ponential and non-exponential decay processes is beyond the scope of this paper. It might bepossible to use some other existing work (Kuntz et al., 2019; Backenköhler, Bortolussi andWolf, 2019).
5. Dataset and Experiments for Interacting Nucleic Acid Strands.
We implementpathway elaboration for interacting nucleic acid strands on top of Multistrand (Schaeffer,2013; Schaeffer, Thachuk and Winfree, 2015) (see Section 3.1). Our framework and datasetare available at https://github.com/DNA-and-Natural-Algorithms-Group/PathwayElaboration .In this section, first, we describe our dataset of nucleic acid kinetics. Then, we describe ourexperimental setup that is common in our experiments. Afterwards, we use pathway elabo-ration in a case study to gain insight on the kinetics of two contrasting reactions. Then, wecompare the performance of pathway elaboration with SSA. After that, we show the effec-tiveness of the δ -pruning step. Finally, we use pathway elaboration for the rapid evaluationof perturbed parameters in parameter estimation.5.1. Dataset of Interacting Nucleic Acid Strands.
We conduct computational experi-ments on interacting nucleic acid strands (see Section 3.1). The speed at which nucleic acidstrands interact is difficult to predict and depends on reaction topology, strands’ sequences,and experimental conditions. The number of secondary structures interacting nucleic strandsmay form is exponentially large in the length of the strands. Typical to these reactions arehigh energy barriers that prevent the reaction from completing, meaning that long periodsof simulation time are required before successful reactions occur. Consider reactions thatoccur with rates lower than M − s − such as three-way strand displacement at roomtemperature (see Table 1). These types of reactions are slow to simulate not because the sim-ulator takes longer to generate trajectories for larger molecules, but the slowness is insteada result of the energy landscape: at low temperatures, duplexes simply are more stable, andrequire longer simulated time until their dissociation is observed. Predicting the kinetics ofinteracting nucleic acid strands is also difficult with classical machine learning methods andneural network models. For example, Zhang et al. (2018) successfully predict hybridizationrates with a weighted neighbour voting prediction algorithm and Angenent-Mari et al. (2020) T ABLE Summary of the dataset of 267 nucleic acid kinetics. The initial concentration of the reactants is denoted as u and k is the experimental reaction rate constant. DatasetNo. Reaction type & source † [ Na + ] (M) T ( ◦ C) u (M) log k D t r a i n Hairpinopening (Bonnet,Krichevsky andLibchaber, 1998)
63 25 0 . – . –
49 1 × − . – . Hairpin closing (Bonnet,Krichevsky andLibchaber, 1998)
62 25 0 . – . –
49 1 × − . – . Helixdissociation (Cisse, Kimand Ha, 2012)
39 18 0 . – . –
37 1 × − − . – . Helix association (Hata,Kitajima and Suyama,2018)
43 46 0 .
195 25 5 × − . – . Helix association (Zhanget al., 2018)
20 72 0 .
75 37 –
55 1 × − . – . Toehold-mediatedthree-way stranddisplacement (Machineket al., 2014)
10 102 0 . ††
23 5 × − – × − . – . D t e s t Helix association (Hata,Kitajima and Suyama,2018) .
195 25 5 × − . – . Toehold-mediatedthree-way stranddisplacement (Machineket al., 2014)
26 100 0 . †† × − – × − . – . † See Figure 2 for example figures of these reactions. †† The experiment was performed without Na + in the buffer. successfully predict toehold switch function with neural networks. However, despite the ac-curate and fast computational prediction of these methods, to treat different type of reactionsor to treat different initial and target states the models have to be adapted. On the other hand,the CTMC model of Multistrand can readily be applied to unimolecular and bimolecularreactions. Moreover, the CTMC model could provide an unlimited number of unexpectedintermediate states. In contrast, with the neural networks models of Angenent-Mari et al.(2020), which use attention maps to interpret intermediate states, the number is limited.We curate a dataset of 267 interacting DNA strands from the published literature, sum-marized in Table 1. The reactions are annotated with the temperature, the buffer condition,and the experimentally determined reaction rate constant. The dataset covers a wide range ofslow and fast unimolecular and bimolecular reactions where the reaction rate constants varyover 8.6 orders of magnitude. For unimolecular reactions, we consider hairpin opening (Bon-net, Krichevsky and Libchaber, 1998), hairpin closing (Bonnet, Krichevsky and Libchaber,1998), and helix dissociation (Cisse, Kim and Ha, 2012). For bimolecular reactions, we con-sider helix association (Hata, Kitajima and Suyama, 2018; Zhang et al., 2018) and toehold-mediated three-way strand displacement (Machinek et al., 2014). The reactions from Cisse,Kim and Ha (2012) and Machinek et al. (2014) may have mismatches between the bases ofthe strands. The type of reactions in Table 1 are widely used in nanotechnology, such as inmolecular beacon probes (Chen et al., 2015). HE PATHWAY ELABORATION METHOD For bimolecular reactions, we Boltzmann sample initial reacting complexes. For reactionsin which we define only one target state, in the pathway construction step, we bias the pathstowards that state. In this work, for reactions in which we define a set of target states, we biaspaths towards only one target state, so that π (cid:48) ( s b ) = 1 for one state and π (cid:48) ( s ) = 0 for all otherstates. Next we describe these states. Hairpin closing and hairpin opening.
For a hairpin opening reaction, we define the initialstate to be the system microstate in which a strand has fully formed a duplex and a loop (seeFigure 2a). We define the target state to be the system microstate in which the strand has nobase pairs. Hairpin closing is the reverse reaction, where a strand with no base pair forms afully formed duplex and a loop.
Helix dissociation and helix association.
For a helix dissociation reaction, we specifythe initial state to be the system microstate in which two strands have fully formed a helix(see Figure 2b). We define the set of target states to be the set of system microstates inwhich the strands have detached and there are no base pairs within one of the strands. Webias paths towards the target state in which there are no base pairs formed within any of thestrands. Helix association is the reverse reaction. We Boltzmann sample the initial reactingcomplexes in which the strands have not formed base pairs with each other. We define thetarget state to be the system microstate in which the duplex has fully formed.
Toehold-mediated three-way strand displacement.
In this reaction, an invader stranddisplaces an incumbent strand in a duplex, where a toehold domain facilitates the reaction(see Figure 2c and Figure 4). We Boltzmann sample initial reacting complexes in which theincumbent and substrate form a complex through base pairing and the invader forms anothercomplex. We define the set of target states to be the set of microstates where the incumbent isdetached from the substrate and there are no base pairs within the incumbent. We bias pathstowards the target state in which the substrate and invader have fully formed base pairs andthere are no base pairs within the incumbent.In datasets No. 1-6 from Table 1, we consider reactions that are feasible with SSA withour parameterization of Multistrand, given two weeks computation time, since we compareSSA results with pathway elaboration results. We indicate these reactions as D train since wealso use them as training set in Section 5.6. We indicate datasets No. 7-8 as D test since weuse them as testing set in Section 5.6.5.2. Experimental Setup.
Experiments are performed on a system with 64 2.13GHz IntelXeon processors and 128GB RAM in total, running openSUSE Leap 15.1. An experiment fora reaction is conducted on one processor. Our framework is implemented in Python, on top ofthe Multistrand kinetic simulator (Schaeffer, 2013; Schaeffer, Thachuk and Winfree, 2015).To solve the matrix equations in Eq. 6, we use the sparse direct solver from SciPy (Virta-nen et al., 2020) when possible . Otherwise we use the sparse iterative biconjugate gradientalgorithm (Fletcher, 1976) from SciPy.In all of our experiments, the thermodynamic parameters for predicting the energy ofthe states are fixed and the energies are calculated with Multistrand. Each reaction uses itsown experimental condition as provided in the dataset. In all our experiments, we use theMetropolis kinetic model from Multistrand. For all experiments except for Section 5.6, wefix the kinetic parameters to the Metropolis Mode parameter set (Zolaktaf et al., 2017), thatis θ = { k uni ≈ . × s − , k bi ≈ . × M − s − } . To obtain MFPTs with SSA, weuse 1000 samples, except for three-way strand displacement reactions in which we use 100samples, since the simulations take a longer time to complete. The implementation we used allowed the sparse direct solver to use only up to 2GB of RAM IncumbentInvaderSubstrate InvaderSubstrateIncumbent (a)
Incumbent SubstrateSubstrate InvaderInvader Incumbent (b)
Invader base pair I n c u m b e n t b a s e p a i r |S x,y | (c) Invader base pair I n c u m b e n t b a s e p a i r − − G x,y (d) Invader base pair I n c u m b e n t b a s e p a i r . . . . . . δ x,y (e) ∆ G bo x ( kc a l / m o l ) (f) Invader base pair I n c u m b e n t b a s e p a i r |S x,y | (g) Invader base pair I n c u m b e n t b a s e p a i r − − G x,y (h) Invader base pair I n c u m b e n t b a s e p a i r . . . . . . δ x,y (i) ∆ G bo x ( kc a l / m o l ) (j) Fig 4: Results of truncated CTMCs built with pathway elaboration ( N = 128 , β = 0 . , K = 1024 , κ = 16 ns) for two toehold-mediated three-way strand displacement reactionsfrom Machinek et al. (2014). ( a ) A toehold-mediated three-way strand displacement reac-tion that has a 6-nt toehold and a 17-nt displacement domain (Machinek et al., 2014). ( b ) Atoehold-mediated three-way strand displacement reaction that has a 6-nt toehold, a 17-nt dis-placement domain, and a mismatch exists between the invader and the substrate at position 6of the displacement domain (Machinek et al., 2014). Figures 4c, 4d, 4e, and 4f correspond toFigure 4a. Figures 4g,4h,4i, and 4j correspond to Figure 4b. In Figures 4c, 4d, 4e, 4g, 4h,and 4i, the x-axis corresponds to the number of base pairs between the invader and thesubstrate, and the y-axis corresponds to the number of base pairs between the incumbentand the substrate. ( c, g ) At coordinate ( x, y ) , |S x,y | is shown, where S x,y is a systemmacrostate (a nonempty set of system microstates) equal to the set of states with coordi-nate ( x, y ) . ( d, h ) At coordinate ( x, y ) , the free energy ∆ G x,y is shown, which is definedas ∆ G x,y = − RT ln (cid:80) s ∈S x,y e − ∆ G ( s ) RT (Schaeffer, 2013). The free energy of the paths in Fig-ures 4f and 4j are also shown with the ◦ marker in Figures 4d and 4h, respectively. ( e, i ) Atcoordinate ( x, y ) , the value of δ x,y = (cid:80) s ∈S x,y w s δ ( s ) (cid:80) s ∈S x,y w s is shown, where δ ( s ) = E [ τ s ] /τ π and w s = e − ∆ G ( s ) RT . For ease of understanding, the green “halfway line” separates coordinateswhere δ x,y is greater than 0.5 from coordinates where δ x,y is less than 0.5. ( f, j ) The free en-ergy landscape of a random path built with pathway elaboration ( N = 1 , β = 0 , K = 0 , κ = 0 ns) and the initial and the final states and some states near the local extrema are illustrated. HE PATHWAY ELABORATION METHOD Case Study.
Here we illustrate the use of pathway elaboration to gain insight on thekinetics of two contrasting reactions from Machinek et al. (2014), one being a rare event.Figures 4a and 4b show the two teohold-mediated three-way strand displacement reactionsthat we consider (Machinek et al., 2014). In the reaction in Figure 4a, the invader and sub-strate are complementary strands in the displacement domain. In the reaction in Figure 4b,there is a mismatch between the invader and the substrate in the displacement domain. Therate of toehold-mediated strand displacement is usually determined by the time to completethe first bimolecular transition, in which the invader forms a base pair with the substrate forthe first time. However, the rate could be controlled by several orders of magnitude by alter-ing positions across the strand, such as using mismatch bases (Machinek et al., 2014). Thereaction in Figure 4b is approximately 3 orders of magnitude slower than the reaction in Fig-ure 4a. For the reaction in Figure 4a, log k = 6 . , log ˆ k PE = 6 . , log ˆ k SSA = 6 . , | ˆ S| = 4 . × , the computation time of pathway elaboration is . × s, and thecomputation time of SSA is . × s. For the reaction in Figure 4b, log k = 3 . , log ˆ k PE = 3 . , | ˆ S| = 7 × , the computation time of pathway elaboration is . × s,and SSA is not feasible within × s.In Figures 4c-4e and 4g-4i, we illustrate different properties of the truncated CTMCs forthe reactions in Figures 4a and 4b, respectively. Comparing Figure 4c with Figure 4g, we seethat many states are sampled midway in Figure 4g due to the mismatch. In Figures 4d and 4h,we compare the energy barrier (increase in free energy) while moving from the beginning ofthe x-axis towards the end of the x-axis. In Figure 4d, we can see a noticeable energy barrierin the beginning. However, in Figure 4h, we can see two noticeable energy barriers, one in thebeginning and one midway. Figures 4e and 4i show states that are δ -close to the target states.These figures show that with δ -pruning, states that are further from the initial states and closerto the target states will be pruned with smaller values of δ , compared to states that are closerto the initial states and further from the target states. Comparing Figure 4e with Figure 4i, thestates quickly reach the target states after the first several transitions in Figure 4e (after theenergy barrier). However, in Figure 4i, the states do not quickly reach the target states untilafter the second energy barrier. Figure 4f and 4j show the free energy landscape and someof the secondary structures for a random path from an initial state to a target state for thereactions in Figures 4a and 4b, respectively. For the reaction in Figure 4a, the barrier is nearthe first transition. For the reaction in Figure 4b, there is a noticeable barrier after severalbase pairs form between the invader and the substrate, presumably near the mismatch.5.4. Mean First Passage Time and Reaction Rate Constant Estimation.
To evaluate theestimations of pathway elaboration, we compare its estimations with estimations obtainedfrom SSA for the reactions in D train . Note that for many of these reactions the size of the statespace is exponentially large in the length of the strands. Therefore, exact matrix equations isnot possible for them. Instead we use SSA since it can generate statistically correct trajecto-ries. We also compare the wall-clock computation time of pathway elaboration with SSA forthese reactions.We evaluate the estimations of pathway elaboration based on the mean absolute error(MAE) with SSA, which is defined over a dataset D as(25) MAE = 1 |D| (cid:88) r ∈D | log ˆ τ r SSA − log ˆ τ r PE | = 1 |D| (cid:88) r ∈D | log ˆ k r SSA − log ˆ k r PE | , where ˆ τ r PE and ˆ τ r SSA are the estimated MFPTs of SSA and pathway elaboration for reaction r ,respectively, and ˆ k r SSA and ˆ k r PE are the estimated reaction rate constants of SSA and pathwayelaboration for reaction r , respectively. The equality follows from Eqs. 17 and 18. We use log differences since the reactions rate constants cover many orders of magnitude. We use | ˆ S| . . . . . . . M A E ( )( )( ) ( ) N = 32 N = 128 N = 256 N = 512 N = 4096 β = 0 β = 0 . β = 0 . β = 0 . K = 0 K = 16 K = 256 K = 1024 κ = 0 ns κ = 1 ns κ = 16 ns κ = 64 ns (a) | ˆ S| . . . . . . . . M A E ( ) ( )( )( ) (b) | ˆ S| . . . . . . . . M A E ( )( )( ) ( ) (c) | ˆ S| . . . . . . . . M A E ( )( )( )( ) (d) Fig 5: MAE vs | ˆ S| for different values of N , β , K and κ . ( a ) datasets No. 1,2, and 3, ( b )dataset No. 4, ( c ) dataset No. 5, and ( d ) dataset No. 6. The annotated values on the figurescorrespond to N , β , K , and κ , respectively.the MAE as our evaluation metric since it is conceptually easy to understand. For example,here, an MAE of means on average the predictions are off by a factor of 10. In the rest ofthis subsection, we first look at the trade-off between the MAEs and the size of the truncatedstate space set ˆ S , with regards to different parameter settings of the pathway elaborationmethod. Then we look at the trade-off between the MAE and the computation time. MAE versus | ˆ S| . Figure 5 shows the MAE versus | ˆ S| of pathway elaboration for differentconfigurations of the N , β , K , and κ parameters. Figure A1 and A2 from the Appendixrepresent Figure 5 by varying only two parameters at a time. The figures show that generallyas N and β increase, the MAE decreases. This is because for a fixed N as β → the ensembleof paths will be generated by SSA. As N → ∞ , the truncated state space becomes larger andis more likely to contain the most probable paths from the initial states to the target states.Comparing the MAE of configurations where K = 0 and κ = 0 with other settings where K > and κ > , shows that the elaboration step helps reduce the MAE (in the Appendix,compare Figures A1a-A1d with Figures A1i-A1l). Particularly, the elaboration step is usefulfor dataset No. 4, helix association from Zhang et al. (2018) where intra-strand base pairscan form before completing hybridization. The plots show that the elaboration step is moreuseful when β is small (in the Appendix, compare Figures A2a-A2d with Figures A2i-A2l). HE PATHWAY ELABORATION METHOD Table 2: Pathway elaboration ( N = 128 , β = 0 . , K = 256 , κ = 16 ns) versus SSA. The mean statistics are averaged over the ’ | ˆ S| is the size of the truncated state space. DatasetNo. | ˆ S| forpathwayelaboration Mean matrixcomputationtime (s) forpathwayelaboration Meancomputationtime (s) forpathwayelaboration Meancomputationtime (s) for SSA .
04 5 . × . × − . × . × .
03 1 . × . × − . × . × .
04 5 . × . × − . × . × .
29 8 . × . × . × . × .
51 3 . × . × . × . × .
31 3 . × . × . × . × Alldatasets
237 0 .
13 6 . × . × . × . × (a) (b) (c) (d) Fig 6: The log ˆ k SSA and log ˆ k PE ( N = 128 , β = 0 . , K = 256 , κ = 16 ns) for ( a ) datasetsNo. 1,2, and 3, and ( b ) dataset No. 4, ( c ) dataset No. 5, and ( d ) dataset No. 6. The reactions areordered along the x-axis by their predicted log ˆ k SSA . The pathway elaboration experimentsare repeated three times. For each reaction, log ˆ k PE is calculated by the average of thethree experiments. The error bars for pathway elaboration indicate the range (minimum tomaximum) of the three experiments. The error bars for SSA indicate the 95% percentilebootstrap of the log ˆ k SSA .This could be because elaboration helps find rate determining states that were not exploreddue to the biased sampling. When β → the pathway elaboration method will perform asSSA and rate determining states can be found without elaboration.Furthermore, the figures show that as K increases, the MAE decreases. However, with alarge value for κ and a small value of K the performance could be diminished (such as inFigure A2c of the Appendix). In particular, consider that K and κ might involve simulationsthat go on excursions outside the “main” densely-visited parts of the enumerated state space,and they might even terminate out there. Such excursions might very well introduce signif-icant local minima into the enumerated state space - even when no significant local minimaexist in the original full state space. For example, consider an excursion that goes off-pathdown a wide slope, perhaps toward the target state. If it terminates before reaching a targetstate, then a hypothetical simulation in the enumerated state space could get stuck, needingto climb back up the slope to the point where the excursion began. The expected hitting timein the enumerated state space will account for such wasted time, thus leading to an overestimation of the MFPT. Therefore, κ should be tuned with respect to K . MAE versus computation time.
Table 2 illustrates the MAE and computation time ofpathway elaboration for when N = 128 , β = 0 . , K = 256 , and κ = 16 ns compared with (a) (b) (c) Fig 7: The effect of δ -pruning with different values of δ on truncated CTMCs that are builtwith pathway elaboration ( N = 128 , β = 0 . , K = 1024 , κ = 16 ns) for dataset No. 6. δ = 0 indicates δ -pruning is not used. ( a ) The log ˆ k . ( b) The size of the truncated state space | ˆ S| .( c ) The computation time for solving Eq. 6.SSA. We illustrate this parameter setting because it provides a good trade-off between ac-curacy and computational time for the larger reactions. For the smaller reactions, we couldachieve the same MAE with less computational time (by using smaller values for the param-eter setting). Figure 6 further shows the prediction of pathway elaboration for this parametersetting compared to the prediction of SSA for individual reactions. In Table 2, the MAE forunimolecular reactions is smaller than . , whereas for bimolecular reactions it is largerthan . . This is because the CTMCs for the bimolecular reactions in our dataset are natu-rally bigger than the CTMCs for the unimolecular reactions in our dataset, and require largertruncated CTMCs. The MAE can be further reduced by changing the parameters (as shownin Figure 5). With our implementation of pathway elaboration, the computation time of path-way elaboration for datasets No. 3, No. 4, and No. 6 are 2 times, 20 times, 3 times smallerthan SSA, respectively. The computation time of SSA for datasets No. 1, No. 2, and No. 5 issmaller than the computation time of pathway elaboration. This is because pathway elabora-tion has some overhead, and in cases where SSA is already fast it can be slow. However, as weshow in Section 5.6, even for these reactions, pathway elaboration could still be useful for therapid evaluation of perturbed parameters. Also, the computation time for pathway elaborationcould be significantly improved with more efficient implementations of the method.5.5. δ -Pruning. Figure 7 shows how δ -pruning affects the quality of the log reactionrate constant estimates, the size of the state spaces, and the computation time of solving thematrix equations, for dataset No. 6. The MFPT estimates satisfy the bound given by Eq. 22whilst δ -pruning reduces the computation time for solving the matrix equations by an orderof magnitude for δ = 0 . . Using larger values of δ we can further decrease the computationtime. If we reuse the CTMCs many times, such as in parameter estimation, δ -pruning couldhelp reduce computation time significantly.5.6. Parameter Estimation.
In the previous subsections the underlying parameters of theCTMCs were fixed. Here we assume the parameters of the kinetic model of the CTMCs arenot calibrated and we use pathway elaboration to build truncated CTMCs to rapidly evaluateperturbed parameter sets during parameter estimation. We use the 237 reactions indicated as D train in Table 1 as our training set. We use the 30 rare event reactions indicated as D test inTable 1 to show that given a well-calibrated parameter set for the CTMC model, the pathway HE PATHWAY ELABORATION METHOD k uni10 k b i θ θ ∗ (a) (b) Fig 8: Results of parameter estimation using pathway elaboration ( N = 128 , β = 0 . , K = 256 , κ = 16 ns). ( a ) The parameters are optimized from an initial simplex of θ andits perturbations to θ ∗ = { k uni ≈ . × s − , k bi ≈ . × M − s − } . ( b ) The pa-rameters are optimized using D train , shown with a line graph, and evaluated on dataset No.7 and No. 8. The markers are annotated with the MSE and MAE of the datasets when thetruncated CTMCs are built from scratch using θ and θ ∗ .elaboration method can estimate MFPTs and reaction rate constants of reactions close to theirexperimental measurement.We seek the parameter set that minimizes the mean squared error (MSE) as θ ∗ = argmin θ |D train | (cid:88) r ∈D train (log τ r − log ˆ τ r PE ( θ )) = argmin θ |D train | (cid:88) r ∈D train (log k r − log ˆ k r PE ( θ )) , (26)which is a common cost function for regression problems. The equality follows from Eqs. 17and 18. We use the Nelder-Mead optimization algorithm (Nelder and Mead, 1965; Virta-nen et al., 2020) to minimize the MSE. We initialize the simplex in the algorithm with θ = { k uni = 5 × s − , k bi = 5 × M − s − } in which we choose arbitrarily and twoperturbed parameter sets. Each perturbed parameter set is obtained from θ by multiplyingone of the parameters by . , which is the default implementation of the optimization soft-ware (Virtanen et al., 2020). For every reaction, we also initialize the Multistrand kineticmodel with θ . We build truncated CTMCs with pathway elaboration ( N = 128 , β = 0 . , K = 256 , κ = 16 ns). Whenever the matrix equation solving time is large (here we considera time of s large), we use δ -pruning (here we use δ values of . − . ) to reduce thetime. During the optimization, for a new parameter set we update the parameters in the kineticmodel of the truncated CTMCs and we reuse the truncated CTMC to evaluate the parameterset. Similar to our previous work (Zolaktaf et al., 2019), to reduce the bias and to ensure thatthe truncated CTMCs are fair with respect to the optimized parameters, we can occasionallyrebuild truncated CTMCs from scratch.Although we use the MSE of pathway elaboration with experimental measurements asour cost function in the optimization procedure, the MAE of pathway elaboration with ex-perimental measurements also decreases. Figure 8 shows how the parameters, the MSE, andthe MAE change during optimization. The markers are annotated with the MSE and theMAE of D train and datasets No. 7-8 when truncated CTMCs are built from scratch. TheMAE of D train with the initial parameter set θ is . . The optimization finds θ ∗ = { k uni ≈ . × s − , k bi ≈ . × M − s − } and reduces the MAE of D train to . . The MAE of dataset No. 7 and dataset No. 8, which are not used in the optimization, reduce from . to . and from . to . , respectively.Overall, the experiment in this subsection shows that pathway elaboration enables MFPTestimation of rare events. It predicts their MFPTs close to their experimental measurementsgiven an accurately calibrated model for their CTMCs. Moreover, it shows that pathway elab-oration enables the rapid evaluation of perturbed parameters and makes feasible tasks suchas parameter estimation which benefit from such methods. On average for the 30 reactionsin the testing set, pathway elaboration takes less than two days, whereas SSA is not feasiblewithin two weeks. The entire experiment in Figure 8 takes less than five days parallelized on40 processors. Note that clearly our optimization procedure could be improved, for exam-ple by using a larger dataset or a more flexible kinetic model. However, this experiment is apreliminary study; we leave a rigorous study on calibrating nucleic acid kinetic models withpathway elaboration and possible improvements to future studies.
6. Discussion.
In this paper, we address the problem of estimating MFPTs of rare eventsin large CTMCs and also the rapid evaluation of perturbed parameters. We propose the path-way elaboration method, which is a time-efficient probabilistic truncation-based approachfor MFPT estimation in CTMCs. We conduct computational experiments on a wide range ofexperimental measurements to show pathway elaboration is suitable for predicting nucleicacid kinetics. In summary, our results are promising, but there is still room for improvement.Using pathway elaboration, in the best possible case, the sampled region of states andtransitions is obtained faster than SSA, but without significant bias in the collected states andtransitions. The sampled region may however qualitatively differ from what would be ob-tained from SSA, which may compromise the MFPT estimates. Moreover, reusing truncatedCTMCs for significantly perturbed parameters could lead to inaccurate estimation of theMFPT in the original CTMC. In Section 4, for exponential decay processes, we introduceda method that could help us quantify the error of the MFPT estimate. However, it might beslow in practice. So how can we efficiently tune these parameters? Similar to SSA, for afixed β and when K = 0 and κ = 0 , we could increase N until the estimated MFPT stopschanging significantly (based on the law of large numbers it will converge). Note that for K = 0 and κ = 0 we could compute the MFPT by computing the average of the biased pathswithout solving matrix equations. As shown in Proposition 4.1, if we set β to less than / ,then biased paths will reach target states in expected time that is linear in the distance frominitial to target states. For setting K , one possibility is to consider the number of neighborsof each state. A reaction where states have a lot of neighbors requires a larger K comparedto a reaction where states have a smaller number of neighbors. κ should be set with respectto K . As stated in Section 5, a large value of κ along with a small value of K could resultin excursions that do not reach any target state and lead to overestimates of the MFPT. Onecould set κ to a small value and then increase K until the MFPT estimate stops changing,and could repeat this process while feasible.In the pathway elaboration method, we estimate MFPTs by solving matrix equations.Thus, its performance depends on the accuracy and speed of matrix equation solvers. Forexample, applying matrix equation solvers may not be suitable if the initial states lie very farfrom the target states, since the size of the truncated CTMCs depends on the shortest-pathdistance between these states. Although solving matrix equations through direct and iterativemethods has progressed, both theoretically and practically (Fletcher, 1976; Virtanen et al.,2020; Cohen et al., 2018), solving stiff (multiple time scales) or very large equations couldstill be problematic in practice. More stable and faster solvers would allow us to estimateMFPTs for stiffer and larger truncated CTMCs. Moreover, it might be possible to use fastupdates for solving the matrix equations (Brand, 2006; Parks et al., 2006). Therefore, if we HE PATHWAY ELABORATION METHOD require to compute MFPT estimates with matrix equations as we monotonically grow thesize of the state space or for a perturbed parameter set, the total cost for solving all the linearsystems would be the same cost as solving the final linear system from scratch.We might be able to improve the pathway elaboration method to relieve the limitationsdiscussed above. For example, it might be possible to use an ensemble of truncated CTMCsto obtain an unbiased estimate of the MFPT (Georgoulas, Hillston and Sanguinetti, 2017).To avoid excursions that lead to overestimation of the MFPT in the state elaboration step, wecould run the pathway construction step from the last states visited in the state elaborationstep. This would also relax the constraint of having reversible or detailed balance transitions.Presumably, an alternating approach of the two steps would make the approach more flexible.Moreover, currently we run the state elaboration step from every state of the pathway withthe same setting, which might not be necessary. Efficiently running the state elaboration stepas necessary, could reduce the time to construct the truncated CTMC in addition to the matrixequation solving time.Finally, we evaluated the pathway elaboration method for predicting the MFPT of nucleicacid kinetics. However, the method is generally applicable to detailed-balance CTMC mod-els, such as protein folding, chemical reaction networks, and molecular evolution.APPENDIX A: APPENDIX: THE MEAN ABSOLUTE ERROR OF THE PATHWAYELABORATION METHOD FOR NUCLEIC ACID KINETICSFigures A1 and A2 represent Figure 5 from the main paper by varying only two parametersat a time. See the main paper for the explanation of these figures. In ( a - h ), K = 0 and κ = 0 ns are fixed (the state elaboration step is not used). M A E N =32 N =128 N =512 N =4096 (a) M A E N =32 N =128 N =512 N =4096 (b) M A E N =32 N =128 N =512 N =4096 (c) M A E N =32 N =128 N =512 N =4096 (d) || N =32 N =128 N =512 N =4096 (e) || N =32 N =128 N =512 N =4096 (f) || N =32 N =128 N =512 N =4096 (g) || N =32 N =128 N =512 N =4096 (h)In ( a - h ), K = 256 and κ = 16 ns are fixed. . . . . . . β . . . . . . . M A E N = 32 N = 128 N = 256 N = 512 (i) . . . . . . β . . . . . . . . M A E N = 32 N = 128 N = 256 N = 512 (j) . . . . . . β . . . . . . . . M A E N = 32 N = 128 N = 256 N = 512 (k) . . . . . . β . . . . . . . . M A E N = 32 N = 128 N = 256 N = 512 (l) . . . . . . β | ˆ S | N = 32 N = 128 N = 256 N = 512 (m) . . . . . . β | ˆ S | N = 32 N = 128 N = 256 N = 512 (n) . . . . . . β | ˆ S | N = 32 N = 128 N = 256 N = 512 (o) . . . . . . β | ˆ S | N = 32 N = 128 N = 256 N = 512 (p) Fig A1: The effect of pathway construction with different values of N and β and fixed valuesof K and κ . MAE for ( a ) datasets No. 1,2, and 3, ( b ) dataset No. 4, ( c ) dataset No. 5, and ( d )dataset No. 6. |S| for ( e ) datasets No. 1,2, and 3, ( f ) dataset No. 4, ( g ) dataset No. 5, and ( h )dataset No. 6. For the missing settings, pathway elaboration did not finish within two weekscomputation time. HE PATHWAY ELABORATION METHOD In ( a - h ), N = 128 and β = 0 . are fixed. M A E K =0 K =16 K =256 K =1024 (a) M A E K =0 K =16 K =256 K =1024 (b) M A E K =0 K =16 K =256 K =1024 (c) M A E K =0 K =16 K =256 K =1024 (d) || K =0 K =16 K =256 K =1024 (e) || K =0 K =16 K =256 K =1024 (f) || K =0 K =16 K =256 K =1024 (g) || K =0 K =16 K =256 K =1024 (h)In ( i - p ), N = 128 and β = 0 . are fixed. M A E K =0 K =16 K =256 K =1024 (i) M A E K =0 K =16 K =256 K =1024 (j) M A E K =0 K =16 K =256 K =1024 (k) M A E K =0 K =16 K =256 K =1024 (l) || K =0 K =16 K =256 K =1024 (m) || K =0 K =16 K =256 K =1024 (n) || K =0 K =16 K =256 K =1024 (o) || K =0 K =16 K =256 K =1024 (p) Fig A2: The effect of state elaboration, with different values of K and κ and fixed values of N and β . N = 0 indicates that the states of the pathway are not elaborated. ( a ), ( e ), ( i ), and( m ) correspond to datasets No. 1,2, and 3. ( b ), ( f ), ( j ), and ( n ) correspond to dataset No. 4. ( c )( g ), ( k ), and ( o ) correspond to dataset No. 5. ( d ), ( h ), ( l ), and ( p ) correspond to dataset No.6. For the missing settings, pathway elaboration did not finish within two weeks computationtime. REFERENCES A LLEN , R. J., F
RENKEL , D. and
TEN W OLDE , P. R. (2006). Simulating rare events in equilibrium or nonequi-librium stochastic systems.
The Journal of Chemical Physics
LLEN , R. J., V
ALERIANI , C. and
TEN W OLDE , P. R. (2009). Forward flux sampling for rare event simulations.
Journal of Physics: Condensed Matter MATO , N. M. and S
ONG , G. (2002). Using motion planning to study protein folding pathways.
Journal ofComputational Biology NDERSON , D. F. and K
URTZ , T. G. (2011). Continuous time Markov chain models for chemical reactionnetworks. In
Design and Analysis of Biomolecular Circuits
NDRIEU , C., D E F REITAS , N., D
OUCET , A. and J
ORDAN , M. I. (2003). An introduction to MCMC for machinelearning.
Machine learning NGENENT -M ARI , N. M., G
ARRUSS , A. S., S
OENKSEN , L. R., C
HURCH , G. and C
OLLINS , J. J. (2020). Adeep learning approach to programmable RNA switches.
Nature Communications SMUSSEN , S. and G
LYNN , P. W. (2007).
Stochastic simulation: algorithms and analysis . Springer Science& Business Media.A ZIMZADEH , P. and F
ORSYTH , P. A. (2016). Weakly chained matrices, policy iteration, and impulse control.
SIAM Journal on Numerical Analysis ACKENKÖHLER , M., B
ORTOLUSSI , L. and W
OLF , V. (2019). Bounding Mean First Passage Times in Popula-tion Continuous-Time Markov Chains. arXiv preprint arXiv:1910.12562 .B OLHUIS , P. G., C
HANDLER , D., D
ELLAGO , C. and G
EISSLER , P. L. (2002). Transition path sampling: Throw-ing ropes over rough mountain passes, in the dark.
Annual Review of Physical Chemistry ONNET , G., K
RICHEVSKY , O. and L
IBCHABER , A. (1998). Kinetics of conformational fluctuations in DNAhairpin-loops.
Proceedings of the National Academy of Sciences RAND , M. (2006). Fast low-rank modifications of the thin singular value decomposition.
Linear Algebra andits Applications
ABRIOLU , R., S
KJELBRED R EFSNES , K. M., B
OLHUIS , P. G. and
VAN E RP , T. S. (2017). Foundations and lat-est advances in replica exchange transition interface sampling. The Journal of Chemical Physics AO , Y., G ILLESPIE , D. T. and P
ETZOLD , L. R. (2007). Adaptive explicit-implicit tau-leaping method withautomatic tau selection.
The Journal of Chemical Physics
HEN , Y.-J., G
ROVES , B., M
USCAT , R. A. and S
EELIG , G. (2015). DNA nanotechnology from the test tube tothe cell.
Nature nanotechnology ISSE , I. I., K IM , H. and H A , T. (2012). A rule of seven in Watson-Crick base-pairing of mismatched sequences. Nature Structural & Moleuclar Biology OHEN , M. B., K
ELNER , J., K
YNG , R., P
EEBLES , J., P
ENG , R., R AO , A. B. and S IDFORD , A. (2018). Solvingdirected laplacian systems in nearly-linear time through sparse LU factorizations. In
OHEN -T ANNOUDJI , C., D
AVIES , P. C., D IU , B., L ALOE , F., D UI , B. et al. (1977). Quantum mechanics .John Wiley & Sons.D INH , K. N. and S
IDJE , R. B. (2016). Understanding the finite state projection and related methods for solvingthe chemical master equation.
Physical Biology INH , K. N. and S
IDJE , R. B. (2017). An application of the Krylov-FSP-SSA method to parameter fitting withmaximum likelihood.
Physical Biology OOB , J. L. (1942). Topics in the theory of Markoff chains.
Transactions of the American Mathematical Society OUCET , A. and J
OHANSEN , A. M. (2009). A tutorial on particle filtering and smoothing: Fifteen years later.
Handbook of Nonlinear Filtering YKEMAN , E. C. (2015). An implementation of the Gillespie algorithm for RNA kinetics with logarithmic timeupdate.
Nucleic Acids Research IDELSON , N. and P
ETERS , B. (2012). Transition path sampling for discrete master equations with absorbingstates.
The Journal of Chemical Physics
LBER , R. (2017). A new paradigm for atomically detailed simulations of kinetics in biophysical systems.
Quar-terly Reviews of Biophysics .E SCOBEDO , F. A., B
ORRERO , E. E. and A
RAQUE , J. C. (2009). Transition path sampling and forward fluxsampling. Applications to biological systems.
Journal of Physics: Condensed Matter ELLER , W. (1968).
An Introduction to Probability Theory and its Applications , 3rd ed. Wiley, New York.F LAMM , C., F
ONTANA , W., H
OFACKER , I. L. and S
CHUSTER , P. (2000). RNA folding at elementary stepresolution.
RNA LETCHER , R. (1976). Conjugate gradient methods for indefinite systems. In
Numerical Analysis G EORGOULAS , A., H
ILLSTON , J. and S
ANGUINETTI , G. (2017). Unbiased Bayesian inference for populationMarkov jump processes via random truncations.
Statistics and Computing ILLESPIE , D. T. (1977). Exact stochastic simulation of coupled chemical reactions.
The Journal of PhysicalChemistry ILLESPIE , D. T. (2001). Approximate accelerated stochastic simulation of chemically reacting systems.
TheJournal of Chemical Physics
ILLESPIE , D. T. (2007). Stochastic simulation of chemical kinetics.
Annu. Rev. Phys. Chem. AJIAGHAYI , M., K
IRKPATRICK , B., W
ANG , L. and B
OUCHARD -C ÔTÉ , A. (2014). Efficient continuous-timeMarkov chain estimation. In
International Conference on Machine Learning
ATA , H., K
ITAJIMA , T. and S
UYAMA , A. (2018). Influence of thermodynamically unfavorable secondary struc-tures on DNA hybridization kinetics.
Nucleic Acids Research OFACKER , I. L. (2003). Vienna RNA secondary structure server.
Nucleic Acids Research AVRAKI , L. E., S
VESTKA , P., L
ATOMBE , J.-C. and O
VERMARS , M. H. (1996). Probabilistic roadmaps forpath planning in high-dimensional configuration spaces.
IEEE transactions on Robotics and Automation UEHLMANN , A., M C M ILLAN , K. L. and B
RAYTON , R. K. (1999). Probabilistic state space search. In
UNTZ , J., T
HOMAS , P., S
TAN , G.-B. and B
ARAHONA , M. (2019). The exit time finite state projection scheme:bounding exit distributions and occupation measures of continuous-time Markov chains.
SIAM Journal onScientific Computing A748–A769.K
UWAHARA , H. and M
URA , I. (2008). An efficient and exact stochastic simulation method to analyze rare eventsin biochemical systems.
The Journal of Chemical Physics IÒ , P. and G OLDMAN , N. (1998). Models of molecular evolution and phylogeny.
Genome Research ACHINEK , R. R., O
ULDRIDGE , T. E., H
ALEY , N. E., B
ATH , J. and T
URBERFIELD , A. J. (2014). Pro-grammable energy landscapes for kinetic control of DNA strand displacement.
Nature Communications .M ADRAS , N. N. (2002).
Lectures on Monte Carlo methods . American Mathematical Soc.M C G IBBON , R. T. and P
ANDE , V. S. (2015). Efficient maximum likelihood parameterization of continuous-timeMarkov processes.
The Journal of Chemical Physics
ETROPOLIS , N., R
OSENBLUTH , A. W., R
OSENBLUTH , M. N., T
ELLER , A. H. and T
ELLER , E. (1953).Equation of state calculations by fast computing machines.
The Journal of Chemical Physics ETZNER , P., S
CHÜTTE , C. and V
ANDEN -E IJNDEN , E. (2009). Transition path theory for Markov jump pro-cesses.
Multiscale Modeling & Simulation ORRISON , L. E. and S
TOLS , L. M. (1993). Sensitive fluorescence-based thermodynamic and kinetic measure-ments of DNA hybridization in solution.
Biochemistry UNSKY , B. and K
HAMMASH , M. (2006). The finite state projection algorithm for the solution of the chemicalmaster equation.
The Journal of Chemical Physics
ELDER , J. A. and M
EAD , R. (1965). A simplex method for function minimization.
The Computer Journal ARKS , M. L., D E S TURLER , E., M
ACKEY , G., J
OHNSON , D. D. and M
AITI , S. (2006). Recycling Krylovsubspaces for sequences of linear systems.
SIAM Journal on Scientific Computing EIMANN , P., S
CHMID , G. and H
ÄNGGI , P. (1999). Universal equivalence of mean first-passage time andKramers rate.
Physical Review E R1.R
IPLEY , B. D. (2009).
Stochastic simulation . John Wiley & Sons.R
UBINO , G. and T
UFFIN , B. (2009).
Rare event simulation using Monte Carlo methods . John Wiley & Sons.S
ARICH , M., B
ANISCH , R., H
ARTMANN , C. and S
CHÜTTE , C. (2014). Markov state models for rare events inmolecular dynamics.
Entropy CHAEFFER , J. M. (2013). Stochastic simulation of the kinetics of multiple interacting nucleic acid strands, PhDthesis, California Institute of Technology.S
CHAEFFER , J. M., T
HACHUK , C. and W
INFREE , E. (2015). Stochastic Simulation of the Kinetics of MultipleInteracting Nucleic Acid Strands. In
DNA Computing and Molecular Programming, Lecture Notes in ComputerScience
HAHABUDDIN , P. (1994). Importance sampling for the simulation of highly reliable Markovian systems.
Man-agement Science IDJE , R. B. and V O , H. D. (2015). Solving the chemical master equation by a fast adaptive finite state projectionbased on the stochastic simulation algorithm. Mathematical Biosciences
IMMONS , G. F. (1972).
Differential equations with applications and historical notes . CRC Press. S INGHAL , N., S
NOW , C. D. and P
ANDE , V. S. (2004). Using path sampling to build better Markovian state mod-els: predicting the folding rate and mechanism of a tryptophan zipper beta hairpin.
The Journal of ChemicalPhysics
UHOV , Y. and K
ELBERT , M. (2008).
Probability and Statistics by Example: Volume 2, Markov Chains: A Primerin Random Processes and Their Applications . Cambridge University Press.S UTTON , R. S. and B
ARTO , A. G. (2018).
Reinforcement learning: An introduction . MIT Press.T
ANG , X. (2010). Techniques for modeling and analyzing RNA and protein folding energy landscapes, PhDthesis, Texas A & M University.T
ANG , X., K
IRKPATRICK , B., T
HOMAS , S., S
ONG , G. and A
MATO , N. M. (2005). Using motion planning tostudy RNA folding kinetics.
Journal of Computational Biology HOMAS , S., T
APIA , L., E
KENNA , C., Y EH , H.-Y. C. and A MATO , N. M. (2013). Rigidity analysis for proteinmotion and folding core identification. In
Workshops at the Twenty-Seventh AAAI Conference on ArtificialIntelligence .T HUBAGERE , A. J., L I , W., J OHNSON , R. F., C
HEN , Z., D
OROUDI , S., L EE , Y. L., I ZATT , G., W
ITTMAN , S.,S
RINIVAS , N., W
OODS , D., W
INFREE , E. and Q
IAN , L. (2017). A cargo-sorting DNA robot.
Science .V AN K AMPEN , N. G. (1992).
Stochastic processes in physics and chemistry . Elsevier.V IRTANEN , P., G
OMMERS , R., O
LIPHANT , T. E., H
ABERLAND , M., R
EDDY , T., C
OURNAPEAU , D.,B
UROVSKI , E., P
ETERSON , P., W
ECKESSER , W., B
RIGHT , J. et al. (2020). SciPy 1.0: fundamental algo-rithms for scientific computing in Python.
Nature Methods
EINAN , E., R EN , W. and V ANDEN -E IJNDEN , E. (2002). String method for the study of rare events.
PhysicalReview B ETMUR , J. G. and D
AVIDSON , N. (1968). Kinetics of renaturation of DNA.
Journal of Molecular Biology HITT , W. (2006). Continuous-time Markov chains.
Dept. of Industrial Engineering and Operations Research,Columbia University, New York .Z ADEH , J. N., S
TEENBERG , C. D., B
OIS , J. S., W
OLFE , B. R., P
IERCE , M. B., K
HAN , A. R., D
IRKS , R. M.and P
IERCE , N. A. (2011). NUPACK: analysis and design of nucleic acid systems.
Journal of ComputationalChemistry HANG , J. X., F
ANG , J. Z., D
UAN , W., W U , L. R., Z HANG , A. W., D
ALCHAU , N., Y
ORDANOV , B., P E - TERSEN , R., P
HILLIPS , A. and Z
HANG , D. Y. (2018). Predicting DNA hybridization kinetics from sequence.
Nature Chemistry OLAKTAF , S., D
ANNENBERG , F., R
UDELIS , X., C
ONDON , A., S
CHAEFFER , J. M., S
CHMIDT , M.,T
HACHUK , C. and W
INFREE , E. (2017). Inferring Parameters for an Elementary Step Model of DNA Struc-ture Kinetics with Locally Context-Dependent Arrhenius Rates. In
DNA Computing and Molecular Program-ming, Lecture Notes in Computer Science
OLAKTAF , S., D
ANNENBERG , F., W
INFREE , E., B
OUCHARD -C ÔTÉ , A., S
CHMIDT , M. and C
ONDON , A.(2019). Efficient Parameter Estimation for DNA Kinetics Modeled as Continuous-Time Markov Chains. In
DNA Computing and Molecular Programming, Lecture Notes in Computer Science
UCKERMAN , D. M. and C
HONG , L. T. (2017). Weighted ensemble simulation: review of methodology, appli-cations, and software.
Annual Review of Biophysics46