We consider the problem of computing first-passage time distributions for reaction processes modeled by master equations. We show that this generally intractable class of problems is equivalent to a sequential Bayesian inference problem for an auxiliary observation process. The solution can be approximated efficiently by solving a closed set of coupled ordinary differential equations (for the low-order moments of the process) whose size scales with the number of species. We apply it to an epidemic model and a trimerization process and show good agreement with stochastic simulations. DOI: 10.1103/PhysRevLett.119.210601 Many systems in nature consist of stochastically interacting agents or particles. Such systems are frequently modeled as reaction processes whose dynamics are described by master equations [1] . There are several examples of the stochastic modeling of reaction processes in the fields of systems biology [2, 3] , ecology [4] , epidemiology [5] , social sciences [6] , and neuroscience [7] . The mathematical analysis of such stochastic processes, however, is highly nontrivial. A particularly important quantity of interest is the firstpassage time (FPT), that is, the random time it takes the process to first cross a certain threshold [8, 9] . FPT distributions play a crucial role both in the theory of stochastic processes and in their applications across various disciplines, as they allow us to investigate quantitatively the uncertainty in the emergence of system properties within a finite time horizon. For example, the time it takes cells to respond to external signals by expressing certain genes may be modeled as a FPT problem. Different characteristics of this first time distribution, for example, the variance of the FPT, may represent evolutionarily different strategies that organisms adopt to filter fluctuations in the environment [10] [11] [12] . Examples from other disciplines include the extinction time of diseases in epidemic models or the time it takes to form a certain number of polymers in polymerization processes. FPTs for stochastic processes have been of interest in statistical physics for many decades [13] . For certain random walk or spatial diffusion processes, analytic solutions have been derived [14] [15] [16] . Recently, analytic results have been found for effective one-dimensional diffusion processes to a target [17] [18] [19] . For multidimensional diffusion processes to small targets, approximate solutions have been derived using singular perturbation methods and matched asymptotic expansions [20] [21] [22] [23] . The problem of computing FPT distributions for reaction processes modeled by master equations, however, is much less explored. Generally, no tractable evolution equations exist except for one-variable, one-step processes [1, 13] or certain linear and/or catalytic processes [24] [25] [26] . For singletime properties of the underlying master equation, efficient approximation methods exist relying on continuous state spaces [27] , but it is not clear how to extend them for the computation of FPTs. Spectral methods constitute efficient approximations for small systems [28, 29] . Since these methods typically scale with the size of the state space, they are not applicable to large systems. Some existing FPT approaches for master equations consider rare events in single-species systems and/or mean FPTs only [30] [31] [32] [33] . In this Letter, we approach the problem of computing FPTs from a novel perspective. We show that the FPT problem can be formulated exactly as a Bayesian inference problem. We achieve this by introducing an auxiliary observation process that determines whether the process has crossed the threshold up to a given time. This novel formulation allows us to derive an efficient approximation scheme that relies on the solution of a small set of ordinary differential equations. We will use this approximation to analyze the FPT distributions in several nontrivial examples. We focus on reaction networks with discrete state spaces modeled by master equations, but the developed method can also be applied to processes with continuous state spaces modeled by Fokker-Planck equations. The standard approach to compute the FPT of a process x t to leave a certain region C is to compute the survival probability Z ½0;t , that is, the probability that the process remains in C on the time interval ½0; t [13]. The FPT distribution is then given by the negative time derivative of Z ½0;t . The latter can be written as a path integral over the process with an absorbing boundary of C [13] . Equivalently, one can reweigh the unconstrained process by an indicator
Schnoerr et al. REF also employ a moment closure approximation and further approximate threshold functions to derive an approximate first passage time distribution.
Efficient low-order approximation of first-passage time distributions
Physical review letters
Physical review letters
Physics Medicine