Exact first passage time distribution for second-order reactions in chemical networks
Changqian Rao David Waxman Wei Lin Zhuoyi Song




             

§ ABSTRACT 

The first passage time (FPT) is a generic measure that quantifies when a random quantity reaches a specific state. We consider the FTP distribution in nonlinear stochastic biochemical networks, where obtaining exact solutions of the distribution is a challenging problem. Even simple two-particle collisions cause strong nonlinearities that hinder the theoretical determination of the full FPT distribution. Previous research has either focused on analyzing the mean FPT, which provides limited information about a system, or has considered time-consuming stochastic simulations that do not clearly expose causal relationships between parameters and the system's dynamics. This paper presents the first exact theoretical solution of the full FPT distribution in a broad class of chemical reaction networks involving A + B → C type of second-order reactions. Our exact theoretical method outperforms stochastic simulations, in terms of computational efficiency, and deviates from approximate analytical solutions. Given the prevalence of bimolecular reactions in biochemical systems, our approach has the potential to enhance the understanding of real-world biochemical processes.



 

A Data Selection Approach for Enhancing Low Resource Machine Translation Using Cross-Lingual Sentence Representations



    
Nidhi Kowtal * 

SCTR's Pune Institute of Computer Technology

Pune, India 

kowtalnidhi@gmail.com

Tejas Deshpande * 

SCTR's Pune Institute of Computer Technology

Pune, India 

tejasdeshpande1112@gmail.com
 

Raviraj Joshi

Indian Institute of Technology Madras, India

L3Cube Labs, Pune 

Pune, India 

ravirajoshi@gmail.com

    September 9, 2024
==============================================================================================================================================================================================================================================================================================================================================







§ INTRODUCTION


The first passage time (FPT) is a fundamental concept that is used to analyze the behavior and dynamics of stochastic processes<cit.>. In biochemical reaction networks, the FPT is a key quantity that refers to the time when a specific event or state first occurs within the network<cit.>.
Examples of the specific event include reaction completion<cit.>, binding or unbinding events<cit.>, protein translocation<cit.>, and state transitions<cit.>. Analyzing the FPT for these events is not just a theoretical exercise; the FPT can provide detailed insight into the timing, efficiency, and reliability of the underlying biochemical reaction system<cit.>. This insight is crucial because it enables not only the understanding of regulatory mechanisms<cit.>, but also their manipulation<cit.>, thereby opening up new possibilities in biochemistry and chemical kinetics.  

When counts of the constituent molecules are low, stochasticity and discreteness are inescapable features of chemical kinetics<cit.>. In this context, the stochastic properties of the FPT require a characterization in terms of its probability distribution. However, past theoretical work has primarily focused on deriving the mean FPT and the global FPT<cit.>. This is because the FPT distribution is hard to measure experimentally; the stochastic timings are disguised in cell population measurements due to cell-cell variabilities, and precisions in single-cell measurements are limited by experimental technologies<cit.>. Attention has shifted recently to focus on obtaining the full FPT distribution beyond the mean<cit.>. This distribution provides much more information about the underlying biochemical system. For example, Thorneywork et al. demonstrate that a purely dynamic measurement of the full FPT distribution uncovers that a short-time, power-law regime of the distribution, rather than the mean FPT, reflects the number of intermediate states in an underlying potential energy landscape<cit.>. Therefore, mathematical methodologies for estimating and analyzing the FPT distribution are highly desirable, and complement advances in measurement technologies of the FPT distribution in biochemical systems.

Traditionally, estimating the FPT distributions depends on solving the underlying chemical master equations (CMEs)<cit.>, which are the primary modelling approach of stochastic biochemical systems<cit.>. Typically, CME solutions can be simulated<cit.>, solved approximately<cit.>, or solved exactly<cit.>. Simulation approaches, such as the Gillespie algorithm<cit.>, approximate the solutions of the CME by generating many realizations of the associated Markov process. Such simulations trade time efficiency for accuracy, since producing tens of thousands of sample paths, to expose the underlying distribution, takes time<cit.>. 

Approximate methods were developed to solve CMEs, trading off estimation accuracy for time efficiency<cit.>. There are two main classes of such methods: (i) closure schemes<cit.> and (ii) linear mapping methods<cit.>. Under closure schemes, the solutions to CMEs are obtained by approximating higher-order moments of the solution by nonlinear functions of lower-order moments, thereby leading to tractable equations. Under linear mapping, bimolecular reactions are approximated by first-order reactions, allowing simpler, solvable systems. However, even with such approximations, it is a formidable task to elucidate the causal mechanisms regulating FPT distributions, as exhaustive parameter searches are usually unrealistic, and the complex compensatory effects of parameter variations are difficult to clarify<cit.>. 

Exact theoretical expressions for the FPT distributions are highly desired, as they may be essential to identifying and quantifying the causal regulatory mechanisms<cit.>. However, this depends on determining time-dependent solutions for the CMEs, which are only known for specific cases<cit.>.To the best of our knowledge, general, exact time-dependent solutions of CMEs are known only for simple reaction systems - with zero and first-order reactions, i.e., linear systems<cit.>.

Time-dependent solutions of CMEs are still challenging for general, non-linear, biochemical reaction networks<cit.>. Even widely presented bimolecular reactions, which constitute one of the simplest core building blocks of a biochemical reaction network, involve highly nonlinear models, making exact CME solutions non-tractable<cit.>. As a result, for biochemical networks with second-order reactions, time-varying CME solutions are only available for specific cases, such as highly simplified systems with only a few states<cit.>, or reversible bimolecular reaction in which the transition matrix has a tridiagonal format<cit.>.  

General time-dependent CME solutions are still unknown for systems with bimolecular reactions<cit.>. Matters are even more challenging if there is a large state space and if non-constant reaction rates are involved<cit.>. Generally, second-order reactions can be grouped into two types: A + A → C and A + B → C, and mathematical analysis treats these two types of second-order reactions separately<cit.>. In earlier work, we derived an exact FPT distribution for a general class of biochemical networks involving an A + A → C second-order reaction downstream of two zero/first-order reactions<cit.>. In the present work, we report the first exact distribution of the FPT for a general class of chemical reaction networks with an A + B → C second-order reaction that is downstream of various zero/first-order reactions.

We note that Anderson et al. have recently provided the first exact time-dependent distribution for a general class of reaction networks with higher-order complexes<cit.>. They find that the time-varying solutions of the CMEs will maintain a Poisson-product form. However, this result requires that: (i) the chemical reaction system initially has such a Poisson-product form, and (ii) the system has to satisfy a dynamical and restricted (DR) condition. This DR condition requires that, for any higher-order reactant pairs, the production and consumption are balanced. Generally, the DR condition is very restrictive, and only applies to specific reaction formats with specific kinetic rates and initial conditions. In this work, we present results that are not subject to the DR condition - by allowing the mean of molecular numbers to follow a stochastic process, as described by a stochastic differential equation (SDE) 
rather than an ordinary differential equation (ODE).

Our exact analytical results are not only novel and efficient, but also highly practical. They are much more time efficient than stochastic simulations, and are much more accurate than traditional approximation methods, such as linear mapping methods (LMA)<cit.>. As we demonstrate, our methods have applications in diverse chemical reaction networks from different areas. These include gene regulatory networks<cit.> and multi-step transition models<cit.>. Our results thus have broad applicability and real-world significance. 



§ PROBLEM FORMULATION


We consider a biochemical system with N chemical species S = [S_1,S_2,...,S_N]^⊤ whose molecules can undergo M+1 chemical reactions, R_m, with m=0,1,2,3,....,M. There are M zero or first-order reactions followed by one second-order reaction, which is of the type  A + B → C. 

We label the second-order reaction by m=0 and the other reactions (zero- or first-order) by m=1,2,...M. The molecules A and B are the reactants of the second-order reaction, which can be any two different molecular species. We denote these two reactants as S_1 and S_2 without loss of generality. 

We shall use the following general notation for such a biochemical system: 


    y_m·Sa_m(t)→   y'_m·S,
    
    S_1+S_2a_0(t)→   *

where m= 1,…,M label the m'th zero or first-order reaction, while y_m and y_m' are the stoichiometric coefficient vectors of the m'th reaction. 
The reaction rate 'constants' are written a_m(t) and a_0(t), and we have assumed these vary with the time, t. The resulting analysis thus includes a broad class of non-linear time-varying biochemical systems with second-order reactions.


Because, in the above biochemical system, the first m reactions are zero or first-order, we require that the sum of the absolute values of all components of y_m, and separately of y_m', should not exceed one; we write these conditions as y_m_1≤ 1 and y'_m_1≤ 1, respectively.

To begin a derivation of the exact FPT distribution of the second-order reaction, we first define it mathematically. We adopt a standard approach by adding an additional species S_0, the number of which counts the time of occurrence of the second-order reaction<cit.>. We denote the number of S_0 molecules by x_0 ≡ x_0(t), then x_0 is a non-decreasing function of t. Thus, we modify the system in (<ref>) to: 

    y_m·Sa_m(t)→   y'_m·S,
    
    S_1+S_2a_0(t)→   S_0.


If at any time t the second order reaction has not occurred, we have x_0(t)=0, and
equivalently the FPT exceeds t. Therefore, the complementary cumulative probability distribution of the FPT can be represented as: 


    P(FPT>t) = P(x_0(t) = 0).




§ RESULT



 §.§ Analytical representation for the auxiliary chemical master equation


For the system given in (<ref>), we adopt the following notation:


    
  * 𝐗=[[ x_0;   𝐱 ]] denotes the complete state vector of the system, representing the numbers of different species of S^*=[[ S_0;   S ]]. The quantity 𝐱=[x_1,⋯,x_N]^⊤ denotes the numbers of different species of S.


  * Y_m=[
         y_m,0,
         y_m 
    ] and Y'_m=[
         y'_m,0,  
         y'_m 
    ] are stoichiometric coefficient vectors of the m'th reaction in system (<ref>), and y_m=[y_m,1,y_m,2,…,y_m,N], y'_m=[y'_m,1,y'_m,2,…,y'_m,N] denote parts of stoichiometric coefficient vectors, excluding y_m,0 and y'_m,0 for the auxiliary species S_0.


  * a_0(t) denotes the rate `constant' at time t of the second-order reaction,
while a_m(t) (m=1,…,M) denotes the rate `constant' at time t of the m'th reaction
 (which is either a zero or first-order reaction).



The CMEs for the system given in (<ref>) correspond to a differential equation for the probability distribution of X at time t<cit.>: 



    d P(X, t )/ d t =∑_k=0^M[ P(X-v_k, t ) c_k(X-v_k,t)- P(X, t ) c_k(X,t)],


where v_k is the transition vector for the k'th reaction, and c_k(X,t) is the propensity function for the k'th reaction, i.e., the probability that the k'th reaction occurs in state X, at time t. The 
quantity c_k(X,t) equals a_k(t)X^Y_k, where, for any vectors with d components, a vector to the power of a
vector is defined by 𝐮^𝐯 def=∏^d_i=1u_i^v_i and we adopt the convention that 0^0=1. Therefore, for k=0, the reaction is given in (<ref>), thus, c_0(X,t)=a_0(t)x_1x_2. For k=1,…, M, the reactions are first or lower-order, and thus c_k(X,t), k=1,…, M are linear functions of states.

We next present a theorem that gives the exact solution of (<ref>) 
in terms of the following variables/notation:




  * λ(t)=[λ_1(t),…
,λ_N(t)]^⊤ is a column vector containing the mean of all numbers of the species in S at
any time t (t≥0), and Λ(t)= [[ λ_0(t);   λ(t) ]] is a column vector containing the mean of all numbers of the species in S^* at
any time t (t≥0).


  * 

Λ(0)= [[ λ_0(0);   λ(0) ]]denotes the initial mean number of each species. 


  * ℳ_1 denotes an N× N matrix with components
η_i,j^1, and ℳ_2 denotes an N×1 vector with
components η_i^2, where the components are given by:

    η_i,j^1=∑_m:y_m=e_ia_m(t)(y_m,j^'
    -y_m,j),  η_i^2=∑_m:y_m=0a_m
    (t)(y_m,i^'-y_m,i),

in which e_i is a vector where only element i is 1, and all other elements
are 0, while y_m,i^' and y_m,i are the i'th component
of y_m^' and y_m, respectively,


  * 𝒩_1 and 𝒩_2 denote N× N matrices with
all elements zero, except the upper left 2×2 block, with

    𝒩_1=(
    [  a_S    0    ⋯;    0 -a_S    ⋯;    ⋮    ⋮    ⋱ ])  ,  𝒩_2=(
    [ ia_S    0    ⋯;    0 ia_S    ⋯;    ⋮    ⋮    ⋱ ]),

where 
a_S=√(2a_0)/2. 





The theorem reads as follows.



For the system in (<ref>), providing:

    
  * the variables λ and λ_S obey the stochastic 
differential equations (SDEs): 

    dλ   =(ℳ_1λ+ℳ_2) dt + 𝒩_1λ dW^1_t+𝒩_2λ dW^2_t,
     dλ_S   = ( a_S λ_1-a_S λ_2 ) dW^1_t + (  ia_S λ_1+ ia_S λ_2 ) dW^2_t,

subject to λ_S(0)=0, and some initial condition λ(0),
    
  * the new variable λ_0 obeys the SDE: 

    dλ_0=a_0λ_1λ_2 dt,

subject to λ_0(0)=0,


then if the initial condition of (<ref>) is a distribution of the Poisson-product form: 

    P(X,0)   =Λ(0)^X/X !exp(-λ(0))=∏_i=0^N λ_i(0)^x_i/x_i !exp(-λ_i(0)),

the solution of (<ref>) is

    P(X,t)=<Λ(t)^X/X!exp(-Λ(t))exp(λ_0(t)+λ_S(t))>,

 where <...> denotes an expectation operation.



Theorem <ref> is proved in Appendix A. 

Theorem <ref> gives an exact theoretical expression for the CME solution of a general class of nonlinear biochemical reaction networks with A + B → C type of second-order reactions. Compared to the previous analytical solutions of CMEs by Anderson et al.<cit.>, the key feature is that Theorem <ref> does not require the underlying system to satisfy a dynamically restricted (DR) condition, which considerably constrains the system's structure and parameters. However, the systems we analyze here, as represented by (<ref>), are not subject to the DR condition. This indicates a broader applicability of Theorem <ref> than previous exact theoretical results in the literature. 

The key reason why Theorem <ref> can break the DR condition is that we allow the mean of each molecular species (contained in λ and λ_S), to follow stochastic processes, as described by the SDEs of (<ref>)), rather than being smoothly changing continuous variables, with no randomness, whose dynamics are governed by ordinary differential equations (ODEs). As a result, the time-dependent CME solution can be written as the average of a distribution over all λ paths. More details are given in Appendix <ref>.

Another significant aspect of Theorem <ref> is its applicability to time-varying systems where the reaction rates are not constant. This is a particularly complex problem, as solving CMEs with time-varying rates is much more challenging than those with constant reaction rates<cit.>. To date, except for our earlier work that provides an exact FPT distribution for a specific type of nonlinear biochemical reaction network with an A + A → C type of second-order reactions<cit.>, we have not encountered any work that offers exact theoretical results for FPTs with A + B → C type of second-order reactions and time-varying reaction rates. This suggests that Theorem <ref> will lead to new avenues of research in this area. 



 §.§ Analytical expression for the FPT distribution


One approach to obtain the solution of Eq. (<ref>) is to calculate the expectation by simulation, i.e., by independently solving the SDEs (for λ, λ_S and λ_0) many times and then averaging. This may be very time-consuming. However, Theorem <ref> can lead to an exact 
form of the FPT distribution, as given in (<ref>), which can then be numerically evaluated.


The complementary cumulative probability distribution of the FPT

equals the probability of occurrence of states with x_0=0:

    P(FPT>t)   = ∑_x|x_0=0P(X,t)
    =<exp(λ_S)>.


Corollory <ref> follows by substituting P(x,x_0,t) into (<ref>). 

Equation (<ref>) is a highly compact exact theoretical result. Next, we present numerical methods for calculating it. 






 §.§ Numerical approximations for the FPT distribution
 

A direct approach to calculating <exp(λ_S)> in (<ref>) would require solving an infinite set of coupled ODEs (see Appendix <ref> for details). To avoid such an issue, we shall introduce a numerical approximation method based on a Padé approximant of  <exp(λ_S)>. This requires the calculation of moments of λ_S(t). As we can show, the n'th moment, <λ_S^n>, has a closed-form expression. This follows because <λ_S^n> is governed by a finite set of coupled ODEs, which results because  λ ( (<ref>)) and λ_S ( (<ref>)) follow linear SDEs.

We proceed by constructing a function H(s,t)=<exp(sλ_S(t))>; the quantity 
required for Eq. (<ref>) is given by H(1,t). 
To approximate H(1,t), we first obtained a Padé approximant of the function H(s,t), denoted as T(s,t), and then set s=1 to calculate T(1,t).  

We determined a Padé approximant of H(s,t) using the following procedure<cit.>: 


  * 
Construct the Maclaurin series of H(s,t) in s, truncated at the Ñ'th
term, i.e. T_Ñ(s,t)= ∑_n=0^Ñs^n/n!×∂^n/∂ s^nH(s,t)|_s=0 = ∑_n=0^Ñb_n(t)s^n. The larger Ñ is, the better the Padé approximate is, but at the cost of a more time-consuming algorithm. 


  * Calculate b_n(t) by determining <λ_S(t)^n>; since ∂^n/∂ s^nH(s,t)|_s=0= <λ_S(t)^n>, b_n(t) is known if <λ_S(t)^n> is calculated. See 
below for details on how to calculate <λ_S(t)^n>.


  * Find the Padé approximant by determining two polynomials P^*_L̃(s,t) and Q^*_Ñ-L̃(s,t), such that the Maclaurin series of P^*_L̃(s,t)/Q^*_Ñ-L̃(s,t), truncated at the Ñ'th term, equals T_Ñ(s,t). 


  * Equate the coefficients of P^*_L̃(s,t)/Q^*_Ñ-L̃(s,t) with that of the corresponding polynomial terms of T_Ñ(s,t), then P^*_L̃(s,t) and Q^*_Ñ-L̃(s,t)  can be uniquely set by solving Ñ+1 algebraic equations. For example, if P^*_L̃(s,t)=p_0(t)+p_1(t)s+p_2(t)s^2+…+p_L̃(t)s^L̃, and 
if
Q^*_Ñ-L̃(s,t) = 1+q_1(t)s+q_2(t)s^2+…+q_Ñ-L̃(t)s^Ñ-L̃, 
then the algebraic equations are:

    b_0=p_0   
    
    b_1+b_0q_1=p_1   
    
    b_2+b_1q_1+b_0q_2=p_2   
    ⋮
    
    b_L̃+b_L̃-1q_1+…+b_0q_L̃=p_L̃   
    
    b_L̃+1+b_L̃q_1+…+b_0q_L̃-1=0   
    
    b_Ñ+b_L̃q_1+…+b_0q_2L̃-Ñ=0   ,

where q_n=0 for n<0 or n>Ñ-L̃.






The key to obtaining a Padé approximant of H(s,t) is the determination of <λ_S(t)^n>. To achieve this, we adopt the following procedure<cit.>. 


    
  * Differentiate λ_S(t)^n, using Ito's rule: 
d(λ_S^n)=nλ_S^n-1 dλ_S+ n(n-1)/2λ_S^n-2( dλ_S)^2, and substitute dλ_S and ( dλ_S)^2 from the SDEs 
in (<ref>). The resulting right-hand side is a polynomial in λ_S^l_0λ^l. 


  * Keep differentiating all new terms of the form λ_S^l_0λ^l and 
substitute dλ_S along with the SDEs in (<ref>), until all terms in the right-hand side of the equations are known. 


  * Average 
the equations for d(λ_S^n) and d(λ_S^l_0λ^l); a set of ODEs are obtained, whose solution yields <λ_S(t)^n>. This set of ODEs is of finite dimension as a consequence of linearity of the SDEs in (<ref>).







 §.§ Theory validation


We first applied our theory to an exemplar biochemical network that composes four zero/first-order reactions upstream of a second-order reaction:

    [                     ∅ [   a_2(t) ]   a_1(t)                S_1 
    ,;                     ∅ [   a_4(t) ]   a_3(t)                S_2 
    ,;             S_1 + S_2                                        S_0, ]


We shall determine the FPT distribution via Eq. (<ref>)

We validated our theoretical results by comparing them to stochastic simulation algorithm (SSA) 
results, which follow from the Gillespie algorithm<cit.>. To the best of our knowledge, there are no other exact results for comparison  (our method is the first that provides an exact FPT distribution for biochemical systems with the network structure of (<ref>)). Thus, to illustrate the method's effectiveness, we compared our method with another recently developed approximation method, namely, the linear-mapping method (LMA)<cit.>, which derives CME solutions of a linearly approximated biochemical network using a mean-field assumption. We chose the LMA as a benchmark because: 1) our method and the LMA method require similar computational time, while moment closure schemes are significantly slower; 2) both our method and the LMA method can be used for biochemical systems with time-varying reaction rates.  



Our method is more accurate and robust than the LMA method, across a broader range of parameters for the network in (<ref>) with time-constant reaction rates. The LMA and our method can provide accurate FPT distributions in some parameter ranges (Fig.<ref> A&B). However, as the second-order reaction goes faster or with more reactant molecules, our method can be more accurate (Fig.<ref> D vs. Fig.<ref> C and Fig.<ref> F vs. Fig.<ref> E ). We used the normalized Wasserstein distance (W-distance) to measure the error between the SSA simulated and the analytical FPT distributions gained from the LMA method or our method. The W-distance is normalized against the standard deviation of the SSA-simulated FPT distribution so that the normalized W-distance is dimensionless, allowing comparisons through various time scales. The heat map in Fig. <ref> E&F were fitted upon 49 error points, obtained from the results of a (7 × 7) set of parameters. 

We further applied our method to the network in (<ref>) with time-varying reaction rates. We restricted consideration to a second-order reaction rate (a_0) that varies with time in a sinusoidal fashion. We tested four cases, where a_0 can be small (Fig. <ref> A-B) or large (Fig. <ref> C-D), or the rate of a_0 can be fast (Fig. <ref> A&C) or slow (Fig. <ref> B&D). Again, with larger a_0 values, our method is more accurate (Fig. <ref> C-D), whereas the rate at which a_0 changes has a slight influences on accuracy (Fig. <ref> A&C vs. Fig. <ref> B&D). 





 §.§ Applications to real biochemical networks


Our method can be applied to various biochemical systems in different fields. Here, we demonstrate two applications in genetic regulation networks (GRNs)<cit.> and a multistep reaction pathway in Ras activation by a protein, called Son Of Sevenless (SOS)<cit.>. 



FPT distribution of a GRN: Gene expression is a fundamental process allowing organisms to create life machinery<cit.>. It is a highly regulated process that involves 
the coordinated action of regulatory proteins that bind to specific DNA sequences to activate or repress gene transcription<cit.>. These genetic regulation networks are crucial in cell differentiation, development, and disease<cit.>. Two decades of research have shown that gene-gene or protein-gene interactions are inherently stochastic, leading to cell-cell variations in mRNA and protein levels<cit.>. Thus, analyzing the stochastic timings of the GRNs can be important for understanding cellular phenomena and their function<cit.>.  

Here, we analyze a simple GRN system (Fig.<ref> A), which has only one gene that can be in two states, inactive (G) and active (G^*)<cit.>. A protein P, which is generated and degraded dynamically, 
can bind to G to activate it, and the activated gene, G^*, may become deactivated after some time. The GRN involves three chemical reactions, including a second-order reaction, as shown in (<ref>). Starting with some substances that can generate protein P and hence G^*, we ask at what time does the protein P activate the gene? 


    [      ∅ [a2]a1      P;    G^*            G,;   G +P          G^*,;        ]


To answer this question, we must derive the FPT distribution of the second-order reaction,  in (<ref>). We applied both the LMA method and our method to solve this problem. Our method is accurate across a wide range of parameters (Fig.<ref> B), whereas the LMA can lead to significant errors when the reaction rate of the second-order reaction is high (Fig.<ref> C). 






FPT distribution of a multistep reaction pathway: A second application is a multistep reaction pathway, which is present in many biological and chemical processes, such as enzymatic reactions, the 
folding and unfolding of RNA molecules, and the conformation changes of ion channels<cit.>. 

We chose a specific multistep reaction pathway<cit.>, the Ras activation by SOS (Fig. <ref>), to illustrate the effectiveness of our method. SOS is a Ras guanine nucleotide exchange factor (GEF) that plays a central role in numerous cellular signaling pathways, such as the epidermal growth factor receptor and T-cell receptor signaling. 

SOS is autoinhibited in cytosol and activates only after recruitment to the membrane. There are two phases in the activation (Fig. <ref>A) the release of autoinhibition at the membrane through several membrane-mediated intermediates, by a sequence of first-order chemical reactions and 2) the binding of Ras at the allosteric site of SOS by a second-order reaction, which enables the activation of Ras<cit.>. What is the time when Ras is activated, i.e., what is the FPT distribution of the second-order reaction of (<ref>)?  

A recent study analyzed the FPT distribution of the first phase of Ras activation<cit.>. It revealed how a faster Ras activation timescale is possible by using much slower activating SOS molecules through multistep reactions. The study resolves the odd discrepancy between the long timescale of individual SOS molecules and the much shorter timescale of Ras activation. More importantly, it demonstrates how rare and early SOS activation events dominate the macroscopic reaction dynamics, implying that the full FPT distribution is required for understanding of this phenomenon, rather than just the mean.

Nevertheless, the analysis presented in<cit.> only focused on first-order multistep reactions. In contrast, our method can include the second-order Ras binding phase of the activation process, pushing the analysis one step further. Without loss of generality, we used two steps in the multistep reaction phase of the SOS activation. The system of reactions is represented in (<ref>), where protein S_0 first needs to change to S_1 and then S_* before binding with a protein R to activate Ras. The S_0 and S_1 proteins degrade, with some time constants, and the activated R_* can deactivate to R. 




    [      ∅ [a2]a1    S_0;    S_1           ∅
    ,;    S_0        S_1 
    ,;    S_1          S_*,;    R^*            R,;  R+S_2            *,;        ]


We applied the LMA, along with our methods, to derive the FPT distribution of the binding reaction between S_* and R in (<ref>) across a wide range of parameters. We compared the FPT distributions with the SSA-simulated ones and measured the errors by a scaled W-distance. The W-distance is scaled by the standard variable of the SSA-simulated FPT distribution, denoted by σ. The scale is to make the W-distance a dimensionless quantity, enabling comparisons across systems with different timescales. The errors are all minimal for all parameter sets tested, with the maximum error at a scale of 10% of σ (Fig. <ref> B). However, the LMA can induce a largely inaccurate FPT distribution with a W-distance of 90% of σ (Fig. <ref> C), 9 times of error compared to our results. 



§ DISCUSSION



In this paper we have derived an exact theoretical result for a general class of biochemical networks that involve nonlinearities caused by two-particle collisions. The networks include a second-order reaction of the A+B → C type that is downstream of a series of zero or first-order reactions. We have derived the exact theoretical distribution 
for the earliest time the second-order reaction occurs - the first passage time (FPT) distribution of the second-order reaction.

Exact theoretical first passage time distributions have not, previously, been derived for such a system 
with a broad range of applications. This is due to a lack of time-dependent solutions 
of the distribution for the system's states, as described by a chemical master equation (CME). It is known that solutions of the CME, for general biochemical systems, 
only hold for time-constant, 
linear systems composed only of zero and first-order chemical reactions. Even the simplest two-particle collisions cause strong nonlinearities that hinder development of theoretical solutions of the CME for a general system. Additionally, there are further difficulties 
caused by the time-varying reaction rates. 

However, theoretical results for the solution of the CME have advanced to deal with time-varying reaction rates and second-order reactions for specific systems. The most recent theoretical results state that for systems that satisfy a dynamically restricted complex balance condition (DR condition), the solution of the CME maintains a Poisson-Product for all times, from the 
initial time. The solution achieves this by deriving time-varying dynamics of the mean number of each molecular species, 
as described by a set of deterministic ordinary differential equations (ODEs). However, the DR condition requires that all reactant pairs are balanced between reactants and products. This is strongly restrictive, and limits
not only to the form of the reactions but also the reaction rates and initial conditions. 

By contrast, the results we have presented in this work represent a solution of the CME in its entirety, and are free of the 
restrictions of the DR condition. This is achieved by 
allowing the mean molecular numbers to be stochastic processes,
rather than having deterministic dynamics; the solution of the CME is obtained by averaging over the stochasticity. 
Based on the exact solution of the CME, we derived the full distribution of the FPT, and developed a numerical scheme, 
based on Pade-approximants, that allow its rapid computation. 

Apart from demonstrating the theoretical value of our results, we have shown its practical value. 
Our method can accurately compute the FPT distributions of nonlinear systems with A+B → C type of second-order reactions. The results are much more accurate than state-of-the-art linear mapping approximation methods (LMA), which we specifically chose
as a comparison, because it is a recent method that can be applied to time-varying systems, as can our results. 
We further demonstrated our method's applicability to various networks, including genetic regulation networks and multi-step reaction pathways. This exhibits the potential of our approach to real-world applications. 

While our analysis represents theoretical progress, it is limited to systems that have only one second-order reaction, and what is derived is the FPT of this particular second-order reaction. Therefore, our analysis cannot be applied to biochemical systems that 
do not conform to the description given in (<ref>). One such system is the widely-used enzymatic process given by
the Michaelis Menten reaction, where a first-order reaction follows a reversible second-order reaction. 
We plan to extend our analysis to such systems. 

Overall, our work is a step forward in the theoretical derivation of exact solutions of full FPT distributions for biochemical networks with second-order reactions. Treating the mean molecular numbers as stochastic, and then averaging the result is, we believe, a new approach to theoretical derivations of the solution of CMEs. Indeed, this approach is the key 
that enabled us to derive theoretical solutions of the CME for more general systems, and this may represent
a new way of analyzing more general nonlinear biochemical networks.

plain

 


§ 
0.3em section.


  §.§.§ 
0.3em  


§ PROOF OF THEOREM 1

In this appendix we give a proof of Theorem 1.



 §.§ General procedure to solve the CME, assuming a solution of Poisson-product form


For a general second-order biochemical system, as defined in (<ref>), the CME can be written as (<ref>):  


    d P(X, t )/ d t =∑_k=0^M[ P(X-v_k, t ) c_k(X-v_k,t)- P(X, t ) c_k(X,t)],

where P(X, t ) is the probability of the system being in state X and time t. 
The quantity c_k(X,t) equals the propensity of the k'th reaction at time t, a_k(t)X^Y_k,  where a_k is the rate constant and, for any vectors with d components, a vector to the power of a
vector is defined by 𝐮^𝐯 def=∏^d_i=1u_i^v_i, with the convention that 0^0=1.
The quantity v_k is the transition vector of the k'th reaction.



A key question is: if the initial condition of (<ref>) has a Poisson product form:

    P(X,0)   =∏_i=0^N λ_0(0)^x_i/x_i!exp(-λ_i(0)),

then under what condition does the solution of (<ref>) remain of Poisson product form, i.e. 

    P(X,t)= ∏_i=0^N λ_i(t)^x_i/x_i!exp(-λ_i(t)),



We define 
    λ(t)=[λ_1(t),…
    ,λ_N(t)]^⊤
 be a column vector containing the mean of all numbers of the species in S at
any time t (t≥0), and 
    Λ(t)= [[ λ_0(t);   λ(t) ]]
 be a column vector containing the mean of all numbers of the species in S^* at
any time t (t≥0).





Proceeding by substituting the Poisson product form of (<ref>) into 
 (<ref>), the left hand side follows by the Chain rule: 


    d P(X, t )/ d t=P(X,t) ×∑_i=0^N λ_i'(t) (X!/(X-e_i)!Λ(t)^-e_i-1 ),

where e_i is a vector within which only element i is 1, and all other elements are zero.

The right hand side of (<ref>) is given by:


    ∑_m=0^M[ P(X-v_m, t ) c_m(X-v_m,t)- P(X, t ) c_m(X,t)]
    
       = P(X,t) [ ∑_i=0^N K_i(t) (X!/(x-e_i)!Λ(t)^-e_i-1 ).
       .+ ∑_i=0^N∑_j=0^N K_ij(t) ( X!/(X-e_i-e_j)!Λ(t)^-e_i-e_j -1 ) ],

where 

    K_i(t)   =∑_m:Y'_m=e_ia_m(t)Λ^Y_m-∑_m:Y_m=e_ia_m(t)Λ^Y_m,
    
    K_ij(t)   =∑_m:Y'_m=e_i+e_ja_m(t)Λ^Y_m-∑_m:Y_m=e_i+e_ja_m(t)Λ^Y_m.



For the solution of (<ref>) to remain Poisson-product form over time, we need to equate the right-hand sides of (<ref>) and (<ref>). The result is:


    ∑_i=0^N λ_i'(t) (X!/(X-e_i)!Λ(t)^-e_i-1 )
       =∑_i=0^N K_i(t) (X!/(X-e_i)!Λ(t)^-e_i-1 ) 
       +∑_i=0^N∑_j=0^N K_ij(t) ( X!/(X-e_i-e_j)!Λ(t)^-e_i-e_j -1 ).


In general, this is a complicated nonlinear problem to solve for the λ_i. In particular, the left-hand side contains only first-order terms of Λ, whereas the right-hand side contains second-order terms of Λ. 



 §.§ Solution that leads to the dynamical and restricted complex balance (DR) condition


One form of solution of (<ref>) occurs when all of the K_ij coefficients are zero (K_ij = 0 for all
i,j).

Anderson et al. assumed this and found the condition under which this holds (namely the DR condition)<cit.>.

 

DR condition: For any higher-order reaction complexes z·S, where the stoichiometric vector z satisfies z∈ℕ^N+1_≥ 0 and z_1≥ 2, then the DR condition is:

    ∑_m:y_m=za_mΛ^z(t)=∑_m:y'_m=za_mΛ^Y_m(t),

where the sum on the left is over those reactions where z·S are reactants, and the right is over those reactions where z·S are products. 


When K_ij = 0 for all i and j we can write Eq. (<ref>) as ∑_i=1^N(λ_i^'(t) -K_i(t))(  x!/(x-e_i
)!λ(t)^-e_i-1) =0 and a solution arises from 
λ_i^'= K_i(t).

The dynamics of Λ is then governed by a set of deterministic ODEs:

    d/ d tΛ   =∑_m=0^Ma_mΛ^Y_m(Y_m'-Y_m).


The DR condition says that any higher-order reactant pairs within the system should be balanced. It is very restrictive for the system to satisfy: it requires that the higher-order reactant pairs remain the same from the reactants to the products and restricts the system's reaction rate coefficients. 

Next, we present our method of solution that does not constrain the system by the DR condition. 



 §.§ General method of solution  - without requiring the DR condition


In this work we solve (<ref>) in its entirety, without being restricted by the DR condition. Thus, instead of requiring all K_ij to be zero, as applies under a DR condition, we allow 
the sum, that contains the K_ij, to be present
in (<ref>). We incorporate the effects of the sum by 'upgrading' Λ to be a 
stochastic process, which we shall subsequently average over. 

To proceed, we use Ito's rule: 
    df(λ)=f'(λ) dλ+1/2f”(λ)( dλ)^2,
 when we differentiate the Poisson-product form of (<ref>), and obtain

    d P(X, t )=P(X,t) ×[∑_i=0^N  (X!/(X-e_i)!Λ(t)^-e_i-1 ) dλ_i(t).
       +.1/2∑_i=0^N∑_j=0^N  (X!/(X-e_i)!Λ(t)^-e_i-1 )(X!/(X-e_j)!Λ(t)^-e_j-1 )  dλ_i(t) dλ_j(t)].




  §.§.§ Determination of the stochastic process, Λ
  

Ultimately, we will equate the averaged right-hand side of (<ref>) with 
the right hand side of (<ref>) to determine an SDE that governs the dynamics of the stochastic process, Λ. 



We proceed as follows:  




  * We assume the following form of the SDE that governs the dynamics of Λ
(for reasons that will be made clear, shortly):


    d([   λ; λ_0 ])=([   b(Λ,t); b_0(Λ,t) ]) dt + ([   σ(Λ,t); σ_0(Λ,t) ]) dW_t,

where b(Λ,t), b_0(Λ,t) are drift terms, while σ(Λ,t) and σ_0(Λ,t) are diffusion terms. In the above equation: b(Λ,t) is a N component column vector; b_0(Λ,t) is a scalar; W_t is an 2 component column vector 
of independent Wiener processes; σ(Λ,t) is an N × 2 matrix; 
σ_0(Λ,t) is a row vector with 2 elements.


  * We determine b(Λ,t) and b_0(Λ,t) by equating them to K_i(t) in (<ref>), as they are both determined by the zero- and first-order reactions. We obtain:


    b(Λ,t)=∑_m: Y_m'_1≤ 1,Y_m_1≤ 1a_mλ^y_m(y_m'-y_m)=ℳ_1λ+ℳ_2,
       b_0(Λ,t)=a_0λ^y_0(y'_0,0-y_0,0)=a_0λ_1λ_2.


Here ℳ_1 is an N× N matrix with components
η_i,j^1, ℳ_2 is an N component column vector with
components η_i^2, while y_m,i^' and y_m,i are the i'th
component of y_m^' and y_m, respectively, and


    η_i,j^1=∑_m:y_m=e_ia_m(t)(y_m,j^'
    -y_m,j),
       η_i^2=∑_m:y_m=0a_m
    (t)(y_m,i^'-y_m,i).



  * We now explicitly write W_t as W_t=[W^1_t,W^2_t]^⊤, then
σ(Λ,t) · dW_t can then be written as: 

    σ(Λ,t) · dW_t    = σ_1(Λ,t) dW^1_t+ σ_2(Λ,t) dW^2_t

where σ_1(Λ,t) and σ_2(Λ,t) are N component column vectors.

We adopted this particular format because the only second-order reactions that occur
within the system are of the type
A + B→ C, with two different reactants. The form of (<ref>) ensures that K_ij in (<ref>) satisfies: 

 

  * K_ij = 0 if i>2 or j>2, which means that species S_k for k=3,4… are not involved in the second-order reaction.


  * K_ij = 0 if i=0 or j=0, which means that S_0 is not a reactant in any one reaction. It follows that 
    σ_0(Λ,t)=0.



  * K_ij = 0 if i=j, which means that there are no square terms of the form x!/(x-2e_i)!λ(t)^-2e_i in the right-hand side of (<ref>).



When  σ_1(Λ,t) and σ_2(Λ,t) satisfy certain relationships, the non-zero square terms of λ can be eliminated in the right-hand side of (<ref>) (see below).



  * K_12 = K_21 = -a_0/2λ_1λ_2, which means that the matrix K = {K_ij} is symmetric, and its components are all products of different components of Λ. To satisfy this constraint, σ_1(Λ,t) and σ_2(Λ,t) are set as linear functions of λ:


    σ_1(Λ,t)= 𝒩_1λ
       σ_2(Λ,t)=𝒩_2λ,

where 𝒩_1 and 𝒩_2 are N × N square matrices. 







  * We determine the matrices 𝒩_1 and 𝒩_2, by equating terms of the form X!/(X-e_i-e_j)!Λ(t)^-e_i-e_j in (<ref>) with the K_ij terms in (<ref>). This leads to an algebraic equation, the solution of which sets the forms of 𝒩_1 and 𝒩_2: 

    1/2(σ_1(Λ,t)·σ_1(Λ,t)^⊤+σ_2(Λ,t)·σ_2(Λ,t)^⊤)=K,

where K is the matrix with elements K_ij (i≥ 1, j≥ 1) in (<ref>). 𝒩_1 and 𝒩_2 are determined as that shown in (<ref>).



At this point, the dynamics of λ_0 and λ are determined by: 


    dλ=(ℳ_1λ+ℳ_2) dt + 𝒩_1λ dW^1_t+𝒩_2λ dW^2_t,
        dλ_0=λ_1λ_2 dt.

However, on substituting (<ref>) into (<ref>), terms are generated in the right-hand side of (<ref>) but are absent in (<ref>). These terms are λ_1λ_2x!/(x-e_1)!λ(t)^-e_1 and λ_1λ_2x!/(x-e_2)!λ(t)^-e_2. 

To eliminate these terms in the right-hand side of (<ref>). We modified the Poisson-product form to P̃(X,t),
as given by


    P̃(X,t)=Λ(t)^X/X!exp(-Λ(t))exp(λ_0(t)+λ_S(t)),

where λ_S(t) is a stochastic process to be determined.

Note that P̃(X,t) is not normalized to unity, and thus is not a probability distribution. 
However, we will prove later that its average, <P̃(X,t)>, is a distribution; it is the solution of the CME 
and obeys the Poisson-product form initial condition of (<ref>).



  §.§.§ Determination of λ_S
  



Changing the Poisson product form from P(X,t) to P̃(X,t) does not change the format of (<ref>):


    ∑_m=0^M[ P̃(X-v_m, t ) c_m(X-v_m,t)-P̃(X, t ) c_m(X,t)] dt
       =P̃(X,t) [ ∑_i=0^N K_i(t) (X!/(X-e_i)!Λ(t)^-e_i-1 ). 
       .+ ∑_i=0^N∑_j=0^N K_ij(t) ( X!/(X-e_i-e_j)!Λ(t)^-e_i-e_j -1 ) ] dt

where K_i(t) is the drift term of dλ_i.

However, when we differentiate P̃(X,t) using Ito's rule, (<ref>) is changed to a more complex format, as in (<ref>).

Setting I_e=exp(-Λ) and I_p=Λ^X/X! allows us to write
P̃(X,t)= I_eI_pexp(λ_0(t)+λ_S(t)). The quantity dP̃(X, t ) then becomes:


    dP̃(X, t )=[I_p d(I_eexp(λ_0(t)+λ_S(t)))+I_eexp(λ_0(t)+λ_S(t)) dI_p]+exp(λ_0(t)) dI_p d(I_eexp(λ_S(t)))
       =P̃(X,t) ×[∑_i=0^N  K_i(t)(X!/(X-e_i)!Λ(t)^-e_i-1) dt.
       +∑_i=0^N∑_j=0^N K_ij(t) ( X!/(X-e_i-e_j)!Λ(t)^-e_i-e_j -1 ) dt
    
       +a_S[x_1-x_2] dW^1_t+ ia_S[x_1+x_2] dW^2_t
       -( a_S λ_1-a_S λ_2 ) dW^1_t - (  ia_S λ_1+ ia_S λ_2 ) dW^2_t
       .+ dλ_S+1/2( dλ_S)^2+∑_i=1^2 dλ_i dλ_S-a_0λ_1λ_2 dt]
       +exp(λ_0(t)) dI_p d(I_eexp(λ_S(t))),

where

    dI_p d(I_eexp(λ_S(t)))
    =P̃(X,t) ×[a_0  (λ_2x_1+λ_1x_2 )  dt.
       .
    +1/2(x_1/λ_1 dλ_1(t) dλ_S(t)+x_2/λ_2 dλ_2(t) dλ_S(t) ) 
    ].


We then determine λ_S by equating the coefficients of the dt terms on the right-hand sides of (<ref>) and (<ref>). This leads to constraints that λ_S has to satisfy: 


    
  * dλ_S is a function of dW_t, as (<ref>) and (<ref>)  already have 
    corresponding d t terms

    
  * setting dλ_S = ξ_1(λ,t) dW^1_t + ξ_2(λ,t) dW^2_t 
    allows dI_p d(I_eexp(λ_S(t))) to be zero, i.e., the right hand side of (<ref>) is zero: 
    

    a_0  (λ_2x_1+λ_1x_2 )  dt
    +1/2(x_1/λ_1 dλ_1(t) dλ_S(t)+x_2/λ_2 dλ_2(t) dλ_S(t) ) =0


    
  * dλ_S should be a quadratic function of dλ_S:


    1/2( dλ_S)^2+∑_i=1^2 dλ_i dλ_S-a_0λ_1λ_2 dt=0.
 



 (<ref>) and (<ref>) determine that the solution of dλ_S is the sum of the diffusion terms of dλ_1 and dλ_2:


    dλ_S   = ( a_S λ_1-a_S λ_2 ) dW^1_t + (  ia_S λ_1+ ia_S λ_2 ) dW^2_t.


The introduction of λ_S(t) allows the elimination of the cross-product terms in (<ref>), 
corresponding to terms of the form λ_1(t)λ_2(t)X!/(X-e_1)!Λ(t)^-e_1 and λ_1(t)λ_2(t)X!/(X-e_2)!Λ(t)^-e_2 in (<ref>). 

We now substitute λ_S into (<ref>) to obtain 


    dP̃(X, t )=P̃(X,t) ×[∑_i=0^N  K_i(t)(X!/(X-e_i)!Λ(t)^-e_i-1) dt.
       +∑_i=0^N∑_j=0^N K_ij(t) ( X!/(X-e_i-e_j)!Λ(t)^-e_i-e_j -1 ) dt
       +.a_S[x_1-x_2] dW^1_t+ ia_S[x_1+x_2] dW^2_t]


To prove Theorem <ref>, we needed to equate the right-hand sides of (<ref>) and (<ref>), which means that the stochastic diffusion terms involving dW_t^1 and dW_t^2 should be eliminated. Next, we prove that the average of P̃(X,t) eliminates these stochastic terms, and <P̃(X,t)> satisfies (<ref>), and hence is the solution of the CME in (<ref>).



  §.§.§ Eliminating stochastic diffusion terms in (<ref>)
 

Averaging both sides of (<ref>) yields


    d<P̃(X, t )> =<P̃(X,t) ×[∑_i=0^N  K_i(t)(X!/(X-e_i)!Λ(t)^-e_i-1) dt..
       +..∑_i=0^N∑_j=0^N K_ij(t) ( X!/(X-e_i-e_j)!Λ(t)^-e_i-e_j -1 ) dt ]> 
       +<a_SP̃(X,t)[x_1-x_2] dW^1_t+ ia_SP̃(X,t)[x_1+x_2] dW^2_t>


Set σ̃_1(X,λ,t)  = a_SP̃(X,t)[x_1-x_2]  and σ̃_2(X,λ,t) = ia_SP̃(X,t)[x_1+x_2], we then prove that: <σ̃_1(X,λ,τ) dW^1_τ+σ̃_2(X,λ,τ) dW^2_τ> = 0 (see Proposition <ref>), and hence the stochastic diffusion terms in (<ref>) are eliminated.  


Define t_n be a series of stopping times t_n= inf{τ | λ(τ)≥ n} for n=1,2,….

Then, for any t∈[0,t_n], the following expectation is zero:
    
    <σ̃_1(X,λ,t) dW^1_t+σ̃_2(X,λ,t) dW^2_t>=0,



 <ref>:

From the definition of t_n, when t∈[0,t_n], λ(τ) is bounded, and hence   |σ̃_1(X,λ,t)| and |σ̃_2(X,λ,t)| are bounded.

Thus, Z(t)=∫^t_0σ̃_1(X,λ,τ) dW^1_τ+σ̃_2(X,λ,τ) dW^2_τ is a martingale, and according to the properties of martingales, E(Z(t)) = Z(0)= 0, which is in equivalent to (<ref>).



Proposition <ref> guarantees that the stochastic diffusion term is eliminated when t<t_n. Next, we prove that <P̃(X, t )> is the solution of the CME in (<ref>).






  §.§.§ <P̃(X, t )> is the solution of the CME in (<ref>)


We first prove that <P̃(X, t )> is the solution of the CME in (<ref>) when t<t_n, then we generalise the result as t_n→∞.

When t< t_n, (<ref>) holds, and <P̃(X, t) > is the solution of the CME (<ref>), because (<ref>) becomes:


    d<P̃(X,t) > = 
          <P̃(X,t) ×[∑_i=0^N  K_i(t)(X!/(X-e_i)!Λ(t)^-e_i-1) dt..
       +..∑_i=0^N∑_j=0^N K_ij(t) ( X!/(X-e_i-e_j)!Λ(t)^-e_i-e_j -1 ) dt ]> 
       =∑_m=0^M[ <P̃(X-v_m, t )> c_m(X-v_m,t)-<P̃(X, t )> c_m(X,t)] dt


Define t∧ t_n=min(t,t_n), and denote 𝒫(X,t) as the solution of the CME in (<ref>). (<ref>) means that 𝒫(X,t∧ t_n)=<P̃(X, t∧ t_n ) > at time t∧ t_n. 

Next, we prove that when n→∞, we have the limit lim_n→∞ t∧ t_n =t, leading to


    <P̃(X,t)>=lim_n→∞<P̃(X,t∧ t_n)>=lim_n→∞𝒫(X,t∧ t_n)=𝒫(X,t),

where the last equal sign stems from 𝒫(X,t), as the solution of CME (<ref>), being a continuous function.

By probability measure, the limit lim_n→∞t∧ t_n=t is equivalent to 

    lim_n→∞ P(t_n≤ t)→ 0.


From the definition t_n= inf{t | λ(t)≥ n}, (<ref>) is equivalent to 

    lim_n→∞ P(max_0≤τ≤ tλ(τ)≥ n)→ 0.


Next, we use Lemma <ref> to prove (<ref>).

The linear SDE in (<ref>), that  governs λ, satisfies the condition of Lemma <ref>, i.e, (<ref>) because:

1) (<ref>) can be rewritten as

    dλ=b(Λ,t) dt + σ(Λ,t) dW_t,

where b is an N column vector of linear functions defined in (<ref>), and σ  is an 
N× 2 matrix of linear functions defined in (<ref>).

2) Linearity of b and σ lead to b(t,y)^2+σ(t,y)^2 being quadratic functions of the vector 
y∈ℝ^N. Thus, (<ref>) holds, which means


    <max_0≤τ≤ tλ(τ)^2> ≤ K(λ(0)) e^Ct,



where C is a constant, and K(λ(0)) is a constant that depends on λ(0).

Then by Chebyshev's inequality

    P(max_0≤τ≤ tλ(τ)≥ n)≤ P(|max_0≤τ≤ tλ(τ)- μ|≥ n-μ)
       ≤<max_0≤τ≤ tλ(t)^2>/(n-μ)^2≤K(λ(0)) e^Ct/(n-μ)^2,

where μ=<λ(t_m)> and λ(t_m)=max_0≤τ≤ tλ(τ).

Because 

    μ^2=<λ(t_m)>^2≤<λ(t_m)^2>≤ K(λ(0)) e^Ct,

we have the limit

    lim_n→∞K(λ(0)) e^Ct/(n-μ)^2=0.

Thus, (<ref>) leads to (<ref>)


 
(Problem 5.3.15 of<cit.>)
    Suppose b_i(t,y) and σ_ij(t,y); 1≤ i ≤ d, 1≤ j ≤ r, are progressively measurable functionals from [0,∞)× C[0,∞)^d into ℝ satisfying
    
    b(t,y)^2+σ(t,y)^2≤ K( 1+y^2);   
    ∀ 0≤ t< ∞, y∈ℝ^d,

    where K is a positive constant. If (X,W), (Ω,ℱ,P), {ℱ_t} is a weak solution to the SDE 
    
    dX=b(t,X) dt+σ(t,X) dW_t,

    with <X_0^2m><∞ for some m>1, then for any finite time T>0, we have
    
    <max_0≤ s≤ tX_t^2m>≤ C( 1+<X_0^2m>) e^Ct;0≤ t ≤ T,

    where C is a positive constant depend only on m, T, K and d.




§  WHY CALCULATE INSTEAD OF ?


In this appendix we answer why we calculated <λ_S(t)^n> instead of <exp(λ_S).

In the main text, we calculated <λ_S(t)^n> because it can be determined from solution of a finite number of ODEs, 
whereas calculating <exp(λ_S)> involves solving an infinite-dimensional system of coupled ODEs.

The details are as follows.  

A natural idea to calculate <exp(λ_S)> is to derive its governing ODE by differentiating exp(λ_S):

    dexp(λ_S)=exp(λ_S) dλ_S+1/2exp(λ_S)( dλ_S)^2.
 

Substituting dλ_S and dλ_S^2 in (<ref>) with the SDEs in (<ref>), some polynomial functions of λ show up in the right-hand side of (<ref>). The new polynomial functions have a general format of  λ^lexp(λ_S), where l=[l_1,…,l_N] is a vector of positive integer numbers.

To determine terms like λ^lexp(λ_S), we kept deriving their ODEs by differentiating them. This process results in terms with higher exponents of λ, i.e. λ^l'exp(λ_S), where l'_1>l_1: 


    d(exp(λ_S)λ^l)=exp(λ_S)λ^l dλ_S+l_1exp(λ_S)λ^l-e_1 dλ_1+l_2exp(λ_S)λ^l-e_2 dλ_2+…
       +1/2exp(λ_S)λ^l( dλ_S)^2+l_1(l_1-1)/2exp(λ_S)λ^l-2e_1( dλ_1)^2+…
       + l_1exp(λ_S)λ^l-e_1 dλ_S dλ_1+l_2exp(λ_S)λ^l-e_2 dλ_S dλ_2+l_1l_2exp(λ_S)λ^l-e_2-e_2 dλ_1 dλ_2+…,


where e_j (j=1,2,…,N) is a N dimensional vector, whose j'th component equals one, and all other components equal zero. 

It can be seen from the right-hand side of (<ref>), that the differentiation process always results 
in nonlinear terms with higher exponents of λ, whose expectations are not zero, e.g. 
the terms like exp(λ_S)λ^l dλ_S^2. 

As a result, the expectations of (<ref>) and (<ref>) compose a set of non-closed-form equations, and directly solving <exp(λ_S)> is dogged by infinite-dimensional ODEs.

To avoid solving infinite number of ODEs for <exp(λ_S)>, we approximated it by calculating <λ_S(t)^n>. This is achieved by introducing a new function H(s,t) = <exp(s*λ_S)>, and one may approximate H(s,t) by its Taylor expansions: 


    H(s,t)≈ T_Ñ(s,t)=s^n/n!×∂^n/∂ s^nH(s,t)|_s=0
    =∑_n=0^Ñs^n/n! <λ_S(t)^n>.


 (<ref>) means that we only need to determine <λ_S(t)^n> for every n to determine H(s,t), and H(1,t) is the required FPT distribution. 

For practical reasons, we can calculate <λ_S(t)^n>, with n less than a predefined integer, i.e. n=0,1,2…,Ñ (Ñ is the highest order of calculated moments). Since exp(sλ_S(t)) is a transcendental function, we used a Padé approximant to control
the approximation errors. We calculated the Padé approximant via the extended Euclidean algorithm. 

It turns out that calculating <λ_S(t)^n> is much easier than calculating <exp(λ_S)>, because the linearity of the SDEs in (<ref>) guarantes that <λ_S(t)^n> is governed by a set of finite dimensional ODEs. 

We derived the ODEs governing <λ_S(t)^n> through a similar process as that of <exp(λ_S)>. We
started differentiating <λ_S(t)^n> by Ito's rule: 

    d(λ_S^n)   =nλ_S^n-1 dλ_S+ n(n-1)/2λ_S^n-2( dλ_S)^2.


Substituting dλ_S and dλ_S^2 in (<ref>) with the SDEs in (<ref>), the right-hand side of (<ref>) would only be polynomial functions of λ_S^l_0λ^l. The linearities of the SDEs guarantee that
dλ_S can be represented by polynomial functions of λ with orders no more than one. Likewise, dλ_S^2 can be represented with polynomial functions of λ with orders no more than 2. As a result, the polynomial order of λ_S^l_0λ^l is not higher than n, i.e. l_0+l_1 ≤ n. 

Furthermore, the differentiation of <λ_S^l_0λ^l> does not increase its polynomial order. Because of the Itô's rule, 


    d(λ_S^l_0λ^l)=l_0λ_S^l_0-1λ^l dλ_S+l_1λ_S^l_0λ^l-e_1 dλ_1+l_2λ_S^l_0λ^l-e_2 dλ_2…
       +l_0(l_0-1)/2λ_S^l_0-2λ^l( dλ_S)^2+l_1(l_1-1)/2λ_S^l_0λ^l-2e_1( dλ_1)^2
       + l_0l_1λ_S^l_0-1λ^l-e_1 dλ_S dλ_1+l_0l_2λ_S^l_0-1λ^l-e_2 dλ_S dλ_2+l_1l_2λ_S^l_0λ^l-e_2-e_2 dλ_1 dλ_2…,


where e_j (j=1,2,…,N) is a N dimensional vector, whose j'th component equals one, and all other components equal zero. Similarly, because of the linearity of (<ref>), the polynomial orders of λ_S and λ of all terms in (<ref>) do not increase.

Therefore, (<ref>) and (<ref>) compose a set of closed-form equations, and <λ_S(t)^n> can be solved by a set of finite-dimensional ODEs. 




















































    






































  






























































   




   


















































































































    






























































