Step 0: Inputs: Initial node N0, aμ, vμ, tol τm, tf, tstep Initialize: $${\boldsymbol{Bound}}_{\boldsymbol{lower}}={\boldsymbol{X}}_{\boldsymbol{K}},{\boldsymbol{b}}_{{\boldsymbol{N}}^i,\boldsymbol{N}}^{\hbox{'}}{\left(\boldsymbol{X}-{\boldsymbol{v}}_{\mu}\right)}^{\hbox{'}}={\boldsymbol{P}}^{\left(\boldsymbol{t}\right)}\left({\boldsymbol{X}}_{\mathbf{0}}\right),\boldsymbol{A}=\left[\right]$$ Step 1: Start from parent node Ni = (X0, ƌl) ← Current State of the system at td, Step 2: Flag the current node as explored, update A and add the state Xi in the domain so that; if 1 − IT exp (t. Aj). P(t)(X0) ≥ τm(leak) holds true go to Step 3; else stop the algorithm Step 3: Sort exp(t. Aj). P(t)(X0) and shift the set of states in $${\mathbf{X}}_K^{\hbox{'}}\;\mathrm{at}\;{t}^{\hbox{'}}$$ having smallest probabilities, if P(t)(XK) ≥ τm(leak) > $${\mathrm{P}}^{(t)}\left({\mathbf{X}}_K^{\hbox{'}}\right)$$ and at td update $${\mathbf{X}}_K\leftarrow {\mathbf{X}}_K-{\mathbf{X}}_K^{\hbox{'}}$$ Step 4a: Extend the graph dictionary Dict by vμ(Xi(t)) by 1 level to check all the nodes nj = (Xj, ƌl , Ͼ$${}_{N_i,{N}_i^{\hbox{'}}}\left.\left(\mathit{\min}\right)\right)$$ adjacent to Ni: Boundupper ← RM(Boundlower) reachable by exactly RM reactions (from fast to slow) having Ͼ$${}_{N_i,{N}_i^{\hbox{'}}}\left(\mathit{\min}\right)$$. If nK = (XK, ƌl, Ͼ$${}_{N_i,{N}_i^{\hbox{'}}}\left.\left(\mathit{\min}\right)\right)$$ be the set of adjacent nodes such that nK ∈ nJ then go to next Step, Step 4b: Compute the BLNP function for nK∈ Boundupper: $$b\left(\left.{N}_{N1,.. NM}\right|{b}_{1,N.\dots {N}^{\hbox{'}}N}^{\hbox{'}}\right)={P}_{N,1.\dots N,N\hbox{'}}\left(\omega \right)\ast {b}_{N^{\hbox{'}},N}^{\hbox{'}}{\left(X-{v}_i\right)}^{\hbox{'}}$$ Step 5a: If nK = (XK, ƌl, Ͼ$${}_{N_i,{N}_i^{\hbox{'}}}\left.\left(\mathit{\min}\right)\right)\in domain$$, then update the values of the set of states XK present in domain and take domain ← domainprevious ∪ domain and go back to Step 1; else If nK = (XK, ƌl, Ͼ$${}_{N_i,{N}_i^{\hbox{'}}}\left.\left(\mathit{\min}\right)\right)\notin domain$$, then add it to the queue in order, according to reachability and go to the next Step, Step 5b: sort $$b\left({N}_{N1,.. NM}\left|{b}_{N1,.. NM}^{\hbox{'}}\right.\right)$$ in descending order and update $$queue\leftarrow \left( queue;b\left({N}_{N1,.. NM}\left|{b}_{1,N.\dots {N}^{\hbox{'}},N}^{\hbox{'}}\right.\right)\right)$$ Step 6: Pull out the nodes nK = (XK, ƌl, Ͼ$${}_{N_i,{N}_i^{\hbox{'}}}\left.\left(\mathit{\min}\right)\right)$$ from the queue in order and add the set of states XK in the domain as domain ← domain + XK and take domainprevious ∪ domain, then go back to Step 1, Output: domain with probable states 